4.3. Stability Functions for Implicit Runge-Kutta Methods#

To determine the stability function for an implicit Runge-Kutta method we consider the stability function of a general Runge-Kutta method given in equations (4.3) and (4.4) which are repeated below

\[\begin{split} \begin{align*} y_{n+1} & = y_n + z \mathbf{b}^\mathsf{T} Y, \\ Y &= \mathbf{e} y_n + z A Y, \end{align*} \end{split}\]

where \(Y = (Y_1, Y_2, \ldots, Y_s)\). Rewriting these gives

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

We can write these as a matrix equation

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

Using Cramer’s rule to solve this system for \(y_{n+1}\) we have

\[\begin{split} \begin{align*} y_{n+1} = \frac{\det \begin{pmatrix} y_n & -zb_1 & -zb_2 & \cdots & -zb_s \\ y_n & 1-za_{11} & -za_{12} & \cdots & -za_{1s} \\ y_n & -za_{21} & 1-za_{22} & \cdots & -za_{2s} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ y_n & -za_{s1} & -za_{s2} & \cdots & 1-za_{ss} \end{pmatrix}}{\det (I-zA)}. \end{align*} \end{split}\]

Performing a row operation of subtracting the first row of matrix from all other rows gives

\[\begin{split} \begin{align*} y_{n+1} =\frac{\det \begin{pmatrix} y_n & -zb_1 & -zb_2 & \cdots & -zb_s \\ 0 & 1-za_{11} +zb_1 & -za_{12} +zb_2 & \cdots & -za_{1s} +zb_s \\ 0 & -za_{21} +zb_1 & 1-za_{22} +zb_2 & \cdots & -za_{2s} +zb_s \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & -za_{s1} +zb_1 & -za_{s2} +zb_2 & \cdots & 1-za_{ss} +zb_s \end{pmatrix}}{\det(I-zA)}. \end{align*} \end{split}\]

This provides the following definition of the stability function for an implicit Runge-Kutta method.

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

The stability function of an implicit Runge-Kutta method is

(4.7)#\[R(z) = \frac{\det (I - zA + z\mathbf{e}\mathbf{b}^\mathsf{T})}{\det(I - zA)}.\]

where \(\mathbf{e}\mathbf{b}^\mathsf{T}\) is a diagonal matrix with the elements of \(\mathbf{b}\) on the main diagonal

\[\begin{split} \begin{align*} \mathbf{e}\mathbf{b}^\mathsf{T} = \begin{pmatrix} b_1 & 0 & \cdots & 0 \\ 0 & b_2 & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & b_s \end{pmatrix}. \end{align*} \end{split}\]

Example 4.2

The Radau IA IRK method is defined by the following Butcher tableau

\[\begin{split} \begin{array}{c|cc} 0 & \frac{1}{4} & -\frac{1}{4} \\ \frac{2}{3} & \frac{1}{4} & \frac{5}{12} \\ \hline & \frac{1}{4} & \frac{3}{4} \end{array} \end{split}\]

Determine the stability function of this method

Solution

Using equation (4.7)

\[\begin{split} \begin{align*} R(z) &= \frac{\det \left( \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} - z \begin{pmatrix} \frac{1}{4} & -\frac{1}{4} \\ \frac{1}{4} & \frac{5}{12} \end{pmatrix} + z \begin{pmatrix} \frac{1}{4} & 0 \\ 0 & \frac{3}{4} \end{pmatrix} \right) } { \det \left( \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} - z \begin{pmatrix} \frac{1}{4} & -\frac{1}{4} \\ \frac{1}{4} & \frac{5}{12} \end{pmatrix} \right) } \\ &= \frac{ \det \begin{pmatrix} 1 & \frac{1}{4}z \\ -\frac{1}{4}z & 1 + \frac{1}{3}z \end{pmatrix} } { \det \begin{pmatrix} 1 - \frac{1}{4}z & \frac{1}{4}z \\ -\frac{1}{4}z & 1 - \frac{5}{12}z \end{pmatrix} } = \frac{1 + \frac{1}{3}z + \frac{1}{16}z^2}{1 - \frac{2}{3}z + \frac{1}{6}z^2}. \end{align*} \end{split}\]

4.3.1. Code#

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

import sympy as sp

# Define IRK method
A = sp.Matrix([[sp.Rational(1,4), -sp.Rational(1,4)],
            [sp.Rational(1,4), sp.Rational(5,12)]])
ebT = sp.Matrix([[sp.Rational(1,4), 0], [0, sp.Rational(3,4)]])

# Calculate R(z)
def P(z):
    return (sp.eye(2) - z * A + z * ebT).det()


def Q(z):
    return (sp.eye(2) - z * A).det()


z = sp.symbols('z')
sp.pprint(P(z) / Q(z))

% Define IRK method
A = [1/4, -1/4 ; 1/4, 5/12];
ebT = [1/4, 0 ; 0, 3/4];

% Calculate R(z)
P = @(z) det(eye(2) - z * A + z * ebT);
Q = @(z) det(eye(2) - z * A);

syms z
Rz = P(z) / Q(z)