1. DiffEqBenchmarks.jl
  2. NonStiffODE
In [1]:
using DifferentialEquations, ParameterizedFunctions, ODE, ODEInterfaceDiffEq, LSODA

f(t) = 0.25*sin(t)^2

g = @ode_def RigidBody begin
  dy1  = I₁*y2*y3
  dy2  = I₂*y1*y3
  dy3  = I₃*y1*y2 + f(t)
end I₁ I₂ I₃

p = [-2.0,1.25,-0.5]
prob = ODEProblem(g,[1.0;0.0;0.9],(0.0,10.0),p)

abstols = 1./10.^(6:13)
reltols = 1./10.^(3:10);
sol = solve(prob,Vern7(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)
using Plots; gr()
Out[1]:
Plots.GRBackend()
In [2]:
sol = solve(prob,Vern7())
plot(sol)
Out[2]:
0 2 4 6 8 10 0 1 2 t y1(t) y2(t) y3(t)
In [3]:
setups = [Dict(:alg=>DP5())
          #Dict(:alg=>ode45()) # fails
          Dict(:alg=>dopri5())
          Dict(:alg=>Tsit5())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=true,numruns=1000,maxiters=10000)
plot(wp)
Out[3]:
10 - 10 10 - 8 10 - 6 10 - 4 10 - 2 10 - 3.5 10 - 3.0 10 - 2.5 Error Time (s) OrdinaryDiffEq.DP5 ODEInterfaceDiffEq.dopri5 OrdinaryDiffEq.Tsit5 OrdinaryDiffEq.Vern6

The DifferentialEquations.jl algorithms once again pull ahead. This is the first benchmark we've ran where ode45 doesn't fail. However, it still doesn't do as well as Tsit5. One reason why it does so well is that the maximum norm that ODE.jl uses (as opposed to the L2 norm of Sundials, DifferentialEquations, and ODEInterface) seems to do really well on this problem. dopri5 does surprisingly bad in this test.

Higher Order

In [4]:
setups = [Dict(:alg=>DP8())
          #Dict(:alg=>ode78()) # fails
          Dict(:alg=>Vern7())
          Dict(:alg=>Vern8())
          Dict(:alg=>dop853())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=1000,maxiters=1000)
plot(wp)
Out[4]:
10 - 10 10 - 5 10 - 3.75 10 - 3.50 10 - 3.25 Error Time (s) OrdinaryDiffEq.DP8 OrdinaryDiffEq.Vern7 OrdinaryDiffEq.Vern8 ODEInterfaceDiffEq.dop853 OrdinaryDiffEq.Vern6
In [2]:
setups = [Dict(:alg=>Vern7())
          Dict(:alg=>Vern8())
          Dict(:alg=>odex())
          Dict(:alg=>CVODE_Adams())
          Dict(:alg=>lsoda())
          Dict(:alg=>ddeabm())
          Dict(:alg=>ARKODE(Sundials.Explicit(),order=6))
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=1000,maxiters=1000)
plot(wp)
Out[2]:
10 - 12 10 - 10 10 - 8 10 - 6 10 - 4 10 - 3.5 10 - 3.4 10 - 3.3 10 - 3.2 10 - 3.1 10 - 3.0 10 - 2.9 Error Time (s) OrdinaryDiffEq.Vern7 OrdinaryDiffEq.Vern8 ODEInterfaceDiffEq.odex Sundials.CVODE_Adams LSODA.lsoda ODEInterfaceDiffEq.ddeabm Sundials.ARKODE

Conclusion

Once again, the OrdinaryDiffEq.jl pull far ahead in terms of speed and accuracy.