Accuracy for discontinuous right hand side

[Stivanification]

I am trying to solve a system with discontinous derivatives numerically, but can't get the accuracy I had hoped for. This is my code:

using DifferentialEquations

function periodic(t)
    C_max = 36.8

    if t < 0
        throw(BoundsError)
    elseif t%150 <= 25
        C = C_max
    elseif t%150 > 25 && t%150 <= 50
        C = -C_max
    else
        C = 0
    end
    C
end

function inaccurateperiodic(du,u,p,t)
du[1] = 4.477e-7 * u[2] * abs(periodic(t))
du[2] = (-1/8280)*(periodic(t))
end

u0 = [0.0;0.55]
tspan = (0.0,150.0)
prob = ODEProblem(inaccurateperiodic,u0,tspan)
sol = solve(prob,abstol=e-12,reltol=e-12)

periodic(t) is a step function so u[2] is a linear function as depicted below:

I would like my code to be able to work with periodic(t) as a step function, but it seems that this doesn't work out as intended, because what I'm getting is

sol(150)
2-element Array{Float64,1}:
 -0.000112729
  0.52394

Now sol(150)[1] should be positive and sol(150)[2] should be 0.55, so obviously this is very inaccurate despite the relatively small tolerance values. Part of the issue seems to be that there are very few interpolation points in the solution, and it seems that even if I prescribe something like dt = 0.0001, only 10 interpolation points are used (at least that is how I interpretate the output).

[ChrisRackauckas]

Hey, no numerical method for solving differential equations can exactly resolve discontinuities. You have to make sure that you exactly hit the time point of discontinuity in order to not have order loss. However, there are two tools which are helpful here. If you know in advance where a discontinuity will be, use tstops to make sure the integrator steps to that time. If you don't know where a discontinuity will be but have an implicit equation for it, use a ContinuousCallback.

Your case is the first, so adding a tstop at the time points of discontinuity will make them resolve accurately:

sol = solve(prob,Tsit5(),tstops = [25,50])
using Plots; gr()
plot(sol,denseplot=false)

and then at high accuracy:

sol = solve(prob,Tsit5(),abstol=1e-12,reltol=1e-12,tstops = [25,50])
using Plots; gr()
plot(sol)

julia> print(sol[end])
[0.000407308, 0.55]

Another thing you can do is help it resolve on its own by lowering dtmax. This makes the integrator steps not grow too fast when the derivative is constant, and makes it more readily catch the discontinuities:

sol = solve(prob,Tsit5(),abstol=1e-12,reltol=1e-12,dtmax=1)
using Plots; gr()
plot(sol)

julia> sol[end]
2-element Array{Float64,1}:
 0.000407308
 0.55

This should probably get an FAQ listing. Your example here is pretty good for that. Note that MATLAB's strategy here is to always set dtmax = (tspan[2] - tspan[1]) / 100, i.e. always require at least 100 steps. That's adhoc and tends to catch it sometimes, but not all the time. I think that's too imposing on efficiency in standard use cases and still not reliable enough though.

[Stivanification]

Very nice! For my particular application the dtmax solution is perfect, because my "step function" is data measured at regular intervals, but I agree with you that the MATLAB strategy is rather inelegant in general. I did look for something like that in the docs, but I didn't expect it to be in the section "Output Control" although that makes sense in a way. After messing around with values like dtmax = 0.001 I have to say that I'm seriously impressed by the performance of this package. Thanks for your help and keep up the good work!