2.3. Derivation of Explicit Runge-Kutta Methods#

The derivation of an explicit Runge-Kutta (ERK) method is achieved by comparing the Taylor series for the first-order ODE \(y' = f(t,y)\) to that of the general Runge-Kutta method and ensuring the coefficients \(a_{ij}\), \(b_i\) and \(c_i\) match.

We will consider the derivation of a second-order ERK method. The second-order Taylor series expansion of \(y(t_n + h)\) is

(2.4)#\[ y(t_n + h) = y(t_n) + h y'(t_n) + \frac{h^2}{2} y''(t_n) + O(h^3).\]

Since we are solving the ODE \(y' = f(t, y)\), we know that \(y'(t_n) = f(t_n, y_n)\), and we can use the chain rule to determine \(y''(t_n)\)

\[ \begin{align*} y''(t_n) &= f_t(t_n, y_n) + f_y(t_n, y_n) y_n' = f_t(t_n, y_n) + f_y(t_n, y_n)f(t_n, y_n). \end{align*} \]

Substituting \(y'(t_n)\) and \(y''(t_n)\) into equation (2.4) and using the notation \(y_{n+1} = y(t_n + h)\) and \(y_n = y(t_n)\) gives

(2.5)#\[ y_{n+1} = y_n + h f(t_n, y_n) + \frac{h^2}{2} [ f_t(t_n, y_n) + f_y(t_n, y_n)f(t_n, y_n) ] + O(h^3) \]

This is the Taylor expansion of the ODE we are attempting to solve. What we now need to do is determine an equivalent expansion for the 2-stage ERK method

(2.6)#\[\begin{split} \begin{align*} y_{n+1} &= y_n + h (b_1 k_1 + b_2 k_2), \\ k_1 &= f(t_n, y_n), \\ k_2 &= f(t_n + c_2h, y_n + a_{21} h k_1). \end{align*} \end{split}\]

For a second-order method, we need to ensure that that values of \(a_{21}\), \(b_1\), \(b_2\) and \(c_2\) result in a method that has truncation errors \(O(h^3)\). Substituting \(k_1\) and \(k_2\) into \(y_{n+1}\)

(2.7)#\[\begin{split} \begin{align*} y_{n+1} &= y_n + h b_1 f(t_n, y_n) + h b_2 f(t_n + c_2h, y_n + a_{21} h k_1) + O(h^3) \\ \end{align*} \end{split}\]

The bivariate Taylor expansion of the two-variable function \(f(t_n + c_2h, y_n + a_{21} h k_1)\) correct to \(O(h^3)\) is

\[ \begin{align*} f(t_n + c_2 h, y_n + a_{21} h k_1) = f(t_n, y_n) + c_2 h f_t(t_n, y_n) + a_{21} h f_y(t_n, y_n) f(t_n, y_n) + O(h^3). \end{align*} \]

Substituting this into equation (2.7)

(2.8)#\[\begin{split} \begin{align*} y_{n+1} &= y_n + h b_1 f(t_n, y_n) + h b_2 [f(t_n, y_n) + c_2 h f_t(t_n, y_n) + a_{21} h f(t_n, y_n)] + O(h^3) \\ &= y_n + h (b_1 + b_2) f(t_n, y_n) + h^2 b_2 [c_2f_t(t_n, y_n) + a_{21} f(t_n, y_n)] + O(h^3) \end{align*} \end{split}\]

This is the Taylor expansion of a general 2-stage ERK method. We need equations(2.5) and (2.8) to be equal so that the local truncation errors in the ERK method are \(O(h^3)\).

Equating the coefficients of \(f(t_n, y_n)\)

(2.9)#\[\begin{split} \begin{align*} h( b_1 + b_2) &= h, \\ \therefore b_1 + b_2 &= 1. \end{align*} \end{split}\]

Equating the cofficients of \(f_t(t_n, y_n)\)

(2.10)#\[\begin{split} \begin{align*} h^2 b_2 c_2 &= \frac{h^2}{2} \\ \therefore b_2 c_2 &= \frac{1}{2}. \end{align*} \end{split}\]

Equating the coefficients of \(f_y(t_n, y_n)f(t_n, y_n)\)

\[\begin{split} \begin{align*} h^2 b_2 a_{21} &= \frac{h^2}{2} \\ \therefore b_2 a_{21} &= \frac{1}{2}, \end{align*} \end{split}\]

since we know from equation (2.10) that \(b_2 = \dfrac{1}{2c_2}\) then

(2.11)#\[ a_{21} = c_2 .\]

So we have the three equations (2.9), (2.10) and (2.11) expressed in four unknowns \(a_{21}\), \(b_1\), \(b_2\) and \(c_2\). Any set of values that satisfy these this system of equations give a ERK method with local truncation errors \(O(h^3)\), so the global truncation errors are \(O(h^2)\) and we have a second-order method.

These conditions are known as the order conditions for a method. Since we have an underdetermined system to get a unique solution we choose a value for one of the unknowns and solve for the others.

Definition 2.1 (Order conditions for a second-order explicit Runge-Kutta method)

(2.12)#\[\begin{split} \begin{align*} b_1 +b_2 &=1,\\ c_2 b_2 &=\frac{1}{2},\\ a_{21} &= c_2. \end{align*} \end{split}\]

Example 2.2

Derive two second-order Runge-Kutta ERK methods where

(i)   \(c_2 = 1\);

Solution (click to show)

Substituting \(c_2 = 1\) into the order conditions gives \(b_2 = \frac{1}{2}\), \(b_1 = \frac{1}{2}\) and \(a_{21} = 1\). So this second-order ERK method is

\[\begin{align*} y_{n+1} &=y_n +\frac{h}{2}(k_1 +k_2 ),\\ k_1 &=f(t_n ,y_n ),\\ k_2 &=f(t_n +h,y_n +hk_1 ), \end{align*}\]

or expressed using a Butcher tableau

\[\begin{align*} \begin{array}{c|cc} 0 & & \\ 1 & 1 & \\ \hline & 1/2 & 1/2 \end{array} \end{align*}\]

This method is known as Heun’s method or the RK2 method.

(ii)   \(b_2 = 1\).

Solution (click to show)

Substituting \(b_2=1\) into the order conditions gives \(b_1=0\), \(c_2 = \frac{1}{2}\) and \(a_{21} = \frac{1}{2}\) so this second-order ERK method is

\[\begin{align*} y_{n+1} &= y_n + h k_2, \\ k_1 &= f(t_n, y_n), \\ k_2 &= f(t_n + \tfrac{1}{2}h, y_n + \tfrac{1}{2}hk_1), \end{align*}\]

or expressed using a Butcher tableau

\[\begin{align*} \begin{array}{c|cc} 0 & \\ 1/2 & 1/2 \\ \hline & 0 & 1 \end{array} \end{align*}\]

This method is known as the midpoint method.

2.3.1. Using Python and MATLAB to solve the order conditions#

The algebra used to solve the order conditions in Example 2.2 is quite simple but for higher order methods it can soon get more complicated (see the derivation of a fourth-order explicit Runge-Kutta method for an example) and it is useful to use software to help with the algebra. The following code shows how Python and MATLAB can be used to solve the order conditions for Example 2.2.

There is a Python library SymPy (short for Symbolic Python) that has functions that can solve algebraic equations. After importing SymPy we need to declare each of the coefficients \(a_{21}\), \(b_1\), \(b_2\) and \(c_2\) as symbolic variables using the sp.symbols() command before setting \(c_2=1\). Each order condition is then defined using these symbolic variables. Note that SymPy assumes that equations are equal to zero which is why we have subtracted the right-hand side. We have also used the sp.Rational(1,2) command for the fraction \(\frac{1}{2}\) so that SymPy will output any fractional values as fractions rather than decimals. The system of three equations eq1, eq2 and eq3 is then solved using the sp.solve command.

import sympy as sp

# Declare symbolic variables
a21, b1, b2, c2 = sp.symbols("a21, b1, b2, c2")
c2 = 1

# Define order conditions
eq1 = b1 + b2 - 1
eq2 = b2 * c2 - sp.Rational(1,2)
eq3 = a21 - c2

# Solve order conditions
sp.solve((eq1, eq2, eq3))

We first declare each of the coefficients \(a_{21}\), \(b_1\), \(b_2\) and \(c_2\) as symbolic variables using the syms command before setting \(c_2=1\). Each order condition is then defined using these symbolic variables. Note that the == command is used for the equals sign in each order condition. The system of three equations eq1, eq2 and eq3 is then solved using the solve command.

% Declare symbolic variables
syms a21 b1 b2 c2
c2 = 1;

% Define order conditions
eq1 = b1 + b2 == 1;
eq2 = b2 * c2 == 1/2;
eq3 = a21 == c2;

% Solve order conditions
solve(eq1, eq2, eq3)