Feagin's Order 10, 12, and 14 Methods

Chris Rackauckas

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.

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) = (linear_bigα*u)

# Add analytical solution so that errors are checked
f_analytic(u0,p,t) = u0*exp(linear_bigα*t)
ff = ODEFunction(f,analytic=f_analytic)
prob = ODEProblem(ff,big(0.5),(0.0,1.0))
sol = solve(prob,Feagin14(),dt=1//16,adaptive=false);
println(sol.errors)
Dict(:l∞=>2.19751e-23,:final=>2.19751e-23,:l2=>1.0615e-23)

Compare that to machine $\epsilon$ for Float64:

eps(Float64)
2.220446049250313e-16

The error for Feagin's method when the stepsize is 1/16 is 8 orders of magnitude below machine $\epsilon$! 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); print("The length was $(length(sol))")
Dict(:l∞=>1.54574e-09,:final=>1.54574e-09,:l2=>8.92507e-10)
The length was 3

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 DiffEqDevTools
dts = 1.0 ./ 2.0 .^(10:-1:4)
sim = test_convergence(dts,prob,Feagin14())
DiffEqDevTools.ConvergenceSimulation{DiffEqBase.ODESolution{BigFloat,1,Arra
y{BigFloat,1},Array{BigFloat,1},Dict{Symbol,BigFloat},Array{Float64,1},Arra
y{Array{BigFloat,1},1},DiffEqBase.ODEProblem{BigFloat,Tuple{Float64,Float64
},false,Nothing,DiffEqBase.ODEFunction{false,typeof(Main.WeaveSandBox22.f),
LinearAlgebra.UniformScaling{Bool},typeof(Main.WeaveSandBox22.f_analytic),N
othing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Dif
fEqBase.StandardODEProblem},OrdinaryDiffEq.Feagin14,OrdinaryDiffEq.Interpol
ationData{DiffEqBase.ODEFunction{false,typeof(Main.WeaveSandBox22.f),Linear
Algebra.UniformScaling{Bool},typeof(Main.WeaveSandBox22.f_analytic),Nothing
,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{BigFloat,1}
,Array{Float64,1},Array{Array{BigFloat,1},1},OrdinaryDiffEq.Feagin14Constan
tCache{BigFloat,Float64}},DiffEqBase.DEStats}}(DiffEqBase.ODESolution{BigFl
oat,1,Array{BigFloat,1},Array{BigFloat,1},Dict{Symbol,BigFloat},Array{Float
64,1},Array{Array{BigFloat,1},1},DiffEqBase.ODEProblem{BigFloat,Tuple{Float
64,Float64},false,Nothing,DiffEqBase.ODEFunction{false,typeof(Main.WeaveSan
dBox22.f),LinearAlgebra.UniformScaling{Bool},typeof(Main.WeaveSandBox22.f_a
nalytic),Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},N
othing,DiffEqBase.StandardODEProblem},OrdinaryDiffEq.Feagin14,OrdinaryDiffE
q.InterpolationData{DiffEqBase.ODEFunction{false,typeof(Main.WeaveSandBox22
.f),LinearAlgebra.UniformScaling{Bool},typeof(Main.WeaveSandBox22.f_analyti
c),Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{B
igFloat,1},Array{Float64,1},Array{Array{BigFloat,1},1},OrdinaryDiffEq.Feagi
n14ConstantCache{BigFloat,Float64}},DiffEqBase.DEStats}[retcode: Success
Interpolation: 3rd order Hermite
t: [0.0, 0.000976563, 0.00195313, 0.00292969, 0.00390625, 0.00488281, 0.005
85938, 0.00683594, 0.0078125, 0.00878906  …  0.991211, 0.992188, 0.993164, 
0.994141, 0.995117, 0.996094, 0.99707, 0.998047, 0.999023, 1.0]
u: BigFloat[0.50, 0.500493, 0.500987, 0.501482, 0.501977, 0.502472, 0.50296
8, 0.503464, 0.503961, 0.504458  …  1.36067, 1.36201, 1.36335, 1.3647, 1.36
605, 1.3674, 1.36874, 1.3701, 1.37145, 1.3728], retcode: Success
Interpolation: 3rd order Hermite
t: [0.0, 0.00195313, 0.00390625, 0.00585938, 0.0078125, 0.00976563, 0.01171
88, 0.0136719, 0.015625, 0.0175781  …  0.982422, 0.984375, 0.986328, 0.9882
81, 0.990234, 0.992188, 0.994141, 0.996094, 0.998047, 1.0]
u: BigFloat[0.50, 0.500987, 0.501977, 0.502968, 0.503961, 0.504956, 0.50595
3, 0.506952, 0.507953, 0.508956  …  1.34864, 1.35131, 1.35397, 1.35665, 1.3
5933, 1.36201, 1.3647, 1.3674, 1.3701, 1.3728], retcode: Success
Interpolation: 3rd order Hermite
t: [0.0, 0.00390625, 0.0078125, 0.0117188, 0.015625, 0.0195313, 0.0234375, 
0.0273438, 0.03125, 0.0351563  …  0.964844, 0.96875, 0.972656, 0.976563, 0.
980469, 0.984375, 0.988281, 0.992188, 0.996094, 1.0]
u: BigFloat[0.50, 0.501977, 0.503961, 0.505953, 0.507953, 0.509961, 0.51197
7, 0.514001, 0.516033, 0.518073  …  1.32491, 1.33015, 1.33541, 1.34069, 1.3
4599, 1.35131, 1.35665, 1.36201, 1.3674, 1.3728], retcode: Success
Interpolation: 3rd order Hermite
t: [0.0, 0.0078125, 0.015625, 0.0234375, 0.03125, 0.0390625, 0.046875, 0.05
46875, 0.0625, 0.0703125  …  0.929688, 0.9375, 0.945313, 0.953125, 0.960938
, 0.96875, 0.976563, 0.984375, 0.992188, 1.0]
u: BigFloat[0.50, 0.503961, 0.507953, 0.511977, 0.516033, 0.520121, 0.52424
1, 0.528394, 0.53258, 0.536799  …  1.27869, 1.28882, 1.29903, 1.30932, 1.31
969, 1.33015, 1.34069, 1.35131, 1.36201, 1.3728], retcode: Success
Interpolation: 3rd order Hermite
t: [0.0, 0.015625, 0.03125, 0.046875, 0.0625, 0.078125, 0.09375, 0.109375, 
0.125, 0.140625  …  0.859375, 0.875, 0.890625, 0.90625, 0.921875, 0.9375, 0
.953125, 0.96875, 0.984375, 1.0]
u: BigFloat[0.50, 0.507953, 0.516033, 0.524241, 0.53258, 0.541051, 0.549658
, 0.558401, 0.567283, 0.576306  …  1.19103, 1.20998, 1.22923, 1.24878, 1.26
864, 1.28882, 1.30932, 1.33015, 1.35131, 1.3728], retcode: Success
Interpolation: 3rd order Hermite
t: [0.0, 0.03125, 0.0625, 0.09375, 0.125, 0.15625, 0.1875, 0.21875, 0.25, 0
.28125  …  0.71875, 0.75, 0.78125, 0.8125, 0.84375, 0.875, 0.90625, 0.9375,
 0.96875, 1.0]
u: BigFloat[0.50, 0.516033, 0.53258, 0.549658, 0.567283, 0.585473, 0.604247
, 0.623623, 0.64362, 0.664258  …  1.03333, 1.06647, 1.10067, 1.13596, 1.172
39, 1.20998, 1.24878, 1.28882, 1.33015, 1.3728], retcode: Success
Interpolation: 3rd order Hermite
t: [0.0, 0.0625, 0.125, 0.1875, 0.25, 0.3125, 0.375, 0.4375, 0.5, 0.5625, 0
.625, 0.6875, 0.75, 0.8125, 0.875, 0.9375, 1.0]
u: BigFloat[0.50, 0.53258, 0.567283, 0.604247, 0.64362, 0.685558, 0.730229,
 0.777811, 0.828493, 0.882477, 0.93998, 1.00123, 1.06647, 1.13596, 1.20998,
 1.28882, 1.3728]], Dict{Any,Any}(:l∞=>BigFloat[3.35435e-49, 5.07978e-45, 6
.96505e-41, 6.99856e-37, 2.7616e-33, 4.96506e-28, 2.19751e-23],:final=>BigF
loat[3.35435e-49, 5.07978e-45, 6.96505e-41, 6.99856e-37, 2.7616e-33, 4.9650
6e-28, 2.19751e-23],:l2=>BigFloat[1.55766e-49, 2.36041e-45, 3.24061e-41, 3.
26457e-37, 1.29478e-33, 2.35149e-28, 1.0615e-23]), 7, Dict(:dts=>[0.0009765
63, 0.00195313, 0.00390625, 0.0078125, 0.015625, 0.03125, 0.0625]), Dict{An
y,Any}(:l∞=>14.2933,:final=>14.2933,:l2=>14.3028), [0.000976563, 0.00195313
, 0.00390625, 0.0078125, 0.015625, 0.03125, 0.0625])

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

using Plots
gr()
plot(sim)

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

Appendix

This tutorial is part of the DiffEqTutorials.jl repository, found at: https://github.com/JuliaDiffEq/DiffEqTutorials.jl

To locally run this tutorial, do the following commands:

using DiffEqTutorials
DiffEqTutorials.weave_file("ode_extras","02-feagin.jmd")

Computer Information:

Julia Version 1.1.1
Commit 55e36cc308 (2019-05-16 04:10 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-3770 CPU @ 3.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, ivybridge)

Package Information:

Status `~/.julia/environments/v1.1/Project.toml`
[7e558dbc-694d-5a72-987c-6f4ebed21442] ArbNumerics 0.5.4
[6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf] BenchmarkTools 0.4.2
[be33ccc6-a3ff-5ff2-a52e-74243cff1e17] CUDAnative 2.2.0
[3a865a2d-5b23-5a0f-bc46-62713ec82fae] CuArrays 1.0.2
[55939f99-70c6-5e9b-8bb0-5071ed7d61fd] DecFP 0.4.8
[abce61dc-4473-55a0-ba07-351d65e31d42] Decimals 0.4.0
[ebbdde9d-f333-5424-9be2-dbf1e9acfb5e] DiffEqBayes 1.1.0
[eb300fae-53e8-50a0-950c-e21f52c2b7e0] DiffEqBiological 3.8.2
[459566f4-90b8-5000-8ac3-15dfb0a30def] DiffEqCallbacks 2.5.2
[f3b72e0c-5b89-59e1-b016-84e28bfd966d] DiffEqDevTools 2.9.0
[1130ab10-4a5a-5621-a13d-e4788d82bd4c] DiffEqParamEstim 1.6.0
[055956cb-9e8b-5191-98cc-73ae4a59e68a] DiffEqPhysics 3.1.0
[6d1b261a-3be8-11e9-3f2f-0b112a9a8436] DiffEqTutorials 0.1.0
[0c46a032-eb83-5123-abaf-570d42b7fbaa] DifferentialEquations 6.4.0
[31c24e10-a181-5473-b8eb-7969acd0382f] Distributions 0.20.0
[497a8b3b-efae-58df-a0af-a86822472b78] DoubleFloats 0.9.1
[f6369f11-7733-5829-9624-2563aa707210] ForwardDiff 0.10.3
[c91e804a-d5a3-530f-b6f0-dfbca275c004] Gadfly 1.0.1
[7073ff75-c697-5162-941a-fcdaad2a7d2a] IJulia 1.18.1
[4138dd39-2aa7-5051-a626-17a0bb65d9c8] JLD 0.9.1
[23fbe1c1-3f47-55db-b15f-69d7ec21a316] Latexify 0.8.2
[eff96d63-e80a-5855-80a2-b1b0885c5ab7] Measurements 2.0.0
[961ee093-0014-501f-94e3-6117800e7a78] ModelingToolkit 0.2.0
[76087f3c-5699-56af-9a33-bf431cd00edd] NLopt 0.5.1
[2774e3e8-f4cf-5e23-947b-6d7e65073b56] NLsolve 4.0.0
[429524aa-4258-5aef-a3af-852621145aeb] Optim 0.18.1
[1dea7af3-3e70-54e6-95c3-0bf5283fa5ed] OrdinaryDiffEq 5.8.1
[65888b18-ceab-5e60-b2b9-181511a3b968] ParameterizedFunctions 4.1.1
[91a5bcdd-55d7-5caf-9e0b-520d859cae80] Plots 0.25.1
[d330b81b-6aea-500a-939a-2ce795aea3ee] PyPlot 2.8.1
[731186ca-8d62-57ce-b412-fbd966d074cd] RecursiveArrayTools 0.20.0
[90137ffa-7385-5640-81b9-e52037218182] StaticArrays 0.11.0
[f3b207a7-027a-5e70-b257-86293d7955fd] StatsPlots 0.11.0
[c3572dad-4567-51f8-b174-8c6c989267f4] Sundials 3.6.1
[1986cc42-f94f-5a68-af5c-568840ba703d] Unitful 0.15.0
[44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9] Weave 0.9.0
[b77e0a4c-d291-57a0-90e8-8db25a27a240] InteractiveUtils
[37e2e46d-f89d-539d-b4ee-838fcccc9c8e] LinearAlgebra
[44cfe95a-1eb2-52ea-b672-e2afdf69b78f] Pkg