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.16)\[\Phi(\tau) = \dfrac{1}{\gamma(\tau)}.\]

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.

\(\tau\)
\(\Phi(\tau)\)	\(\displaystyle\sum_i b_i\)	\(\displaystyle\sum_i b_ic_i\)	\(\displaystyle\sum_i b_ic_i^2\)	\(\displaystyle\sum_{i,j} b_ia_{ij}c_j\)	\(\displaystyle\sum_i b_ic_i^3\)	\(\displaystyle\sum_{i,j} b_ia_{ij}c_ic_j\)	\(\displaystyle\sum_{i,j} b_ia_{ij}c_j^2\)	\(\displaystyle\sum_{i,j,k} b_ia_{ij}a_{jk}c_k\)
\(\gamma(\tau)\)	1	2	3	6	4	8	12	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. 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.2 (Order conditions for a fourth-order explicit Runge-Kutta method)

(2.17)\[\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}, \\ c_2 &= a_{21}, \\ c_3 &= a_{31} + a_{32}, \\ c_4 &= a_{41} + a_{42} + a_{43}. \end{align*} \end{split}\]

A fourth-order explicit Runge-Kutta method has 11 order conditions expressed in 14 unknowns (6 \(a_{ij}\) coefficients and 4 each of the \(b_i\) and \(c_i\) coefficients). To derive a method we need to select values for at least 3 of the unknowns and solve for the remaining.

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.

Python

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

MATLAB

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

Python

# 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))

MATLAB

% 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.3 (The fourth-order explicit Runge-Kutta method (RK4))

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