100 Independent Linear Work-Precision Diagrams

Chris Rackauckas

For these tests we will solve a diagonal 100 independent linear differential equations. This will demonstrate the efficiency of the implementation of the methods for handling large systems, since the system is both large enough that array handling matters, but f is cheap enough that it is not simply a game of calculating f as few times as possible. We will be mostly looking at the efficiency of the work-horse Dormand-Prince Order 4/5 Pairs: one from DifferentialEquations.jl (DP5), one from ODE.jl rk45, one from ODEInterface (Hairer's famous dopri5, and one from SUNDIALS' ARKODE suite.

Also included is Tsit5. While all other ODE programs have gone with the traditional choice of using the Dormand-Prince 4/5 pair as the default, DifferentialEquations.jl uses Tsit5 as one of the default algorithms. It's a very new (2011) and not widely known, but the theory and the implimentation shows it's more efficient than DP5. Thus we include it just to show off how re-designing a library from the ground up in a language for rapid code and rapid development has its advantages.

Setup

using OrdinaryDiffEq, Sundials, DiffEqDevTools, Plots, ODEInterfaceDiffEq, ODE, LSODA
using Random
Random.seed!(123)
gr()
# 2D Linear ODE
function f(du,u,p,t)
@inbounds for i in eachindex(u)
du[i] = 1.01*u[i]
end
end
function f_analytic(u₀,p,t)
u₀*exp(1.01*t)
end
tspan = (0.0,10.0)
prob = ODEProblem(ODEFunction(f,analytic=f_analytic),rand(100,100),tspan)

abstols = 1.0 ./ 10.0 .^ (3:13)
reltols = 1.0 ./ 10.0 .^ (0:10);


Speed Baseline

First a baseline. These are all testing the same Dormand-Prince order 5/4 algorithm of each package. While all the same Runge-Kutta tableau, they exhibit different behavior due to different choices of adaptive timestepping algorithms and tuning. First we will test with all extra saving features are turned off to put DifferentialEquations.jl in "speed mode".

setups = [Dict(:alg=>DP5())
Dict(:alg=>ode45())
Dict(:alg=>dopri5())
Dict(:alg=>ARKODE(Sundials.Explicit(),etable=Sundials.DORMAND_PRINCE_7_4_5))
Dict(:alg=>Tsit5())]
solnames = ["OrdinaryDiffEq";"ODE";"ODEInterface";"Sundials ARKODE";"OrdinaryDiffEq Tsit5"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=solnames,save_everystep=false,numruns=100)
plot(wp)


Full Saving

setups = [Dict(:alg=>DP5(),:dense=>false)
Dict(:alg=>ode45(),:dense=>false)
Dict(:alg=>dopri5()) # dense=false by default: no nonlinear interpolation
Dict(:alg=>ARKODE(Sundials.Explicit(),etable=Sundials.DORMAND_PRINCE_7_4_5),:dense=>false)
Dict(:alg=>Tsit5(),:dense=>false)]
solnames = ["OrdinaryDiffEq";"ODE";"ODEInterface";"Sundials ARKODE";"OrdinaryDiffEq Tsit5"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=solnames,numruns=100)
plot(wp)


Continuous Output

Now we include continuous output. This has a large overhead because at every timepoint the matrix of rates k has to be deep copied.

setups = [Dict(:alg=>DP5())
Dict(:alg=>ode45())
Dict(:alg=>dopri5())
Dict(:alg=>ARKODE(Sundials.Explicit(),etable=Sundials.DORMAND_PRINCE_7_4_5))
Dict(:alg=>Tsit5())]
solnames = ["OrdinaryDiffEq";"ODE";"ODEInterface";"Sundials ARKODE";"OrdinaryDiffEq Tsit5"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=solnames,numruns=100)
plot(wp)


Other Runge-Kutta Algorithms

Now let's test it against a smattering of other Runge-Kutta algorithms. First we will test it with all overheads off. Let's do the Order 5 (and the 2/3 pair) algorithms:

setups = [Dict(:alg=>DP5())
Dict(:alg=>BS3())
Dict(:alg=>BS5())
Dict(:alg=>Tsit5())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;save_everystep=false,numruns=100)
plot(wp)


Higher Order

Now let's see how OrdinaryDiffEq.jl fairs with some higher order algorithms:

setups = [Dict(:alg=>DP5())
Dict(:alg=>Vern6())
Dict(:alg=>TanYam7())
Dict(:alg=>Vern7())
Dict(:alg=>Vern8())
Dict(:alg=>DP8())
Dict(:alg=>Vern9())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;save_everystep=false,numruns=100)
plot(wp)


Higher Order With Many Packages

Now we test OrdinaryDiffEq against the high order methods of the other packages:

setups = [Dict(:alg=>DP5())
Dict(:alg=>Vern7())
Dict(:alg=>dop853())
Dict(:alg=>ode78())
Dict(:alg=>odex())
Dict(:alg=>lsoda())
Dict(:alg=>ddeabm())
Dict(:alg=>ARKODE(Sundials.Explicit(),order=8))
wp = WorkPrecisionSet(prob,abstols,reltols,setups;save_everystep=false,numruns=100)
plot(wp)


Interpolation Error

Now we will look at the error using an interpolation measurement instead of at the timestepping points. Since the DifferentialEquations.jl algorithms have higher order interpolants than the ODE.jl algorithms, one would expect this would magnify the difference. First the order 4/5 comparison:

setups = [Dict(:alg=>DP5())
#Dict(:alg=>ode45())
Dict(:alg=>Tsit5())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;error_estimate=:L2,dense_errors=true,numruns=100)
plot(wp)


Note that all of ODE.jl uses a 3rd order Hermite interpolation, while the DifferentialEquations algorithms interpolations which are specialized to the algorithm. For example, DP5 and Tsit5 both use "free" order 4 interpolations, which are both as fast as the Hermite interpolation while achieving far less error. At higher order:

setups = [Dict(:alg=>DP5())
Dict(:alg=>Vern7())
#Dict(:alg=>ode78())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;error_estimate=:L2,dense_errors=true,numruns=100)
plot(wp)


Comparison with Fixed Timestep RK4

Let's run the first benchmark but add some fixed timestep RK4 methods to see the difference:

abstols = 1.0 ./ 10.0 .^ (3:13)
reltols = 1.0 ./ 10.0 .^ (0:10);
dts = [1,1/2,1/4,1/10,1/20,1/40,1/60,1/80,1/100,1/140,1/240]
setups = [Dict(:alg=>DP5())
Dict(:alg=>ode45())
Dict(:alg=>dopri5())
Dict(:alg=>RK4(),:dts=>dts)
Dict(:alg=>Tsit5())]
solnames = ["DifferentialEquations";"ODE";"ODEInterface";"DifferentialEquations RK4";"DifferentialEquations Tsit5"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=solnames,
save_everystep=false,verbose=false,numruns=100)
plot(wp)


Comparison with Non-RK methods

Now let's test Tsit5 and Vern9 against parallel extrapolation methods and an Adams-Bashforth-Moulton:

setups = [Dict(:alg=>Tsit5())
Dict(:alg=>Vern9())
Dict(:alg=>VCABM())
solnames = ["Tsit5","Vern9","VCABM","AitkenNeville","Midpoint Deuflhard","Midpoint Hairer Wanner"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=solnames,
save_everystep=false,verbose=false,numruns=100)
plot(wp)

setups = [Dict(:alg=>ExtrapolationMidpointDeuflhard(min_order=1, max_order=9, init_order=9, threading=false))

setups = [Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=10, threading=true))
Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=15, init_order