Which Runge-Kutta method is more accurate: Dormand-Prince or Cash-Karp?

Question

I simply want to know whether the Dormand-Prince Numerical Method or the Cash-Karp Numerical Method is more accurate.

Answer 1

Since I just finished optimizing a lot of them in software, DifferentialEquations.jl, I decided to just lay out a comparison of the main Order 4/5 methods. The Fehlberg method was left out because it's commonly known to be less efficient than the DP5 method.

Backstories

Dormand-Prince 4/5

The Dormand-Prince method was developed to be accurate as a 4/5 pair with local extrapolation usage (i.e. step with the order 5 pair. This is because it was designed to be close to optimal (i.e. minimal) principle truncation error coefficient (under the restraint of also having the minimal number of steps to achieve order 5). It has an order 4 interpolation which is free but needs extra steps for an order 5 interpolation.

Cash-Karp 4/5

The Cash-Karp method was developed to satisfy different constraints, namely to deal with non-smooth problems better. They chose the $c_i$, the percentage of the timestep in the $i$th step (i.e. $t+c_i \Delta t$ is the time the $i$th step is calculated at) to be as uniform as possible, yet still achieve order 5. Then it also was derived to have embedded 1st, 2nd, 3rd, and 4th order methods with this uniformity of the $c_i$. They are spaced in such a manner that you can find out where a stiff part starts by which difference is large. Moreover, note that the more stiff the equation, the worse a higher order method does (because it needs bounds on higher derivatives). So they develop a strategy that uses the 5 embedded methods to "quit early": i.e. if you detect stiffness, stop at stage $i<6$ to decrease the number of function calls and save time. So in the end, this "pair" was developed with a lot of other constraints in mind, and so there's no reason to expect it would be "more accurate", at least as a 4/5 pair. If you add all of this other machinery then, on (semi-)stiff problems, it will be more accurate (but in that case you may want to use a different method like a W-Rosenbrock method). This is one reason why this pair hasn't become standard over the DP5 pair, but it still can be useful (maybe it would be good for a hybrid method that switches to a stiff solver when stiffness is encountered?).

Bogacki & Shampine 4/5

To round out the answer, let's discuss the Bogacki & Shampine pair that was mentioned in the comment. The BS5 method drops the constraint of "using the least function calls" (it uses 8 instead of 6) in order to do 2 things:

• Get really low principle truncation error coefficients.

• Produce an order 5 interpolation with lower error coefficients.

These coefficients are so low that for many problems with tolerances that users likely use, it measures as though it's 6th order. Their paper shows that for cheap function calls, this can be more efficient than DP5 by about the same amount as DP5 was over RKF5 (the Fehlberg method).

You might put two-and-two together and see: wait a second, Shampine is the same person who developed the MATLAB ODE suite, this was after the BS5 pair paper was published, why doesn't MATLAB's ode45 use the BS5 pair? One reason is, it was mostly done before the BS5 pair was released. The other reason is that the ode45 function was developed to minimize time. While the BS5 pair is more efficient (i.e. gets lower accuracy), the purpose of ode45 is to have a good enough error to make a good enough plot. This means that, in order to deal with the large steps, it also produces two extra interpolated solutions between every step. For the DP5 method, there is a "free" order 4 interpolation, and so this is much faster than using BS5. Since it is also "accurate enough" at moderate tolerances, this method is set as the standard because it gives a better standard user experience than BS5 when doing interactive computing (so this choice was context-specific).

Tsitouras 4/5

Here's one fewer people know about. It's derived in this paper. It's derived using fewer assumptions than the DP5 method and tries to get a pair with lower principle truncation error coefficients. In its tests, it states that it achieves this. It also has a free order 4 interpolation like the DP5 method.

Numerical tests

I wrote the numerical package DifferentialEquations.jl to be a pretty comprehensive set of solvers for Julia. Along the warpath, I implemented over 100 Runge-Kutta methods, and hand-optimized plenty. Three of the hand-optimized integrators are the DP5, BS5, and Tsit5 methods (I did not do CK5 because, as noted in the backstory, its main case is for problems that are kind of stiff. I think the better way to handle them is to use DP5/BS5 and switch to stiff solvers as necessary in a manner like LSODE, but that's a story for a different time) (one way to see they are close to optimal is that these methods are faster than the Hairer dopri5 implementations, so they are at least decent implementations). Tests between a lot of Runge-Kutta methods on nonstiff equations can be found on this benchmarks page. I am adding more as I go along, but you can see from the linear ODE and the Three-Body problem work precision diagrams, I measure the DP5 and Tsit5 methods to have almost identical efficiency, beating out the BS5 method in the linear ODE, while it's DP5 and BS5 that are almost identical on the Three-Body problem with Tsit5 behind. From this information, at least for now, I have settled on the DP5 method as the default, matching previous recommendations. That may change with future tests (or you could add benchmarks! Feel free to contribute, or Star the repo to give this effort more support).

---

Conclusion

In conclusion, the Order 5 pairs go like this:

• The Dormand-Prince 4/5 pair is a good go-to pair since it's well-optimized in terms of principle truncation error coefficient and has a cheap order 4 interpolation, which makes it fast for producing decent plots.

• The Cash-Karp pair has more constraints on it to better handle stiff equations. However, to get the full benefit you'll want to use the full algorithm with the 5 embedded methods.

• Bogacki & Shampine Order 5 method may be the most efficient in terms of producing error per function calls (it has a double error estimator, so in harder problems, it probably does better), but that allows it to take larger timesteps. However, if you just want to produce a smooth plot, you then have to counter-act this method: use a lower tolerance (so it will take longer than DP5 but with less error) or use more interpolated steps. In the end, this meant that it might not be better for interactive applications, although it might be better for some scientific computing applications.

• The Tsitouras 4/5. Was developed fairly recently (2011) to beat out the DP5 in a head-to-head comparison. My tests don't give me a reason to believe that it's so much better than DP5 that it should now be considered the new standard method, but future tests may begin to side in its favor.

Edit

---

I did improve the Tsit5 implementation. It now does better than DP5 on most tests, both the DifferentialEquations.jl and the Hairer dopri implementations (though one might be surprised that the DifferentialEquations.jl implementations are actually faster, which of course helps the Tsit5 implementation). I now recommend it as the default order 4/5 method.

Answer 2

If you are interested in comparing two integrators solve

$$\frac{dq}{dt} = p \\ \frac{dp}{dt} = -q$$

with initial values $(q,p) = (1,0)$ and $dt = 0.1$ for both of them. Then plot the errors

$$dq = q - \cos(t) ; dp = p + \sin(t)$$

for a reasonable range of $t$ (0 to 100, one thousand time steps in all) and choose the integrator with the smaller error. This harmonic oscillator test will show you the phase and amplitude errors with very little effort.

• P J Channell and C Scovel. Symplectic integration of Hamiltonian systems. Journal of Computational Physics 73, 468 (1973).