Week 2 – Explicit Runge-Kutta Methods
6G6Z3017 Computational Methods in ODEs
Pop Quiz
- What is an ODE?
An Ordinary Differential Equation which is an equation expressed in terms of an independent variable, a function of the independent variable and derivatives of that function.
- What three things a required for an IVP?
An initial value problem requires an ODE, a range for the value of the independent variable and the initial values of the solution.
- What are the local and global truncation errors of a method?
The local truncation error is the error accrued over one step of a method assuming the previous step was exact.
The global truncation error is the accumulation of the local truncation errors.
- What does it mean when we say the error of a method is \(O(h^n)\)?
As the step length \(h\) tends to zero, the value of the error will tend to zero at least as fast as \(h^n\).
Runge-Kutta Methods
Definition: General form of a Runge-Kutta Method
The general form of a Runge-Kutta method for solving the initial value problem \(y' =f(t,y)\), \(t \in [t_0, t_{\max}]\), \(y(0) = y_0\) is
\[ \begin{align} y_{n+1} &=y_n + h \sum_{i=1}^s b_i k_i,\\ k_i &= f(t_n +c_i h,y_n + h \sum_{j=1}^s a_{ij} k_j ), \end{align} \]
Runge-Kutta methods are known as single step methods because they only need information from a single step, e.g., \(t_n\) and \(y_n\) to compute the solution at the next step.
The \(k_i\) values are intermediate calculations known as stage values
Heun’s Method
An example of a Runge-Kutta method is Heun’s method which 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 + h k_1), \\ \end{align*} \qquad(1)\]
The general form of a 2-stage method is
\[ \begin{align*} y_{n+1} &= y_n + h (b_1k_1 + b_2 k_2 ), \\ k_1 &= f(t_n + c_1 h, y_n + h(a_{11} k_1 + a_{12} k_2)), \\ k_2 &= f(t_n + c_2 h, y_n + h(a_{21} k_1 + a_{22} k_2)). \end{align*} \qquad(2)\]
Comparing Equation 1 and Equation 2 we can see that
\[ \begin{align*} b_1 &= \frac{1}{2}, & b_2 &= \frac{1}{2}, & c_1 &= 0, & c_2 &= 1, \\ a_{11} &= 0, & a_{12} &= 0, & a_{21} &= 1, & a_{22} &= 0. \end{align*} \]
Butcher Tableau
Runge-Kutta methods are often summarised in a Butcher tableau named after the New Zealand mathematician John Butcher.
A Butcher tableau is a table of values containing the coefficients \(a_{ij}\), \(b_i\) and \(c_i\) for a Runge-Kutta method
\[ \begin{align*} \begin{array}{c|cccc} c_1 & a_{11} & a_{12} & \cdots & a_{1s} \\ c_2 & a_{21} & a_{22} & \cdots & a_{2s} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ c_s & a_{s1} & a_{s2} & \cdots & a_{ss} \\ \hline & b_1 & b_2 & \cdots & b_s \end{array} \end{align*} \]
For example, recall Heun’s method
\[ \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 + h k_1), \\ \end{align*} \]
where \(a_{11} = a_{12} = a_{22} = 0\), \(b_1 = b_2 = \frac{1}{2}\), \(c_1 = 0\) and \(c_2 = 1\). So this can be expressed as the Butcher tableau
\[ \begin{align*} \begin{array}{c|cccc} 0 & \\ 1 & 1 \\ \hline & \frac{1}{2} & \frac{1}{2} \end{array} \end{align*} \]
Activity
Write down the Butcher tableau for the 4-stage Runge-Kutta method given below as a Butcher tableau
\[ \begin{align*} y_{n+1} &= y_n + \frac{h}{6}(k_1 + 2k_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 + h k_3). \end{align*} \]
Solution \[ \begin{align*} \begin{array}{c|cc} 0 & 0 & 0 & 0 & 0 \\ \frac{1}{2} & \frac{1}{2} & 0 & 0 & 0 \\ \frac{1}{2} & 0 & \frac{1}{2} & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ \hline & \frac{1}{6} & \frac{1}{3} & \frac{1}{3} & \frac{1}{6} \end{array} \end{align*} \]
Implicit Runge-Kutta Methods
The stage values for an \(s\)-stage Runge-Kutta method are
\[ \begin{align*} k_1 &=f(t_n +c_1 h,y_n +h(a_{11} k_1 +a_{12} k_2 +\cdots +a_{1s} k_s )),\\ k_2 &=f(t_n +c_2 h,y_n +h(a_{21} k_1 +a_{22} k_2 +\cdots +a_{2s} k_s )),\\ &\vdots \\ k_s &=f(t_n +c_s h,y_n +h(a_{s1} k_1 +a_{s2} k_s +\cdots +a_{ss} k_s )). \end{align*} \]
These are implicit functions, e.g., the value \(k_1\) appears on the right-hand side of the equation for \(k_1\).
Runge-Kutta methods where the stage values are expressed using implicit functions are known as Implicit Runge-Kutta (IRK) methods.
Explicit Runge-Kutta Methods
If the upper limit of the sum in the expression for \(k_i\) is changed to \(i-1\), i.e.,
\[ \begin{align*} k_i = f(t_n + c_i h,y_n + h\sum_{j=1}^{i-1} a_{ij} k_j) \end{align*} \]
and let \(c_1 = 0\), then we have the following equations for calculating the stage values
\[ \begin{align*} k_1 &=f(t_n ,y_n),\\ k_2 &=f(t_n +c_2 h,y_n +ha_{21} k_1 ),\\ k_3 &=f(t_n +c_3 h,y_n +h(a_{31} k_1 +a_{32} k_2 )),\\ &\vdots \\ k_s &=f(t_n +c_s h,y_n +h(a_{s1} k_1 +a_{s2} k_s +\cdots +a_{s,s-1} k_{s-1} )). \end{align*} \]
These are explicit functions so Runge-Kutta methods and are the stage values for an explicit Runge-Kutta method.
Butcher Tableau of Explicit and Implicit Runge-Kutta Methods
Explicit Method
\[ \begin{array}{c|ccccc} 0 & 0 & & & & \\ c_2 & a_{21} & & & & \\ c_3 & a_{31} & a_{32} & & & \\ \vdots & \vdots & \vdots & \ddots & & \\ c_s & a_{s1} & a_{s2} & \cdots & a_{s,s-1} & \\ \hline & b_1 & b_2 & \cdots & b_{s-1} & b_s \end{array} \]
Implicit Method
\[ \begin{array}{c|ccccc} c_1 & a_{11} & a_{12} & a_{13} & \cdots & a_{1s} \\ c_2 & a_{21} & a_{22} & a_{23} & \cdots & a_{2s} \\ c_3 & a_{31} & a_{32} & a_{33} & \cdots & a_{3s} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ c_s & a_{s1} & a_{s2} & a_{s3} & \cdots & a_{ss} \\ \hline & b_1 & b_2 & b_3 & \cdots & b_s \end{array} \]
- Explicit methods are quicker and more straight forward to compute but do not perform well on stiff problems.
- Implicit methods are more complicated to compute but can be used to solve stiff problems
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
The second-order Taylor series is
\[ \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*} \qquad(3)\]
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 3 and using the notation \(y_{n+1} = y(t_n + h)\) and \(y_n = y(t_n)\) we have
\[ \begin{align*} 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) \end{align*} \]
This is the Taylor expansion of the ODE \(y' = f(t, y)\) which are attempting to solve.
What we now need to do is determine an equivalent expansion for the 2-stage ERK method
\[ \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*} \]
Substituting \(k_1\) and \(k_2\) into \(y_{n+1}\)
\[ \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*} \qquad(4)\]
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 4
\[ \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_y(t_n, y_n) 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_y(t_n, y_n) f(t_n, y_n)] + O(h^3). \\ \end{align*} \]
This is the Taylor expansion of a general 2-stage ERK method.
We need the following two equations to be equivalent.
\[ \begin{align*} y_{n+1} &= y_n + h f(t_n, y_n) + \frac{h^2}{2} f_t(t_n, y_n) + \frac{h^2}{2} f_y(t_n, y_n)f(t_n, y_n) + O(h^3), \\ y_{n+1} &= y_n + h (b_1 + b_2) f(t_n, y_n) + h^2 b_2 c_2 f_t(t_n, y_n) + h^2 b_2a_{21} f_y(t_n, y_n) f(t_n, y_n) + O(h^3). \end{align*} \]
Equating the coefficients of \(hf(t_n, y_n)\)
\[ \begin{align*} b_1 + b_2 &= 1 \end{align*} \]
Equating the coefficients of \(h^2 f_t(t_n, y_n)\)
\[ \begin{align*} b_2 c_2 & = \frac{1}{2}. \end{align*} \qquad(5)\]
Equating the coefficients of \(h^2 f_y(t_n, y_n)f(t_n, y_n)\)
\[ \begin{align*} b_2 a_{21} &= \frac{1}{2}, \end{align*} \]
since we know from Equation 5 that \(b_2 = \dfrac{1}{2c_2}\) so \(a_{21} = c_2\).
Order Conditions for a Second-order ERK Method
Order conditions for a second-order explicit Runge-Kutta method
\[ \begin{align*} b_1 +b_2 &=1,\\ c_2 b_2 &=\frac{1}{2},\\ a_{21} &= c_2. \end{align*} \]
We have three equations in four unknowns so to derive a second-order ERK method we need to pick a value for one of \(b_1\), \(b_2\), \(c_2\) or \(a_{21}\) and solve for the rest.
Activity
- Derive a second-order Runge-Kutta ERK methods where \(c_2 = 1\)
Solution Substituting \(c_2 = 1\) into the order conditions gives
\[ \begin{align*} b_1 +b_2 &=1,\\ b_2 &=\frac{1}{2},\\ a_{21} &= c_2. \end{align*} \]
so \(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*}\]
\[\begin{align*} \begin{array}{c|cc} 0 & & \\ 1 & 1 & \\ \hline & \frac{1}{2} & \frac{1}{2} \end{array} \end{align*}\]
This is Heun’s method seen earlier.
- Derive a second-order Runge-Kutta ERK methods where \(b_2 = 1\)
Solution Substituting \(b_2 = 1\) into the order conditions gives
\[ \begin{align*} b_1 + 1 &= 1,\\ c_2 &=\frac{1}{2}, \\ a_{21} &= c_2. \end{align*} \]
so \(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*}\]
\[\begin{align*} \begin{array}{c|cc} 0 & \\ \frac{1}{2} & \frac{1}{2} \\ \hline & 0 & 1 \end{array} \end{align*}\]
This method is known as the midpoint method.
Using Python 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.
The following code shows how Python and MATLAB can be used to solve the order conditions.
Using MATLAB to solve the order conditions
Rooted Trees
A rooted tree is a tree where a single node is designated as the root. Rooted trees are usually drawn with the root node at the bottom of the tree.
The order of a rooted tree \(\tau\), denoted by \(r(\tau)\), is the total number of nodes in the tree.
The density of a rooted tree \(\tau\), denoted by \(\gamma(\tau)\), is a measure related to the number of nodes and the children of the notes in the tree. The density can be calculated using the following:
- assign a value of 1 to each leaf node in the tree;
- for all non-leaf nodes assign a value of 1 more than the sum of the values of their child nodes;
- calculate the product the values of all nodes in the tree which gives \(\gamma(\tau)\).
Activity
Determine the order and density of this rooted tree.
Solution
There are 7 nodes in this tree so \(r(\tau)=7\).
Assigning the weights to the nodes gives
So \(\gamma(\tau) = 7 \cdot 1 \cdot 3 \cdot 2 \cdot 1 \cdot 1 \cdot 1 = 42\)
Determining derivatives of \(f(y)\) using rooted trees
Consider the derivatives of the function \(y = f(y)\) using the product and chain rules
\[ \begin{align*} y' &= f(y), \\ y'' &= f'(y)y' \\ &= (f'f)(y), \\ y''' &= f''(y)(f(y), f'(y)) + f'(y)f'(y)y' \\ &= (f''(f,f) + f'f'f)(y) \\ &\vdots \end{align*}\]
As you can see this soon becomes complicated
English mathematician Robert Merson noticed that the structure of rooted trees matches that of the derivatives of \(f(y)\).
Defining a rooted tree notation where:
- \(\mathbf{f}\) denotes a node of the tree and the number of children the node has is denoted by a prime
- a child node is written immediately following the parent node
- sibling nodes are separated by a comma (using brackets if there are more than one child).
| \(\tau\) |  |
 |
 |
 |
|---|---|---|---|---|
| \(r(\tau)\) | 1 | 2 | 3 | 3 |
| \(\gamma(\tau)\) | \(\mathbf{f}\) | \(\mathbf{f}'\mathbf{f}\) | \(\mathbf{f}''(\mathbf{f},\mathbf{f})\) | \(\mathbf{f}'\mathbf{f}'\mathbf{f}\) |
If we sum the two order 3 trees then
\[ \mathbf{f}''(\mathbf{f},\mathbf{f}) + \mathbf{f}'\mathbf{f}'\mathbf{f} = (f''(f,f) + f'f'f)(y) = y'''\]
Elementary Weights
To determine the values of \(a_{ij}\), \(b_i\) and \(c_i\) in the Butcher tableau we define a set of elementary weights, denoted by \(\Phi(\tau)\) which are generated using the following steps:
- associate a label \(i\) to the root node and for all non-leaf nodes associate unique labels \(j\), \(k\), \(\ell\) etc;
- write down a sequence of factors for which the first is \(b_i\);
- for each edge that does not terminate in a leaf node, write down another factor \(a_{pq}\) where \(p\) and \(q\) are the labels associated with the parent and child node respectively;
- for each leaf node write down a factor \(c_p\) where \(p\) is the label associated with the parent node;
- sum the product of factors for all possible choices of the labels.
Activity
Determine the elementary weight of this rooted tree.
Solution Assigning the labels to the non-leaf nodes gives
\[ \begin{align*} \Phi(\tau) &= \sum_{i,j,k} b_ia_{ij} a_{ik} c_ic_jc_jc_k = \sum_{i,j,k} b_ia_{ij} a_{ik} c_ic_j^2c_k. \end{align*} \]
Deriving Order Conditions using Rooted Trees
John Butcher noticed that comparing the Taylor series expansions for the autonomous ODE \(y'=f(y)\) and the general form of a Runge-Kutta method gives the order condition
\[\Phi(\tau) = \dfrac{1}{\gamma(\tau)}. \qquad(6)\]
To derive the order conditions for a \(k\)-th order \(s\)-stage Runge-Kutta method we determine the set \(T\) of all rooted trees of order up to and include \(k\) and use Equation 6 for each tree in \(T\).
In addition to these order conditions there is also the following row-sum condition placed on the value of \(c_i\)
\[c_i = \sum_{j=1}^sa_{ij},\]
i.e., \(c_i\) is equal to the sum of row \(i\) of \(A\)
Activity
Use trees to derive the order conditions for a 2-stage second-order explicit Runge-Kutta method.
Solution
| Tree |  |
 |
|---|---|---|
| \(\Phi(\tau)\) | \(\sum b_i\) | \(\sum b_ic_i\) |
| \(\gamma(\tau)\) | 1 | 2 |
So the order conditions are (including the row-sum condition) \[ \begin{align*} b_1 + b_2 &= 1, \\ b_1c_1 + b_2c_2 &= \frac{1}{2}, \\ c_1 &= a_{11} + a_{12}, \\ c_2 &= a_{21} + a_{22}. \end{align*} \]
For an ERK method, \(c_1 = 0\) and \(a_{11} = a_{12} = a_{22} = 0\) so the order conditions reduce to
\[ \begin{align*} b_1 + b_2 &= 1, \\ b_2c_2 &= \frac{1}{2}, \\ c_2 &= a_{21}. \end{align*} \]
Derivation of a Fourth-order ERK Method
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
| \(\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 fewest number of stages required for a fourth-order explicit Runge-Kutta method is 4, so \(s = 4\)
\[ \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*} \]
Since \(c_1 = 0\), then for the first 4 order conditions on the top row
\[ \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} \]
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{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*} \]
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*} \]
Order conditions for a 4-stage fourth-order ERK method
\[ \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*} \]
Here we have 11 order conditions and 14 coefficients
We choose values for 3 coefficients and solve for the rest
Derivation of the RK4 Method (Python)
Derive a fourth-order Runge-Kutta method where \(c_2 = c_3 = \frac{1}{2}\) and \(c_4 = 1\).
Using Python to solve the order conditions and the row sum conditions. Solving the first four order conditions for the \(b_i\) coefficients.
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))Which outputs \(b_1 = b_4 = \frac{1}{6}\) and \(b_2 = \frac{2}{3} - b_3\). Lets choose \(b_3 = \frac{1}{3}\) so that \(b_2 = \frac{1}{3}\).
Solve the remaining order conditions to determine the other coefficients
# 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))Which outputs \(a_{31} = a_{41} = a_{42} = 0\), \(a_{21} = a_{32} = \frac{1}{2}\) and \(a_{43} = 1\).
So the RK4 method is
\[ \begin{align*} \begin{array}{c|cccc} 0 & & & & \\ \frac{1}{2} & \frac{1}{2} & & & \\ \frac{1}{2} & 0 & \frac{1}{2} & & \\ 1 & 0 & 0 & 1 & \\ \hline & \frac{1}{6} & \frac{1}{3} & \frac{1}{3} & \frac{1}{6} \end{array} \end{align*} \]
\[ \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*} \]
Derivation of the RK4 Method (MATLAB)
Solve the first four order conditions
% 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)Which outputs \(b_1 = b_4 = \frac{1}{6}\) and \(b_2 = \frac{2}{3} - b_3\). Lets choose \(b_3 = \frac{1}{3}\) so that \(b_2 = \frac{1}{3}\).
Solve the remaining order conditions to determine the other coefficients
% 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)Which outputs \(a_{31} = a_{41} = a_{42} = 0\), \(a_{21} = a_{32} = \frac{1}{2}\) and \(a_{43} = 1\).
Summary
- Runge-Kutta methods are single step methods since they only require values from 1 step to compute the next step
- Runge-Kutta methods are commonly presented in tabular form known as a Butcher tableau
- Explicit Runge-Kutta (ERK) methods use explicit functions for the stage values (straight forward to calculate)
- Implicit Runge-Kutta (IRK) methods use implicit functions for the stage values (more complex to calculate but can be applied to stiff problems)
- The Taylor series expansion of the ODE \(y' = f(t, y)\) and the general Runge-Kutta method are compared to give a set of order conditions the govern the coefficients of the method
- Rooted trees present an easier way of deriving the order conditions
- Choosing values for some of the coefficients allows us to solve the order conditions (using software) and derive Runge-Kutta methods