For these tests we will solve an 100x100 system of linear differential equations. This will demonstrate the efficiency of the methods for handling large systems. 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, and one from ODEInterface: Hairer's famous dopri5).

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.

In [1]:
using DifferentialEquations, DiffEqDevTools, Plots, ODEInterfaceDiffEq, ODE
gr()
# 2D Linear ODE
f = (t,u,du) -> begin
  for i in 1:length(u)
    du[i] = 1.01*u[i]
  end
end
function (p::typeof(f))(::Type{Val{:analytic}},t,u₀)
  u₀*exp(1.01*t)
end
tspan = (0.0,10.0)
prob = ODEProblem(f,rand(100,100),tspan)

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

Speed Baseline

First a baseline. These settings make DifferentialEquations match the Hairer algorithm. Extra saving features are turned off to put DifferentialEquations.jl in "speed mode".

In [2]:
setups = [Dict(:alg=>DP5())
          Dict(:alg=>ode45())
          Dict(:alg=>dopri5())
          Dict(:alg=>Tsit5())]
names = ["DifferentialEquations";"ODE";"ODEInterface";"DifferentialEquations Tsit5"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=names,save_everystep=false)
plot(wp)
Out[2]:
10 - 6 10 - 4 10 - 2 10 0 10 2 10 - 2 10 - 1 Error Time (s) DifferentialEquations ODE ODEInterface DifferentialEquations Tsit5

DifferentialEquations.jl has is clearly far in the lead, being more than an order of magnitude faster for the same amount of error.

Full Saving

Now if we include full saving at every timestep, we receive the following:

In [3]:
setups = [Dict(:alg=>DP5(),:dense=>false)
          Dict(:alg=>ode45(),:dense=>false)
          Dict(:alg=>dopri5()) # dense=false by default: no nonlinear interpolation
          Dict(:alg=>Tsit5(),:dense=>false)]
names = ["DifferentialEquations";"ODE";"ODEInterface";"DifferentialEquations Tsit5"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=names)
plot(wp)
Out[3]:
10 - 6 10 - 4 10 - 2 10 0 10 2 10 - 2 10 - 1 Error Time (s) DifferentialEquations ODE ODEInterface DifferentialEquations Tsit5

While not as dramatic as before, DifferentialEquations.jl is still far in the lead. Since the times are log scaled, this comes out to be almost a 10x lead over ODEInterface, and about a 50x lead over ODE.jl at default tolerances.

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. This feature is also only in DifferentialEquations.jl, but leads to the nice sol(t) output.

In [4]:
setups = [Dict(:alg=>DP5())
          Dict(:alg=>ode45())
          Dict(:alg=>dopri5())
          Dict(:alg=>Tsit5())]
names = ["DifferentialEquations";"ODE";"ODEInterface";"DifferentialEquations Tsit5"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=names)
plot(wp)
Out[4]:
10 - 6 10 - 4 10 - 2 10 0 10 2 10 - 2 10 - 1 Error Time (s) DifferentialEquations ODE ODEInterface DifferentialEquations Tsit5

As you can see, even with this large overhead, DifferentialEquations.jl essentially ties with ODEInterface. This shows that the fully featured DP5 solver holds its own with even the classic "great" methods.

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:

In [5]:
setups = [Dict(:alg=>DP5())
          Dict(:alg=>BS3())
          Dict(:alg=>BS5())
          Dict(:alg=>Tsit5())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;save_everystep=false)
plot(wp)
Out[5]:
10 - 6 10 - 4 10 - 2 10 0 10 2 10 - 2 10 - 1 Error Time (s) OrdinaryDiffEq.DP5 OrdinaryDiffEq.BS3 OrdinaryDiffEq.BS5 OrdinaryDiffEq.Tsit5

As you can see, the Tsit5 algorithm is the most efficient, beating DP5 which is more efficient than the Bogacki-Shampine algorithms. However, you can see that when the tolerance is high, BS3 could be of use since its slope is so steep.

Higher Order

Now let's see how it fairs against some higher order algorithms:

In [6]:
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)
plot(wp)
Out[6]:
10 - 5 10 0 10 - 2.5 10 - 2.0 10 - 1.5 Error Time (s) OrdinaryDiffEq.DP5 OrdinaryDiffEq.Vern6 OrdinaryDiffEq.TanYam7 OrdinaryDiffEq.Vern7 OrdinaryDiffEq.Vern8 OrdinaryDiffEq.DP8 OrdinaryDiffEq.Vern9

Vern7 looks to be the winner here, with DP5 doing well at higher tolerances but trailing of when it gets lower as one would expect with lower order algorithms. Some of the higher order methods, such as Vern9, would do better at lower tolerances than what's tested (outside of floating point range).

Higher Order Against Other Packages

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

In [7]:
setups = [Dict(:alg=>DP5())       
          Dict(:alg=>Vern7())
          Dict(:alg=>dop853())
          Dict(:alg=>ode78())
          Dict(:alg=>CVODE_Adams())]    
wp = WorkPrecisionSet(prob,abstols,reltols,setups;save_everystep=false)
plot(wp)
Out[7]:
10 - 6 10 - 4 10 - 2 10 0 10 2 10 - 2 10 - 1 Error Time (s) OrdinaryDiffEq.DP5 OrdinaryDiffEq.Vern7 ODEInterfaceDiffEq.dop853 ODE.ode78 Sundials.CVODE Adams :Functional,:None

Here you can once again see the DifferentialEquations algorithms far in the lead. It's well known that for cheap function costs Adams methods are inefficient. ODE.jl one again has a bad showing with feh78 an order of magnitude slower than the DiffererentialEquations algorithm.

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:

In [13]:
setups = [Dict(:alg=>DP5())
          #Dict(:alg=>ode45())
          Dict(:alg=>Tsit5())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;error_estimate=:L2,dense_errors=true)
plot(wp)
Out[13]:
10 - 6 10 - 4 10 - 2 10 0 10 2 10 - 1.0 10 - 0.8 10 - 0.6 10 - 0.4 Error Time (s) OrdinaryDiffEq.DP5 OrdinaryDiffEq.Tsit5

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:

In [11]:
setups = [Dict(:alg=>DP5())       
          Dict(:alg=>Vern7())
          #Dict(:alg=>ode78())
          ]    
wp = WorkPrecisionSet(prob,abstols,reltols,setups;error_estimate=:L2,dense_errors=true)
plot(wp)
Out[11]:
10 - 6 10 - 3 10 0 10 - 1.0 10 - 0.8 10 - 0.6 10 - 0.4 Error Time (s) OrdinaryDiffEq.DP5 OrdinaryDiffEq.Vern7

Conclusion

DifferentialEquations's default choice of Tsit5 does well for quick and easy solving at normal tolerances. However, at low tolerances the higher order algorithms are faster. In every case, the DifferentialEquations algorithms are far in the lead, many times an order of magnitude faster than the competitors. Vern7 with its included 7th order interpolation looks to be a good workhorse for scientific computing in floating point range. ODE.jl is far far behind. These along with many other benchmarks are why these algorithms were chosen as part of the defaults.