Choosing a solver

ContinuousDynamicalSystems are evolved using solvers from DifferentialEquations.jl. In this page we discuss the importance of which solver to choose.

Default Solver

The default solver is:

using DynamicalSystems
DynamicalSystemsBase.DEFAULT_SOLVER

SimpleDiffEq.SimpleATsit5()

which is a Runge-Kutta-like solver. The number in the solver's name is the "order" of the solver.

Speed of a solver

Estimating a given solver's performance for a particular problem is not trivial. The following are general rules of thumb:

1. Higher order solvers call the equations of motion function more times per step.
2. Higher order solvers can cover larger timespans per step.
3. Higher order solvers do better at small tolerances.

This means that there is a delicate balance between how expensive is your function and how large of a step a solver can take while it is still efficient. In general you want to strike a point of taking large steps but also not calling the function exceedingly often.

How do I pick?

The answer to this question is easy: benchmarks!

Here is a simple case: let's compute the Lyapunov spectrum of the Lorenz system using lyapunovs:

ds = Systems.lorenz()
tols = (abstol = 1e-6, reltol = 1e-6)
lyapunovs(ds, 2000; Ttr = 100.0, tols...)

3-element Array{Float64,1}:
   0.9061277794060383
  -0.0003445080325231465
 -14.572426217477933

The above uses the default solver. Let's now benchmark using two different solvers, SimpleATsit5 and Vern9. Since the SimpleATsit5 case is of lower order, naively one might think it is faster because it makes less function calls. This argument is not necessarily true though.

It is important to understand that when calling lyapunovs(ds, 2000) you want the system (and the tangent space) to be evolved so that it reaches a total time of 2000*dt, which by default is 2000.0 units of time. Even though SimpleATsit5 requires less function calls per step, Vern9 can cover larger timespans per step.

Here are the numbers:

using BenchmarkTools, OrdinaryDiffEq, SimpleDiffEq, Statistics
b1 = @benchmark lyapunovs(ds, 2000; alg = SimpleATsit5(), Ttr = 100.0, tols...);
b2 = @benchmark lyapunovs(ds, 2000; alg = Vern9(),        Ttr = 100.0, tols...);
println("Timing for SimpleATsit5:")
println(mean(b1))
println("Timing for Vern9:")
println(mean(b2))

Timing for SimpleATsit5:
TrialEstimate(64.804 ms)
Timing for Vern9:
TrialEstimate(33.897 ms)

As you can see Vern9 is faster in doing the entire computation! Of course this does not have to be universally true. It is true for the Lorenz system, but for your specific system you should do dedicated benchmarks!

DifferentialEquations.jl

For more info about the possible solvers be sure to head over to the documentation of DifferentialEquations.jl!