Chapter 8 · Section 8.6

8.6 Numerical Methods III: Runge–Kutta

Heun's method (Section 8.5) improved on Euler's by sampling the slope twice per step instead of once, and the reward was a jump from first to second order. A natural question follows: what if we sample the slope at even more places inside the step and combine the samples cleverly? Around 1900 the German mathematicians Carl Runge and Wilhelm Kutta answered it. The family of methods that carries their names contains the most widely used general-purpose ODE solver of all time, the classical fourth-order Runge–Kutta method, usually called simply RK4.

As before, we solve

y^{'} = f (x, y), y (x_{0}) = y_{0},

on the grid $x_{n} = x_{0} + nh$ , writing $y (x_{n})$ for the exact value and $y_{n}$ for its approximation.

The classical RK4 formulas

Definition 8.6.1 (Classical fourth-order Runge–Kutta method)

Given $y_{n}$ at $x_{n}$ , compute four slopes

k_{1} k_{2} k_{3} k_{4} = f (x_{n}, y_{n}), = f (x_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{1}), = f (x_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{2}), = f (x_{n} + h, y_{n} + h k_{3}),

and then step with their weighted average:

y_{n + 1} = y_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}) .

Each step costs four evaluations of $f$ . Notice the structure: every slope after the first is evaluated at a point reached by moving from $(x_{n}, y_{n})$ with the previous slope. The method is explicit: each line uses only quantities already computed.

Another common notation

Many books and programs put the factor $h$ inside the stages: $K_{1} = h f (x_{n}, y_{n})$ , $K_{2} = h f (x_{n} + \frac{h}{2}, y_{n} + \frac{1}{2} K_{1})$ , and so on, with $y_{n + 1} = y_{n} + \frac{1}{6} (K_{1} + 2 K_{2} + 2 K_{3} + K_{4})$ . This is the same method; only the bookkeeping differs ( $K_{i} = h k_{i}$ ). In this course the $k_{i}$ are always slopes.

What the four slopes mean

It helps to read the formulas as a small exploration of the direction field across one step.

$k_{1}$ : slope at the start. This is the Euler slope at $(x_{n}, y_{n})$ .
$k_{2}$ : first look at the midpoint. Follow $k_{1}$ for half a step to the point $(x_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{1})$ and read the slope there.
$k_{3}$ : second, better look at the midpoint. The point used for $k_{2}$ was only an Euler guess. Go back to the start, follow the improved slope $k_{2}$ for half a step, and read the slope at this new midpoint estimate.
$k_{4}$ : slope at the end. Follow $k_{3}$ for a full step to $(x_{n} + h, y_{n} + h k_{3})$ and read the slope there.
Weighted average. The end slopes get weight $1$ each and the two midpoint slopes weight $2$ each; the weights add up to $6$ , hence the $h /6$ .

Figure 8.6.1 shows the four samples for $y^{'} = y$ , $y (0) = 1$ , with the large step $h = 1$ . Here $k_{1} = 1$ , $k_{2} = f (0.5, 1.5) = 1.5$ , $k_{3} = f (0.5, 1.75) = 1.75$ , $k_{4} = f (1, 2.75) = 2.75$ , and

y_{1} = 1 + \frac{1}{6} (1 + 3 + 3.5 + 2.75) = 1 + \frac{10.25}{6} = 2.708333,

only $0.010$ below $e = 2.718282$ . With the same step, Euler gives $2$ and Heun gives $2.5$ .

Figure 8.6.1. The four slopes of one RK4 step for

y^{'} = y

y (0) = 1

h = 1

. Dotted rays show how each sample point is reached from the start. The weighted step (orange) ends at

2.7083

, almost on the exact curve.

Why the weights 1, 2, 2, 1?

Look at the special case where $f$ depends on $x$ only. Then the $y$ -arguments play no role, $k_{2} = k_{3} = f (x_{n} + \frac{h}{2})$ , and the RK4 step becomes

y_{n + 1} = y_{n} + \frac{h}{6} [f (x_{n}) + 2 f (x_{n} + \frac{h}{2}) + 2 f (x_{n} + \frac{h}{2}) + f (x_{n} + h)] = y_{n} + \frac{h}{6} [f (x_{n}) + 4 f (x_{n} + \frac{h}{2}) + f (x_{n} + h)] .

The bracket is Simpson's rule for $\int_{x_{n}}^{x_{n} + h} f (x) d x$ , using the two ends and the midpoint. So RK4 is to Simpson's rule what Heun's method is to the trapezoid rule. Splitting the midpoint weight $4$ into $2 + 2$ between two different midpoint estimates is what lets the method handle a slope that also depends on $y$ .

Method in a nutshell: one RK4 step

$k_{1} = f (x_{n}, y_{n})$ .
$k_{2} = f (x_{n} + h /2, y_{n} + (h /2) k_{1})$ .
$k_{3} = f (x_{n} + h /2, y_{n} + (h /2) k_{2})$ (note: $k_{2}$ , not $k_{1}$ ).
$k_{4} = f (x_{n} + h, y_{n} + h k_{3})$ (full step, slope $k_{3}$ ).
$y_{n + 1} = y_{n} + (h /6) (k_{1} + 2 k_{2} + 2 k_{3} + k_{4})$ ; set $x_{n + 1} = x_{n} + h$ and repeat.

Common mistakes

(1) Using $h$ instead of $h /2$ in $k_{2}$ or $k_{3}$ . (2) Building $k_{3}$ from $k_{1}$ again instead of from $k_{2}$ . (3) Forgetting to advance $x$ in $k_{2}$ , $k_{3}$ , $k_{4}$ . (4) Dividing by $4$ instead of $6$ . A quick check: if $f$ is a constant $c$ , all four slopes equal $c$ and the step must be $y_{n} + h c$ .

Worked examples

We return to the test problem of Section 8.5,

y^{'} = y - x^{2}, y (0) = 1, exact solution y (x) = x^{2} + 2 x + 2 - e^{x} .

As before, we carry full precision and round only what we display.

Example 8.6.1

Use RK4 with $h = 0.2$ to approximate the solution of $y^{'} = y - x^{2}$ , $y (0) = 1$ , at $x = 0.2, 0.4, \dots, 1.0$ , and compare with the exact values.

Solution

First step in full. $x_{0} = 0$ , $y_{0} = 1$ , $h /2 = 0.1$ . $k_{1} k_{2} k_{3} k_{4} = f (0, 1) = 1 - 0 = 1, = f (0.1, 1 + 0.1 \cdot 1) = 1.1 - 0.01 = 1.09, = f (0.1, 1 + 0.1 \cdot 1.09) = 1.109 - 0.01 = 1.099, = f (0.2, 1 + 0.2 \cdot 1.099) = 1.2198 - 0.04 = 1.1798.$
Combine. $y_{1} = 1 + \frac{0.2}{6} (1 + 2 (1.09) + 2 (1.099) + 1.1798) = 1 + \frac{0.2}{6} (6.5578) = 1.218593.$ The exact value is $y (0.2) = 1.218597$ : the error is only $3.9 \times 1 0^{- 6}$ , compared with $0.002597$ for Heun and $0.018597$ for Euler with the same $h$ .

All steps.

$x_{n}$	$y_{n}$	$k_{1}$	$k_{2}$	$k_{3}$	$k_{4}$	$y_{n + 1}$
0.0	1.000000	1.000000	1.090000	1.099000	1.179800	1.218593
0.2	1.218593	1.178593	1.246453	1.253239	1.309241	1.468167
0.4	1.468167	1.308167	1.348984	1.353066	1.378780	1.737869
0.6	1.737869	1.377869	1.385656	1.386434	1.375156	2.014442
0.8	2.014442	1.374442	1.341886	1.338631	1.282168	2.281697

Errors.

$x_{n}$	exact	RK4	error
0.2	1.218597	1.218593	$3.9 \times 1 0^{- 6}$
0.4	1.468175	1.468167	$8.1 \times 1 0^{- 6}$
0.6	1.737881	1.737869	$1.2 \times 1 0^{- 5}$
0.8	2.014459	2.014442	$1.7 \times 1 0^{- 5}$
1.0	2.281718	2.281697	$2.1 \times 1 0^{- 5}$

With only five steps, RK4 is correct to about five significant digits on the whole interval.

The solver below shows every $k_{i}$ of every step. Start with the preset (Example 8.6.1), then try $h = 0.1$ : the error at $x = 1$ should drop by a factor close to $16$ . You can also enter the problem of Example 8.6.2.

Example 8.6.2 (A nonlinear equation)

Take one RK4 step with $h = 0.2$ for $y^{'} = - 2 x y^{2}$ , $y (0) = 1$ . Compare with the exact solution.

Solution

Exact solution. Separate the variables: $- y^{- 2} d y = 2 x d x$ , so $y^{- 1} = x^{2} + c$ . From $y (0) = 1$ , $c = 1$ , and $y = \frac{1}{1 + x ^{2}}$ . Thus $y (0.2) = 1/1.04 = 0.961538$ .
The slopes. With $f (x, y) = - 2 x y^{2}$ : $k_{1} k_{2} k_{3} k_{4} = f (0, 1) = 0, = f (0.1, 1 + 0.1 \cdot 0) = - 2 (0.1) (1)^{2} = - 0.2, = f (0.1, 1 + 0.1 (- 0.2)) = - 0.2 (0.98)^{2} = - 0.192080, = f (0.2, 1 + 0.2 (- 0.192080)) = - 0.4 (0.961584)^{2} = - 0.369858.$
Combine. $y_{1} = 1 + \frac{0.2}{6} (0 - 0.4 - 0.384160 - 0.369858) = 1 + \frac{0.2}{6} (- 1.154018) = 0.961533.$ The error is about $6 \times 1 0^{- 6}$ . Note that $k_{1} = 0$ : Euler's method would have predicted no change at all, $y_{1} = 1$ .

How accurate is RK4?

The analysis follows the pattern of Section 8.5. We expand the RK4 step and the exact solution in powers of $h$ and compare. For Heun's method they agreed up to $h^{2}$ . For RK4 they agree up to $h^{4}$ . The algebra is long, because every $k_{i}$ must be expanded with the two-variable Taylor formula up to fourth order, so we state the result and then check it on a case where the computation is short.

Theorem 8.6.1 (Order of the classical Runge–Kutta method)

Suppose $f$ has continuous partial derivatives up to order four near the solution curve, so that the solution has a continuous fifth derivative on $[x_{0}, b]$ . Then RK4 has local truncation error $O (h^{5})$ and its global error satisfies

∣ y (x_{n}) - y_{n} ∣ \leq C h^{4} (x_{0} \leq x_{n} \leq b)

with a constant $C$ independent of $h$ . RK4 is a fourth-order method: halving $h$ divides the error by about $2^{4} = 16$ .

A check on the test equation $y^{'} = λ y$

For $f (x, y) = λ y$ ( $λ$ constant) put $z = hλ$ . The slopes are

k_{1} k_{2} k_{3} k_{4} = λ y_{n}, = λ (y_{n} + \frac{h}{2} k_{1}) = λ y_{n} (1 + \frac{z}{2}), = λ (y_{n} + \frac{h}{2} k_{2}) = λ y_{n} (1 + \frac{z}{2} + \frac{z ^{2}}{4}), = λ (y_{n} + h k_{3}) = λ y_{n} (1 + z + \frac{z ^{2}}{2} + \frac{z ^{3}}{4}) .

Add them with weights $1, 2, 2, 1$ :

k_{1} + 2 k_{2} + 2 k_{3} + k_{4} = λ y_{n} (6 + 3 z + z^{2} + \frac{z ^{3}}{4}) .

Multiply by $h /6$ (remember $hλ = z$ ) and add $y_{n}$ :

y_{n + 1} = (1 + z + \frac{z ^{2}}{2} + \frac{z ^{3}}{6} + \frac{z ^{4}}{24}) y_{n} .

The exact solution satisfies $y (x_{n + 1}) = e^{z} y (x_{n})$ , and $e^{z} = 1 + z + \frac{z ^{2}}{2} + \frac{z ^{3}}{6} + \frac{z ^{4}}{24} + \frac{z ^{5}}{120} + \dots$ . RK4 reproduces the series exactly through $z^{4}$ and first misses at $z^{5} = h^{5} λ^{5}$ : a local error of order $h^{5}$ , as the theorem says. For comparison, Euler's multiplier is $1 + z$ and Heun's is $1 + z + \frac{z ^{2}}{2}$ (Exercise 4 of Section 8.5).

The Runge–Kutta family

Euler, Heun, the midpoint method and RK4 are all explicit Runge–Kutta methods: each step evaluates $f$ at a few trial points built from earlier slopes and takes a weighted average. With $s$ slope evaluations ("stages") one can reach order $s$ for $s = 1, 2, 3, 4$ . Beyond that the pattern breaks: order $5$ needs at least six stages. This is one reason why four stages and fourth order are such a popular compromise.

A fair comparison: Euler, Heun and RK4

Counting steps is unfair to Euler, since an RK4 step does four times as much work. The honest measure of cost is the number of evaluations of $f$ , because in real problems that is where the computing time goes. So let us give every method the same budget.

Example 8.6.3 (Equal cost)

For $y^{'} = y - x^{2}$ , $y (0) = 1$ , allow each method exactly $20$ evaluations of $f$ to reach $x = 1$ . Which method is most accurate?

Solution

Step sizes. Euler uses $1$ evaluation per step, so $20$ steps: $h = 0.05$ . Heun uses $2$ per step, so $10$ steps: $h = 0.1$ . RK4 uses $4$ per step, so $5$ steps: $h = 0.2$ .
Results at $x = 1$ (exact $y (1) = 2.281718$ ):

method $h$ $y_{n}$ at $x = 1$ error

Euler 0.05 2.264037 0.017681

Heun 0.1 2.277757 0.003961

RK4 0.2 2.281697 0.000021
Conclusion. For the same work, RK4 is about $840$ times more accurate than Euler and about $190$ times more accurate than Heun. The extra effort per step pays off many times over.

method	$h$	$y_{n}$ at $x = 1$	error
Euler	0.05	2.264037	0.017681
Heun	0.1	2.277757	0.003961
RK4	0.2	2.281697	0.000021

Figure 8.6.2. Errors along

[0, 1]

for the three methods of Example 8.6.3, each using 20 slope evaluations. Note the logarithmic vertical axis: RK4 is almost three orders of magnitude better than Euler.

The advantage grows as we ask for more accuracy. The table lists the error at $x = 1$ for several step sizes; the last column shows how much the RK4 error shrinks each time $h$ is halved.

$h$	Euler	Heun	RK4	RK4 ratio
0.2	$6.8 \times 1 0^{- 2}$	$1.5 \times 1 0^{- 2}$	$2.1 \times 1 0^{- 5}$	–
0.1	$3.5 \times 1 0^{- 2}$	$4.0 \times 1 0^{- 3}$	$1.3 \times 1 0^{- 6}$	15.95
0.05	$1.8 \times 1 0^{- 2}$	$1.0 \times 1 0^{- 3}$	$8.2 \times 1 0^{- 8}$	16.02
0.025	$8.9 \times 1 0^{- 3}$	$2.5 \times 1 0^{- 4}$	$5.1 \times 1 0^{- 9}$	16.02

Figure 8.6.3. Error at

x = 1

against

h

on logarithmic axes. The three methods follow lines of slope

1

2

and

4

: the slope of each line is the order of the method.

Estimating the error without the exact solution

In practice we do not know the exact solution, so how can we judge the accuracy of an RK4 result? A simple and reliable trick is step doubling: compute the same value twice, with steps $h$ and $h /2$ .

Call the results $y_{h}$ and $y_{h /2}$ , and the true value $y$ . Since the method is fourth order,

y - y_{h} \approx C h^{4}, y - y_{h /2} \approx C (\frac{h}{2})^{4} = \frac{C h ^{4}}{16} .

Subtracting the second relation from the first eliminates the unknown $y$ :

y_{h /2} - y_{h} so y - y_{h /2} \approx C h^{4} - \frac{C h ^{4}}{16} = \frac{15}{16} C h^{4}, \approx \frac{C h ^{4}}{16} \approx \frac{y _{h /2} - y _{h}}{15} .

Example 8.6.4 (Step doubling)

RK4 for $y^{'} = y - x^{2}$ , $y (0) = 1$ gives $y_{h} = 2.28169712$ at $x = 1$ with $h = 0.2$ , and $y_{h /2} = 2.28171685$ with $h = 0.1$ . Estimate the error of the better value, pretending we do not know the exact solution.

Solution

Estimate. $y - y_{h /2} \approx \frac{2.28171685 - 2.28169712}{15} = \frac{1.973 \times 1 0 ^{- 5}}{15} \approx 1.3 \times 1 0^{- 6} .$
Check. The true error is $y (1) - y_{h /2} = 5 - e - 2.28171685 = 1.3 \times 1 0^{- 6}$ . The estimate is excellent. Adaptive solvers in engineering software use ideas of exactly this kind to choose $h$ automatically.

Example 8.6.5 (Computing

π

with RK4)

Since $\int_{0}^{1} \frac{d x}{1 + x ^{2}} = arctan 1 = \frac{π}{4}$ , the solution of $y^{'} = \frac{1}{1 + x ^{2}}$ , $y (0) = 0$ , satisfies $y (1) = π /4$ . Use RK4 with $h = 0.5$ to estimate $π$ .

Solution

First step ( $x_{0} = 0$ ). Since $f$ does not involve $y$ , $k_{2} = k_{3} = f (0.25)$ : $k_{1} y_{1} = 1, k_{2} = k_{3} = \frac{1}{1.0625} = 0.941176, k_{4} = \frac{1}{1.25} = 0.8, = 0 + \frac{0.5}{6} (1 + 4 (0.941176) + 0.8) = 0.463725.$
Second step ( $x_{1} = 0.5$ ). $k_{1} y_{2} = 0.8, k_{2} = k_{3} = \frac{1}{1.5625} = 0.64, k_{4} = \frac{1}{2} = 0.5, = 0.463725 + \frac{0.5}{6} (0.8 + 4 (0.64) + 0.5) = 0.785392.$
Result. $π \approx 4 y_{2} = 3.141569$ , compared with $π = 3.141593$ . Four slope evaluations at five distinct points give five correct significant digits. (This is Simpson's rule with two panels.)

Summary

RK4 samples the slope at the start, twice at the midpoint, and at the end of each step, and averages with weights $1, 2, 2, 1$ .
It is fourth order: global error $O (h^{4})$ , local error $O (h^{5})$ . Halving $h$ divides the error by about $16$ .
For $y^{'} = f (x)$ it reduces to Simpson's rule.
At equal cost (same number of $f$ -evaluations) it is far more accurate than Euler or Heun, which is why it is the default choice for smooth problems.
Step doubling gives a practical error estimate: $y - y_{h /2} \approx (y_{h /2} - y_{h}) /15$ .

Exercises

Take one RK4 step with $h = 0.1$ for $y^{'} = 3 x - y$ , $y (0) = 2$ , listing $k_{1}, \dots, k_{4}$ . Compare with the exact value $y (0.1) = 1.824187$ and with the Heun result of Exercise 1 in Section 8.5.
Answer
$k_{1} = - 2$ , $k_{2} = - 1.75$ , $k_{3} = - 1.7625$ , $k_{4} = - 1.52375$ ; $y_{1} = 1.824188$ (error $\approx - 4 \times 1 0^{- 7}$ , versus $- 8.1 \times 1 0^{- 4}$ for Heun).
Use RK4 with $h = 0.2$ to approximate $y (0.4)$ for $y^{'} = 1 + y^{2}$ , $y (0) = 0$ . Compare with $tan 0.4$ .
Answer
$y_{1} = 0.202707$ , $y_{2} = 0.422789$ ; $tan 0.4 = 0.422793$ .
Use RK4 with $h = 0.25$ to approximate $y (0.5)$ for $y^{'} = x^{2} - y^{2}$ , $y (0) = 1$ . (Compare with the Heun result of Exercise 3, Section 8.5, and with the accurate value $0.701768$ .)
Answer
$y_{1} = 0.804726$ , $y_{2} = 0.701792$ (Heun gave $0.714083$ ).
For the draining tank $H^{'} = - 0.6 H$ , $H (0) = 4$ of Example 8.5.3, take one RK4 step with $h = 0.5$ .
Answer
$k_{1} = - 1.2$ , $k_{2} = - 1.154123$ , $k_{3} = - 1.155910$ , $k_{4} = - 1.109926$ ; $H_{1} = 3.422501$ (exact $3.4225$ ).
Show directly from Definition 8.6.1 that RK4 gives the exact answer, for any $h$ , when $f (x, y) = a + b x + c x^{2} + d x^{3}$ (a cubic in $x$ alone).
Answer
Here $k_{2} = k_{3} = f (x_{n} + h /2)$ , so the step is Simpson's rule on $[x_{n}, x_{n} + h]$ , which integrates cubics exactly. Hence $y_{n + 1} - y_{n} = \int_{x_{n}}^{x_{n + 1}} f d x$ exactly.
For $y^{'} = λ y$ the RK4 multiplier is $R (z) = 1 + z + \frac{z ^{2}}{2} + \frac{z ^{3}}{6} + \frac{z ^{4}}{24}$ with $z = hλ$ . Evaluate $R (- 0.5)$ and compare with $e^{- 0.5}$ . Do the same for Euler's multiplier $1 + z$ .
Answer
$R (- 0.5) = 0.606771$ , $e^{- 0.5} = 0.606531$ (difference $2.4 \times 1 0^{- 4}$ ); Euler gives $0.5$ (difference $0.107$ ).
Use RK4 with $h = 0.5$ on $y^{'} = 1/ x$ , $y (1) = 0$ , to estimate $ln 2$ .
Answer
$y_{1} = 0.405556$ , $y_{2} = 0.693254$ ; $ln 2 = 0.693147$ .
In Example 8.6.1 the RK4 error at $x = 1$ was $2.1 \times 1 0^{- 5}$ with $h = 0.2$ . Using the fourth-order law, estimate the step size needed for an error below $1 0^{- 8}$ . How many evaluations of $f$ is that?
Answer
$C \approx 2.1 \times 1 0^{- 5} /0. 2^{4} \approx 0.013$ , so $h \leq (1 0^{- 8} /0.013)^{1/4} \approx 0.03$ ; e.g. $h = 0.025$ (40 steps, 160 evaluations), which gives $5.1 \times 1 0^{- 9}$ .
For $y^{'} = 1 + y^{2}$ , $y (0) = 0$ , RK4 gives $y_{h} = 0.42278899$ at $x = 0.4$ with $h = 0.2$ and $y_{h /2} = 0.42279299$ with $h = 0.1$ . Estimate the error of the second value by step doubling and compare with the true error.
Answer
$(0.42279299 - 0.42278899) /15 \approx 2.7 \times 1 0^{- 7}$ ; the true error is $tan 0.4 - 0.42279299 = 2.3 \times 1 0^{- 7}$ .
A population model reads $P^{'} = 0.8 P (1 - P /10)$ , $P (0) = 1$ (in thousands, $t$ in years). Use RK4 with $h = 0.5$ to estimate $P (1)$ , and compare with the exact logistic solution $P = 10/ (1 + 9 e^{- 0.8 t})$ .
Answer
$P_{1} = 1.421850$ , $P_{2} = 1.982469$ ; exact $P (1) = 1.982569$ .

The classical RK4 formulas

What the four slopes mean

Why the weights 1, 2, 2, 1?

Worked examples

How accurate is RK4?

A check on the test equation y′=λy

A fair comparison: Euler, Heun and RK4

Estimating the error without the exact solution

Summary

Exercises

A check on the test equation $y^{'} = λ y$