# Why are higher-order Runge–Kutta methods not used more often?

I was just curious as to why high-order (i.e. greater than 4) Runge–Kutta methods are almost never discussed/employed (at least to my knowledge). I understand it requires greater computational time per step (e.g. RK14 with 12th-order embedded step), but are there any other downsides of using higher order Runge–Kutta methods (e.g. stability issues)? When applied to equations with highly oscillating solutions on extreme time scales, wouldn't such higher-order methods be typically preferred?

• I think this is a very subjective question. The biggest downside as you have already noted is the cost of computation. Generally we try to balance between accuracy and computational time. In PDE's when people talk about higher order they are generally thinking about 3rd or 4th order. And the time stepping is also kept at the same order. – Vikram Nov 17 '16 at 14:13
• In PDE, a high order accuracy scheme for temporal dependence makes no sense if spatial accuracy is worse. In fact, the accuracy of spatial dependence is mostly about 2nd or 3rd order, especially when working on unstructured meshes. People needs to control global error truncation with the least cost, hence, considers Runge-Kutta with high enough order of accuracy in particular cases. – tqviet Nov 17 '16 at 17:34
• @tqviet If using backward or central difference approximations of up to order 8 for the spatial derivatives, RK8 would be suitable, no? In general, is there any accuracy or stability issues with using such high order finite difference approximations of the spatial derivatives? – Mathews24 Nov 18 '16 at 0:30
• @Mathews24: I did not mention to stability, which depends strongly on the equation. When a highly accurate scheme is applied to spatial dependence, we adopt RK to temporal dependence with at least the same order of accuracy, but the stability condition may require a smaller value of $\Delta t$. – tqviet Nov 18 '16 at 12:40
• I want to know why Taylor methods + autodiff'ing the rhs isn't used more often. – user14717 Feb 19 at 10:28

There are thousands of papers and hundreds of codes out there using Runge-Kutta methods of fifth order or higher. Note that the most commonly used explicit integrator in MATLAB is ODE45, which advances the solution using a 5th-order Runge-Kutta method.

# Examples of widely-used high-order Runge-Kutta methods

The paper of Dormand & Prince giving a 5th-order method has over 1700 citations according to Google Scholar. Most of those are papers using their method to solve some problem. The Cash-Karp method paper has over 400 citations. Perhaps the most widely-used method of order higher than 5 is the 8th-order method of Prince-Dormand which has over 400 citations on Google Scholar. I could give many other examples; and keep in mind that many (if not most) of the people using these methods never cite the papers.

# High-order methods and rounding error

If your accuracy is limited by rounding errors then you should use a higher-order method. This is because higher-order methods require fewer steps (and fewer function evaluations, even though there are more evaluations per step), so they commit fewer rounding errors. You can easily verify this yourself with simple experiments; it is a good homework problem for a first course in numerical analysis.

Tenth-order methods are extremely useful in double-precision arithmetic. On the contrary, if all we had was Euler's method, then rounding error would be a major issue and we would need very high-precision floating point numbers for many problems where high-order solvers do just fine.

# High order methods can be just as stable

@RichardZhang has referenced the second Dahlquist barrier, but that applies only to multistep methods. The question posted here is about Runge-Kutta methods, and there are Runge-Kutta methods of every order that are not only $$A$$-stable, but also $$B$$-stable (a stability property useful for some nonlinear problems). To learn about these methods, see for instance the text of Hairer & Wanner.

# High-order methods in celestial mechanics

When applied to equations with highly oscillating solutions on extreme time scales, wouldn't such higher-order methods be typically preferred?

You're exactly right! A prime example of this is celestial mechanics. I'm not an expert in that area. But this paper, for instance, compares methods for celestial mechanics and doesn't even consider order lower than 5. It concludes that methods of order 11 or 12 are often the most efficient (with the Prince-Dormand method of order 8 also often very efficient).

• Ketchson: could you please provide some evidence or explanation about this statement: "high-order extrapolation and deferred correction methods are Runge-Kutta methods"? Especially the "deferred correction methods". Thanks. – tqviet Nov 18 '16 at 6:43
• @David Ketcheson Can you discuss how your answer would change if using validated (verified) computing techniques, such as outward rounded interval or radial arithmetic? How about if higher than double precision outward rounded interval or radial arithmetic were used? What will happen with wrapping and dependency as Runge-Kutta order is increased, and just for fun, let's say that the ODE is very stiff.? – Mark L. Stone Nov 18 '16 at 12:57
• @MarkL.Stone That's a completely different set of questions. If you want to ask them, please post them as separate questions. However, I'm not an expert on those things and won't be able to answer. – David Ketcheson Nov 18 '16 at 19:01
• @tqviet Take a look at this paper for an explanation. – David Ketcheson Nov 18 '16 at 19:03

As long as you're using standard double precision floating point arithmetic, very high order methods aren't needed to get a solution with high accuracy in a reasonable number of steps. In practice I find that the accuracy of the solution is normally limited to a relative error of 1.0e-16 by the double precision floating point representation rather than the number/length of the steps that are taken with RKF45.

If you switch to some higher than double precision floating point arithmetic scheme, then using a 10th order method might well be worth while.

• I think this answer is misleading. High-order methods lead to much less rounding error, whereas low-order methods suffer from rounding error being dominant when the required accuracy is great or the time interval is long; see my answer below. – David Ketcheson Nov 18 '16 at 5:09
• The point is that in double precision floating point you can't even represent a solution with better than 1.0e-16 relative accuracy. In many practical situations good old RKF45 will get you to that accuracy level over the period you're interested in without requiring tiny steps. It might not be a good choice for stiff systems or situations where a symplectic integrator is called for, but a higher order Runge Kutta method isn't a great solution in those situations either. I agree that for very long time periods higher order Runge Kutta methods can make some sense. – Brian Borchers Nov 18 '16 at 5:43

## The Benchmark Setup

In the Julia software DifferentialEquations.jl we implemented plenty of higher order methods, including the Feagin methods. You can see it in our list of methods, and then there are tons of others you can use as supplied tableaus. Because all of these methods are put together, you can easily benchmark between them. You can see the benchmarks I have online here, and see that it's very simple to benchmark many different algorithms. So if you want to take a few minutes to run the benchmarks, go for it. Here's a summary of what comes out.

First off it's important to note that, if you look at each of the benchmarks, you will see that our DP5 (Dormand-Prince Order 5) and DP8 methods are faster than the Hairer Fortran codes (dopri5 and dop853), and so these implementations are very well optimized. These show that as noted in another thread the over-use of the Dormand-Prince methods is because the methods are already written, not because they are still the best. So the real comparison between the most optimized implementations is between the Tsitorous methods, the Verner methods, and the Feagin methods from DifferentialEquations.jl.

## The Results

In general, the methods of an order higher than 7 have an added computational cost which is usually not outweighed by order, given the tolerances chosen. One reason for this is that the coefficient choices for lower order methods are more optimized (they have small "principle truncation error coefficients", which matter more when you're not asymtopically small). You can see that in many problems like here the Verner Efficient 7 and 8 methods do extremely well, but methods like the Verner Efficient 9 can have a lower slope. This is because the "gains" of higher order are compounding at lower tolerances, so there is always a tolerance where the higher order methods will be more efficient.

However, the question is then, how low? In a well-optimized implementation, that gets pretty low for two reasons. The first reason is because lower order methods implement something called FSAL (first same as last). This property means that the lower order methods re-use a function evaluation from the previous step in the next step, and thus have effectively one less function evaluation. If this is used properly, then something like a 5th order method (Tsitorous or Dormand-Prince) is actually taking 5 function evaluations instead of the 6 that the tableaus would suggest. This is also true for the Verner 6 method.

The other reason is due to interpolations. One reason to use a very high order method is to take less steps and simply interpolate intermediate values. However, in order to get the intermediate values, the interpolating function may need more function evaluations than were used for taking the step. If you look at the Verner methods, it takes 8 extra function evaluations for the Order 8 method to get an Order 8 interpolant. Many times the low order methods provide a "free" interpolant, for example most 5th order methods have a free 4th order interpolation (no extra function evaluations). So this means that if you need intermediate values (which you will need for a good plot if you're using a high order method), there are some extra hidden costs. Factor in the fact that these interpolated values are really important for event handling and solving delay differential equations and you see why the extra interpolation cost factors in.

## So What About the Feagin Methods?

So you will see that the Feagin methods are suspiciously missing from the benchmarks. They are fine, convergence tests work on arbitrary precision numbers, etc., but to actually get them to do well you need to be asking for some pretty absurdly low tolerances. For example, I found in unpublished benchmarks that the Feagin14 outperforms Vern9 (the 9th order Verner Efficient Method) at tolerances like 1e-30. For applications with chaotic dynamics (like in the Pleides or 3-body astrophysics problems), you may want this amount of accuracy due to sensitive dependence (errors in chaotic systems compound fast). However, Most people are probably computing with double-precision floating point numbers, and I haven't found a benchmark where they outperform in this domain of tolerance.

Here's a work-precision efficiency plot on this 100 dimensional ODE which showcases this: And here's a work-precision efficiency plot for the Pleiades ODE: With these we do not find good results when testing the Feagin methods with adaptivity in tolerances within the Float64-range. Possibly the issue could be the stability of the error estimate. Note that this is done with PI-adaptivity for stabilizing the stepping of each of the methods.

In addition, there is no interpolant to go along with the Feagin methods (Edit 2020: See Feagin's remark that one now exists, but is not available online. But note that this has no effect on the benchmarks since the benchmarks do not use the interpolant). So what I do is simply put a third order Hermite interpolation on them so that way one exists (and it works surprisingly well). However, if there's no standard interpolating function, you can do the recursive Hermite method (use this interpolation to get the midpoint, then do a 5th order interpolation, etc.) to get a high order interpolation, but this is very costly and the resulting interpolation doesn't necessarily have a low principle truncation error term (so it's only good when dt is really small, which is the exact opposite of the case we want!). So if you ever need a really good interpolation to match your accuracy, you need to at least go back to something like Vern9.

Note that extrapolation methods are simply algorithms for generating arbitrary order Runge-Kutta methods. However, for their order they take more steps than necessary and have high principle truncation error coefficients, and so they are not as efficient as a well-optimized RK method at a given order. But given the previous analysis, this means that there is a domain of extremely low tolerance where these methods will do better than the "known" RK methods. But in every benchmark I've ran, it seems I haven't gotten that low.

The choice really has nothing to do with stability issues. In fact, if you go through the DifferentialEquations.jl tableaus (you can just plot(tab) for the stability regions) you will see that most of the methods have suspiciously similar stability regions. This is actually a choice. Usually when deriving the methods, the author usually does the following:

1. Find the lowest principle truncation error coefficients (that is, the coefficients for the next order terms)
2. Subject to the order constraints
3. And make the stability region close to that of the Dormand-Prince Order 5 method.

Why the last condition? Well, because that method tends to be always stable with the way PI-controlled adaptive stepsize choices are done, so it's a good bar for "good enough" stability regions. So it's no coincidence that the stability regions all tend to be similar.

## Conclusion

There's tradeoffs in every choice of method. The highest order RK methods are simply not as efficient at lower tolerances both because it's harder to optimize the choice of coefficients, and because the number of function evaluations compounds (and grows even faster when interpolations are involved). However, if the tolerance gets low enough they win out, but the tolerances that are required can be far below "standard" applications (i.e. really only applicable to chaotic systems).

### Appendix: Benchmarking Code

using OrdinaryDiffEq, DiffEqDevTools, Plots
using Random
Random.seed!(123)
gr()
# 2D Linear ODE
function f(du,u,p,t)
@inbounds for i in eachindex(u)
du[i] = 1.01*u[i]
end
end
function f_analytic(u₀,p,t)
u₀*exp(1.01*t)
end
tspan = (0.0,10.0)
prob = ODEProblem(ODEFunction(f,analytic=f_analytic),rand(100,100),tspan)

abstols = 1.0 ./ 10.0 .^ (3:13)
reltols = 1.0 ./ 10.0 .^ (0:10);

setups = [Dict(:alg=>Tsit5())
Dict(:alg=>Vern7())
Dict(:alg=>Vern8())
Dict(:alg=>Feagin10())
Dict(:alg=>Feagin12())
Dict(:alg=>Feagin14())
Dict(:alg=>Vern9())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;save_everystep=false,numruns=100)
plot(wp)
savefig("feagin1.png")

f = (du,u,p,t) -> begin
@inbounds begin
x = view(u,1:7)   # x
y = view(u,8:14)  # y
v = view(u,15:21) # x′
w = view(u,22:28) # y′
du[1:7] .= v
du[8:14].= w
for i in 14:28
du[i] = zero(u)
end
for i=1:7,j=1:7
if i != j
r = ((x[i]-x[j])^2 + (y[i] - y[j])^2)^(3/2)
du[14+i] += j*(x[j] - x[i])/r
du[21+i] += j*(y[j] - y[i])/r
end
end
end
end

prob = ODEProblem(f,[3.0,3.0,-1.0,-3.0,2.0,-2.0,2.0,3.0,-3.0,2.0,0,0,-4.0,4.0,0,0,0,0,0,1.75,-1.5,0,0,0,-1.25,1,0,0],(0.0,3.0))

abstols = 1.0 ./ 10.0 .^ (9:12)
reltols = 1.0 ./ 10.0 .^ (6:9);

sol = solve(prob,Vern8(),abstol=1/10^12,reltol=1/10^10,maxiters=1000000)
test_sol = TestSolution(sol);

setups = [Dict(:alg=>Tsit5())
Dict(:alg=>Vern7())
Dict(:alg=>Vern8())
Dict(:alg=>Feagin10())
Dict(:alg=>Feagin12())
Dict(:alg=>Feagin14())
Dict(:alg=>Vern9())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=100,maxiters=1000)
plot(wp)
savefig("feagin2.png")


Just to add to Brian Borcher's excellent answer, many real-life applications admit highly stiff ODEs or DAEs. Intuitively, these problems experience nonsmooth, abrupt changes over time, so are better modeled using low-order polynomials spread finely over short step-sizes, as opposed to high-order polynomials stretched over long step-sizes. Also, stability often necessitates the use of implicit methods, for which the computational penalty of higher order methods is much steeper.

More rigorously, higher-order methods are less stable than lower-order methods for stiff problems. We have for example, the Dahlquist barriers for linear multistep methods.

Theorem (Second Dahlquist barrier). An A-stable multistep method must be of order $r\le2$. Among all multistep methods of order 2, the trapezoidal rule has the smallest error constant.

Similar (but far more complicated) statements can be made for L-stability in the RK formulas. In all cases, the increase in order often does not always lead to more accurate solutions. The following is an excerpt from Prothero and Robinson's seminal 1974 paper:

In using A-stable one-step methods to solve large systems of stiff nonlinear differential equations, we have found that
(a) some A-stable methods give highly unstable solutions, and
(b) the accuracy of the solutions obtained when the equations are stiff often appears to be unrelated to the order of the method used.

For even more rigorous treatments of this topic, see the classic text by Hairer & Wanner, "Solving ordinary differential equations II: Stiff and Differential - Algebraic Problems", 1991.

In practice, stiff equations are almost always solved using the trapezoidal rule or the TR-BDF2 formula (ode23t and ode23tb functions in MATLAB). Both of these are implicit second-order methods. Of course, where stability isn't an issue (i.e. in nonstiff equations) we are free to chose from a number of options; RK45 is the most common choice.

• Very interesting. Is there any (intuitive) explanation for why the order has to be less than or equal to 2 for it to be an A-stable multistep method? And just to clarify, when you say similar statements can be made for RK formulas, is it of order 2 once again? – Mathews24 Nov 18 '16 at 3:58
• But for Runge-Kutta methods, there are A-stable methods of arbitrary order. – David Ketcheson Nov 18 '16 at 4:37
• @DavidKetcheson Yes, but they are not strongly A-stable (i.e. L-stable). They have a lot of issues when used to solve DAEs, e.g. simulate simple transistor circuits. Indeed, TR is infamous for causing artificial ringing in SPICE, which is what motivated the development of TR-BDF2. – Richard Zhang Nov 18 '16 at 8:07
• @DavidKetcheson For reference, see doi.org/10.1090/S0025-5718-1974-0331793-2. The notion of A-stability isn't quite strong enough for DAEs, and high-order A-stable methods often produce weird results when used to solve DAEs. – Richard Zhang Nov 18 '16 at 8:14
• Sure, but the question isn't about DAEs or about multistep methods. – David Ketcheson Nov 18 '16 at 18:59

Well, in answer to the question about the use of higher-order RK methods, I have a few things to add to the discussion:

1. As an example, it is often helpful to plot an "efficiency diagram" for several different methods being applied to a given initial value problem (using a variety of stepsizes). On such a diagram, one can plot (for each method) how the accuracy (usually measured according to the number of significant digits obtained) varies with log10 of the number function/derivative calls used to integrate over a particular interval.

On such a plot, one can plot several different methods simultaneously in order to obtain a pretty good idea of which method performs the best (or provides the most efficient integration for a given desired accuracy). On such a plot, the slope of the line for each method is approximately equal to the order of the method (for those stepsizes for which the truncation error is dominant).

Eventually, of course, as the stepsize becomes very small and the number of function/derivative calls becomes very large, the rounding errors dominate and the curves bend more-or-less downward. On the efficiency diagram the rounding error is not visible for the fairly large stepsizes taken. For this diagram, the problem being solved is the two-body problem with the eccentricity of 0.4 and the methods being compared are some Explicit RK methods of orders 4, 6, 8, 10, 12. One may observe that if you only want 4-6 digits of accuracy, then it really doesn't make a lot of difference which method you might wish to utilize. If you want more accuracy, then it becomes more important which method you employ.

The methods shown on the diagram include RK4, Rk6B (A 6th order method due to Butcher), RK8CV (an 8th order method of Cooper and Verner), RK10H (a 10th order method due to Hairer), and RK12 (my 12th order method, which happens to have an embedded 10th order methods so that you can estimate the local truncation errors). My 14th order method is not shown, and it is not particularly efficient unless you need extremely high accuracy), as noted by others above.

1. My RK12(10), as it is now known, DOES have the ability to provide interpolated results of order 7. You can still obtain 12th order results at the end of each step and in addition, if you need so-called "dense output", you can use the 7th order results at any point you wish within the step. The best way to describe the method would be for me to send you a copy of the code. I would rather not publish the code on a website. That way, I can make sure that any updates, etc. can be distributed to those who are interested.

Just send any request to feagin@uhcl.edu (I teach Computer Science at the Univ. of Houston - Clear Lake).

• Hey, thank you for commenting! I will be asking for the interpolation to start testing it in the interpolation-based tests as well. I noticed some of the links from my response needed an update (post from 2016!), so they are now updated and I put some plots (and code) for the work-precision efficiency graphs in Float64-tolerances. Hopefully that clarifies a few things. If there have been changes lately to the coefficients that we should know about, please let us know! – Chris Rackauckas Feb 19 at 1:17