Chapter 8 · Section 8.5

8.5 Numerical Methods II: The Improved Euler Method

Euler's method from Section 8.4 is easy to understand and easy to program, but it is slow to converge: halving the step size only halves the error. To get three more correct digits you need a thousand times more steps. In this section we make one small, clever change to Euler's method and obtain a method whose error drops by a factor of four each time the step is halved. The price is one extra evaluation of the slope function per step, which is almost always worth paying.

Throughout, we work with the initial-value problem

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

on equally spaced points $x_{n} = x_{0} + nh$ . Section 8.4 called the exact solution $φ$ ; from here on, to keep formulas short, we write it simply as $y (x)$ . So $y (x_{n})$ is the exact value and $y_{n}$ is our approximation of it.

Where Euler's method loses accuracy

Start from the exact relation between two neighbouring values of the solution. Integrating $y^{'} = f (x, y (x))$ from $x_{n}$ to $x_{n + 1}$ gives

y (x_{n + 1}) = y (x_{n}) + \int_{x_{n}}^{x_{n + 1}} f (x, y (x)) d x .

This formula is exact. Every one-step method is simply a way to estimate the integral on the right.

Euler's method replaces the integrand by its value at the left end of the interval. Geometrically, it approximates the area under the curve $f (x, y (x))$ by a rectangle of height $f (x_{n}, y_{n})$ . If the slope of the solution changes across the step (and it almost always does), the rectangle misses a whole sliver of area. That sliver is the source of the $O (h^{2})$ local error you met in Section 8.4.

From calculus we know a better estimate of an integral: the trapezoid rule, which averages the integrand at both ends,

\int_{x_{n}}^{x_{n + 1}} g (x) d x \approx \frac{h}{2} [g (x_{n}) + g (x_{n + 1})] .

Applied to our integral, it suggests

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

There is a catch. The unknown $y_{n + 1}$ appears on both sides, buried inside $f$ . For a nonlinear $f$ we cannot simply solve for it. A formula of this kind is called implicit.

Predict, then correct

The way out is wonderfully simple. We need $y_{n + 1}$ only to evaluate the slope at the right end of the interval, and for that purpose a rough value is good enough. So we first predict $y_{n + 1}$ with an ordinary Euler step, and then use this prediction inside the trapezoid formula to correct it.

Definition 8.5.1 (Improved Euler method, or Heun's method)

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

k_{1} y_{n + 1}^{*} k_{2} y_{n + 1} = f (x_{n}, y_{n}) = y_{n} + h k_{1} = f (x_{n + 1}, y_{n + 1}^{*}) = y_{n} + \frac{h}{2} (k_{1} + k_{2}) (slope at the left end) (predictor: one Euler step) (slope at the predicted right end) (corrector: trapezoid rule)

Methods built from a rough explicit guess followed by an improving formula are called predictor–corrector methods. Heun's method is the simplest of them.

Figure 8.5.1 shows one step. The Euler line (slope $k_{1}$ ) undershoots the curve, because the solution bends upward. At the predicted point we sample the slope again and get $k_{2}$ , which is too steep for the beginning of the step but about right for its end. Their average is a much better description of the whole step.

Figure 8.5.1. One step of Heun's method for

y^{'} = y

y (0) = 1

, with the deliberately large step

h = 1

. Euler predicts

2

; averaging the slopes

1

and

2

gives

2.5

, much closer to

e \approx 2.718

Method in a nutshell: one step of Heun's method

Evaluate the slope at the current point: $k_{1} = f (x_{n}, y_{n})$ .
Take a trial Euler step: $y_{n + 1}^{*} = y_{n} + h k_{1}$ .
Evaluate the slope at the trial point: $k_{2} = f (x_{n} + h, y_{n + 1}^{*})$ .
Step with the average slope: $y_{n + 1} = y_{n} + \frac{h}{2} (k_{1} + k_{2})$ . Discard $y_{n + 1}^{*}$ .
Repeat from the new point $(x_{n + 1}, y_{n + 1})$ .

Common mistake

The predicted value $y_{n + 1}^{*}$ is only scaffolding. The next step starts from the corrected value $y_{n + 1}$ , never from $y_{n + 1}^{*}$ . Also note that $k_{2}$ uses the new abscissa $x_{n + 1} = x_{n} + h$ , not $x_{n}$ .

Worked examples

We will use the following problem as a running test case in this section and the next two:

y^{'} = y - x^{2}, y (0) = 1.

It is linear, so Chapter 2 gives its exact solution. With the integrating factor $e^{- x}$ one finds $y = x^{2} + 2 x + 2 + C e^{x}$ , and $y (0) = 1$ forces $C = - 1$ :

y (x) = x^{2} + 2 x + 2 - e^{x} .

Knowing the exact answer lets us measure the true error of every method. In all tables we carry full precision in the calculation (errors included) and round only the displayed numbers to six decimals, so a last digit may occasionally differ by one from what you get by subtracting two rounded entries.

Example 8.5.1

Use Heun's method 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$ . Compare with Euler's method and with the exact solution.

Solution

First step ( $n = 0$ ). Here $x_{0} = 0$ , $y_{0} = 1$ , $x_{1} = 0.2$ . $k_{1} y_{1}^{*} k_{2} y_{1} = f (0, 1) = 1 - 0 = 1, = 1 + 0.2 (1) = 1.2, = f (0.2, 1.2) = 1.2 - 0.04 = 1.16, = 1 + 0.1 (1 + 1.16) = 1.216.$ The exact value is $y (0.2) = 0.04 + 0.4 + 2 - e^{0.2} = 1.218597$ , so the error is $0.002597$ .
Second step ( $n = 1$ ). Now $x_{1} = 0.2$ , $y_{1} = 1.216$ , $x_{2} = 0.4$ . $k_{1} y_{2}^{*} k_{2} y_{2} = f (0.2, 1.216) = 1.216 - 0.04 = 1.176, = 1.216 + 0.2 (1.176) = 1.4512, = f (0.4, 1.4512) = 1.4512 - 0.16 = 1.2912, = 1.216 + 0.1 (1.176 + 1.2912) = 1.46272.$

Remaining steps. Continuing in the same way gives the table below.

$x_{n}$	$y_{n}$	$k_{1}$	$y_{n + 1}^{*}$	$k_{2}$	$y_{n + 1}$
0.0	1.000000	1.000000	1.200000	1.160000	1.216000
0.2	1.216000	1.176000	1.451200	1.291200	1.462720
0.4	1.462720	1.302720	1.723264	1.363264	1.729318
0.6	1.729318	1.369318	2.003182	1.363182	2.002568
0.8	2.002568	1.362568	2.275082	1.275082	2.266334

Comparison. Euler's method with the same

h

gives

y_{n + 1} = y_{n} + 0.2 (y_{n} - x_{n}^{2})

$x_{n}$	exact $y (x_{n})$	Euler	Euler error	Heun	Heun error
0.2	1.218597	1.200000	0.018597	1.216000	0.002597
0.4	1.468175	1.432000	0.036175	1.462720	0.005455
0.6	1.737881	1.686400	0.051481	1.729318	0.008563
0.8	2.014459	1.951680	0.062779	2.002568	0.011891
1.0	2.281718	2.214016	0.067702	2.266334	0.015385

x = 1

, Heun's error is about

4.4

times smaller than Euler's, although each Heun step costs only two slope evaluations instead of one.

Figure 8.5.2. Example 8.5.1. Left: the Heun points (squares) sit almost on the exact curve, while the Euler points (circles) drift below it. Right: both errors grow with

x

, but Heun's grows far more slowly.

Try it yourself. The solver below reproduces the table of Example 8.5.1 with every intermediate quantity. Change $h$ to $0.1$ and watch the error at $x = 1$ fall by about a factor of four; then type your own $f$ , for instance 1 + y^2 with exact solution tan(x).

How accurate is Heun's method?

A single table does not tell us how the error behaves as $h \to 0$ . For that we compare the method with the Taylor expansion of the exact solution, exactly as Section 8.4 did for Euler's method.

The exact solution over one step

Suppose we start a step from the exact value, $y_{n} = y (x_{n})$ . Taylor's theorem gives

y (x_{n + 1}) = y (x_{n}) + h y^{'} (x_{n}) + \frac{h ^{2}}{2} y^{''} (x_{n}) + O (h^{3}) .

We know $y^{'} = f (x, y)$ . To find $y^{''}$ , differentiate $f (x, y (x))$ with the chain rule (write $f_{x}$ , $f_{y}$ for the partial derivatives):

y^{''} = \frac{d}{d x} f (x, y (x)) = f_{x} + f_{y} y^{'} = f_{x} + f f_{y} .

So, with every function evaluated at $(x_{n}, y_{n})$ ,

y (x_{n + 1}) = y_{n} + h f + \frac{h ^{2}}{2} (f_{x} + f f_{y}) + O (h^{3}) . (A)

Heun's step over the same interval

Now expand $k_{2} = f (x_{n} + h, y_{n} + h k_{1})$ , where $k_{1} = f$ . Taylor's theorem in two variables, with increments $h$ in $x$ and $h f$ in $y$ , gives

k_{2} = f + h f_{x} + (h f) f_{y} + O (h^{2}) = f + h (f_{x} + f f_{y}) + O (h^{2}) .

Insert this into the corrector:

y_{n + 1} = y_{n} + \frac{h}{2} (k_{1} + k_{2}) = y_{n} + \frac{h}{2} (2 f + h (f_{x} + f f_{y}) + O (h^{2})) = y_{n} + h f + \frac{h ^{2}}{2} (f_{x} + f f_{y}) + O (h^{3}) . (B)

Compare (A) and (B): they agree in the terms of order $1$ , $h$ and $h^{2}$ . They first differ at order $h^{3}$ . Euler's method, $y_{n} + h f$ , already disagrees with (A) at order $h^{2}$ .

Theorem 8.5.1 (Order of Heun's method)

Suppose $f$ and its partial derivatives up to second order are continuous near the solution curve, and that the solution has a continuous third derivative on $[x_{0}, b]$ . Then Heun's method has local truncation error $O (h^{3})$ , and its global error at any fixed $x$ in $[x_{0}, b]$ satisfies

∣ y (x_{n}) - y_{n} ∣ \leq C h^{2}

for some constant $C$ that does not depend on $h$ . In short, Heun's method is a second-order method.

The step from "local $O (h^{3})$ " to "global $O (h^{2})$ " uses the same reasoning as for Euler's method in Section 8.4. To reach a fixed point $x$ we take $(x - x_{0}) / h$ steps; each adds a local error of size about $h^{3}$ ; and, under the stated smoothness, earlier errors are not amplified by more than a bounded factor. The total is roughly $(x - x_{0}) / h \times h^{3}$ , which is of order $h^{2}$ .

What "second order" means in practice

If the error is about $C h^{2}$ , then halving $h$ divides it by $2^{2} = 4$ , and dividing $h$ by $10$ divides it by $100$ . For Euler's method (first order) the same changes divide the error only by $2$ and by $10$ .

Example 8.5.2 (Seeing the order in the numbers)

Solve $y^{'} = y - x^{2}$ , $y (0) = 1$ by Euler's method and by Heun's method with $h = 0.2, 0.1, 0.05, 0.025$ . Tabulate the error at $x = 1$ and check the predicted orders.

Solution

Errors. The exact value is $y (1) = 5 - e = 2.281718$ . Running both methods to $x = 1$ (that is, $1/ h$ steps) gives:

$h$ Euler error ratio Heun error ratio

0.2 0.067702 – 0.015385 –

0.1 0.034835 1.94 0.003961 3.88

0.05 0.017681 1.97 0.001003 3.95

0.025 0.008909 1.98 0.000252 3.98

The "ratio" column is the previous error divided by the current one.
Reading the ratios. For Euler the ratios approach $2 = 2^{1}$ ; for Heun they approach $4 = 2^{2}$ . This is exactly the behaviour predicted for first- and second-order methods.
Estimating the order. If error $\approx C h^{p}$ , then the ratio for halving is $2^{p}$ , so $p \approx lo g_{2} (ratio)$ . For the last Heun ratio, $p \approx lo g_{2} 3.98 \approx 1.99$ .

$h$	Euler error	ratio	Heun error	ratio
0.2	0.067702	–	0.015385	–
0.1	0.034835	1.94	0.003961	3.88
0.05	0.017681	1.97	0.001003	3.95
0.025	0.008909	1.98	0.000252	3.98

A log–log plot makes this visible at a glance. If the error is $C h^{p}$ , then $lo g (error) = lo g C + p lo g h$ , a straight line of slope $p$ .

Figure 8.5.3. Error at

x = 1

against

h

on logarithmic axes for

y^{'} = y - x^{2}

y (0) = 1

. The Euler data follow a line of slope

1

, the Heun data a line of slope

2

(dotted reference lines).

Example 8.5.3 (A draining tank)

Water drains from a cylindrical tank through a hole in its bottom. By Torricelli's law the water depth $H$ (in metres) obeys $H^{'} = - 0.6 H$ , with $t$ in minutes. The initial depth is $H (0) = 4$ . Use Heun's method with $h = 0.5$ to estimate $H (0.5)$ and $H (1)$ , and compare with the exact solution.

Solution

Exact solution. The equation is separable (Section 2.2): $H^{- 1/2} d H 2 H = - 0.6 d t = - 0.6 t + c .$ With $H (0) = 4$ we get $c = 4$ , so $H = 2 - 0.3 t$ and $H (t) = (2 - 0.3 t)^{2}$ for $0 \leq t \leq 20/3$ (the tank is empty at $t = 20/3$ min).
First step. $f (t, H) = - 0.6 H$ , $t_{0} = 0$ , $H_{0} = 4$ . $k_{1} H_{1}^{*} k_{2} H_{1} = - 0.6 4 = - 1.2, = 4 + 0.5 (- 1.2) = 3.4, = - 0.6 3.4 = - 1.106345, = 4 + 0.25 (- 1.2 - 1.106345) = 3.423414.$ Exact: $H (0.5) = 1.8 5^{2} = 3.4225$ . Error: $- 0.000914$ .
Second step. $k_{1} H_{2}^{*} k_{2} H_{2} = - 0.6 3.423414 = - 1.110148, = 3.423414 + 0.5 (- 1.110148) = 2.868340, = - 0.6 2.868340 = - 1.016170, = 3.423414 + 0.25 (- 1.110148 - 1.016170) = 2.891834.$ Exact: $H (1) = 1. 7^{2} = 2.89$ . Error: $- 0.001834$ , less than $0.07%$ of the depth, even with a step of half a minute. (Euler's method with the same step, Example 8.4.3, was off by $0.043$ here.)

A special case: pure integration

When $f$ depends on $x$ only, the differential equation $y^{'} = f (x)$ is just an integration problem: $y (b) = y (a) + \int_{a}^{b} f (x) d x$ . In that case $k_{2} = f (x_{n + 1})$ does not depend on the predictor at all, and Heun's step becomes

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

which is exactly the trapezoid rule of calculus. So Heun's method is the natural generalisation of the trapezoid rule to equations whose slope depends on $y$ .

Example 8.5.4 (An integral with no elementary antiderivative)

The function $e^{- x^{2}}$ has no elementary antiderivative, yet its integral appears constantly in heat conduction and probability. Use Heun's method with $h = 0.25$ on $y^{'} = e^{- x^{2}}$ , $y (0) = 0$ to estimate $y (1) = \int_{0}^{1} e^{- x^{2}} d x$ .

Solution

Slopes needed. Since $k_{2}$ at one step equals $k_{1}$ at the next, we need $f$ only at the five nodes: $f (0) f (0.75) = 1, = 0.569783, f (0.25) f (1) = 0.939413, = 0.367879. f (0.5) = 0.778801,$
Four Heun steps. $y_{1} y_{2} y_{3} y_{4} = 0 + 0.125 (1 + 0.939413) = 0.242427, = 0.242427 + 0.125 (0.939413 + 0.778801) = 0.457203, = 0.457203 + 0.125 (0.778801 + 0.569783) = 0.625776, = 0.625776 + 0.125 (0.569783 + 0.367879) = 0.742984.$
Check. The true value is $\frac{π}{2} erf (1) = 0.746824$ . The error, $0.003840$ , is the familiar error of the trapezoid rule with four panels.

Why not iterate the corrector?

One could feed the corrected value back into $k_{2}$ and correct again, repeating until the numbers stop changing. The limit is the solution of the implicit trapezoid formula from the start of this section. This does not raise the order: the trapezoid rule itself is second order, so the iterations only shrink the constant $C$ a little, at the cost of extra slope evaluations. To gain real accuracy it is better to sample the slope at new places inside the step. That is the idea behind the Runge–Kutta methods of Section 8.6.

Summary

Heun's method estimates $\int_{x_{n}}^{x_{n + 1}} f d x$ with the trapezoid rule, using an Euler predictor to supply the unknown right-end value.
It costs two evaluations of $f$ per step and is second order: global error $O (h^{2})$ , local error $O (h^{3})$ .
Halving $h$ divides the error by about $4$ , a fact you can use to check a computation.
For $y^{'} = f (x)$ it reduces to the trapezoid rule.

Exercises

Use Heun's method with $h = 0.1$ to approximate $y (0.1)$ and $y (0.2)$ for $y^{'} = 3 x - y$ , $y (0) = 2$ . Compare with the exact solution $y = 3 x - 3 + 5 e^{- x}$ .
Answer
$y_{1} = 1.825000$ , $y_{2} = 1.695125$ ; exact $1.824187$ and $1.693654$ (errors $- 0.000813$ , $- 0.001471$ ).
Use Heun's method with $h = 0.1$ to approximate the solution of $y^{'} = 1 + y^{2}$ , $y (0) = 0$ , at $x = 0.1, 0.2, 0.3$ . The exact solution is $y = tan x$ .
Answer
$0.100500$ , $0.203035$ , $0.309814$ ; exact $0.100335$ , $0.202710$ , $0.309336$ .
Use Heun's method with $h = 0.25$ to approximate $y (0.5)$ for $y^{'} = x^{2} - y^{2}$ , $y (0) = 1$ (this equation has no elementary solution).
Answer
$y_{1} = 0.812500$ , $y_{2} = 0.714083$ . (A very accurate computation gives $y (0.5) = 0.701768$ .)
Apply one Heun step to the test equation $y^{'} = λ y$ ( $λ$ a constant). Show that $y_{n + 1} = (1 + hλ + \frac{1}{2} h^{2} λ^{2}) y_{n}$ , and compare the multiplier with the exact factor $e^{hλ}$ .
Answer
$k_{1} = λ y_{n}$ , $k_{2} = λ (1 + hλ) y_{n}$ , so $y_{n + 1} = y_{n} + \frac{h}{2} λ (2 + hλ) y_{n}$ . The multiplier matches $e^{hλ} = 1 + hλ + \frac{1}{2} h^{2} λ^{2} + \frac{1}{6} h^{3} λ^{3} + \dots$ up to the $h^{2}$ term; the first difference is $\frac{1}{6} h^{3} λ^{3}$ , a local error of order $h^{3}$ .
The midpoint method uses $y_{n + 1} = y_{n} + h f (x_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} f (x_{n}, y_{n}))$ . (a) Take one step with $h = 0.2$ for $y^{'} = y - x^{2}$ , $y (0) = 1$ . (b) Using the Taylor expansion from this section, show that the midpoint method is also second order.
Answer
(a) $y_{1} = 1 + 0.2 f (0.1, 1.1) = 1 + 0.2 (1.09) = 1.218$ (exact $1.218597$ ). (b) $f (x_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} f) = f + \frac{h}{2} (f_{x} + f f_{y}) + O (h^{2})$ , so $y_{n + 1} = y_{n} + h f + \frac{h ^{2}}{2} (f_{x} + f f_{y}) + O (h^{3})$ , which matches (A).
For the first step of Example 8.5.1, iterate the corrector: start from $y^{(0)} = 1.2$ and repeat $y^{(k + 1)} = 1 + 0.1 [f (0, 1) + f (0.2, y^{(k)})]$ . Find the first three iterates and the limit. Is the limit closer to $y (0.2) = 1.218597$ than Heun's value $1.216$ ?
Answer
$1.216$ , $1.2176$ , $1.21776$ ; the limit solves $y = 1.096 + 0.1 y$ , so $y = 274/225 = 1.217778$ . Yes, a little closer (error $0.000819$ instead of $0.002597$ ), but the order is unchanged.
A second-order method gives an error of $2.4 \times 1 0^{- 3}$ at some point when $h = 0.1$ . Estimate the error with $h = 0.05$ and with $h = 0.01$ . If a computation instead gives errors $2.4 \times 1 0^{- 3}$ and $6.1 \times 1 0^{- 4}$ for $h = 0.1$ and $h = 0.05$ , what order does this suggest?
Answer
About $6 \times 1 0^{- 4}$ and $2.4 \times 1 0^{- 5}$ . The ratio $2.4/0.61 \approx 3.93$ gives $p \approx lo g_{2} 3.93 \approx 1.98$ , i.e. order 2.
Use Heun's method with $h = π /4$ on $y^{'} = cos x$ , $y (0) = 0$ , to estimate $y (π /2)$ . Which classical rule have you used, and what is the error?
Answer
The trapezoid rule: $y_{2} = \frac{π}{8} (1 + 2 cos \frac{π}{4} + 0) = \frac{π}{8} (1 + 2) = 0.948059$ . Exact $1$ ; error $0.051941$ .
For $y^{'} = y - x^{2}$ , $y (0) = 1$ , Example 8.5.2 shows Heun's error at $x = 1$ is about $0.396 h^{2}$ . Roughly what step size makes the error smaller than $1 0^{- 4}$ ? How many slope evaluations is that, and how many would Euler's method need (its error is about $0.35 h$ )?
Answer
$h \leq 1 0^{- 4} /0.396 \approx 0.016$ ; e.g. $h = 0.0125$ (80 steps, 160 evaluations; actual error $6.3 \times 1 0^{- 5}$ ). Euler needs $h ≲ 2.9 \times 1 0^{- 4}$ , about $3500$ steps and evaluations.
Use Heun's method with $h = 0.5$ for the tank of Example 8.5.3 to estimate $H (1.5)$ , continuing from $H_{2} = 2.891834$ .
Answer
$k_{1} = - 1.020324$ , $H_{3}^{*} = 2.381672$ , $k_{2} = - 0.925960$ , $H_{3} = 2.405263$ ; exact $1.5 5^{2} = 2.4025$ (error $- 0.002763$ ).