Derivation of Explicit Runge-Kutta Methods

2.3. Derivation of Explicit Runge-Kutta Methods#

Recall the general form of an explicit Runge-Kutta methods

\[\begin{split} \begin{align*} y_{n+1} &= y_n + h \sum_{i = 1}^s b_i k_i, \\ k_i &= f(t_n + c_ih, y_n + h \sum_{j = 1}^{i - 1} a_{ij} k_j). \end{align*} \end{split}\]

We need to address a fundamental question: how are the coefficients, \(a_{ij}\), \(b_i\) and \(c_i\), of a Runge-Kutta method chosen? To construct a method of a given order, we compare the exact solution with the Taylor series expansion of the numerical method. By matching the coefficients of corresponding powers of \(h\), we obtain equations that the Runge-Kutta coefficients must satisfy. These equations are known as order conditions.

Here we derive a second-order ERK method. A second-order method has global error \(O(h^2)\), which corresponds to a local truncation error of \(O(h^3)\). The Taylor series expansion of \(y(t_n + h)\) to \(O(h^3)\) is

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

We want to solve \(y' = f(t, y)\), for convenience let \(f = f(t_n, y_n)\) so the series expansion is

\[ y(t + h) = y(t_n) + hf + \frac{h^2}{2} f' + O(h^3). \]

The total derivative of \(f\) is

\[ f' = f_x + f_y y', \]

and substituting \(y' = f\) and \(y'' = f'\) into equation (2.4)

(2.5)#\[ \begin{align*} y(t+h) = y(t) + h f + \frac{h^2}{2} ( f_t + f_yf ) + O(h^3) \end{align*} \]

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

\[\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 + h a_{21} k_1). \end{align*} \end{split}\]

Substituting \(k_1\) and \(k_2\) into \(y_{n+1}\)

(2.6)#\[ y_{n+1} = y_n + h (b_1 f(t_n, y_n) + b_2 f(t_n + c_2h, y_n + ha_{21}k_1)). \]

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 error \(O(h^3)\). We need a Taylor series expansion for \(k_2\) for which we can use the bivariate Taylor series expansion up to \(O(h^3)\)

\[ \begin{align*} f(t + h_1, y + h_2) = f + h_1 f_t + h_2 f_y + \frac{h_1^2}{2}f_{tt} + h_1h_2 f_{ty} + \frac{h_2^2}{2}f_{yy} + O((|h_1| + |h_2|)^3). \end{align*} \]

Let \(h_1 = c_2h\) and \(h_2 = h a_{21} k_1\) then

\[ \begin{align*} f(t + c_2h, y + ha_{21}k_1) = f + c_2h f_t + ha_{21}k_1 f_y + \frac{c_2^2h^2}{2}f_{tt} + c_2a_{21}k_1h^2 f_{ty} + \frac{h^2 a_{21}^2 k_1^2}{2}f_{yy} + O(h^2). \end{align*} \]

Since this expansion will subsequently be multiplied by \(h\) in equation (2.6), we only require terms up to \(O(h)\) in the expansion of \(f\). Higher-order terms contribute only to \(O(h^3)\) and can therefore be absorbed into the remainder.

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

Since \(k_1 = f(t_n, y_n) = f\), we may replace \(k_1\) by \(f\) int he first-order terms

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

Substituting this into equation (2.6)

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

Equation (2.5) represents the exact solution expanded about \(t_n\)​, whereas equation (2.7) represents the numerical Runge-Kutta method​. To achieve second-order accuracy, these expansions must agree through all terms of order \(h^2\).

2.3.1. Order conditions#

To derive a second-order ERK method to solve the ODE \(y' = f(t, y)\), we need equations(2.5) and (2.7) to be equal so that the local truncation errors in the ERK method are \(O(h^3)\). We have four unknowns in (2.7), \(b_1\), \(b_2\), \(c_2\) and \(a_{21}\), so that the two expressions are equal.

Since \(f\), \(f_t\) and \(f_yf\) are independent quantities, the coefficients multiplying each term must be equal in the two expansions. For the coefficents of \(f\)

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

for the coefficients of \(f_t\)

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

and for the coefficients of \(f_yf\)

\[\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.9) that \(b_2 = \dfrac{1}{2c_2}\) then

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

So we have the three equations (2.8), (2.9) and (2.10) expressed in the four unknowns. Any set of values that satisfy these this system of equations give an 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 under determined system to get a unique solution we choose a value for one of the unknowns and solve for the others.

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

By matching the Taylor series expansion of the exact solution with that of the Runge-Kutta method, we obtain the following order conditions:

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

Any coefficients satisfying these conditions define a second-order explicit Runge-Kutta method.

2.3.2. Deriving a second-order ERK method#

Let’s derive a second-order ERK where \(c_2 = 1\), so the order conditions are

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

So here we already have the solution \(b_2 = \frac{1}{2}\) and \(a_{21} = 1\) and we just need to solve the first order condition to give \(b_1 = 1 - \frac{1}{2} = \frac{1}{2}\). Now we have values for the four unknowns we can write the ERK method

\[\begin{split} \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*} \end{split}\]

or alternatively as the Butcher tableau

\[\begin{split} \begin{array}{c|cc} 0 & 0 & 0 \\ 1 & 1 & 0 \\ \hline & \frac{1}{2} & \frac{1}{2} \end{array} \end{split}\]

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

Example 2.2

Derive a second-order ERK method where \(b_2 = 1\).

Solution

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 & \\ \frac{1}{2} & \frac{1}{2} \\ \hline & 0 & 1 \end{array} \end{align*}\]

Since \(c_2 = \frac{1}{2}\), the derivative used to advance the solution is evaluated at the midpoint of the interval, hence this method is known as the midpoint method.

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

The algebra used to solve the order conditions for a second-order method is quite simple but for higher order methods it can soon get more complicated. Therefore 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.

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

This produces the output

{a21: 1, b1: 1/2, b2: 1/2}

Here we have used the Python library SymPy (short for Symbolic Python) that has functions that can solve algebraic equations. After importing SymPy we have declared symbolic variables for each of the coefficients \(a_{21}\), \(b_1\), \(b_2\) and \(c_2\) 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.

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

This produces the output

ans = struct with fields:
    a21: 1
     b1: 1/2
     b2: 1/2

Here we have declared symbolic variables for each of the coefficients \(a_{21}\), \(b_1\), \(b_2\) and \(c_2\) 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.