Stability Functions for Implicit Runge-Kutta Methods

Contents

4.3. Stability Functions for Implicit Runge-Kutta Methods#

In the previous section, we derived the stability function for a general Runge-Kutta method solution to the text equation \(y' = \lambda y\) which is

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

where \(z = h\lambda\) and \(\mathbf{e} = (1, \ldots, 1)^\mathsf{T}\).

For implicit methods, \(A\) is not lower triangular so we cannot use a power series summation to determine the stability function. Instead, we can use a result from linear algebra to to express the stability function as a quotient of determinants. The matrix determinant lemma is

\[ \det(M + \mathbf{uv}^{\mathsf{T}}) = (1 + \mathbf{v}^\mathsf{T}M^{-1}\mathbf{u}) \det(M), \]

where \(M\) is an invertible square matrix and \(\mathbf{u}, \mathbf{v}\) are column vectors. Let \(M = I - zA\), \(\mathbf{u} = z \mathbf{e}\) and \(\mathbf{v} = \mathbf{b}\) then

\[ \det(I - zA + z \mathbf{eb}^\mathsf{T}) = (1 + z \mathbf{b}^\mathsf{T}(I - zA)^{-1} \mathbf{e}) \det(I - zA), \]

therefore

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

This is the stability function for an implicit Runge-Kutta method.

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

The stability function of an implicit Runge-Kutta method is

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

Note that \(\mathbf{eb}^\mathsf{T}\) is the outer product of the column vector \(\mathbf{e} = (1, \ldots, 1)^\mathsf{T}\) and the row vector \(\mathbf{b}^\mathsf{T} = (b_1, b_2, \ldots, b_s)\)

\[\begin{split} \mathbf{eb}^\mathsf{T} = \begin{pmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{pmatrix} \begin{pmatrix} b_1 & b_2 & \cdots & b_s \end{pmatrix} = \begin{pmatrix} b_1 & b_2 & \cdots & b_s \\ b_1 & b_2 & \cdots & b_s \\ \vdots & \vdots & \ddots & \vdots \\ b_1 & b_2 & \cdots & b_s \end{pmatrix}. \end{split}\]

Since each entry of \(I - zA\) is a polynomial in \(z\), both the numerator and denominator are polynomials. Therefore the stability function of an implicit Runge-Kutta method is a rational function

\[R(z) = \frac{P(z)}{Q(z)}. \]

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.8)

\[\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} & \frac{3}{4} \\ \frac{1}{4} & \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 & z \\ 0 & 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}{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)]])
b = sp.Matrix([[sp.Rational(1,4)], 
               [sp.Rational(3, 4)]])
s = len(b)
e = sp.ones(s, 1)

# Define P(z) and Q(z) functions
def P(z):
    return (sp.eye(s) - z * A + z * e * b.T).det()


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


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

% Define IRK method
A = [1/4, -1/4 ; 1/4, 5/12];
b = [1/4 ; 3/4];
s = length(b);
e = ones(s, 1);

% Define P(z) and Q(z) functions
P = @(z) det(eye(s) - z * A + z * e * b');
Q = @(z) det(eye(s) - z * A);

% Determine R(z)
syms z
Rz = P(z) / Q(z)
Contents