Change specific heaviside() interface to Maxima

[rwst]

Sage interfaces its Heaviside function with a Maxima function named hstep which seems to implement the same function. However, there is no documentation on it, and it is not supported in Runge-Kutta DE computations, nor in integration. This leads to mathematically wrong results when calling desolve_rk4() and integrate() with expressions containing heaviside (but it's working with unit_step).

The ticket should either replace hstep with unit_step in the Maxima interfaces, or remove it altogether with a warning. Alternatively, use substitute_function in desolve_rk4() and integrate(). Additionally, the issue with hstep should be reported upstream.

Upstream: Not yet reported upstream; Will do shortly.

CC: @tscrim @slel @kliem

Component: calculus

Author: Frédéric Chapoton

Branch/Commit: u/chapoton/22850 @ 15df4b5

Issue created by migration from https://trac.sagemath.org/ticket/22850

[mforets]

comment:1

i am curious about the 2nd alternative: remove it altogether with a warning. does it mean that we can write a custom integrate method that fixes integration etc and symbolics?

because if you compare with the similar problematic behaviour with the dirac delta:

sage: integrate(dirac_delta(x-1)*sin(x), x, 0, 2, algorithm='maxima')  # ??
integrate(dirac_delta(x - 1)*sin(x), x, 0, 2)
sage: integrate(dirac_delta(x-1)*sin(x), x, 0, 2, algorithm='sympy')  # ok
sin(1)

then in this other case there is no other "dirac delta" that can be used as backup, in analogy to the unit step. hence in this case i have the impression that the 2nd alternative is the only choice, isn't it?

[rwst]

comment:2

Replying to @mforets:

i am curious about the 2nd alternative: remove it altogether with a warning. does it mean that we can write a custom integrate method that fixes integration etc and symbolics?

because if you compare with the similar problematic behaviour with the dirac delta:

sage: integrate(dirac_delta(x-1)*sin(x), x, 0, 2, algorithm='maxima')  # ??
integrate(dirac_delta(x - 1)*sin(x), x, 0, 2)
sage: integrate(dirac_delta(x-1)*sin(x), x, 0, 2, algorithm='sympy')  # ok
sin(1)

then in this other case there is no other "dirac delta" that can be used as backup, in analogy to the unit step. hence in this case i have the impression that the 2nd alternative is the only choice, isn't it?

It would be if the result (the unevaluated integral) was mathematically wrong, but it's not, so it's not a bug.

[rwst]

comment:3

However, another alternative for the enhancement of dirac_delta integration would then be to always call SymPy.

[mforets]

comment:4

However, another alternative for the enhancement of dirac_delta integration would then be to always call SymPy.

is it possible to dispatch a particular integrator when the symbolic expression has the presence of a given function (such as dirac delta)? like a notion of "weak default"

[rwst]

comment:5

Replying to @mforets:

is it possible to dispatch a particular integrator when the symbolic expression has the presence of a given function (such as dirac delta)? like a notion of "weak default"

Recognition works even in Python:

sage: (1 + sin(dirac_delta(x))).has(dirac_delta(SR.wild()))
True

Changing the integrator would be done in src/sage/symbolic/integration/integral.py. Why not?

[mforets]

comment:6

ok, and this would require some decorator technology?

(just in case, in general i'm -1 to change the default integrator to other different than maxima)

think that a similar trick could be applied to laplace with heaviside for instance

[rwst]

comment:7

I have no idea (I only fiddled once with changing the result after the integrators had a go). Please go ahead.

[mforets]

comment:8

we should also consider extending the available software interfaces to desolvers.py. SymPy has dsolve. Giac has odesolve. one could be interested in choosing some of these, not only for this but for other examples.

for the integral customization, there is some initial effort in the DefiniteIntegral class, at /symbolic/integration/integral.py.

[fchapoton]

Commit: 15df4b5

[fchapoton]

Branch: u/chapoton/22850

[fchapoton]

comment:9

tiny commit

---

New commits:

15df4b5

using unit_step for heaviside

[fchapoton]

Author: Frédéric Chapoton

[fchapoton]

comment:10

here is a minimal fix, that apparently has no wrong side effect. Please review

[tscrim]

comment:11

Interestingly Sage does not define the Heaviside step function at 0. However, one subtle consequence of this is change is that it changes the value for Maxima at 0 Before, we had

sage: maxima(heaviside(0))
1/2

and now we have

sage: maxima(heaviside(0))
0

From reading the wiki page, this can matter for some applications. This fixes how I think most people will use this (i.e., as a distribution as part of an integration problem), but there seems to be a bug in the conversion code that we should fix instead.

I am comparing the following two:

sage: integrate(heaviside(x), x, -1, 1, algorithm='maxima')
sage: integrate(sin(cos(x)), x, -1, 1, algorithm='maxima')

In sr_to_max, both of these add something to the *_op_dict. In the latter case, it is %SIN but in the former it is $HSTEP. When it tries to convert back from %HSTEP, it cannot find it. This leads to an error when integrating.

Interestingly enough, if I do unit_step, then what gets passed back from ECL is $UNIT_STEP.

With your proposed fix, I get an error:

sage: integrate(sin(cos(heaviside(x))), x, -1, 1, algorithm='maxima')
sage: integrate(sin(cos(unit_step(x))), x, -1, 1, algorithm='maxima') # Boom

Very roughly speaking, this is because Sage realizes that something is wrong with the round trip because it ends up in a difference place then where it started. Hence, I don't think we can use this approach.

[tscrim]

comment:12

Something very subtle is going on:

sage: heaviside(x)._maxima_lib_().ecl()
<ECL: ((%HSTEP SIMP) |$_SAGE_VAR_x|)>
sage: integrate(heaviside(x), x, -1, 1, algorithm='maxima')
integrate(heaviside(x), x, -1, 1)

This "works" (the fact it cannot simplify the integral is quite bad on Maxima's part IMO, but that is a separate issue upstream) because it sets the correct conversion name. I am not sure if it is on our side, but this suggests that it is.

[tscrim]

comment:13

I don't quite understand the mechanism of how this works. I also don't know how to verify if this is a bug on Maxima's side or not. At this point I don't know how to fix this. (It really is baffling to me that maxima can handle the unit step function but not the Heaviside in the same way wrt integrals, but I digress.)

[fchapoton]

comment:14

The solution would be to have, both in maxima and sage, one unique function (call it either unit_step or heaviside)

Then we meet the other issue that there is no general agreement on what should be the value at 0. Some say 0, some say 1 and some say 1/2. On my branch, one gets

sage: heaviside(0)._giac_()
1
sage: heaviside(0)._sympy_()
1/2
sage: heaviside(0)._maxima_()
0

[tscrim]

comment:15

Hmmm...probably having the functions heaviside = unit_step would be a possible solution. Especially since unit_step seems to have better behavior in Maxima (for whatever reason...).