Topic 14.2: Heun's Method

Contents Previous Chapter Start of Chapter Previous Topic Introduction Notes Theory HOWTO Examples Engineering Error Questions Matlab Maple Next Topic Next Chapter

The trapezoidal rule ½(f(a) + f(b))(b - a) for integration is much better than using only the left end-point with a formula like f(a)(b - a) because it samples the value of the function at two points. If the function is reasonably smooth, then the interpolation of these two points is appropriate. Similarly, Euler's method only samples the slope at the left end-point, that is, the initial point (t0, y0). Heun's method is an attempt to use information about the slope at both end points to find the average slope.

Like Euler, the eu in Heun sounds like oi, that is, Heun rhymes with coin.

Background

Useful background for this topic includes:

References

Interactive Maplet

A Demonstration of Heun's Method

Theory


Euler's method only samples the slope at the left end-point, that is, the initial point (t0, y0). Let's call this first slope K0 = f(t0, y0). Then if h = t1 - t0, then y1 = y0 + h K0 is Euler's approximation of y(t1).

Let's now sample the slope at this new point (t1, y1), that is, K1 = f(t1, y1). This is only an approximation of the actual slope at y(t1) because y1 is only an approximation of y(t1). Thus, K1 is not as accurate as K0, but let us assume it is still a reasonable approximation.

Therefore, if the change in slope is approximately linear on the interval [t0, t1], then why not use ½(K0 + K1) as a reasonable average value of the slope on the interval?

Thus, use Euler's method to approximate y1 and then use this approximation to approximate K1. Average these two slopes and use this to approximate y1:

Graphical Representation

To view this graphically, let's consider the IVP

y(1)(t) = ½y(t) - t + 1
y(0) = 0.5

Figure 1 shows Euler's approximation of y(1), that is, following the slope K0 = f(t0, y0) = 1.25, we get that y(t1) ≈ 0.5 + 1 ⋅ 1.25 = 1.75.

Figure 1. Euler's approximation to the value of y(1).

If we sample the slope at (1, 1.75), we find that it is K1 = f(t1, y1) = f(1, 1.75) = 0.875. This slope is shown in Figure 2.

Figure 2. The slope at the point (1, 1.75).

The average of these two slopes is ½(K0 + K1) = 1.0625. Therefore, we approximate y(t1) by following this average slope out a distance 1: y1 = y0 + ½(K0 + K1)h = 0.5 + 1.0625 ⋅ 1 = 1.5625. This is shown in Figure 3. Note the significant improvement in the approximation as compared to that shown in Figure 1.

Figure 3. Heun's approximation to y(1).


HOWTO


Problem

Given the IVP

y(1)(t) = f( t, y(t) )
y(t0) = y0

approximate y(t1).

Assumptions

The function f(t, y) should be continuous in both variables.

Tools

We will use Taylor series.

Initial Conditions

Set h = t1t0. Let y1 be the approximation of y(t1).

Process

Define

K0 = f( t0, y0 ) and
K1 = f( t0 + h, y0 + h K0 )

Then let

Note that for K1, t0 + h = t1, however we use the above notation to provide similarity with the next topic.


Examples


Example 1

Given the IVP y(1)(t) = 1 - t y(t) with y(0) = 1, approximate y(1) with one step.

First, let t0 = 0, y0 = 1, and h = 1. Thus, we calculate

K0 = f(0, 1) = 1
K1 = f(1, 2) = -1
½(K0 + K1) = 0

Therefore, approximation is y1 = y0 + h = 1 + 1⋅0 = 1. The actual value is 1.331309118 and therefore the absolute error is approximately 0.331, or approximately half that of what we found using Euler's method.

Example 2

Given the same IVP shown in Example 1, approximate y(0.5).

K0 = f(0, 1) = 1
K1 = f(0.5, 1.5) = 0.25
½(K0 + K1) = 0.625

Therefore, approximation is y1 = y0 + h = 1 + 0.5⋅0.625 = 1.3125. The actual value is 1.342841185 and therefore the absolute error is approximately 0.0303, or approximately 1/8 the error when we used h = 1.

Example 3

Repeat Examples 1 and 2 but with with the initial value y(0.5) = 2.5 and approximating y(1.5) and y(1.0).

First, let t0 = 0.5, y0 = 2.5, and h = 1. Therefore

K0 = f(0.5, 2.5) = -0.25
K1 = f(1.5, 2.5 - 1 ⋅ 0.25) = -2.375
½(K0 + K1) = -1.3125

Therefore, our approximation is y1 = y0 + h(-1.3125) = 2.5 + 1(-1.3125) = 1.1875. The actual value is 1.502483616 and therefore the absolute error is approximately 0.315.

Next, let t0 = 0.5, y0 = 2.5, and h = 0.5. Therefore

K0 = f(0.5, 2.5) = -0.25
K1 = f(1.5, 2.5 - 1 ⋅ 0.25) = -1.375
½(K0 + K1) = -0.8125

Therefore, our approximation is y1 = y0 + h(-0.8125) = 2.5 + 0.5 (-0.8125) = 2.09375. The actual value is 2.126611964 and therefore the absolute error is approximately 0.0329. Again, the error is approximately 1/8th that of the previous calculation where h = 1.

If we tabulate the errors for various values of h, as is shown in Table 1, we note that as h gets smaller, the error drops by almost exactly 1/8 each time we divide h by two.

Table 1. Errors when approximating y(t0 + h) for decreasing values of h.

hApproximation
of y(0.5 + h)
Error
11.18750.315
0.52.093750.0329
0.252.3652343750.00338
0.1252.4504394531250.000366
0.06262.4797668460.0000421
0.031252.4910316470.00000502
0.016252.4958043100.000000612
0.00781252.4979744550.0000000755
0.003906252.4990053250.00000000937

You will note that with each step, the error goes down by approximately a factor of 1/8.


Engineering


To be completed.


Error


To determine the error for Heun's method, we must look at the Taylor series:

Consider the IVP:

y(1)(t) = f(t, y(t))
y(t0) = y0

We would like to approximate y(t0 + h). The Taylor series gives us that:

Note that we may approximate the 2nd derivative by taking the forward-divided difference formula [i.e., f(1)(x) = (f(x + h) − f(x))/h − ½f(2)(ξ)h] for the derivative of the 1st derivative:

First, let us denote K0 = y(1)(t0) = f(t0, y0) and note that y(t0 + h) ≈ y0 + hK0.

We also observe that y(1)(t0 + h) = f(t0 + h, y(t0 + h)), however, by substituting in our approximation, we have y(1)(t0 + h) ≈ f(t0 + h, y0 + hK0)).

By letting K1 = f(t0 + h, y0 + hK0)), we may rewrite the approximation as

Substituting these into our Taylor series, we get that:

and by canceling terms, we get:

Therefore, the error is O(h3).


Questions


Question 1

Given the IVP

y(1)(t) = 1 - 0.25 y(t) + 0.2 t
y(0) = 1

approximate y(0.5), y(1), and y(1.5) using Heun's method.

Answer: 1.3765625, 1.75625, 2.1390625.

Question 2

Given the same ODE as in Question 1, but with the initial condition y(1) = 2, approximate y(1.5) and y(2.0).

Answer: 2.353125 and 2.7125

Question 3

Given an IVP with an initial condition y(0) = y0, if the second derivative is bounded by -8 < y(3)(t) < 8, on how large an interval can we estimate y(t) if we want to ensure that the error is less than 0.0001? Compare this with the range for Question 3 of Euler's method.

Answer: (-0.0531, 0.0531).

Question 4

Consider the following modification to Huen's method:

Instead of defining K1 = f(t1, y0 + hK0), use this as an initial value and iterate:

K1 = f(t1, y0 + ½h(K0 + K1)).

For Question 1, use this modification, iterating the previous step 10 times.


Matlab


To be completed.


Maple


Given the IVP y(1)(t = t y(t) - t2 + 1 with y(0) = 3, to approximate y(1), we can do:

t[0] := 0;
t[1] := 1;
y[0] := 3;
f := (t, y) -> t*y - t^2 + 1;
h := t[1] - t[0];
K[0] := f(t[0],     y[0]         );
K[1] := f(t[0] + h, y[0] + h*K[0]);
y[1] := y[0] + h*(K[0] + K[1])/2;


Copyright ©2005 by Douglas Wilhelm Harder. All rights reserved.