3.1. Deriving implicit Runge-Kutta methods#

Deriving Implicit Runge-Kutta (IRK) methods is done similarly to that for Explicit Runge-Kutta (ERK) methods where the values of the \(a_{ij}\), \(b_i\) and \(c_i\) coefficients are chosen to satisfy a set of order conditions. These order conditions were derived by Butcher [1964] based on the structure of routed trees using a similar method we used to derive the order conditions for ERK methods

Definition 3.1 (The \(B(k)\), \(C(k)\) and \(D(k)\) order conditions)

(3.1)#\[\begin{split} \begin{align*} B(k): && \sum_{i=1}^s b_i c_i^{j-1} = & \frac{1}{j}, & j&=1\ldots k, \\ C(k): && \sum_{j=1}^s a_{ij} c_j^{\ell-1} = & \frac{1}{\ell}c_i^{\ell} , & i&=1 \ldots s, & \ell &=1 \ldots k,\\ D(k): && \sum_{i=1}^s b_i c_i^{\ell-1} a_{ij} = & \frac{1}{\ell}b_j (1-c_j^{\ell}), & j&=1 \ldots s, & \ell &=1 \ldots k. \end{align*}\end{split}\]

3.1.1. Gauss-Legendre methods#

Gauss-Legendre methods are a family of methods that are derived using Gauss-Legendre quadrature. An \(s\)-stage Gauss-Legendre method has order \(2s\). They are derived using Legendre polynomials which are are a system of orthogonal polynomials.

For an order \(k\) Gauss-Legendre IRK method, the values of the \(c_i\) coefficients are chosen to be the roots of the Legendre polynomial

\[ P_k (x)=\sum_{i=0}^n \binom{n}{i} \binom{n+i}{i} \left( x - 1 \right)^i, \]

where \(\displaystyle{n \choose i} = \dfrac{n!}{i!(n - i)!}\) is the binomial coefficient. The values of \(a_{ij}\) and \(b_i\) are then chosen to satisfy the \(B(2k)\) and \(C(k)\) order conditions from equation (3.1).

Example 3.1

Derive a fourth-order Gauss-Legendre method.

Solution

A fourth-order Gauss-Legendre method requires 2 stages. The values of \(c_1\) and \(c_2\) are the roots of \(P_2(x) = 0\)

\[\begin{split} \begin{align*} P_2(x) &= \binom{2}{0} \binom{2}{0} (x - 1)^0 + \binom{2}{1}\binom{3}{1} (x - 1)^1 + \binom{2}{2} \binom{4}{2} (x - 1)^2 \\ &= 1 + 6(x - 1) + 6(x - 1)^2 \\ &= 6x^2 - 6x + 1, \end{align*} \end{split}\]

and the values of \(b_i\) and \(a_{ij}\) satisfy the \(B(4)\) and \(C(2)\) order conditions respectively which are

\[\begin{split} \begin{align*} b_1 + b_2 &= 1, \\ c_1b_1 + c_2b_2 &= \frac{1}{2}, \\ a_{11} + a_{12} &= c_1, \\ a_{21} + a_{22} &= c_2, \\ a_{11}c_1 + a_{12}c_2 &= \frac{1}{2}c_1^2, \\ a_{21}c_1 + a_{22}c_2 &= \frac{1}{2}c_1^2. \end{align*} \end{split}\]

Solving the order conditions using Python or MATLAB:

import sympy as sp

# Define symbolic variables
a11, a12, a21, a22, b1, b2, c1, c2, x = sp.symbols("a11, a12, a21, a22, b1, b2, c1, c2, x")

# Calculate c values 
c1, c2 = sp.solve(6 * x ** 2 - 6 * x + 1)
print(f"c1 = {c1}\nc2 = {c2}")

# Define order conditions
eq1 = b1 + b2 - 1
eq2 = b1 * c1 + b2 * c2 - sp.Rational(1,2)
eq3 = a11 + a12 - c1
eq4 = a21 + a22 - c2
eq5 = a11 * c1 + a12 * c2 - sp.Rational(1,2) * c1 ** 2
eq6 = a21 * c1 + a22 * c2 - sp.Rational(1,2) * c2 ** 2

# Solve order conditions
sp.solve((eq1, eq2, eq3, eq4, eq5, eq6))

syms a11 a12 a21 a22 b1 b2 x

% Calculate c coefficients
c = solve(6 * x ^ 2 - 6 * x + 1)
c1 = c(1);
c2 = c(2);

% Define order conditions
eq1 = b1 + b2 == 1;
eq2 = b1 * c1 + b2 * c2 == 1/2;
eq3 = a11 + a12 == c1;
eq4 = a21 + a22 == c2;
eq5 = a11 * c1 + a12 * c2 == 1/2 * c1^2;
eq6 = a21 * c1 + a22 * c2 == 1/2 * c2^2;

% Solve order conditions
solve(eq1, eq2, eq3, eq4, eq5, eq6)

Output

c1 = 1/2 - sqrt(3)/6
c2 = sqrt(3)/6 + 1/2

{b1: 1/2,
 b2: 1/2,
 a11: 1/4,
 a12: 1/4 - sqrt(3)/6,
 a21: 1/4 + sqrt(3)/6,
 a22: 1/4}

So the fourth-order Gauss-Legendre method is

\[\begin{split} \begin{align*} \begin{array}{c|cc} \frac{1}{2} - \frac{\sqrt{3}}{6} & \frac{1}{4} & \frac{1}{4} - \frac{\sqrt{3}}{6} \\ \frac{1}{2} + \frac{\sqrt{3}}{6} & \frac{1}{4} + \frac{\sqrt{3}}{6} & \frac{1}{4} \\ \hline & \frac{1}{2} & \frac{1}{2} \end{array} \end{align*} \end{split}\]

3.1.2. Radau methods#

Gauss-Legendre methods give us maximal order for the number of stages, however sometimes it is better the sacrifice order to gain better stability properties. An \(s\)-stage Radau method has order \(k=2s-1\) and is A-stable (see A-stability).

There are two types of Radau methods, Radau IA and Radau IIA:

  • Radau IA – the \(c_i\) coefficients are the roots of \(P_s(x) + P_{s-1}(x) = 0\) and the \(a_{ij}\) and \(b_i\) coefficients are chosen to satisfy \(B(2s)\) and \(D(s)\) order conditions respectively.

  • Radau IIA – the \(c_i\) coefficients are the roots of \(P_s(x) - P_{s-1}(x) = 0\) and the \(a_{ij}\) and \(b_i\) coefficients are chosen to satisfy \(B(2s)\) and \(C(s)\) order conditions respectively.

Example 3.2

Derive a third-order Radau IA method.

Solution

A third-order Radau IA method requires \(s=2\) stages. The values of \(c_i\) are the roots of \(0 = P_2(t) + P_1(t)\)

\[\begin{split} \begin{align*} 0 &= 1 + \binom{2}{1} \binom{3}{1} (x - 1) + \binom{2}{2} \binom{4}{2} (x - 1)^2 + 1 + \binom{1}{1} \binom{2}{1} (x - 1) \\ &= 2 + 6(x - 1) + 6(x - 1)^2 + 2(x - 1) \\ &= x(3x - 2) \end{align*} \end{split}\]

so \(c_1 = 0\) and \(c_2 = \frac{2}{3}\). The \(b_i\) and \(a_{ij}\) coefficients need to satisfy the \(B(3)\) and \(D(2)\) order conditions respectively which are

\[\begin{split} \begin{align*} b_1 + b_2 &= 1, \\ b_2c_2 &= \frac{1}{2}, \\ b_2c_2^2 &= \frac{1}{3}, \\ b_1a_{11} + b_2a_{21} &= b_1(1 - c_1), \\ b_1a_{12} + b_2a_{22} &= b_2(1 - c_2), \\ b_1c_1a_{11} + b_2c_2a_{21} &= \frac{1}{2}b_1(1 - c_1^2), \\ b_1c_1a_{12} + b_2c_2a_{22} &= \frac{1}{2}b_2(1 - c_2^2). \end{align*} \end{split}\]

Solving the order conditions using Python or MATLAB

import sympy as sp

# Define symbolic variables
a11, a12, a21, a22, b1, b2, c1, c2 = sp.symbols("a11, a12, a21, a22, b1, b2, c1, c2")
c1, c2 = 0, sp.Rational(2, 3)

# Define order conditions
eq1 = b1 + b2 - 1
eq2 = b1 * c1 + b2 * c2 - sp.Rational(1,2)
eq3 = b1 * c1 ** 2 + b2 * c2 ** 2 - sp.Rational(1,3)
eq4 = b1 * a11 + b2 * a21 - b1 * (1 - c1)
eq5 = b1 * a12 + b2 * a22 - b2 * (1 - c2)
eq6 = b1 * c1 * a11 + b2 * c2 * a21 - sp.Rational(1,2) * b1 * (1 - c1 ** 2)
eq7 = b1 * c1 * a12 + b2 * c2 * a22 - sp.Rational(1,2) * b2 * (1 - c2 ** 2)

# Solve order conditions
sp.solve((eq1, eq2, eq3, eq4, eq5, eq6, eq7))

syms a11 a12 a21 a22 b1 b2 c1 c2
c1 = 0;
c2 = 2/3;

% Define order conditions
eq1 = b1 + b2 == 1;
eq2 = b2 * c2 == 1/2;
eq3 = b2 * c2 ^ 2 == 1/3;
eq4 = b1 * a11 + b2 * a21 == b1 * (1 - c1);
eq5 = b1 * a12 + b2 * a22 == b2 * (1 - c2);
eq6 = b1 * c1 * a11 + b2 * c2 * a21 == 1/2 * b1 * (1 - c1^2);
eq7 = b1 * c1 * a12 + b2 * c2 * a22 == 1/2 * b2 * (1 - c2^2);

% Solve order conditions
solve(eq1, eq2, eq3, eq4, eq5, eq6)

Output:

c1 = 0
c2 = 2/3

[{a11: 1/4, a12: -1/4, a21: 1/4, a22: 5/12, b1: 1/4, b2: 3/4}]

So the third-order Radau IA method is

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

3.1.3. DIRK methods#

Diagonally Implicit Runge-Kutta (DIRK) methods have an \(A\) matrix that is lower triangular with equal elements on the main diagonal.

\[\begin{split} \begin{align*} \begin{array}{c|ccccc} c_1 & \lambda \\ c_2 & a_{21} & \lambda \\ \vdots & \vdots & \ddots & \ddots \\ c_s & a_{s1} & \cdots & a_{s,s-1} & \lambda \\ \hline & b_1 & \cdots & b_{s-1} & b_s \end{array} \end{align*} \end{split}\]

Writing out the stage values for a DIRK method

\[\begin{split} \begin{align*} k_1 &= f(t_n + \lambda h, y_n + h\lambda k_1), \\ k_2 &= f(t_n + c_2h, y_n + h(a_{21}k_1 + \lambda k_2)), \\ & \vdots \\ k_s &= f(t_n + c_sh, y_n + h(a_{s1}k_1 + \cdots + a_{s,s-1}k_{s-1} + \lambda k_s)). \end{align*} \end{split}\]

We see that whilst these are still implicit equations, there is only one unknown in each one, i.e., the equation for \(k_i\) does not include \(k_{i+1} \dots k_s\). So the advantage of DIRK methods is that the solutions to these can be obtained sequentially and is more computationally efficient than non-DIRK methods.

Definition 3.2 (2-stage DIRK method)

\[\begin{split} \begin{array}{c|cc} \frac{1}{4} & \frac{1}{4} \\ c_2 & c_2 - \frac{1}{4} & \frac{1}{4} \\ \hline & \frac{2(2c_2 - 1)}{4c_2 - 1} & \frac{1}{4c_2 - 1} \end{array} \end{split}\]

[Jackiewicz and Tracogna, 1995]

For example, when \(c_2 = \frac{3}{4}\) the 2-stage DIRK method given above is

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

Definition 3.3 (3-stage DIRK method)

\[\begin{split} \begin{array}{c|ccc} c_1 & c_1 \\ \frac{1}{2}(1 + c_1) & \frac{1}{2}(1 - c_1) & c_1 \\ 1 & \frac{1}{4}(-6c_1^2 + 16c_1 - 1) & \frac{1}{4}(6c_1^2 - 20c_1 + 5) & c_1 \\ \hline & \frac{1}{4}(-6c_1^2 + 16c_1 -1) & \frac{1}{4}(6c_1^2 - 20c_1 + 5) & c_1 \end{array} \end{split}\]

This DIRK method is A-stable where \(c_1\) satisfies \(0 = \frac{1}{6} - \frac{3}{2} c_1 + 3c_1^2 - c_1^3\) [Butcher, 2008].