Burden & Faires, 9th Ed. - §5.4–5.5
Chapter 12: Runge-Kutta Methods
Midpoint & Heun's methods · Classical RK4 · Systems of ODEs · A first look at stability
Chapter 11 ended on a specific complaint: Taylor methods buy accuracy with analytic derivatives of $f$ that quickly become unmanageable to compute by hand. Runge-Kutta methods solve exactly that problem - they match a Taylor method's order of accuracy using nothing but extra evaluations of $f$ itself, at cleverly chosen intermediate points within each step. The classical fourth-order Runge-Kutta method (RK4) that closes this chapter is, practically speaking, the default ODE solver most scientific software reaches for first.
Section 5.4
Second-Order Runge-Kutta
Idea - Sample the Slope Twice per Step
Instead of using the slope at just $(t_i,w_i)$ (Euler), evaluate $f$ at a second, well-chosen point within the step and average. Two standard choices:
$$\textbf{Midpoint Method:}\quad w_{i+1} = w_i + h\,f\!\left(t_i+\tfrac h2,\; w_i+\tfrac h2 f(t_i,w_i)\right)$$
$$\textbf{Heun's Method:}\quad k_1=f(t_i,w_i),\quad k_2=f(t_i+h,w_i+hk_1), \quad w_{i+1}=w_i+\frac h2(k_1+k_2)$$
Both match the Taylor series through the $h^2$ term without ever computing $f_t$ or $f_y$ analytically - the extra function evaluation does the work a derivative would have done. Both are $O(h^2)$ global error, same order as Chapter 11's Taylor-2 method.
📘 Example 12.1 - Heun's Method on the Chapter 11 Benchmark
Apply Heun's Method to $y'=y-t^2+1,\;y(0)=0.5$ on $[0,2]$, $h=0.2$ - the same problem used throughout Chapter 11.
Solution
Final error $0.0724$ vs. Euler's $0.4397$ (Chapter 11, Example 11.1) - better by a factor of $6$, using the same $h$ but twice the function evaluations per step (2 vs. 1). No analytic derivative of $f$ was needed anywhere, unlike the Taylor-2 method that achieved similar accuracy in Chapter 11's Example 11.2.
| $t_i$ | $w_i$ (Heun) | true $y(t_i)$ | error |
|---|---|---|---|
| 0.2 | 0.826000 | 0.829299 | 0.003299 |
| 0.6 | 1.637242 | 1.648941 | 0.011698 |
| 1.0 | 2.617688 | 2.640859 | 0.023171 |
| 1.4 | 3.693686 | 3.732400 | 0.038714 |
| 2.0 | 5.233055 | 5.305472 | 0.072417 |
Section 5.4 continued
Classical Fourth-Order Runge-Kutta (RK4)
The RK4 Formula
$$k_1 = f(t_i,w_i)$$
$$k_2 = f\!\left(t_i+\tfrac h2,\; w_i+\tfrac h2 k_1\right)$$
$$k_3 = f\!\left(t_i+\tfrac h2,\; w_i+\tfrac h2 k_2\right)$$
$$k_4 = f(t_i+h,\; w_i+hk_3)$$
$$w_{i+1} = w_i + \frac h6(k_1+2k_2+2k_3+k_4)$$
Four function evaluations per step, matching the Taylor series through the $h^4$ term: local truncation error $O(h^5)$, global error $O(h^4)$. $k_1$ is the slope at the start; $k_2,k_3$ are two independent estimates of the midpoint slope; $k_4$ is the slope at the endpoint - the $\frac16(1,2,2,1)$ weighting is a Simpson's-Rule-like weighted average (Chapter 9) of these four slope estimates.
Same IVP ($y'=y-t^2+1,\,y(0)=0.5$, true solution gray), same $h=0.2$, three different methods. RK4's curve visually overlaps the true solution almost perfectly at this scale - its error is roughly 4000$\times$ smaller than Euler's at identical cost per step in code complexity (though 4$\times$ the function evaluations).
📘 Example 12.2 - RK4's Accuracy, Quantified
Apply RK4 to the same benchmark IVP, $h=0.2$, and compare convergence rate to Euler's as $h$ shrinks.
Solution
Halving $h$ divides RK4's error by $\approx16=2^4$ (confirming $O(h^4)$) versus Euler's factor of $\approx2$ ($O(h)$). At $h=0.2$, RK4 is already $4000\times$ more accurate than Euler - and the gap widens every time $h$ shrinks further.
| $t_i$ | $w_i$ (RK4) | error |
|---|---|---|
| 0.2 | 0.829293 | $5\times10^{-6}$ |
| 1.0 | 2.640823 | $3.6\times10^{-5}$ |
| 2.0 | 5.305363 | $1.09\times10^{-4}$ |
| $h$ | RK4 error at $t=2$ | Euler error at $t=2$ (Ch. 11) |
|---|---|---|
| 0.2 | $1.09\times10^{-4}$ | $4.40\times10^{-1}$ |
| 0.1 | $6.99\times10^{-6}$ | $2.42\times10^{-1}$ |
| 0.05 | $4.42\times10^{-7}$ | $1.27\times10^{-1}$ |
Section 5.4–5.5
Systems of ODEs & Stability
Extending to Systems
A system $\mathbf{y}' = \mathbf{f}(t,\mathbf{y})$, $\mathbf{y}(a)=\boldsymbol\alpha$, with $\mathbf y=(y^{(1)},\ldots,y^{(m)})$, uses exactly the RK4 formula above with every quantity replaced by a vector - compute $\mathbf k_1,\mathbf k_2,\mathbf k_3,\mathbf k_4$ as vectors, all components advancing together with the same $h$. This is also how a higher-order ODE $y^{(m)}=g(t,y,y',\ldots,y^{(m-1)})$ is solved: introduce $u_1=y,\,u_2=y',\ldots,u_m=y^{(m-1)}$ to rewrite it as a first-order system.
📘 Example 12.3 - RK4 on a System: the Simple Harmonic Oscillator
Solve $y''=-y,\;y(0)=0,\,y'(0)=1$ (true solution $y=\sin t$) by converting to the system $u_1'=u_2,\;u_2'=-u_1$, $u_1(0)=0,\,u_2(0)=1$, using RK4 with $h=0.1$.
Solution
Agreement to 5–6 decimal places at every step - RK4 tracks both the position and the derivative of a genuinely oscillating system accurately, something Euler's Method (Chapter 11) would visibly fail at over any appreciable number of periods, spiraling outward due to its systematic energy drift.
| $t$ | $u_1$ (RK4) | true $\sin t$ | $u_2$ (RK4) | true $\cos t$ |
|---|---|---|---|---|
| 0.1 | 0.099833 | 0.099833 | 0.995004 | 0.995004 |
| 0.3 | 0.295520 | 0.295520 | 0.955337 | 0.955336 |
| 0.5 | 0.479425 | 0.479426 | 0.877583 | 0.877583 |
A First Word on Stability
Every method in this chapter has a stability region - a range of $h\lambda$ (where $\lambda$ relates to $\partial f/\partial y$) for which errors stay bounded rather than growing without limit step to step. For "stiff" problems - systems with widely different time scales - even RK4 can require an impractically small $h$ to remain stable, motivating implicit methods that lie beyond this course but build directly on everything developed here: the same $k_i$ structure, solved implicitly rather than explicitly.
Practice
Practice Problem
Try this problem yourself before expanding the solution below.
✏️ Practice Problem 12.1 - Two Steps of RK4 for $y'=x+y$
Solve the IVP $y'=x+y$, $y(0)=1$ using classical RK4 with $h=0.2$, for two steps. The exact solution is $y(x)=-x-1+2e^x$; report the global error at each step.
Solution
Step 1 ($x_0=0,\,w_0=1$):
$$k_1=1.000000,\ k_2=1.200000,\ k_3=1.220000,\ k_4=1.444000$$
$$w_1 = w_0+\frac{h}{6}(k_1+2k_2+2k_3+k_4) = 1.242800.$$
Exact: $y(0.2)=1.242806$, error $\approx5.5\times10^{-6}$.
Step 2 ($x_1=0.2,\,w_1=1.242800$): $$k_1=1.442800,\ k_2=1.687080,\ k_3=1.711508,\ k_4=1.985102$$ $$w_2 = w_1+\frac{h}{6}(k_1+2k_2+2k_3+k_4) = 1.583636.$$ Exact: $y(0.4)=1.583649$, error $\approx1.3\times10^{-5}$. Both steps land within $10^{-5}$ of the true solution using only $h=0.2$ - the same accuracy Euler's Method would need a vastly smaller step size to reach, since RK4's error is $O(h^4)$ per step versus Euler's $O(h)$.
Step 2 ($x_1=0.2,\,w_1=1.242800$): $$k_1=1.442800,\ k_2=1.687080,\ k_3=1.711508,\ k_4=1.985102$$ $$w_2 = w_1+\frac{h}{6}(k_1+2k_2+2k_3+k_4) = 1.583636.$$ Exact: $y(0.4)=1.583649$, error $\approx1.3\times10^{-5}$. Both steps land within $10^{-5}$ of the true solution using only $h=0.2$ - the same accuracy Euler's Method would need a vastly smaller step size to reach, since RK4's error is $O(h^4)$ per step versus Euler's $O(h)$.
Connections Across the Course
- RK4's $\frac16(k_1+2k_2+2k_3+k_4)$ weighting is structurally the same weighted-average idea as Simpson's Rule (Chapter 9) - both average multiple sample points with weights chosen to cancel low-order error terms.
- RK methods trade analytic derivatives for extra function evaluations exactly the way the Secant Method (Chapter 4) traded $f'$ for an extra evaluation of $f$ itself.
- Systems of ODEs solved here are where Chapter 13's vector and matrix norms become essential - bounding error growth for a system needs a norm on $\mathbf y$, not just absolute value.