4.2. Stability Function of a Runge-Kutta Method#

We have seen that the stability function is the function \(R(z)\) such that \(y_{n+1} = R(z) y_n\) where \(y_{n+1}\) is the single step solution of the test ODE \(y' = \lambda y\) and \(z = h\lambda\). The general form of a Runge-Kutta method is

\[\begin{split} \begin{align*} y_{n+1}&=y_n + h \sum_{i=1}^s b_i k_i ,\\ k_i &=f(t_n +c_i h,y_n +h\sum_{j=1}^s a_{ij} k_j ). \end{align*} \end{split}\]

Let \(Y_i = y_n + h \displaystyle \sum_{j=1}^s a_{ij} k_i\) and applying the method to the test equation \(y' = \lambda y\) we have

\[\begin{split} \begin{align*} y_{n+1} &=y_n +h\lambda \sum_{i=1}^s b_i Y_i ,\\ Y_i &=y_n +h\lambda \sum_{j=1}^s a_{ij} Y_j . \end{align*} \end{split}\]

Let \(z = h\lambda\) and expanding out the summations in the stage values

\[\begin{split} \begin{align*} Y_1 &=y_n +z(a_{11} Y_1 + a_{12} Y_2 + \cdots + a_{1s} Y_s),\\ Y_2 &=y_n +z(a_{21} Y_1 + a_{22} Y_2 + \cdots + a_{2s} Y_s),\\ &\vdots \\ Y_s &=y_n +z(a_{s1} Y_1 + a_{s2} Y_2 + \cdots + a_{ss} Y_s). \end{align*} \end{split}\]

We can write these as the matrix equation

\[\begin{split} \begin{align*} \begin{pmatrix} Y_1 \\ Y_2 \\ \vdots \\ Y_s \end{pmatrix} = \begin{pmatrix} y_n \\ y_n \\ \vdots \\ y_n \end{pmatrix} + z \begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1s} \\ a_{21} & a_{22} & \cdots & a_{2s} \\ \vdots & \vdots & \ddots & \vdots \\ a_{s1} & a_{s2} & \cdots & a_{ss} \end{pmatrix} \begin{pmatrix} Y_1 \\ Y_2 \\ \vdots \\ Y_s \end{pmatrix}. \end{align*} \end{split}\]

Let \(Y = (Y_1 ,Y_2 ,\dots ,Y_s )\) and \(\mathbf{e}=(1,1,\dots ,1)\) then

(4.3)#\[ Y = \mathbf{e} y_n + z A Y. \]

Expanding the summation in the equation for updating the solution over a single step and using \(z = h \lambda\) gives

\[ y_{n+1} = y_n + z (b_1 Y_1 + b_2 Y_2 + \cdots + b_s Y_s), \]

which can be written as the matrix equation

(4.4)#\[ y_{n+1} =y_n + z \mathbf{b}^\mathsf{T} Y. \]

The stability functions for explicit and implicit Runge-Kutta methods are derived using different approaches.

4.2.1. Stability function of an explicit Runge-Kutta method#

To derive the stability function of an explicit Runge-Kutta method we rearrange equation (4.3)

\[\begin{split} \begin{align*} Y - zAY &= \mathbf{e}y_n \\ (I - zA)Y &= \mathbf{e}y_n\\ Y &= (I - zA)^{-1} \mathbf{e}y_n, \end{align*} \end{split}\]

and substituting into equation (4.4) gives

\[\begin{split} \begin{align*} y_{n+1} &= y_n +z\mathbf{b}^\mathsf{T} (I - zA)^{-1} \mathbf{e} y_n \\ &= (1 + z \mathbf{b}^\mathsf{T} (I - zA)^{-1} \mathbf{e})y_n , \end{align*} \end{split}\]

so the stability function is

\[ \begin{align*} R(z) = 1 + z\mathbf{b}^\mathsf{T} (I - zA)^{-1} \mathbf{e} \end{align*} \]

Using the geometric series of matrices

\[ \begin{align*} (I - zA)^{-1} = \sum_{k=0}^{\infty} (zA)^k \end{align*} \]

then

\[ \begin{align*} R(z) &= 1 + z \mathbf{b}^\mathsf{T} \sum_{k=0}^{\infty} (zA)^k \mathbf{e} = 1 + \sum_{k=0}^{\infty} \mathbf{b}^\mathsf{T} A^k \mathbf{e}\,z^{k+1}. \end{align*} \]

This is the stability function for an ERK method.

Definition 4.3 (Stability function of an explicit Runge-Kutta method)

(4.5)#\[\begin{split} \begin{align*} R(z) &= 1 + \sum_{k=0}^{\infty} \mathbf{b}^\mathsf{T} A^k \mathbf{e}\,z^{k+1} \\ &= 1 + \mathbf{b}^\mathsf{T} \mathbf{e}\,z+\mathbf{b}^\mathsf{T} A\mathbf{e}\,z^2 +\mathbf{b}^\mathsf{T} A^2 \mathbf{e}\,z^3 + \cdots \end{align*} \end{split}\]

Since the solution to the test equation is \(y = e^{\lambda t}\), over one step of an explicit Runge-Kutta method we would expect the local truncation errors to change at a rate of \(e^z\). The series expansion of \(e^z\) is

(4.6)#\[ e^z = \sum_{k=0}^{\infty} \frac{1}{k!} z^k = 1 + z + \frac{1}{2!}z^2 + \frac{1}{3!}z^3 + \frac{1}{4!}z^4 + \cdots \]

Comparing the coefficients of \(z^k\) in equations (4.5) and (4.6) we have

\[ \begin{align*} \frac{1}{k!}= \mathbf{b}^\mathsf{T} A^{k-1} \mathbf{e}, \end{align*} \]

which must be satisfied up to the \(k\)th term for an order \(k\) explicit method to be stable.

Example 4.1

An explicit Runge-Kutta method is defined by the following Butcher tableau

\[\begin{split} \begin{array}{c|cccc} 0 & & & & \\ \frac{1}{2} & \frac{1}{2} & & & \\ \frac{3}{4} & 0 & \frac{3}{4} & & \\ 1 & \frac{2}{9} & \frac{1}{3} & \frac{4}{9} & \\ \hline & \frac{7}{24} & \frac{1}{4} & \frac{1}{3} & \frac{1}{8} \end{array} \end{split}\]

Determine the stability function for this Runge-Kutta method and hence find its order.

Solution

Calculating \(\mathbf{b}^\mathsf{T}A^{k - 1}\mathbf{e}\) for \(k = 1, \ldots, 4\).

\[\begin{split} \begin{align*} k &= 1: & \mathbf{b}^\mathsf{T}A^{0}\mathbf{e} &= 1, \\ k &= 2: & \mathbf{b}^\mathsf{T}A^{1}\mathbf{e} &= (\tfrac{7}{24}, \tfrac{1}{4}, \tfrac{1}{3}, \tfrac{1}{8}) \begin{pmatrix} 0 & 0 & 0 & 0 \\ \frac{1}{2} & 0 & 0 & 0 \\ 0 & \frac{3}{4} & 0 & 0 \\ \frac{2}{9} & \frac{1}{3} & \frac{4}{9} & 0 \end{pmatrix} \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \end{pmatrix} = \frac{1}{2}, \\ k &= 3: & \mathbf{b}^\mathsf{T}A^{2}\mathbf{e} &= (\tfrac{7}{24}, \tfrac{1}{4}, \tfrac{1}{3}, \tfrac{1}{8}) \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ \frac{3}{8} & 0 & 0 & 0 \\ \frac{1}{6} & \frac{1}{3} & 0 & 0 \end{pmatrix} \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \end{pmatrix} = \frac{3}{16}, \\ k &= 4: & \mathbf{b}^\mathsf{T}A^{3}\mathbf{e} &= (\tfrac{7}{24}, \tfrac{1}{4}, \tfrac{1}{3}, \tfrac{1}{8}) \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ \frac{1}{6} & 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \end{pmatrix} = \frac{1}{48}. \\ \end{align*} \end{split}\]

Therefore the stability function is

\[ \begin{align*} R(z) = 1 + z + \frac{1}{2} z^2 + \frac{3}{16} z^3 + \frac{1}{48} z^4, \end{align*} \]

which agrees to the series expansion of \(e^z\) from equation (4.6) up to and including the \(z^2\) term. Therefore this method is of order 2.

4.2.2. Code#

The Python and MATLAB code used to determine the stability function for the ERK method from Example 4.1 is given below.

import sympy as sp

# Define ERK method
A = sp.Matrix([[0, 0, 0, 0],
               [sp.Rational(1,2), 0, 0, 0],
               [0, sp.Rational(3,4), 0, 0],
               [sp.Rational(2,3), sp.Rational(1,3), sp.Rational(4,9), 0]])
b = sp.Matrix([sp.Rational(7,24), sp.Rational(1,4), sp.Rational(1,3), sp.Rational(1,8)])
e = sp.ones(4,1)

# Determine coefficients of the stability function
for k in range(4):
    sp.pprint(b.T * A ** k * e)

% Define ERK method
A = [0, 0, 0, 0 ;
     1/2, 0, 0, 0 ;
     0, 3/4, 0, 0 ;
     2/3, 1/3, 4/9, 0];
b = [7/24 ; 1/4 ; 1/3 ; 1/8];
e = ones(4, 1);

% Determine coefficients for the stability function
for k = 1 : 4
    sym(b' * A ^ (k - 1) * e)
end