Something is wrong with ode45 integration

[cdsousa]

Check this code where the two integrations differ only in the final time:

using ODE

f(t,y) = 1.0 <= t < 2.0 ? 1.0 : 0.0

_, y = ode45(f, 0.0, 0.0:0.1:16.0);
println(maximum(y))

_, y = ode45(f, 0.0, 0.0:0.1:17.0);
println(maximum(y))

I get 0.999068012862035 for the first, which is correct, but I get 0.0 for the second... I tried this in Julia 0.6 both locally and in JuliaBox.

[ChrisRackauckas]

This is not a problem with ode45 per se (though there are lots of issues that it has, as evidenced by how it diverges for most of the problems in DiffEqBenchmarks.jl). The issue here is that your f is discontinuous. There's no guarantee that the solver will step in the interval [1.0,2.0]. It accidentally does so in the first case because of how the adaptive maximum dt corresponds to the interval length. Because of this, in the second time it never hits a point where the derivative is non-zero.

Theoretically, there is no way to handle a discontinuity without order loss unless you hit the exact timepoints of the discontinuity. ode45 does not have that kind of fine grained control of the solver, so it won't be able to handle discontinuities like this.

This is one of the many reasons why ODE.jl has been deprecated in favor of DifferentialEquations.jl and its ODE solvers in OrdinaryDiffEq.jl. These solvers have ways to declare discontinuities. For example, in your problem, you tell the solver to hit those timepoints via the tstops argument:

using DifferentialEquations # or OrdinaryDiffEq
f(t,y) = 1.0 <= t < 2.0 ? 1.0 : 0.0
prob = ODEProblem(f,0.0,(0.0,16.0))
sol = solve(prob,Tsit5(),tstops=[1.0,2.0])
using Plots; pyplot(); plot(sol)

println(sol[end]) # 0.9976278099151336, gets closer to 1 by adjusting the tolerances
sol = solve(prob,Tsit5(),tstops=[1.0,2.0],abstol=1e-6,reltol=1e-6)
println(sol[end]) # 0.9999984055205092

If the discontinuity is not known in advance, then the timepoint can be hit using a ContinuousCallback. See the docs for details.

To fix this generally and allow user inputs for discontinuities, it would require a pretty substantial rewrite of the solver. That said, OrdinaryDiffEq.jl already includes these features which is why this library is deprecated, and so I am closing this as a "won't fix". On that note I am closing the issue, but please feel free to ask any additional questions.

[cdsousa]

Thank you, Chris, for your very insightful reply! I have little knowledge about ODE solvers and I had no clue about discontinuity issues. I had not noticed that this package was deprecated, I've been using it in an "old" code.

Thank you for your amazing work (I see you have been doing great things with Julia)!