2.5. Derivation of a fourth-order explicit Runge-Kutta method#

We saw in the previous section that we derive a set of order conditions by listing all of the rooted trees up to a certain order and using the following relation between the elementary weight \(\Phi\) and density \(\gamma\) for each tree

(2.15)#\[\Phi = \dfrac{1}{\gamma}.\]

To derive the order conditions for a fourth-order explicit method we consider all of the rooted trees up to and including order 4 which are tabulated below along with their elementary weights and densities.

Tree

\(\Phi\)

\(\gamma\)

Tree

\(\Phi\)

\(\gamma\)

\(\displaystyle\sum_i b_i\)

1

\(\displaystyle\sum_i b_ic_i\)

2

\(\displaystyle\sum_i b_ic_i^2\)

3

\(\displaystyle\sum_i b_i c_i^3\)

4

\(\displaystyle\sum_{i,j} b_i a_{ij} c_j\)

6

\(\displaystyle\sum_{i,j} b_i a_{ij} c_ic_j\)

8

\(\displaystyle\sum_{i,j} b_ia_{ij}c_j^3\)

12

\(\displaystyle\sum_{i,j,k} b_i a_{ij} a_{jk} c_k\)

24

The minimum number of stages for explicit Runge-Kutta method of different orders are shown in the table below [Hairer and Wanner, 1993].

Order

1

2

3

4

5

6

7

Stages

1

2

3

4

6

7

9

So the fewest number of stages required for a fourth-order explicit Runge-Kutta method is \(s=4\), therefore the order conditions are

\[\begin{split} \begin{align*} \sum_{i=1}^4 b_i &= 1, \\ \sum_{i=1}^4 b_ic_i &= \frac{1}{2}, \\ \sum_{i=1}^4 b_ic_i^2 &= \frac{1}{3}, \\ \sum_{i=1}^4 b_ic_i^3 &= \frac{1}{4}, \\ \sum_{i,j=1}^4 b_ia_{ij}c_j &= \frac{1}{6}, \\ \sum_{i,j=1}^4 b_ia_{ij}c_ic_j &= \frac{1}{8}, \\ \sum_{i,j=1}^4 b_ia_{ij}c_j^2 &= \frac{1}{12}, \\ \sum_{i,j,k=1}^4 b_ia_{ij}a_{jk} c_k &= \frac{1}{24}. \end{align*} \end{split}\]

Since we are deriving the order conditions for an explicit method we know that \(c_1=0\) and \(a_{ij} = 0\) where \(i \leq j\) and we have the Butcher tableau

\[\begin{split} \begin{align*} \begin{array}{c|cccc} 0 & 0 \\ c_2 & a_{21} \\ c_3 & a_{31} & a_{32} \\ c_4 & a_{41} & a_{42} & a_{43} \\ \hline & b_1 & b_2 & b_3 & b_4 \end{array} \end{align*} \end{split}\]

So we can omit any terms in the summations where the values of the \(c_1\) and \(a_{ij}\) coefficients are zero. For the first 4 order conditions we have

\[\begin{split} \begin{align} b_1 + b_2 + b_3 + b_4 &= 1, \\ b_2c_2 + b_3c_3 + b_4c_4 &= \frac{1}{2}, \\ b_2c_2^2 + b_3c_3^2 + b_4c_4^2 &= \frac{1}{3}, \\ b_2c_2^3 + b_3c_3^3 + b_4c_4^3 &= \frac{1}{4}, \\ \end{align} \end{split}\]

For the next 3 order conditions, the \(a_{ij}c_j\) term is only non-zero when \(i = 3, 4\) and \(j = 2, 3\) so we have

\[\begin{split} \begin{align*} b_3a_{32}c_2 + b_4a_{42}c_2 + b_4a_{43}c_3 &= \frac{1}{6}, \\ b_3a_{32}c_2c_3 + b_4a_{42}c_2c_4 + b_4a_{43}c_3c_4 &= \frac{1}{8}, \\ b_3a_{32}c_2^2 + b_4a_{42}c_2^2 + b_4a_{43}c_3^2 &= \frac{1}{12}. \end{align*} \end{split}\]

For the final order condition, the \(a_{ij}a_{jk}c_j\) term is only non-zero when \(i=4\), \(j=3\) and \(k=2\) so we have

\[ \begin{align*} b_4a_{43}a_{32}c_2 &= \frac{1}{24}. \end{align*} \]

We have now derived the order conditions for a fourth-order explicit Runge-Kutta method which are combined with the row sum condition equation (2.14). Using rooted trees to do this is much easier than expanding out the Taylor series for the ODE \(y'=f(t, y)\) and the general form of a 4-stage Runge-Kutta method which would require several pages of complicated calculus and algebra.

Definition 2.6 (Order conditions for a fourth-order explicit Runge-Kutta method)

\[\begin{split} \begin{align*} b_1 + b_2 + b_3 + b_4 &= 1, \\ b_2c_2 + b_3c_3 + b_4c_4 &= \frac{1}{2}, \\ b_2c_2^2 + b_3c_3^2 + b_4c_4^2 &= \frac{1}{3}, \\ b_2c_2^3 + b_3c_3^3 + b_4c_4^3 &= \frac{1}{4}, \\ b_3a_{32}c_2 + b_4a_{42}c_2 + b_4a_{43}c_3 &= \frac{1}{6}, \\ b_3a_{32}c_2c_3 + b_4a_{42}c_2c_4 + b_4a_{43}c_3c_4 &= \frac{1}{8}, \\ b_3a_{32}c_2^2 + b_4a_{42}c_2^2 + b_4a_{43}c_3^2 &= \frac{1}{12}, \\ b_4a_{43}a_{32}c_2 &= \frac{1}{24}. \end{align*} \end{split}\]

A fourth-order explicit Runge-Kutta method has 8 order conditions expressed in 14 unknowns (6 \(a_{ij}\) coefficients and 4 each of the \(b_i\) and \(c_i\) coefficients). The row sum condition gives a further 3 order conditions so we have 11 order conditions and 14 unknowns. To derive a method we choose the value of 3 of the unknowns and solve for the rest.

Example 2.6

Derive a fourth-order Runge-Kutta method where \(c_2 = c_3 = \frac{1}{2}\) and \(c_4 = 1\).

Solution (click to show)

Use Python or MATLAB to solve the order conditions and the row sum conditions. First lets solve the first four order conditions for the \(b_i\) coefficients.

import sympy as sp

# Declare symbolic variables
a21, a31, a32, a41, a42, a43, b1, b2, b3, b4 = sp.symbols("a21, a31, a32, a41, a42, a43, b1, b2, b3, b4")
c2, c3, c4 = sp.Rational(1,2), sp.Rational(1,2), 1

# Define order conditions
eq1 = b1 + b2 + b3 + b4 - 1
eq2 = b2 * c2 + b3 * c3 + b4 * c4 - sp.Rational(1,2)
eq3 = b2 * c2 ** 2 + b3 * c3 ** 2 + b4 * c4 ** 2 - sp.Rational(1,3)
eq4 = b2 * c2 ** 3 + b3 * c3 ** 3 + b4 * c4 ** 3 - sp.Rational(1,4)

# Solve for the b coefficients
sp.solve((eq1, eq2, eq3, eq4))

% Define symbolic variables
syms a21 a31 a32 a41 a42 a43 b1 b2 b3 b4

% Set c values
c2 = 1/2;
c3 = 1/2;
c4 = 1;

% Define order conditions
eq1 = b1 + b2 + b3 + b4 == 1;
eq2 = b2 * c2 + b3 * c3 + b4 * c4 == 1/2;
eq3 = b2 * c2 ^ 2 + b3 * c3^2 + b4 * c4^2 == 1/3;
eq4 = b2 * c2^3 + b3 * c3^3 + b4 * c4^3 == 1/4;

% Solve for the b coefficients
solve(eq1, eq2, eq3, eq4)

Output:

{b1: 1/6, b2: 2/3 - b3, b4: 1/6}

So \(b_1 = b_4 = \frac{1}{6}\) and \(b_2 = \frac{2}{3} - b_3\), which means we need to choose a value for \(b_3\) (or \(b_2\)). Lets choose \(b_3 = \frac{1}{3}\) so that \(b_2\) is also \(\frac{1}{3}\) (this gives a nice symmetry to \(b_2\) and \(b_3\) but also simplifies the other coefficients). We can now solve for the remaining order conditions to get the \(a_{ij}\) coefficients.

# Set b values
b1, b2, b3, b4 = sp.Rational(1, 6), sp.Rational(1, 3), sp.Rational(1,3), sp.Rational(1, 6)

# Define order conditions
eq5 = b3 * a32 * c2 + b4 * a42 * c2 + b4 * a43 * c3 - sp.Rational(1,6)
eq6 = b3 * a32 * c2 * c3 + b4 * a42 * c2 * c4 + b4 * a43 * c3 * c4 - sp.Rational(1,8)
eq7 = b3 * a32 * c2 ** 2 + b4 * a42 * c2 ** 2 + b4 * a43 * c3 ** 2 - sp.Rational(1,12)
eq8 = b4 * a43 * a32 * c2 - sp.Rational(1,24)
eq9 = a21 - c2
eq10 = a31 + a32 - c3
eq11 = a41 + a42 + a43 - c4

# Solve order conditions
sp.solve((eq5, eq6, eq7, eq8, eq9, eq10, eq11))

% Set b values
b1 = 1/6;
b2 = 1/3;
b3 = 1/3;
b4 = 1/6;

% Define order conditions
eq5 = b2 * a32 * c2 + b4 * a42 * c2 + b4 * a43 * c3 == 1/6;
eq6 = b3 * a32 * c2 * c3 + b4 * a42 * c2 * c4 + b4 * a43 * c3 * c4== 1/8;
eq7 = b3 * a32 * c2^2 + b4 * a42 * c2^2 + b4 * a43 * c3^2 == 1/12;
eq8 = b4 * a43 * a32 * c2 == 1/24;

% Define row-sum conditions
eq9 = a21 == c2;
eq10 = a31 + a32 == c3;
eq11 = a41 + a42 + a43 == c4;

% Solve order conditions
solve(eq5, eq6, eq7, eq8, eq9, eq10, eq11)

Output:

{a31: 0, a32: 1/2, a41: 0, a42: 0, a43: 1, a21: 1/2}

So \(a_{31} = a_{41} = a_{42} = 0\), \(a_{21} = a_{32} = \frac{1}{2}\) and \(a_{43} = 1\). The Butcher tableau for this fourth-order Runge-Kutta method is

\[\begin{split} \begin{align*} \begin{array}{c|cccc} 0 & & & & \\ 1/2 & 1/2 & & & \\ 1/2 & 0 & 1/2 & & \\ 1 & 0 & 0 & 1 & \\ \hline & 1/6 & 1/3 & 1/3 & 1/6 \end{array} \end{align*} \end{split}\]

or alternatively

\[\begin{split} \begin{align*} y_{n+1} &=y_n +\frac{h}{6}(k_1 +2 k_2 +2k_3 +k_4 ),\\ k_1 &=f(t_n ,y_n ),\\ k_2 &=f(t_n +\tfrac{1}{2}h,y_n +\tfrac{1}{2}h k_1 ),\\ k_3 &=f(t_n +\tfrac{1}{2}h,y_n +\tfrac{1}{2}h k_2 ),\\ k_4 &=f(t_n +h,y_n +hk_3 ), \end{align*} \end{split}\]

This fourth-order explicit Runge-Kutta method we derived in Example 2.6 is often referred to as the Runge-Kutta method or RK4 for short. It is simple, robust and accurate and has become the standard goto method for solving ODEs.

Definition 2.7 (The fourth-order explicit Runge-Kutta method (RK4))

(2.16)#\[\begin{split} \begin{align*} y_{n+1} &= y_n + \frac{h}{6}(k_1 + 2 k_2 + 2 k_3 + k_4), \\ k_1 &= f(t_n, y_n), \\ k_2 &= f(t_n + \tfrac{1}{2} h, y_n + \tfrac{1}{2} h k_1), \\ k_3 &= f(t_n + \tfrac{1}{2} h, y_n + \tfrac{1}{2} h k_2), \\ k_4 &= f(t_n + h, y_n + h k_3). \end{align*} \end{split}\]
# Set b values
b1, b2, b3, b4 = sp.Rational(1, 6), sp.Rational(1, 3), sp.Rational(1,3), sp.Rational(1, 6)

# Define order conditions
eq5 = b3 * a32 * c2 + b4 * a42 * c2 + b4 * a43 * c3 - sp.Rational(1,6)
eq6 = b3 * a32 * c2 * c3 + b4 * a42 * c2 * c4 + b4 * a43 * c3 * c4 - sp.Rational(1,8)
eq7 = b3 * a32 * c2 ** 2 + b4 * a42 * c2 ** 2 + b4 * a43 * c3 ** 2 - sp.Rational(1,12)
eq8 = b4 * a43 * a32 * c2 - sp.Rational(1,24)
eq9 = a21 - c2
eq10 = a31 + a32 - c3
eq11 = a41 + a42 + a43 - c4

# Solve order conditions
sp.solve((eq5, eq6, eq7, eq8, eq9, eq10, eq11))

{a31: 0, a32: 1/2, a41: 0, a42: 0, a43: 1, a21: 1/2}