2.3. Derivation of Explicit Runge-Kutta Methods#

The derivation of an explicit Runge-Kutta 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.

For example, consider the second-order Taylor series expansion of \(y_{n+1}\)

(2.4)#\[ y_{n+1} = y_n + h y_n' + \frac{h^2}{2} y_n'' + O(h^3).\]

We wish to solve the ODE \(y' = f(t,y)\), applying the chain rule to differentiate \(y\) we have

\[\begin{split} \begin{align*} y''(t, y) &= f'(t, y) \\ &= (f_t(t, y) + f_y y')(t, y) \\ &= (f_t + f_yf)(t, y), \end{align*} \end{split}\]

where the subscript notation \(f_y\) denotes the partial derivative of \(f\) with respect to \(y\).

So the Taylor series in equation (2.4) becomes

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

This is the second-order Taylor series expansion for the ODE \(y' = f(t,y)\). To be able to solve this ODE we need the Taylor series expansion of a Runge-Kutta method to be equivalent to equation (2.5). The general form of a 2-stage explicit Runge-Kutta method is

(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_2 h, y_n + h a_{21} k_1) \end{align*} \end{split}\]

Using the bivariate Taylor series which is

\[ f(t + h, y + k) = (f + h f_t + k f_y)(t,y) + O(h^3), \]

then replacing \(h\) and \(k\) with \(c_2h\) and \(h_{a21}k_1\) respectively, then we can rewrite \(k_2\) as

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

Substituting \(k_1\) and \(k_2\) into the equation for \(y_{n+1}\) in equation (2.6) gives

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

We need equation (2.7) to be equal to equation (2.5). Equating the coefficients of \(f\) we have (the \(h\)’s cancel out)

(2.8)#\[ b_1 + b_2 = 1. \]

Equating the coefficients of \(f_t\) we have

(2.9)#\[ b_2c_2 = \frac{1}{2}, \]

and equating the coefficients of \(ff_y\) we have

\[ a_{21}b_2 = \frac{1}{2}. \]

From equation (2.9) we know that \(b_2 = \dfrac{1}{2c_2}\) so

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

So we have the three equations (2.8), (2.9) and (2.10) expressed in four unknowns (\(a_{21}\), \(b_1\), \(b_2\), \(c_2\)). Any set of values that satisfy these this system of equations give a valid second-order explicit Runge-Kutta 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.3 (Order conditions for a second-order explicit Runge-Kutta method)

(2.11)#\[\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.

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

Note

There are an infinite number of combinations for the values of \(a_{ij}\), \(b_i\) and \(c_i\) which satisfy the order conditions for a Runge-Kutta method. All of the possible Runge-Kutta methods of a particular order will the same solutions (computational rounding permitted).

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)