Runge-Kutta Method as a Multi-step Scheme
Previously, we introduced the forward and backward Euler methods. They are quite crude and only operate on a single time step basis, with a relatively large error. The Runge-Kutta (RK4) method is a significant improvement on them, by further splitting a single time step into four smaller sub-steps, hence becoming a multi-step scheme. As before, assume that the first-order ODE in question has the general form of
\begin{equation}
\frac{dy}{dt} = f(t,y) \tag{1}
\end{equation}
For a chosen time step size \( h = \Delta t \), the RK4 scheme then sequentially computes four quantities as follows.
\begin{align}
\left\{\begin{aligned}
k_1 &= f(t, y(t)) \\
k_2 &= f(t+\frac{h}{2}, y(t)+k_1\frac{h}{2}) \\
k_3 &= f(t+\frac{h}{2}, y(t)+k_2\frac{h}{2}) \\
k_4 &= f(t+h, y(t)+k_3h) \\
\end{aligned}\right. \tag{2}
\end{align}
They are the estimated slopes/tendencies \( dy/dt \) of the unknown at different points: \( k_1 \) is the slope at the beginning of the time step given by \(y\); \( k_2, k_3 \) are the slopes at the middle of the time step using \( y \) extrapolated from \( k_1, k_2 \); \(k_4\) is the slope at the end of the time step using \( y \) extrapolated from \(k_3\) instead. Then the value of \(y\) at the next time step is computed by adding a weighted average of \(k_1, k_2, k_3, k_4\):
\begin{equation}
y(t + \Delta t) = y(t) + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) \tag{3}
\end{equation}
This can be compactly encapsulated in the so-called Butcher tableau:
The left column represents the fractional time taken within a major time step (the first argument of \( k_j \) in (2)), and the bottom row indicates the contribution of each \( k_j \) to the total increment of \(y\) over a major time step as indicated by (3). Meanwhile, the square matrix records the coefficients of all \( k_i \) needed to approximate \( y \) for the second argument of \( k_j \) during the sub-steps in (2).
The RK4 method is a generalization of the so-called Simpson’s rule for approximating a definite integral. It is a fourth-order explicit method and has a truncation error on the order of the fifth power of the time step size, i.e. \( O((\Delta t)^5) \). Compared to the two Euler methods, which are of \( O((\Delta t)^2) \), as the resolution of time step gets finer, \( \Delta t \) becomes smaller, the RK4 method will have a much higher accuracy.
Generalizing the RK4 method to a system of first-order ODEs is not hard. Just like the Euler methods, we replace \( f , y \) by vector-valued \( \textbf{f}, \textbf{y} \) in (2) and (3) and work in the same way.
Example
We will redo the same example as before and compare the performance. That is, to numerically integrate
\begin{equation}
\frac{dy}{dt} = -\alpha y \tag{4}
\end{equation}
with \( \alpha = 0.2 \), \( y(0) = 1 \), \( \Delta t = 0.1 \), by RK4.
Since the RK4 method requires four substeps for a single time step, we will just show the calculations for the first three time steps \( t = 0.1, 0.2, 0.3 \) and the expected output at \(t = 1.0\) to keep it short.
For \( t = 0.1 \), (2) and (3) give
\begin{align}
\begin{aligned}
k_1 &= -0.2(1) = -0.2 \\
k_2 &= -0.2(1+(-0.2)(0.1/2)) = -0.198 \\
k_3 &= -0.2(1+(-0.198)(0.1/2)) = -0.19802\\
k_4 &= -0.2(1+(-0.19802)(0.1)) = -0.19604 \\
\end{aligned} \tag{4}
\end{align}
\begin{equation}
y(0.1) = 1 + \frac{0.1}{6}((-0.2) + 2(-0.198) + 2(-0.19802) + (-0.19604)) = 0.98020 \tag{5}
\end{equation}
At \( t = 0.2 \),
\begin{align}
\begin{aligned}
k_1 &= -0.2(0.98020) = -0.19604 \\
k_2 &= -0.2(0.98020+(-0.19604)(0.1/2)) = -0.19408 \\
k_3 &= -0.2(0.98020+(-0.19408)(0.1/2)) = -0.19410 \\
k_4 &= -0.2(0.98020+(-0.19410)(0.1)) = -0.19216 \\
\end{aligned} \tag{6}
\end{align}
\begin{equation}
y(0.2) = 0.98020 + \frac{0.1}{6}((-0.19604) + 2(-0.19408) + 2(-0.19410) + (-0.19216)) = 0.96079 \tag{7}
\end{equation}
At \( t = 0.3 \),
\begin{align}
\begin{aligned}
k_1 &= -0.2(0.96079) = -0.19216 \\
k_2 &= -0.2(0.96079+(-0.19216)(0.1/2)) = -0.19024 \\
k_3 &= -0.2(0.96079+(-0.19024)(0.1/2)) = -0.19026 \\
k_4 &= -0.2(0.96079+(-0.19026)(0.1)) = -0.18835 \\
\end{aligned} \tag{8}
\end{align}
\begin{equation}
y(0.3) = 0.96079 + \frac{0.1}{6}((-0.19216) + 2(-0.19024) + 2(-0.19026) + (-0.18835)) = 0.94176 \tag{9}
\end{equation}
Repeat the procedure until \( t = 1.0 \), and the numerical answer should be approximately \( y(1.0) = 0.81873 \). In fact, compared to the true solution \( e^{-0.2} \), it can be checked that it is correct up to the \(10\)-th decimal place.
Exercise
Numerically integrate
\begin{equation}
\frac{dy}{dt} = -e^{-t}y \tag{10}
\end{equation}
by the RK4 method with a time step size of \( \Delta t = 0.2 \) until \( t = 2.0 \), where the I.C. is \( y(0) = 1 \).
Answer
Again, to avoid congestion, we will just note down the values of \( k_1, k_2, k_3, k_4 \) and \( y \) that should be obtained at each time step for reference.
| \(t\) | \(y\) | \( k_1 \) | \( k_2 \) | \( k_3 \) | \( k_4 \) |
| \( 0.0 \) | \( 1.0 \) | \( -1.0 \) | \( -0.81435 \) | \( -0.83115 \) | \( -0.68263 \) |
| \( 0.2 \) | \( 0.83421 \) | \( -0.68299 \) | \( -0.56740 \) | \( -0.57597 \) | \( -0.48197 \) |
| \( 0.4 \) | \( 0.71916 \) | \( -0.48206 \) | \( -0.40695 \) | \( -0.41151 \) | \( -0.34951 \) |
| \( 0.6 \) | \( 0.63687 \) | \( -0.34952 \) | \( -0.29890 \) | \( -0.30142 \) | \( -0.25908 \) |
| \( 0.8 \) | \( 0.57656 \) | \( -0.25907 \) | \( -0.22388 \) | \( -0.22531 \) | \( -0.19553 \) |
| \( 1.0 \) | \( 0.53146 \) | \( -0.19551 \) | \( -0.17040 \) | \( -0.17124 \) | \( -0.14976 \) |
| \( 1.2 \) | \( 0.49718 \) | \( -0.14975 \) | \( -0.13142 \) | \( -0.13192 \) | \( -0.11610 \) |
| \( 1.4 \) | \( 0.47076 \) | \( -0.11609 \) | \( -0.10245 \) | \( -0.10276 \) | \( -0.09090 \) |
| \( 1.6 \) | \( 0.45018 \) | \( -0.09090 \) | \( -0.08058 \) | \( -0.08077 \) | \( -0.07174 \) |
| \( 1.8 \) | \( 0.43400 \) | \( -0.07174 \) | \( -0.06384 \) | \( -0.06396 \) | \( -0.05701 \) |
At \( t = 2.0 \), the end numerical result will be \( 0.43400 + (0.2/6)((-0.07174)+2(-0.06384)+2(-0.06396)+(-0.05701)) = 0.42119 \).
Leave a Reply