Section 9.5 RK4 Method
A large and important family of numerical methosds for ODEs is the Runge-Kutta one. The basic idea of these methods is to generalize the idea that led us to the Heun method: we evaluate numerically the slope at several points between \(t_n\) and \(t_{n+1}\) and then we use, in the evaluation of \(x_{n+1}\text{,}\) some weighted average of these slopes. If you want to know more details about all these coefficients and their relations, you can get them all from the wikipedia page linked above. Here we illustrate a single but significative example called RK4. As the name suggests, our final slope will be the average of four slopes. Suppose we are at the \(n\)-th step of the algorithm and we reached the poin \((t_n,x_n).\) The first slope is just the slope at the starting point:
\begin{equation*}
k_1 = f(t_n,x_n).
\end{equation*}
The second slope is the slope at the point we reach by moving by half time step with Euler's method:
\begin{equation*}
k_2 = f(t_n+\frac{h}{2},x_n+\frac{h}{2}k_1).
\end{equation*}
The third slope is the slope at the point we reach by moving by half time step with Euler's method, but this time using the slope at the point we reached with the standard Euler's method:
\begin{equation*}
k_3 = f(t_n+\frac{h}{2},x_n+\frac{h}{2}k_2).
\end{equation*}
The fourth slope is the slope at the point we reach by moving by a full time step with Euler's method, this time using the slope at the point we reached in the third step:
\begin{equation*}
k_4 = f(t_n+h,x_n+hk_3).
\end{equation*}
Finally, the RK4 rule is now
\begin{equation*}
x_{n+1} = x_n + h\frac{k_1+2k_2+2k_3+k_4}{6}.
\end{equation*}
If you are puzzled by where do the weights \(1,2,2,1\) of the slopes come from, consider the case when \(f\) does not depend on \(x\text{.}\) In this case clearly \(k_2=k_3\) and RK4 reduces to
\begin{equation*}
x_{n+1} = x_n + h\frac{f(t_n)+4f\left(\frac{t_n+t_{n+1}}{2}\right)+f(t_{n+1})}{6}.
\end{equation*}
Yes, this is exactly Simpson's rule for numerical integration! So RK4 is nothing but the adaptation to ODEs of Simpson's rule. Since in general \(f\) does depend also on \(x\text{,}\) the middle term \(f(\frac{t_n+t_{n+1}}{2})\) splits "evenly" into two terms, both evaluated at the same \(t\) but at two different values of \(x\text{.}\) This particular implementation subdivides the "time" interval \([0,1.5]\) in 10 subintervals, namely \(h=0.15\text{.}\) The code saves all intermediate steps and plots them together with the corresponding values of the exact solution \(x(t)=\sqrt{1.25+\sin(\pi t)}\text{.}\)
One can prove that RK4 method`s local error is \(O(h^5)\) and global error is \(O(h^4)\text{.}\) Below we verify this fact numerically: