# 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

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

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

## 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 6 and 7 methods do extremely well, but methods like the Verner Efficient 8 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.

In addition, there is no interpolant to go along with the Feagin methods. 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: