4.2. Stability Functions

To study the stability properties of a numerical method we consider how it propagates small error from one step to the next. A convenient model problem for this purpose of the linear test equation

(4.2)\[ y' = \lambda y, \]

where \(\lambda \in \mathbb{C}\) is a constant. This equation is simple enough to analyse exactly, having solution

\[ y(t) = e^{\lambda t},\]

yet it captures the essential growth or decay behaviour that numerical methods must reproduce.

Suppose a numerical mtehod is applied to this test equation. Since the method advances the numerical solution from \(y_n\) to \(y_{n+1}\), any perturbation or error is propagated by exactly the same update formula. The resulting factor by which the solution (and hence the error) is multiplied in one step is called the stability function.

Definition 4.2 (Stability function)

The stability function of a method, denoted by \(R(z)\), is the rate of growth over a single step of the method when applied to calculate the solution of an ODE of the form \(y'=\lambda t\) where \(z = h \lambda\) and \(h\) is the step size

\[ y_{n+1} = R(z)y_n. \]

For example, if the Euler method is used to solve the test equation \(y'=\lambda y\) then the solution will be updated over one step using

\[ y_{n+1} = y_n + h f(t_n ,y_n) = y_n + h \lambda y_n, \]

Let \(z = h\lambda\) then

\[ y_{n+1} = y_n + z y_n =(1+z)y_n. \]

So the stability function of the Euler method is \(R(z) = 1 + z\).

---

4.2.1. Stability function of a Runge-Kutta method

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_j\) 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)^\mathsf{T}\) and \(\mathbf{e}=(1,1,\dots ,1)^\mathsf{T}\) 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. \]

Rearranging equation (4.3) gives

\[\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}\]

Substituting into equation (4.4)

\[\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*} \]

This formula is valid for any Runge-Kutta method.

Definition 4.3 (Stability function for a general Runge-Kutta method)

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

---

4.2.2. Stability function of an explicit Runge-Kutta method

Equation (4.5) provides the stability function for a general Runge-Kutta method. For explicit Runge-Kutta methods, we can represent the inverse matrix \((I - zA)^{-1}\) as a power series summation.

The Neumann series for matrices is

\[ (I - M)^{-1} = \sum_{k=0}^\infty M^k, \]

so

\[ (I - zA)^{-1} = \sum_{k=0}^\infty (zA)^k. \]

Since for an explicit method, the matrix \(A\) is strictly lower triangular then

\[A^s = 0,\]

and

\[ (I - zA)^{-1} = \sum_{k=1}^{s-1} (zA)^k = I + zA + z^2A^2 + \cdots + z^{s-1}A^{s-1}.\]

Substituting into equation (4.5)

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

which is the stability function for an ERK method.

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

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

So the stability function of an \(s\)-stage explicit method is a polynomial function of at most degree \(s\)

(4.6)\[ R(z) = 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 + \mathbf{b}^\mathsf{T} A^{s-1} \mathbf{e} z^s. \]

The exact solution satisfies

\[ y_{n+1} = e^z y_n, \]

therefore, for an order \(p\) method, the stability function should approximate \(e^z\) with an error of \(O(z^{p+1})\) as \(z \to 0\). The series expansion of \(e^z\) is

(4.7)\[ e^z = \sum_{k=0}^{\infty} \frac{1}{k!} z^k = 1 + z + \frac{1}{2!}z^2 + \cdots + \frac{z^p}{p!} + O(z^{p+1}). \]

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

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

These are necessary conditions for an ERK method to have order \(k\). For general Runge-Kutta methods additional order conditions are also required.

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, 3\).

\[\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.7) up to and including the \(z^2\) term showing that the method has linear order 2.

---

4.2.3. Code

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

Python

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,9), sp.Rational(1,3), sp.Rational(2,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)

MATLAB

% Define ERK method
A = [0, 0, 0, 0 ;
     1/2, 0, 0, 0 ;
     0, 3/4, 0, 0 ;
     2/9, 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