4 Runge-Kutta Methods
Although the Euler method works adequately in many applications, its accuracy and efficiency diminish for large problems, high dimensions, or stiff equations. For more accurate results, especially when smaller error tolerances are required, we turn to Runge-Kutta methods (see [1],[2])
4.1 Explicit Runge-Kutta Methods
The general form of an explicit Runge-Kutta method is a one-step method that approximates the solution by using multiple intermediate stages within each time step. The method is represented as follows:
\begin{aligned} \bm{z}_{k+1} &= \bm{z}_k + \sum_{i=1}^{s} b_i \bm{K}_i, \\ &\left\{\begin{aligned} \bm{K}_1 &= h \bm{f}\Big(t_k, \bm{z}_k\Big), \\[0.5em] \bm{K}_2 &= h \bm{f}\Big(t_k + h c_2, \bm{z}_k + a_{21} \bm{K}_1\Big), \\[0.5em] \bm{K}_3 &= h \bm{f}\Big(t_k + h c_3, \bm{z}_k + a_{31} \bm{K}_1 + a_{32} \bm{K}_2\Big), \\[0.5em] \vdots & \\[0.5em] \bm{K}_s &= h \bm{f}\Big(t_k + h c_s, \bm{z}_k + \sum_{j=1}^{s-1} a_{sj} \bm{K}_j\Big), \end{aligned} \right. \end{aligned}
where \bm{K}_1, \bm{K}_2, \dots, \bm{K}_s are intermediate values called stages, and s is the number of stages. The constants a_{ij}, b_i, and c_i define the specific Runge-Kutta method and are chosen to achieve higher-order accuracy. The method computes multiple weighted sums of function evaluations (the \bm{K}_i values) to achieve more accurate approximations of the solution.
The values t_k + h c_i represent quadrature points in the time interval [t_k, t_{k+1}], and the quantities K_i / h approximate the function values at these points.
Butcher Tableau
To simplify the notation for Runge-Kutta methods, the coefficients a_{ij}, b_i, and c_i are often summarized in a compact form called the Butcher tableau. The tableau for the explicit Runge-Kutta method is given by:
\begin{array} {c|ccccc} 0 & 0 & 0 & 0 & 0 & 0 \\ c_2 & a_{21} & 0 & 0 & 0 & 0 \\ c_3 & a_{31} & a_{32} & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & 0 & 0 \\ c_s & a_{s1} & a_{s2} & \ldots & a_{s,s-1} & 0 \\ \hline & b_1 & b_2 & \ldots & b_{s-1} & b_s \end{array}
Explicit Euler
The Explicit Euler Method Represented Using a Butcher Tableau:
\def\arraystretch{1.5} \begin{array}{c|c} 0 & 0 \\ \hline & 1 \end{array} \tag{4.1}
Midpoint method
The Midpoint Method, also known as the Collatz Method, is a variant of the explicit Euler Method. This approach is classified as a second-order method and consists of two stages, allowing for improved accuracy in numerical integration.
\def\arraystretch{1.5} \begin{array}{c|cc} 0 & 0 & 0 \\ 1/2 & 1/2 & 0 \\ \hline & 0 & 1 \end{array} \tag{4.2}
Heun’s method
Heun’s Method, often referred to as the explicit trapezoid rule, enhances the standard explicit Euler Method. This technique is a second-order method that employs two stages to achieve greater accuracy in numerical integration.
\def\arraystretch{1.5} \begin{array}{c|cc} 0 & 0 & 0 \\ 1 & 1 & 0 \\ \hline & 1/2 & 1/2 \end{array} \tag{4.3}
Ralston’s method
Ralston’s method is a second-order technique with two stages, specifically designed to achieve a minimal local error bound.
\def\arraystretch{1.5} \begin{array}{c|cc} 0 & 0 & 0 \\ 2/3 & 2/3 & 0 \\ \hline & 1/4 & 3/4 \end{array} \tag{4.4}
Generic second-order method
\def\arraystretch{1.5} \begin{array}{c|cc} 0 & 0 & 0 \\ \alpha & \alpha & 0 \\ \hline\\[-1em] & 1-\dfrac{1}{2\alpha} & \dfrac{1}{2\alpha} \end{array} \qquad \alpha\in(0,1) \tag{4.5}
Runge-Kutta’s third-order method (RK3)
\def\arraystretch{1.5} \begin{array} {c|ccc} 0 & 0 & 0 & 0 \\ 1/2 & 1/2 & 0 & 0 \\ 1 & -1 & 2 & 0 \\ \hline & 1/6 & 2/3 & 1/6 \end{array} \tag{4.6}
Heun’s third-order method
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1/3 & 1/3 & 0 & 0 \\ 1/3 & 0 & 2/3 & 0 \\ \hline & 1/4 & 0 & 3/4 \end{array} \tag{4.7}
Ralston’s third-order method
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1/2 & 1/2 & 0 & 0 \\ 3/4 & 0 & 3/4 & 0 \\ \hline & 2/9 & 1/3 & 4/9 \end{array} \tag{4.8}
Strong Stability Preserving Runge-Kutta (SSPRK3)
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1/2 & 1/4 & 1/4 & 0 \\ \hline & 1/6 & 1/6 & 2/3 \end{array} \tag{4.9}
Generic third-order method
With \alpha \neq 0, 2/3, 1:
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ \alpha & \alpha & 0 & 0 \\ 1 & 1+\dfrac{1-\alpha}{\alpha(3\alpha-2)} & -\dfrac{1-\alpha}{\alpha(3\alpha-2)} & 0 \\[1em] \hline \\[-1em] & \dfrac{1}{2}-\dfrac{1}{6\alpha} & \dfrac{1}{6\alpha(1-\alpha)} & \dfrac{2-3\alpha}{6\alpha(1-\alpha)} \end{array} \tag{4.10}
Classical Runge-Kutta’s method (RK4)
The classical fourth-order Runge-Kutta:
\def\arraystretch{1.5} \begin{array} {c|cccc} 0 & 0 & 0 & 0 & 0 \\ 1/2 & 1/2 & 0 & 0 & 0 \\ 1/2 & 0 & 1/2 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ \hline & 1/6 & 2/6 & 2/6 & 1/6 \end{array} \tag{4.11}
Ralston’s fourth-order method
It is a second-order method with two stages, designed to minimize truncation error.
\def\arraystretch{1.5} \begin{array}{c|ccccc} 0 & 0 & 0 & 0 & 0 \\ 2/5 & 2/5 & 0 & 0 & 0 \\ c_3 & a_{3,1} & a_{3,2} & 0 & 0 \\ 1 & a_{4,1} & a_{4,2} & a_{4,3} & 0 \\ \hline\\ & b_1 & b_2 & b_3 & b_4 \end{array} \tag{4.12}
where
\begin{aligned} b_1 & = \dfrac{263}{1812}+\dfrac{2}{151}\sqrt{5}\\ b_2 & = \dfrac{125}{3828}-\dfrac{250}{957}\sqrt{5} \\ b_3 & = \dfrac{3426304}{5924787}+\dfrac{553984}{1974929}\sqrt{5} \\ b_4 & = \dfrac{10}{41}-\dfrac{4}{123}\sqrt{5} \\ c_3 & = \dfrac{7}{8}-\dfrac{3}{16}\sqrt{5} \\ a_{3,1} & = \dfrac{357}{256}\sqrt{5}-\dfrac{2889}{1024} \\ a_{4,1} & = \dfrac{1047}{3020}\sqrt{5}-\dfrac{673}{1208} \\ a_{3,2} &= \dfrac{3785}{1024}-\dfrac{405}{256}\sqrt{5} \\ a_{4,2} &= -\dfrac{975}{2552}-\dfrac{1523}{1276}\sqrt{5} \\ a_{4,3} &= \dfrac{93408}{48169}+\dfrac{203968}{240845}\sqrt{5} \end{aligned}
3/8-rule fourth-order method
It is a second-order method with two stages that minimizes truncation error.
\def\arraystretch{1.5} \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0 \\ 1/3 & 1/3 & 0 & 0 & 0 \\ 2/3 & -1/3 & 1 & 0 & 0 \\ 1 & 1 & -1 & 1 & 0 \\ \hline & 1/8 & 3/8 & 3/8 & 1/8 \end{array}
4.2 Implicit Runge-Kutta Methods
The methods discussed so far are all explicit Runge-Kutta methods. However, when dealing with stiff equations, it is advisable to use implicit Runge-Kutta methods.
The generic implicit Runge-Kutta method is a one-step method characterized by the following representation:
\begin{aligned} \bm{z}_{k+1} &= \bm{z}_k + \sum_{i=1}^{s}b_i\bm{K}_i \\ &\left\{\begin{aligned} \bm{K}_1 &= h\bm{f}\Big(t_k+hc_1,\bm{z}_k + \displaystyle\sum_{j=1}^{s}a_{1j}\bm{K}_j\Big) \\[0.5em] \bm{K}_2 &= h\bm{f}\Big(t_k+hc_2,\bm{z}_k + \displaystyle\sum_{j=1}^{s}a_{2j}\bm{K}_j\Big) \\[0.5em] \vdots \\[0.5em] \bm{K}_s &= h\bm{f}\Big(t_k+hc_s,\bm{z}_k + \displaystyle\sum_{j=1}^{s}a_{sj}\bm{K}_j\Big) \end{aligned} \right. \end{aligned} \tag{4.13}
The butcher tableau for Equation 4.13 will be then of the form:
\def\arraystretch{1.5} \begin{array} {c|cccc} c_1 & a_{11} & a_{12} & \ldots & a_{1s} \\ c_2 & a_{21} & a_{22} & \ldots & a_{2s} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ c_s & a_{s1} & a_{s2} & \ldots & a_{ss} \\ \hline & b_1 & b_2 & \ldots & b_ss \end{array}
Implicit Euler method
The Implicit Euler Method Represented Using a Butcher Tableau:
\def\arraystretch{1.5} \begin{array}{c|c} 1 & 1 \\ \hline & 1 \end{array} \tag{4.14}
Implicit Midpoint Method
The implicit midpoint method is an adaptation of the explicit Euler method, often referred to as the Collatz method. It is a second-order method characterized by two stages, providing enhanced accuracy for solving ordinary differential equations.
\def\arraystretch{1.5} \begin{array}{c|c} 1/2 & 1/2 \\ \hline & 1 \end{array} \tag{4.15}
Crank-Nicolson method
The Crank-Nicolson method corresponds to the implicit trapezoidal rule.
\def\arraystretch{1.5} \begin{array}{c|cc} 0 & 0 & 0 \\ 1 & 1/2 & 1/2 \\ \hline & 1/2 & 1/2 \end{array} \tag{4.16}
Gauss-Legendre methods
These methods are based on the points of Gauss-Legendre quadrature.
Gauss-Legendre second-order method
\def\arraystretch{3} \begin{array}{c|cc} \dfrac{1}{2}-\dfrac{\sqrt3}{6} & \dfrac{1}{4} & \dfrac{1}{4}-\dfrac{\sqrt3}{6} \\ \dfrac{1}{2}+\dfrac{\sqrt3}{6} & \dfrac{1}{4}+\dfrac{\sqrt3}{6} & \dfrac{1}{4} \\ \hline & \dfrac{1}{2} & \dfrac{1}{2} \\ & \dfrac{1}{2}+\dfrac{\sqrt{3}}{2} & \dfrac{1}{2}-\dfrac{\sqrt{3}}{2} \end{array} \tag{4.17}
Gauss-Legendre sixth-order method
\def\arraystretch{3} \begin{array}{c|ccc} \dfrac{1}{2} - \dfrac{\sqrt{15}}{10} & \dfrac{5}{36} & \dfrac{2}{9}- \dfrac{\sqrt{15}}{15} & \dfrac{5}{36} - \dfrac{\sqrt{15}}{30} \\ \dfrac{1}{2} & \dfrac{5}{36} + \dfrac{\sqrt{15}}{24} & \dfrac{2}{9} & \dfrac{5}{36} - \dfrac{\sqrt{15}}{24} \\ \dfrac{1}{2} + \dfrac{\sqrt{15}}{10} & \dfrac{5}{36} + \dfrac{\sqrt{15}}{30} & \dfrac{2}{9} + \dfrac{\sqrt{15}}{15} & \dfrac{5}{36} \\[1em] \hline & \dfrac{5}{18} & \dfrac{4}{9} & \dfrac{5}{18} \\ & -\dfrac{5}{6} & \dfrac{8}{3} & -\dfrac{5}{6} \end{array} \tag{4.18}
Lobatto methods
Lobatto IIIA fourth-order method
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1/2 & 5/24& 1/3 & -1/24 \\ 1 & 1/6 & 2/3 & 1/6 \\ \hline & 1/6 & 2/3 & 1/6 \end{array} \tag{4.19}
Lobatto IIIB fourth-order method
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 1/6 & -1/6 & 0 \\ 1/2 & 1/6 & 1/3 & 0 \\ 1 & 1/6 & 5/6 & 0 \\ \hline & 1/6 & 2/3 & 1/6 \end{array} \tag{4.20}
Lobatto IIIC fourth-order method
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 1/6 & -1/3 & 1/6 \\ 1/2 & 1/6 & 5/12 & -1/12 \\ 1 & 1/6 & 2/3 & 1/6 \\ \hline & 1/6 & 2/3 & 1/6 \end{array} \tag{4.21}
Lobatto IIIC^\star fourth-order method
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1/2 & 1/4 & 1/4 & 0 \\ 1 & 0 & 1 & 0 \\ \hline & 1/6 & 2/3 & 1/6 \end{array} \tag{4.22}
Radau methods
Radau IA fifth-order method
\def\arraystretch{3} \begin{array}{c|ccc} 0 & \dfrac{1}{9} & \dfrac{-1 - \sqrt{6}}{18} & \dfrac{-1 + \sqrt{6}}{18} \\ \dfrac{3}{5} - \dfrac{\sqrt{6}}{10} & \dfrac{1}{9} & \dfrac{11}{45} + \dfrac{7\sqrt{6}}{360} & \dfrac{11}{45} - \dfrac{43\sqrt{6}}{360}\\ \dfrac{3}{5} + \dfrac{\sqrt{6}}{10} & \dfrac{1}{9} & \dfrac{11}{45} + \dfrac{43\sqrt{6}}{360} & \dfrac{11}{45} - \dfrac{7\sqrt{6}}{360} \\ \hline & \dfrac{1}{9} & \dfrac{4}{9} + \dfrac{\sqrt{6}}{36} & \dfrac{4}{9} - \dfrac{\sqrt{6}}{36} \end{array} \tag{4.23}
Radau II fifth-order method
\def\arraystretch{3} \begin{array}{c|ccc} \dfrac{2}{5} - \dfrac{\sqrt{6}}{10} & \dfrac{11}{45} - \dfrac{7\sqrt{6}}{360} & \dfrac{37}{225} - \dfrac{169\sqrt{6}}{1800} & -\dfrac{2}{225} + \dfrac{\sqrt{6}}{75} \\ \dfrac{2}{5} + \dfrac{\sqrt{6}}{10} & \dfrac{37}{225} + \dfrac{169\sqrt{6}}{1800} & \dfrac{11}{45} + \dfrac{7\sqrt{6}}{360} & -\dfrac{2}{225} - \dfrac{\sqrt{6}}{75}\\ 1 & \dfrac{4}{9} - \dfrac{\sqrt{6}}{36} & \dfrac{4}{9} + \dfrac{\sqrt{6}}{36} & \dfrac{1}{9} \\ \hline & \dfrac{4}{9} - \dfrac{\sqrt{6}}{36} & \dfrac{4}{9} + \dfrac{\sqrt{6}}{36} & \dfrac{1}{9} \end{array} \tag{4.24}