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

Question

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?

Answer 1

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.

Note also that high-order extrapolation and deferred correction methods are Runge-Kutta methods.

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

You ask

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

Answer 2

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 Tsitouras 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 (Tsitouras 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 About Extrapolation

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.

Note About Stability

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[1])
  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")

Answer 3

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.

Answer 4

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.

2. 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 [email protected] (I teach Computer Science at the Univ. of Houston - Clear Lake).

Answer 5

Dr Chasnov gives the reason for popularity of the 4th order RK method here (at 3:50): https://youtu.be/40NvJNrmHhw.

I am curious about why the number of evaluations for f(x,y) jumps at certain orders of the expansion e.g. the 5th order. Very tempted to just go at it by hand to see how this falls out of the calculations.