Puzzling remark about stability region of fifth-order Runge-Kutta method

Question

I came across a puzzling remark in the paper

P. J. van der Houwen, The development of Runge-Kutta methods for partial differential equations, Appl. Num. Math. 20:261, 1996

On lines 8ff on page 264, van der Houwen writes:

"For the Taylor polynomials this implies that the imaginary stability interval is empty for $p = 1, 2, 5, 6, 9, 10, \cdots$"

where Taylor polynomial refers to the stability polynomial (truncated expansion of $\exp(x)$ around $x=0$) of the Runge-Kutta method and p is the order (see page 263). I assume that I misunderstand something because the fifth-order Runge-Kutta method does not have an empty imaginary stability interval as far as I know. From what I remember, the imaginary limit is about 3.4 or so.

What's my misunderstanding?

Answer 1

van der Houwen's statement is correct, but it is not a statement about all fifth-order Runge-Kutta methods. The "Taylor polynomials" he is referring to are (as you seem to know) just the polynomials of degree $p$ that approximate $\exp(z)$ to order $p$:

$$P_p(z) = \sum_{j=1}^p \frac{z^j}{j!}$$

For the fifth-order polynomial, it turns out that $|P_5(i\epsilon)|>1$ for small $\epsilon$, so the stability region of a method having $P_5(z)$ as its stability polynomial does not include any neighborhood of the origin on the imaginary axis. That is, in precise terms, what van der Houwen says.

The most likely source of your confusion is what is meant by "the fifth-order Runge-Kutta method". There are (infinitely) many fifth-order Runge-Kutta methods, but the most well-known ones do not have $P_5(z)$ as their stability polynomial. Why? As John Butcher famously proved, a fifth-order Runge-Kutta method must have at least six stages. Usually, the stability polynomial of a method with six (or more) stages would have degree six (or more). For instance, each of the fifth-order methods listed on this Wikipedia page uses exactly six stages and has a stability polynomial of degree six.

Is it possible for a fifth-order method to have $P_5(z)$ as its stability polynomial? Yes; a fifth-order explicit extrapolation method (like the well-known ones reviewed in this paper of mine) would do so. Note also that a $p$-stage Runge-Kutta method with stability polynomial $P_5(z)$ will be accurate to order 5 for linear ODEs, though not for nonlinear ODEs.

Finally, it is easy to make mistakes when determining the extent of the imaginary stability interval for high-order Runge-Kutta methods. That is because the boundary of the stability region for such methods lies extremely close to the imaginary axis. Therefore, roundoff errors can lead to incorrect conclusions; only exact calculations should be used (of course, the relevance of the stability region boundary for practical purposes under these circumstances could certainly be debated).

For instance, here is a plot of the stability region of the fifth-order method from the Fehlberg 5(4) pair:

The imaginary stability interval is empty, but you can't tell from the picture at this resolution! Note that the region clearly includes part of the imaginary axis, but no interval about the origin.

Meanwhile, here is the plot for the fifth-order method from the Dormand-Prince 5(4) pair:

It has imaginary stability interval approximately $[-1,1]$.

For a precise characterization of the stability region boundary near the imaginary axis for $P_p(z)$ (which is rather fascinating!), see my recent paper.

You may also be interested in the NodePy package, which produced the plots above and which can be used to accurately determine things like the imaginary stability interval of a method (disclaimer: I created NodePy).