DP5() algorithm failing to solve simple SIR problem?

Specific Domains Numerics
,
1

While testing out different ODE solvers, I found that DP5() completely fails when given a simple SIR model (while other solvers handle it perfectly fine):

function SIR!(du,u,pars,t)    
    β, γ = pars
    S, I, R = u
 
    dS = -β*S*I
    dI = β*S*I - γ*I
    dR = γ*I
    du .= [dS, dI, dR]
end

t_obs = collect(0:1.0:218)
prob = ODEProblem(SIR!, [0.99, 0.01, 0.0], (t_obs[1], t_obs[end]), [0.20, 0.15])
sol = solve(prob, DP5(), reltol = 1e-6, abstol = 1e-6, saveat = t_obs)
y = DataFrame(sol',[:S,:I,:R])
plot(t_obs, y."I")

This produces the following nonsensical plot:
DP5_fail

I tried it with extreme low tolerances as well. Here’s what it looks like when rtol = 1e-22 and abstol = 1e-22:
DP5_fail2

Much better, but still not at all good…

Tsit5() handles it fine though:
DP5_Tsit

Any idea what’s going on here? I thought DP5 was quite robust.

1 Like
2

This is more of an issue than a Discourse post. There was a default that was somehow being set incorrectly which caused this. It’s a rather specific case so it wasn’t caught in the other tests. It’s fixed in Fix calck default and add a test for DP5 only saveat interpolation by ChrisRackauckas · Pull Request #1450 · SciML/OrdinaryDiffEq.jl · GitHub and should get merged and tagged by tomorrow, and this has been turned into a test. Thanks for sharing.

8 Likes