Deriving implicit Runge-Kutta methods

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 rooted 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].

3.1.4. Determining the order of an implicit Runge-Kutta method#

One of the differences between implicit and explicit methods is that an implicit method can achieve the same accuracy as an explicit method but using fewer stages. To determine the order of an implicit method we need to consider the order conditions that govern the values of the coefficients \(a_{ij}\), \(b_i\) and \(c_i\)

Let \(G(k)\) represent the fact that a given implicit method has order \(k\) then it can be shown that

\[B(k) \text{ and } C(\left\lfloor \tfrac{1}{2}k \right\rfloor) \text{ and } D (\left\lfloor \tfrac{1}{2}k \right\rfloor) \implies G(2s)\]

In other words, to determine the order of an IRK method we need find the highest value of \(k\) for which \(B(k)\) and \(C(\left\lfloor \frac{1}{2}k \right\rfloor)\) and \(D(\left\lfloor \frac{1}{2}k \right\rfloor)\) are satisfied. If no value of \(k\) satisfies these conditions then it is not a valid implicit Runge-Kutta method.

Example 3.3

Determine the order of the following IRK method

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

Solution

First check the highest value of \(k\) that satisfies \(B(k)\)

\[\begin{split} \begin{align*} k &= 1: & B(1) &= b_1 + b_2 + b_3 = \frac{1}{6} + \frac{2}{3} + \frac{1}{6} = 1, \\ k &= 2: & B(2) &= b_1c_1 + b_2c_2 + b_3c_3 = \frac{1}{6}(0) + \frac{2}{3}\left( \frac{1}{2} \right) + \frac{1}{6} (1) = \frac{1}{2}, \\ k &= 3: & B(3) &= b_1c_1^2 + b_2c_2^2 + b_3c_3^2 = \frac{1}{6}(0)^2 + \frac{2}{3} \left(\frac{1}{2}\right)^2 + \frac{1}{6}(1)^2 = \frac{1}{3}, \\ k &= 4: & B(4) &= b_1c_1^3 + b_2c_2^3 + b_3c_3^3 = \frac{1}{6}(0)^3 + \frac{2}{3} \left(\frac{1}{2}\right)^3 + \frac{1}{6}(1)^3 = \frac{1}{4}, \\ k &= 5: & B(5) &= b_1c_1^4 + b_2c_2^4 + b_3c_3^4 = \frac{1}{6}(0)^4 + \frac{2}{3} \left(\frac{1}{2}\right)^4 + \frac{1}{6}(1)^4 = \frac{5}{24}. \end{align*} \end{split}\]

So the highest value of \(k\) for which \(B(k)\) is satisfied is \(k=4\). We now need to check whether \(C(2)\) and \(D(2)\) are satisfied.

\[\begin{split} \begin{align*} \ell &= 1, & i &= 1: & LHS &= a_{11} + a_{12} + a_{13} = 0 + 0 + 0 = 0, \\ &&&& RHS &= \frac{1}{1}c_1^1 = 0, \\ \ell &= 1, & i &= 2: & LHS &= a_{21} + a_{22} + a_{23} = \frac{1}{4} + \frac{1}{4} + 0 = \frac{1}{2}, \\ &&&& RHS &= \frac{1}{1}c_2^1 = \frac{1}{2}, \\ \ell &= 1, & i &= 3: & LHS &= a_{31} + a_{32} + a_{33} = 0 + 1 + 0 = 1, \\ &&&& RHS &= \frac{1}{1}c_3^1 = 1, \\ \ell &= 2, & i &= 1: & LHS &= a_{11}c_1 + a_{12}c_2 + a_{13}c_3 = 0(0) + 0\left(\frac{1}{2}\right) + 0(1) = 0, \\ &&&& RHS &= \frac{1}{2}c_1^2 = 0, \\ \ell &= 2, & i &= 2: & LHS &= a_{21}c_1 + a_{22}c_2 + a_{23}c_3 = \frac{1}{4}(0) + \frac{1}{4}\left(\frac{1}{2}\right) + 0(1) = \frac{1}{8}, \\ &&&& RHS &= \frac{1}{2}c_2^2 = \frac{1}{2}\left(\frac{1}{2}\right)^2 = \frac{1}{8}, \\ \ell &= 2, & i &= 3: & LHS &= a_{31}c_1 + a_{32}c_2 + a_{33}c_3 = 0(0) + 1\left(\frac{1}{2}\right) + 0(1) = \frac{1}{2}, \\ &&&& RHS &= \frac{1}{2}c_3^3 = \frac{1}{2}(1)^3 = \frac{1}{2}. \end{align*} \end{split}\]

So the \(C(2)\) condition is satisfied. Checking the \(D(2)\) condition

\[\begin{split} \begin{align*} \ell &= 1, & j &= 1: & LHS &= b_1a_{11} + b_2a_{21} + b_3a_{31} = \frac{1}{6}(0) + \frac{2}{3}\left(\frac{1}{4}\right) + \frac{1}{6}(0) = \frac{1}{6}, \\ &&&& RHS &= \frac{1}{1}b_1(1 - c_1^1) = \frac{1}{6}(1 - 0) = \frac{1}{6}, \\ \ell &= 1, & j &= 2: & LHS &= b_1a_{12} + b_2a_{22} + b_3a_{32} = \frac{1}{6}(0) + \frac{2}{3}\left( \frac{1}{4} \right) + \frac{1}{6}(1) = \frac{1}{3}, \\ &&&& RHS &= \frac{1}{1}b_2(1 - c_2^1) = \frac{2}{3}\left(1 - \frac{1}{2}\right) = \frac{1}{3}, \\ \ell &= 1, & j &= 3: & LHS &= b_1a_{13} + b_2a_{23} + b_3a_{33} = \frac{1}{6}(0) + \frac{2}{3}\left(\frac{1}{4}\right) + \frac{1}{6}(0) = \frac{1}{6}, \\ &&&& RHS &= \frac{1}{1}b_3(1 - c_3^1) = \frac{1}{6}(1 - 0) = \frac{1}{6}, \\ \ell &= 2, & j &= 1: & LHS &= b_1c_1a_{11} + b_2c_2a_{21} + b_3c_3a_{31} = \frac{1}{6}(0)(0) + \frac{2}{3}\left(\frac{1}{2}\right)\left(\frac{1}{4}\right) + \frac{1}{6}(1)(0) = \frac{1}{12}, \\ &&&& RHS &= \frac{1}{2}b_1(1 - c_1^2) = \frac{1}{2}\left(\frac{1}{6}\right)(1 - 0) = \frac{1}{12}, \\ \ell &= 2, & j &= 2: & LHS &= b_1c_1a_{12} + b_2c_2a_{22} + b_3c_3a_{23} = \frac{1}{6}(0)(0) + \frac{2}{3}\left(\frac{1}{2}\right)\left(\frac{1}{4}\right) + \frac{1}{6}(1)(1) = \frac{1}{4}, \\ &&&& RHS &= \frac{1}{2}b_2(1 - c_2^2) = \frac{1}{2} \left(\frac{2}{3}\right) \left(1 - \left(\frac{1}{2}\right)^2\right) = \frac{1}{4}, \\ \ell &= 2, & j &= 3: & LHS &= b_1c_1a_{13} + b_2c_2a_{23} + b_3c_2a_{33} = \frac{1}{6}(0)(0) + \frac{2}{3}\left(\frac{1}{2}\right)(0) + \frac{1}{6}(1)(0) = 0, \\ &&&& RHS &= \frac{1}{2}b_3(1 - c_3^3) = \frac{1}{2}\left(\frac{1}{6}\right)(1 - 1^2) = 0. \end{align*} \end{split}\]

So the \(D(2)\) condition is also satisfied, therefore this method is a fourth-order method.