# The Formula of Explicit Runge-Kutta Fourteen order

I need an explicit Runge-Kutta 14th order formula. If you know about some reference that discusses at least 10th order (or higher, since I'm looking for the 14th) of Runge-Kutta and there is information about the whole formula like

$$y=y+\frac{h}{C}\left(c_1 k_1+ c_2 k_2+c_3 k3+\cdots +c_n k_n\right)$$

with defined coefficients $$k_1,\ldots,k_n$$ and maybe a code written in C or C++.

## 1 Answer

The 14th order methods due to Feagin can be found in DifferentialEquations.jl. An example using them with 128-bit floating point arithmetic is as follows:

using OrdinaryDiffEq, DoubleFloats
function lorenz(du,u,p,t)
du[1] = 10.0(u[2]-u[1])
du[2] = u[1]*(28.0-u[3]) - u[2]
du[3] = u[1]*u[2] - (8/3)*u[3]
end
u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz,u0,tspan)
sol = solve(prob,Feagin14(),abstol=1e-12.reltol=1e-12)


and the tolerances can be reduced of course. The implementation is fairly straightforward and can be found here:

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/perform_step/feagin_rk_perform_step.jl#L608-L764

with the tableau defined here:

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/tableaus/feagin_tableaus.jl#L1546-L1980

Notice that you can use the built-in uncertainty quantification tools to see if this accuracy is good enough that the chaotic system didn't diverge from the true trajectory yet, so this is a nice application of a 14th order integrator. See the ProbInts documentation for details. However, it's generally not a good idea to use this, since it's only efficient at tolerances small than machine epsilon for 64-bit floating point numbers, so you need a really weird case like this to justify using it.