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
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=>lsoda())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)

Out[2]:
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]:
In [4]:
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2)
plot(wp)

Out[4]:
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
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)

Out[8]:
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]:
In [10]:
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2)
plot(wp)

Out[10]:
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]:
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]:

### 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=>lsoda())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)

Out[9]:
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]:
In [11]:
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2)
plot(wp)

Out[11]:
In [14]:
setups = [
Dict(:alg=>Rodas5()),
Dict(:alg=>Kvaerno5()),
Dict(:alg=>CVODE_BDF()),
Dict(:alg=>KenCarp4()),
Dict(:alg=>Rodas4()),
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)

Out[14]:
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]:
In [16]:
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2)
plot(wp)

Out[16]:
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]:
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]:

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.