Runge-Kutta 8(5,3)

Question

This is actually three small very related questions about Runge-Kutta methods.

1. I have programmed a RK 7(8) method also RK 4(5). At the beginning I was assuming that the RK 7(8) uses two approximations of different order, one of order 7 an another of order 8. The difference between the two approximations is used to estimate the error of integration, and the algorithm returns the approximation of order 8. But by using an system of ODE's for which I know the exact answer (as a test method), I have seen that the approximation of order 7 gives a smaller error. As when we write RK 7(8) we write first the 7, is it supposed that the method is of order 7 or 8?

2. When we say order $k$, do we mean that the approximation is up to order $k$ or that the error is of order $k$?

3. Python programming language provides a routine called dop853 that performs a Runge Kutta 8(5,3). What does it mean exactly when the method is specified by three numbers.

Thank you very much.

Answer 1

I am still learning but possibly,

(1) Could be accumulation of floating point approximation errors causing this. (2) The approximation is up to order 4, truncation error is then order 5.

Answer 2

Don't know much about python, so can't help there. As for the other two questions:

Each numerical intergration method has a local truncation error $e(h)$. Its meaning is an estimate of the error one makes, when you start a precisely calculated point (with no error) and then make one step of size $h$ with your numerical integrator. For instance, if you start from the initial point at $t = 0$, where both the "true" $f(t)$ and numerical $\hat f(t)$ solutions are identical, and then you calculate the solution value at $t = h$ using, say, $RK_4$, you will have $|f(h) - \hat f(h)| = O(h^5)$.

And no, that's not a typo. We say that a method has order $p$, when the local truncation error is $O(h^{p+1})$. The intuitive reason for that is that the errors will accumulate. So if you solve from $0$ to $1$, and have a step size $h = \frac1n$, then you will make $n$ steps. If each steps error is $O(\frac{1}{n^{p+1}})$ then the total error after $n$ summations will be roughly $O(\frac{1}{n^{p}})$.

Now. Concerning your $RK_{4,5}$ methods. Again, no idea what actually is implemented in python, but I'll make a wild guess.

My guess is that the method is actually not fixed-step-size, but variable step size. (We use those because some ODE systems are "nice", and do not require small step sizes in order to achieve good order of approximation; while others "misbehave", and do require extremely small $h$. You do not want to solve a "nice" system with unnecessarily small step size - that's a waste of time.)

These methods vary the step size in order to keep the local error at a pre-fixed (by the user) order of magnitude. They do that by calculating the solution twice - once with a lower order method ($4$ or $7$) and once with the higher order ($5$ or $8$), and use the difference as an estimate of the local truncation error.

Since there is no longer a fixed step size, the notion of order for these methods is no longer as easily applicable. Both $RK_{45}$ and $RK_{78}$ will produce a solution that will have local truncation error fixed at the user-inputed value at each step. The difference between the two, I think, is that $RK_{78}$ will make fewer steps, but with more function evaluations at each step.

Answer 3

Error-estimating Runge-Kutta methods, as you note, are often written as "RK M(N)", where typically $M$ specifies the order of the method used to obtain the solution, and $N$ specifies the order of the method used to obtain the error estimate.

The purpose of using error control is to increase the accuracy of a method by way of killing two birds with one stone. Since we're already evaluating the function $f(t,x)$ at least $M$ times anyways (and usually more -- it depends on the stage of the method), the idea is that maybe we can clamp down the error by throwing in an extra function call or so and getting an increase in accuracy.

In essence, the higher-order method has the lower-order method embedded within it; they share a Butcher tableau in common except for the final row.

---

DOP853 is a more complicated routine.

At its core, DOP is a Dormand & Prince 8th order method that uses a 5th order error estimator and bootstraps on a 3rd order estimator to obtain a dense output of order 7. Therefore, DOP853 could be thought of ultimately as a 7th order method. Dense output is an interesting concept -- essentially, a dense output of a method provides a solution that can be used as a good linear interpolant for the solution everywhere between output points.

Dormand & Prince methods are constructed in such a way that one or more of the later stages at time $t_i$ are identical to the earlier stages at $t_{i+1}$; for instance, in the classical 7-stage Dormand-Prince 4(5) method, the stage $X_7$ at time $t_i$ is identical to the stage $X_1$ at time $t_{i+1}$.

The 8th order method with the 5th order error estimator used in DOP853, if used alone, would generate a dense output of order 6; by adding on another 3-stage 3rd order method, we can increase the dense output to order 7 by forming a 16-by-16 linear system at each step. This means that not only is our output order 7 accurate at each output point, but can be linearly extrapolated from any point to obtain an order 7 solution, as well.