Heaviside in numerical resolutions

[mforets]

• numerical ode:

sage: y = var('y')
sage: desolve_rk4(x,y,ics=[0,0],end_points=1,step=0.5) # solution is x^2/2
[[0, 0], [0.5, 0.125], [1.0, 0.4999999999999999]]
sage: desolve_rk4(x*unit_step(1+x),y,ics=[0,0],end_points=1,step=0.5) # OK
[[0, 0], [0.5, 0.125], [1.0, 0.4999999999999999]]
sage: desolve_rk4(x*heaviside(1+x),y,ics=[0,0],end_points=1,step=0.5) # ??
[[0, 0]]

• integration:

sage: integral(unit_step(x), (x, 0, 1))
1
sage: integral(heaviside(x), (x, 0, 1)) # ??
integrate(heaviside(x), x, 0, 1)

• custom numerical value at 0:

sage: heaviside(0)
heaviside(0)
sage: f(x) = heaviside(x, 1/2); f(0)  # new 2nd argument (?)
1/2

CC: @rwst @kcrisman @egourgoulhon @tscrim

Component: symbolics

Keywords: heaviside, integrate

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

[rwst]

comment:3

sage: desolve_rk4(x*unit_step(1+x),y,ics=[0,0],end_points=1,step=0.5) # OK
[[0, 0], [0.5, 0.125], [1.0, 0.4999999999999999]]
sage: desolve_rk4(x*heaviside(1+x),y,ics=[0,0],end_points=1,step=0.5) # ??
[[0, 0]]

This is due to Maxima because with Maxima 5.38.1:

(%i3) rk(x*unit_step(x+1),y,0,[x,0,1,0.500000000000000]);
(%o3)               [[0.0, 0.0], [0.5, 0.125], [1.0, 0.5]]
(%i4) rk(x*hstep(x+1),y,0,[x,0,1,0.500000000000000]);
(%o4)                            [[0.0, 0.0]]

[rwst]

comment:4

I cannot even find online documentation on hstep, nor a mention of the Heaviside function. I don't think we can rely on Maxima having implemented it, so we cannot provide DE or integral services with it.

[rwst]

comment:5

This is now #22850.