7
$\begingroup$

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

$\endgroup$
  • 6
    $\begingroup$ They're both 4th/5th order methods that use 6 function evaluations per step. Performance and accuracy are likely to be very similar for most problems but on any particular problem one method might be slightly better than the other. $\endgroup$ – Brian Borchers Sep 17 '16 at 3:17
  • $\begingroup$ I'll post a more thorough answer when I have time, but I'll wager that you'll get more accuracy than either with the same effort using the pair of Bogacki & Shampine. $\endgroup$ – David Ketcheson Sep 17 '16 at 13:59
31
$\begingroup$

Since I just finished optimizing a lot of them in a 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 which 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 which 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 its 6th order. Their paper shows that for cheap function calls, this can be more efficient than DP5 by about the same amount at 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 relaesed. The other reason is because 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 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).

Tsitorous 4/5

Here's one less people know about. It's derived in this paper. It's derived using less 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, it's main case is for problems which 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 in the benchmarks folder. 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're just wanting 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 Tsitorous 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 as 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.

$\endgroup$
2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.