3.2. Deriving implicit Runge-Kutta methods

There are several types of implicit Runge-Kutta methods which are derived using the \(B(k)\), \(C(k)\) and \(D(k)\) order conditions.

\[\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.2.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.

Definition 3.2 (Legendre polynomials)

The \(n\)-th order Legendre polynomial is

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

where \(\displaystyle{n \choose k} = \dfrac{n!}{k!(n - k)!}\) is the binomial coefficient.

Definition 3.3 (Gauss-Legendre method)

A Gauss-Legendre implicit Runge-Kutta method of order \(2s\) has \(c\) coefficients that are the roots of the Legendre polynomial \(P_s(x)=0\) and the \(a_{ij}\) and \(b\) are chosen to satisfy the \(B(2s)\) and \(C(s)\) order conditions.

Example 3.2

Derive a fourth-order Gauss-Legendre method.

Solution (click to show)

A 4th order Gauss-Legendre method requires \(s = 2\) stages. The values of \(c_1\) and \(c_2\) are the roots of \(P_2(x) = 0\) 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:

Python

import sympy as sp

# Define Legendre polynomial
def P(n):
    Pn, x = 0, sp.symbols('x')
    for k in range(n + 1):
        Pn += sp.binomial(n, k) * sp.binomial(n + k, k) * (x - 1) ** k

    return Pn

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

# Calculate c values
c1, c2 = sp.solve(P(2))
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))

MATLAB

syms a11 a12 a21 a22 b1 b2

% Calculate c coefficients
c = solve(P(2))
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)

% -------------------------------------------------------------------
function Pn = P(n)

syms x
Pn = 0;
for k = 0 : n
    Pn = Pn + nchoosek(n, k) * nchoosek(n + k, k) * (x - 1) ^ k;
end

end

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.2.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.

Definition 3.4 (Radau IA methods)

A Radau IA method of order \(2s-1\) has \(c\) coefficients where \(c_1=0\) and the remaining \(c\) 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.

Definition 3.5 (Radau IIA methods)

A Radau IIA method of order \(2s-1\) has \(c\) coefficients where \(c_s=1\) and the remaining \(c\) 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.3

Derive a third-order Radau IA method.

Solution (click to show)

A third-order Radau IA method will have \(s=2\) stages. We know that \(c_1=0\) and the value of \(c_2\) is the non-zero root of \(0 = P_2(t) + P_1(t)\). 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

Python

import sympy as sp

# Define Legendre polynomial
def P(n):
    Pn, x = 0, sp.symbols('x')
    for k in range(n + 1):
        Pn += sp.binomial(n, k) * sp.binomial(n + k, k) * (x - 1) ** k

    return Pn


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

# Calculate c coefficients
c1, c2 = sp.solve(P(2) + P(1))
print(f"c1 = {c1}\nc2 = {c2}")

# Define order conditions
eq1 = b1 + b2 - 1
eq2 = b2 * c2 - sp.Rational(1,2)
eq3 = 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))

MATLAB

syms a11 a12 a21 a22 b1 b2 c1 c2

% Calculate c coefficients
c = solve(P(2) + P(1))
c1 = c(1);
c2 = c(2);

% 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, eq7)

% -------------------------------------------------------------------
function Pn = P(n)

syms x
Pn = 0;
for k = 0 : n
    Pn = Pn + nchoosek(n, k) * nchoosek(n + k, k) * (x - 1) ^ k;
end

end

Output:

c1 = 0
c2 = 2/3

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

3.2.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.

The Butcher tableau for a 2-stage DIRK method is

\[\begin{split} \begin{array}{c|cc} \lambda & \lambda \\ c_2 & c_2 - \lambda & \lambda \\ \hline & b_1 & b_2 \end{array} \end{split}\]

The \(B(k)\) order conditions are

\[\begin{split} \begin{align*} b_1 + b_2 &= 1, \\ \lambda b_1 + b_2c_2 &= \frac{1}{2}, \end{align*} \end{split}\]

which has the solutions \(b_1 = \dfrac{2c_2 - 1}{2(c_2 - \lambda)}\) and \(b_2 = \dfrac{1 - 2\lambda}{2(c_2 - \lambda)}\).

Definition 3.6 (2-stage DIRK method)

A 2-stage DIRK method has the Butcher tableau

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

for some choice of \(\lambda\).

Example 3.4

Derive a 2-stage DIRK method with \(\lambda = 1 + \frac{1}{2}\sqrt{2}\) and \(c_2 = 1\).

Solution (click to show)

The values of \(b_1\) and \(b_2\) are

\[\begin{split} \begin{align*} b_1 &= \frac{2(1) - 1}{2(1 - 1 + \frac{\sqrt{2}}{2})} = \frac{1}{\sqrt{2}} = \frac{\sqrt{2}}{2}, \\ b_2 &= 1 - b_1 = 1 - \frac{\sqrt{2}}{2} \end{align*} \end{split}\]

so this DIRK method is

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