In [1]:
using DifferentialEquations, DiffEqProblemLibrary, ParameterizedFunctions, Plots, ODE, ODEInterfaceDiffEq, LSODA
gr() #gr(fmt=:png)

f = @ode_def Hires begin
  dy1 = -1.71*y1 + 0.43*y2 + 8.32*y3 + 0.0007
  dy2 = 1.71*y1 - 8.75*y2
  dy3 = -10.03*y3 + 0.43*y4 + 0.035*y5
  dy4 = 8.32*y2 + 1.71*y3 - 1.12*y4
  dy5 = -1.745*y5 + 0.43*y6 + 0.43*y7
  dy6 = -280.0*y6*y8 + 0.69*y4 + 1.71*y5 -
           0.43*y6 + 0.69*y7
  dy7 = 280.0*y6*y8 - 1.81*y7
  dy8 = -280.0*y6*y8 + 1.81*y7
end

u0 = zeros(8)
u0[1] = 1
u0[8] = 0.0057
prob = ODEProblem(f,u0,(0.0,321.8122))

sol = solve(prob,Rodas4(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)
abstols = 1./10.^(4:11)
reltols = 1./10.^(1:8);
In [2]:
plot(sol)
Out[2]:
In [3]:
plot(sol,tspan=(0.0,5.0))
Out[3]:

Omissions

The following were omitted from the tests due to convergence failures. ODE.jl's adaptivity is not able to stabilize its algorithms, while GeometricIntegrators.jl's methods either fail to converge at comparable dts (or on some computers errors)

In [4]:
sol = solve(prob,ode23s()); println("Total ODE.jl steps: $(length(sol))")
using GeometricIntegratorsDiffEq
try
    sol = solve(prob,GIRadIIA3(),dt=1/10)
catch e
    println(e)
end
Total ODE.jl steps: 1
MethodError(convert, (GeometricIntegrators.Solvers

High Tolerances

This is the speed when you just want the answer.

In [2]:
abstols = 1./10.^(5:8)
reltols = 1./10.^(1:4);
setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>Rodas3()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>rodas()),
          Dict(:alg=>radau()),
          Dict(:alg=>lsoda())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)
Out[2]:
10 - 8 10 - 7 10 - 6 10 - 5 10 - 4 10 - 3 10 - 3.4 10 - 3.3 10 - 3.2 10 - 3.1 10 - 3.0 10 - 2.9 10 - 2.8 Error Time (s) OrdinaryDiffEq.Rosenbrock23 OrdinaryDiffEq.Rodas3 OrdinaryDiffEq.TRBDF2 Sundials.CVODE_BDF ODEInterfaceDiffEq.rodas ODEInterfaceDiffEq.radau LSODA.lsoda
In [3]:
wp = WorkPrecisionSet(prob,abstols,reltols,setups;dense = false,verbose=false,
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
plot(wp)
Out[3]:
10 - 8 10 - 7 10 - 6 10 - 5 10 - 4 10 - 3 10 - 3.3 10 - 3.2 10 - 3.1 10 - 3.0 10 - 2.9 10 - 2.8 Error Time (s) OrdinaryDiffEq.Rosenbrock23 OrdinaryDiffEq.Rodas3 OrdinaryDiffEq.TRBDF2 Sundials.CVODE_BDF ODEInterfaceDiffEq.rodas ODEInterfaceDiffEq.radau LSODA.lsoda
In [4]:
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2)
plot(wp)
Out[4]:
10 - 6 10 - 5 10 - 4 10 - 3 10 - 3.3 10 - 3.2 10 - 3.1 10 - 3.0 10 - 2.9 10 - 2.8 Error Time (s) OrdinaryDiffEq.Rosenbrock23 OrdinaryDiffEq.Rodas3 OrdinaryDiffEq.TRBDF2 Sundials.CVODE_BDF ODEInterfaceDiffEq.rodas ODEInterfaceDiffEq.radau LSODA.lsoda
In [8]:
setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>Kvaerno3()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>KenCarp3()),
    # Dict(:alg=>SDIRK2()), # Removed because it's bad
          Dict(:alg=>radau())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)
Out[8]:
10^-8 10^-7 10^-6 10^-5 10^-4 10^-3 10^-3.3 10^-3.2 10^-3.1 10^-3.0 10^-2.9 10^-2.8 Error Time (s) OrdinaryDiffEq.Rosenbrock23 OrdinaryDiffEq.Kvaerno3 Sundials.CVODE_BDF OrdinaryDiffEq.KenCarp4 OrdinaryDiffEq.TRBDF2 OrdinaryDiffEq.KenCarp3 ODEInterfaceDiffEq.radau
In [9]:
wp = WorkPrecisionSet(prob,abstols,reltols,setups;dense = false,verbose=false,
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
plot(wp)
Out[9]:
10^-8 10^-7 10^-6 10^-5 10^-4 10^-3 10^-3.3 10^-3.2 10^-3.1 10^-3.0 10^-2.9 10^-2.8 Error Time (s) OrdinaryDiffEq.Rosenbrock23 OrdinaryDiffEq.Kvaerno3 Sundials.CVODE_BDF OrdinaryDiffEq.KenCarp4 OrdinaryDiffEq.TRBDF2 OrdinaryDiffEq.KenCarp3 ODEInterfaceDiffEq.radau
In [10]:
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2)
plot(wp)
Out[10]:
10^-5.5 10^-5.0 10^-4.5 10^-4.0 10^-3.5 10^-3.0 10^-3.3 10^-3.2 10^-3.1 10^-3.0 10^-2.9 10^-2.8 Error Time (s) OrdinaryDiffEq.Rosenbrock23 OrdinaryDiffEq.Kvaerno3 Sundials.CVODE_BDF OrdinaryDiffEq.KenCarp4 OrdinaryDiffEq.TRBDF2 OrdinaryDiffEq.KenCarp3 ODEInterfaceDiffEq.radau
In [13]:
setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>KenCarp5()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>KenCarp3()),
          Dict(:alg=>ARKODE(order=5)),
          Dict(:alg=>ARKODE()),
          Dict(:alg=>ARKODE(order=3))]
names = ["Rosenbrock23" "KenCarp5" "KenCarp4" "KenCarp3" "ARKODE5" "ARKODE4" "ARKODE3"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      names=names,save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)
Out[13]:
10 - 8 10 - 7 10 - 6 10 - 5 10 - 4 10 - 3 10 - 3.25 10 - 3.00 10 - 2.75 10 - 2.50 10 - 2.25 Error Time (s) Rosenbrock23 KenCarp5 KenCarp4 KenCarp3 ARKODE5 ARKODE4 ARKODE3
In [14]:
wp = WorkPrecisionSet(prob,abstols,reltols,setups;dense = false,verbose=false,
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
plot(wp)
Out[14]:
10 - 7 10 - 6 10 - 5 10 - 4 10 - 3 10 - 3.2 10 - 3.0 10 - 2.8 10 - 2.6 10 - 2.4 10 - 2.2 10 - 2.0 Error Time (s) OrdinaryDiffEq.Rosenbrock23 OrdinaryDiffEq.KenCarp5 OrdinaryDiffEq.KenCarp4 OrdinaryDiffEq.KenCarp3 Sundials.ARKODE Sundials.ARKODE Sundials.ARKODE

Low Tolerances

This is the speed at lower tolerances, measuring what's good when accuracy is needed.

In [9]:
abstols = 1./10.^(7:13)
reltols = 1./10.^(4:10)

setups = [Dict(:alg=>GRK4A()),
          Dict(:alg=>Rodas4P()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>ddebdf()),
          Dict(:alg=>Rodas5()),
          Dict(:alg=>rodas()),
          Dict(:alg=>radau()),
          Dict(:alg=>lsoda())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)
Out[9]:
10 - 14 10 - 12 10 - 10 10 - 8 10 - 6 10 - 3.0 10 - 2.5 10 - 2.0 10 - 1.5 Error Time (s) OrdinaryDiffEq.GRK4A OrdinaryDiffEq.Rodas4P Sundials.CVODE_BDF ODEInterfaceDiffEq.ddebdf OrdinaryDiffEq.Rodas5 ODEInterfaceDiffEq.rodas ODEInterfaceDiffEq.radau LSODA.lsoda
In [10]:
wp = WorkPrecisionSet(prob,abstols,reltols,setups;verbose=false,
                      dense=false,appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
plot(wp)
Out[10]:
10 - 14 10 - 12 10 - 10 10 - 8 10 - 6 10 - 3.0 10 - 2.5 10 - 2.0 10 - 1.5 Error Time (s) OrdinaryDiffEq.GRK4A OrdinaryDiffEq.Rodas4P Sundials.CVODE_BDF ODEInterfaceDiffEq.ddebdf OrdinaryDiffEq.Rodas5 ODEInterfaceDiffEq.rodas ODEInterfaceDiffEq.radau LSODA.lsoda
In [11]:
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2)
plot(wp)
Out[11]:
10 - 12 10 - 10 10 - 8 10 - 6 10 - 4 10 - 3.0 10 - 2.5 10 - 2.0 10 - 1.5 Error Time (s) OrdinaryDiffEq.GRK4A OrdinaryDiffEq.Rodas4P Sundials.CVODE_BDF ODEInterfaceDiffEq.ddebdf OrdinaryDiffEq.Rodas5 ODEInterfaceDiffEq.rodas ODEInterfaceDiffEq.radau LSODA.lsoda
In [14]:
setups = [
          Dict(:alg=>Rodas5()),
          Dict(:alg=>Kvaerno5()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>Rodas4()),
          Dict(:alg=>radau())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)
Out[14]:
10^-14 10^-12 10^-10 10^-8 10^-6 10^-3.00 10^-2.75 10^-2.50 10^-2.25 10^-2.00 10^-1.75 Error Time (s) OrdinaryDiffEq.Rodas5 OrdinaryDiffEq.Kvaerno5 Sundials.CVODE_BDF OrdinaryDiffEq.KenCarp4 OrdinaryDiffEq.Rodas4 ODEInterfaceDiffEq.radau
In [15]:
wp = WorkPrecisionSet(prob,abstols,reltols,setups;verbose=false,
                      dense=false,appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
plot(wp)
Out[15]:
10^-12 10^-10 10^-8 10^-6 10^-3.0 10^-2.7 10^-2.4 10^-2.1 10^-1.8 Error Time (s) OrdinaryDiffEq.Rodas5 OrdinaryDiffEq.Kvaerno5 Sundials.CVODE_BDF OrdinaryDiffEq.KenCarp4 OrdinaryDiffEq.Rodas4 ODEInterfaceDiffEq.radau
In [16]:
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2)
plot(wp)
Out[16]:
10^-12 10^-10 10^-8 10^-6 10^-4 10^-3.0 10^-2.7 10^-2.4 10^-2.1 10^-1.8 Error Time (s) OrdinaryDiffEq.Rodas5 OrdinaryDiffEq.Kvaerno5 Sundials.CVODE_BDF OrdinaryDiffEq.KenCarp4 OrdinaryDiffEq.Rodas4 ODEInterfaceDiffEq.radau
In [15]:
setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>KenCarp5()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>KenCarp3()),
          Dict(:alg=>ARKODE(order=5)),
          Dict(:alg=>ARKODE()),
          Dict(:alg=>ARKODE(order=3))]
names = ["Rosenbrock23" "KenCarp5" "KenCarp4" "KenCarp3" "ARKODE5" "ARKODE4" "ARKODE3"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      names=names,save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)
Out[15]:
10 - 14 10 - 12 10 - 10 10 - 8 10 - 6 10 - 3.0 10 - 2.5 10 - 2.0 10 - 1.5 10 - 1.0 Error Time (s) Rosenbrock23 KenCarp5 KenCarp4 KenCarp3 ARKODE5 ARKODE4 ARKODE3
In [16]:
wp = WorkPrecisionSet(prob,abstols,reltols,setups;verbose=false,
                      dense=false,appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
plot(wp)
Out[16]:
10 - 12 10 - 10 10 - 8 10 - 6 10 - 3.0 10 - 2.5 10 - 2.0 10 - 1.5 10 - 1.0 Error Time (s) OrdinaryDiffEq.Rosenbrock23 OrdinaryDiffEq.KenCarp5 OrdinaryDiffEq.KenCarp4 OrdinaryDiffEq.KenCarp3 Sundials.ARKODE Sundials.ARKODE Sundials.ARKODE

The following algorithms were removed since they failed.

In [17]:
#setups = [#Dict(:alg=>Hairer4()),
          #Dict(:alg=>Hairer42()),
          #Dict(:alg=>Rodas3()),
          #Dict(:alg=>Kvaerno4()),
#Dict(:alg=>KenCarp5()),
          #Dict(:alg=>Cash4())
#]
#wp = WorkPrecisionSet(prob,abstols,reltols,setups;
#                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
#plot(wp)

Conclusion

At high tolerances, Rosenbrock23 and lsoda hitsthe the error estimates and are fast. At lower tolerances and normal user tolerances, Rodas5 is extremely fast. When you get down to reltol=1e-10 radau begins to become as efficient as Rodas4, and it continues to do well below that.