Feagin's Order 10, 12, and 14 Methods

DifferentialEquations.jl includes Feagin's explicit Runge-Kutta methods of orders 10/8, 12/10, and 14/12. These methods have such high order that it's pretty much required that one uses numbers with more precision than Float64. As a prerequisite reference on how to use arbitrary number systems (including higher precision) in the numerical solvers, please see the Solving Equations in With Chosen Number Types notebook.

Ju-1.2-DE-Static

Download as Docker image from:

Copy

This image was imported from: docker.nextjournal.com/environment@sha256:32897c89c624757eff78334db434d94a6c945a6406e1f6ff946e71d493183221

Investigation of the Method's Error

We can use Feagin's order 16 method as follows. Let's use a two-dimensional linear ODE. Like in the Solving Equations in With Chosen Number Types notebook, we change the initial condition to BigFloats to tell the solver to use BigFloat types.

using DifferentialEquations
const linear_bigα = big(1.01)
f(u,p,t) = p * u

# Add analytical solution so that errors are checked
f_analytic(u0,p,t) = u0 * exp(p * t)
ff = ODEFunction(f,analytic=f_analytic)
prob = ODEProblem(ff,big(0.5),(0.0,1.0), linear_bigα)
sol = solve(prob,Feagin14(),dt=1//16,adaptive=false);
sol.errors

9.3s

Julia

Dict{Symbol,BigFloat} with 3 entries: :l∞ => 2.19751e-23 :final => 2.19751e-23 :l2 => 1.0615e-23

Compare that to machine for Float64:

eps(Float64)

0.6s

Julia

2.22045e-16

The error for Feagin's method when the stepsize is 1/16 is 8 orders of magnitude below machine ! However, that is dependent on the stepsize. If we instead use adaptive timestepping with the default tolerances, we get

sol = solve(prob,Feagin14())
println(sol.errors)
println("The length was $(length(sol))")

1.9s

Julia

Notice that when the stepsize is much higher, the error goes up quickly as well. These super high order methods are best when used to gain really accurate approximations (using still modest timesteps). Some examples of where such precision is necessary is astrodynamics where the many-body problem is highly chaotic and thus sensitive to small errors.

Convergence Test

The Order 14 method is awesome, but we need to make sure it's really that awesome. The following convergence test is used in the package tests in order to make sure the implementation is correct. Note that all methods have such tests in place.

using Pkg
pkg"add DiffEqDevTools"

22.9s

Julia

using DiffEqDevTools
dts = 1.0 ./ 2.0 .^(10:-1:4)
sim = test_convergence(dts,prob,Feagin14())

11.0s

Julia

For a view of what's going on, let's plot the simulation results.

using Plots
gr()
plot(sim)

36.3s

Julia

This is a clear trend indicating that the convergence is truly Order 14, which is the estimated slope.