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.
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.
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.
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