Derivation of a fourth-order explicit Runge-Kutta method

Contents

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

In the previous section we derived the rooted-tree order conditions for an explicit Runge-Kutta method. We now use these conditions to construct a fourth-order method. Unlike the derivation of RK2, which could be done using Taylor series expansions directly, the rooted-tree approach allows us to derive higher-order methods systematically without lengthy algebraic calculations.

For a fourth-order method we require all rooted trees of order up to 4. There are eight such trees. For each tree we determine its elementary weight \(\Phi(\tau)\), density \(\gamma(\tau)\), and then apply the condition

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

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

\(r(\tau)\)

1

2

3

3

4

4

4

4

\(\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

1st

2nd

3rd

4th

5th

6th

7th

Minimum stages required

1

2

3

4

6

7

9

Therefore we seek a 4-stage ERK method, so \(s=4\) and the order conditions are

(2.15)#\[\begin{split} \begin{align*} \sum_{i=1}^4 b_i &= 1, \\ \sum_{i=1}^4 b_ic_i &= \frac12. \\ \sum_{i=1}^4 b_ic_i^2 &= \frac13, \\ \sum_{i,j=1}^4 b_ia_{ij}c_j &= \frac{1}{6}, \\ \sum_{i,j=1}^4 b_ia_{ij}c_j &= \frac{1}{6}, \\ \sum_{i,j=1}^4 b_ia_{ij}c_ic_j &= \frac18, \\ \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 (i.e., the upper triangular elements of \(A\)). For example

\[\begin{split} \begin{align*} \sum_{i,j=1}^4 b_ia_{ij}c_j &= b_1a_{11}c_1 + b_1a_{12}c_2 + b_1a_{13}c_3 + b_1a_{14}c_4 \\ & \qquad b_2a_{21}c_1 + b_2a_{22}c_2 + b_2a_{23}c_3 + b_2a_{24}c_4 \\ & \qquad b_3a_{31}c_1 + b_3a_{32}c_2 + b_3a_{33}c_3 + b_3a_{34}c_4 \\ & \qquad b_4a_{41}c_1 + b_4a_{42}c_2 + b_4a_{43}c_3 + b_4a_{44}c_4, \end{align*} \end{split}\]

since for an ERK method, \(a_{ij} = 0\) for \(i \le j\) and \(c_1 = 0\), then this simplifies to

\[ \sum_{i,j=1}^4 = b_3a_{32}c_2 + b_4a_{42}c_2 + b_4a_{43}c_3. \]

So the first four order conditions from equation (2.15) simplify to

\[\begin{split} \begin{align*} b_1 + b_2 + b_3 + b_4 &= 1, \\ b_2c_2 + b_3c_3 + b_4c_4 &= \frac12, \\ b_2c_2^2 + b_3c_3^2 + b_4c_4^2 &= \frac13, \\ b_2c_2^3 + b_3c_3^3 + b_4c_4^3 &= \frac14, \\ \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 &= \frac18, \\ 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 from equation (2.15), 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 conditions. 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.1 (Order conditions for a fourth-order explicit Runge-Kutta method)

(2.16)#\[\begin{split} \begin{align*} b_1 + b_2 + b_3 + b_4 &= 1, \\ b_2c_2 + b_3c_3 + b_4c_4 &= \frac12, \\ b_2c_2^2 + b_3c_3^2 + b_4c_4^2 &= \frac13, \\ b_2c_2^3 + b_3c_3^3 + b_4c_4^3 &= \frac14, \\ 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 &= \frac18, \\ 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}\]

The method contains 14 coefficients, but the row-sum conditions relate the \(c_i\) values to the \(a_{ij}\) coefficients. Consequently the system possess three degrees of freedom, allowing infinitely many fourth-order ERK methods. So to derive a method we need to select values for at least three of the unknown coefficients and solve for the remaining ones.

Example 2.6

Derive a fourth-order Runge-Kutta method where \(c_2 = c_3 = \frac12\), \(c_4 = 1\) and \(b_2 = b_3\).

Solution

Substituting \(c_2 = c_3 = \frac12\) and \(c_4 = 1\) into the first four order conditions gives

\[\begin{split} \begin{align*} b_1 + b_2 + b_3 + b_4 &= 1, \\ \frac12b_2 + \frac12b_3 + b_4 &= \frac12, \\ \frac14b_2 + \frac14b_3 + b_4 &= \frac13, \\ \frac18b_2 + \frac18b_3 + b_4 &= \frac14 \end{align*} \end{split}\]

Subtracting \(\frac12\) times the third equation from the second

\[ \frac{1}{2}b_4 = \frac{1}{12}, \]

so \(b_4 = \frac16\) and the second equation is

\[ \frac12 b_2 + \frac12 b_3 = \frac13. \]

Since in this case we want \(b_2\) and \(b_3\) to be the same we have \(b_2 = b_3 = \frac13\). Substituting \(b_2\), \(b_3\) and \(b_4\) into the first equation gives \(b_1 = \frac16\).

Substituting the values of \(b_i\) and \(c_i\) into the fifth, sixth and seventh order conditions gives

(2.17)#\[\begin{split} \begin{align*} \frac{1}{6}a_{32} + \frac{1}{12} a_{42} + \frac{1}{12} a_{43} &= \frac{1}{6}, \\ \frac{1}{12} a_{32} + \frac{1}{12} a_{42} + \frac{1}{12} a_{43} &= \frac{1}{8}, \\ \frac{1}{12} a_{32} + \frac{1}{24} a_{42} + \frac{1}{24} a_{43} &= \frac{1}{12}. \end{align*} \end{split}\]

Subtracting \(\frac{1}{2}\) times the first equation from the second

\[ \frac{1}{24} a_{42} + \frac{1}{24} a_{43} = \frac{1}{24}, \]

so \(a_{42} = 1 - a_{43}\). Substituting \(b_4\) and \(c_2\) into the eighth order condition gives

\[ \frac{1}{12}a_{43}a_{32} = \frac{1}{24},\]

so \(a_{32} = \frac{1}{2a_{43}}\). Substituting \(a_{32}\) and \(a_{42}\) into the first equation from (2.17) gives

\[\begin{split} \begin{align*} \frac{1}{6} \left( \frac{1}{2a_{43}} \right) + \frac{1}{12}(1 - a_{43}) + \frac{1}{12}a_{43} &= \frac{1}{6} \\ 1 + a_{43} - a_{43}^2 + a_{43}^2 &= 2a_{43} \\ a_{43} &= 1. \end{align*} \end{split}\]

So \(a_{43} = 1\) and \(a_{32} = \frac{1}{2}\). Entering in the known values of \(a_{ij}\), \(b_i\) and \(c_i\) into the Butcher tableau we have

\[\begin{split} \begin{align*} \begin{array}{c|cccc} 0 & & & & \\ \frac12 & a_{21} & & & \\ \frac12 & a_{31} & \frac12 & & \\ 1 & a_{41} & 0 & 1 & \\ \hline & \frac{1}{6} & \frac13 & \frac13 & \frac{1}{6} \end{array} \end{align*} \end{split}\]

The remaining values of \(a_{21}\), \(a_{31}\) and \(a_{41}\) are found using the row-sum condition

\[ c_i = \sum_{j = 1}^{s} a_{ij}, \]

so \(a_{21} = \frac{1}{2}\), \(a_{31} = 0\) and \(a_{41} = 0\) and the Butcher tableau is

\[\begin{split} \begin{align*} \begin{array}{c|cccc} 0 & & & & \\ \frac12 & \frac12 & & & \\ \frac12 & 0 & \frac12 & & \\ 1 & 0 & 0 & 1 & \\ \hline & \frac{1}{6} & \frac13 & \frac13 & \frac{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 is often referred to as the Runge-Kutta method or RK4 for short. It is simple, robust and accurate and has become the standard go to method for solving ODEs.

2.5.1. Code#

We saw in the example derivation of the RK4 method that solving the order conditions can be quite an involved process. In practice, we utilise mathematical software to do the algebra for us.

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

# Choose b2 = b3 (this is an additional condition for the RK4 method)
eq12 = b2 - b3

# Solve the order conditions
sp.solve((eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12))

Output

{a31: 0,
 a32: 1/2,
 a41: 0,
 a42: 0,
 a43: 1,
 b1: 1/6,
 b2: 1/3,
 b3: 1/3,
 b4: 1/6,
 a21: 1/2}

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

% Define known coefficients
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;
eq5 = b3 * 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;
eq9 = a21 == c2;
eq10 = a31 + a32 == c3;
eq11 = a41 + a42 + a43 == c4;

% Solve the order conditions
solve(eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11)

Output:

ans = struct with fields:
    a21: 1/2
    a31: 0
    a32: 1/2
    a41: 0
    a42: 0
    a43: 1
     b1: 1/6
     b2: 1/3
     b3: 1/3
     b4: 1/6
Contents