## Mathematics & Science
Learning Center |
|||

## Numerical Methods for Solving Differential Equations## Heun's Method## Theoretical IntroductionIn the last lab you learned to use Euler's Method to generate a numerical solution to an initial value problem of the form:
Now it's time for a confession: In the real-world of using computers to derive numerical solutions to differential equations, no-one actually uses Euler's Method. Its shortcomings, discussed in detail in the last lab, nameley its inaccuracy and its slowness, are just too great. Though it is of some historical interest, the primary reason we bother to discuss it at all is that while it is a fairly simple algorithm to understand, it still enables us to start thinking more clearly about algorithmic approaches to the problem in general. We will now develop a better method than Euler's for numerically solving this same kind of initial value problem, but we'll use Euler's method as a foundation. For this reason Heun's Method is sometimes referred to as the Improved Euler Method. So where should we begin? Let's start by asking what it is about Euler's Method that makes it so poor at its job. Recall that the basic idea is to use the tangent line to the actual solution curve as an estimate of the curve itself, with the thought that provided we don't project too far along the tangent line on a particular step, these two won't drift too far apart. In reality, this turns out to be asking too much. Even when extremely small step sizes are used, over a large number of steps the error starts to accumulate and the two solutions part company! We can even go so far as to theoretically predict what kind of error the method will introduce.
Where the actual solution curve is
Ideally, we would like our "prediction line" to
(Of course, where the actual solution curve is
However, we should remind ourselves that our discussion is purely theoretical. We do Now let's think about how we might begin to fix this problem of the actual solution curve and the numerical solution drifting apart. Can we beat the "concavity problem?" Somehow we need to correct the over- or under-estimate problem that Euler's Method inevitably encounters. The Heun algorithm cleverly addresses this correction requirement. Let's take our concave-up example from above, and consider it more carefully this time. Instead of focussing on the initial, left end-point where we formed our tangent line, consider the interval spanned by the tangent line segment as a whole. As we've already noted, the tangent line from the left underestimates the slope of the curve for the entire width of the interval from the current point to the next predicted point, and as a result of this, the points along the tangent line have vertical coordinates which are all underestimates of those of the points lying along the actual solution curve, including the all-important right end-point of the interval under consideration. (After all, it is this point which forms the next component of our evolving numerical solution.)
To overcome this deficiency we would need to have used a line with a
Heun's clever approach to this problem is to consider the tangent lines to the solution curve at
Obviously its slope is
Now look at the relationship between the errors made by both our left and right tangent line predictions. One is The point we really want, i.e. the "ideal point" discussed earlier, appears to lie approximately halfway between these erroneous estimates based on tangent lines. In order to get this "ideal point" we still need to ride along an "ideal prediction line," but what should we use as this line's slope? If you've really been following the foregoing discussion you can probably guess by now what slope we'll use. We will take the So that's the basic idea behind Heun's method—using a prediction line whose slope is the average of the slopes of the tangent lines at either end of the interval.
Now we have to tackle a problem that I hinted at earlier: We don't actually
(As a side note, numerical analysts would refer to
OK, since we've managed to conceptualize the Heun algorithm graphically, it's now time to develop the method symbolically so that we can more easily see how to write the We already have our iterative formulas from the previous lab on Euler's method to get the ball rolling. These were:
h
h f(x, _{n}y)_{n}where y) represented the known (left) point, and (_{n}x_{n+1}, y_{n+1}) represented the new (right)
point. Remember, the reason we care about these formulas for our new Heun method is that we'll be using Euler's method to make a rough prediction of the location of the predicted next point so that these coordinates may be used for our estimate of the slope of the tangent line at the right end of the interval in question.
So let's get those The y), so the _{n}slope of the tangent line at the left end-point, just as with Euler's method, is the right hand side of the original differential equation evaluated at this point:slope = _{left}f(x, _{n}y)
_{n}
Now for that x_{n+1}, y_{n+1}) = (x + _{n}h, y +
_{n}h f(x, _{n}y))
_{n}
Recalling that one of our fundamental assumptions, based on our interpretation of the original differential equation, is that the quantity slope = _{right}f(x + _{n}h, y +
_{n}h f(x, _{n}y))
_{n}
Next, recall from the discussion and the graphs we considered above that the slope of our "ideal prediction line" is the slope = (1/2) (_{ideal}slope + _{left}slope)
_{right}All that is left to do now is to utilize this slope in creating the "ideal prediction line" itself, and then use this "ideal prediction line" to actually find the coordinates of the right hand end-point. To a great extent the remainder of the work is simply a repetition of what we saw as we developed the formulas for Euler's method in the last lab.
Consider the following illustration. We wish to predict the right hand end-point's coordinates. In order to get them it would be sufficient to find the value of Δ
Using the elementary idea that
y / h
which is easily rearranged to give us a formula for Δ y = h slope_{ideal}
Finally, then, we can predict the coordinates of the x-coordinate plus the step size, h. The next is the current y-coordinatey-coordinate plus Δy. Formulaically, this would be:
x_{n+1} = x_{n} + h
and y_{n+1} = y_{n} + Δy
Replacing Δ y_{n+1} = y_{n} + h slope_{ideal}
And replacing y_{n+1} = y_{n} + (1/2) h (slope + _{left}slope)
_{right}One last substitution, and our formula is complete. We recently found that: slope = _{left}f(x, _{n}y)
_{n}and that slope = _{right}f(x + _{n}h, y +
_{n}h f(x, _{n}y))
_{n}
so if we substitute these values into the above equation for y_{n+1} = y_{n} + (1/2) h (f(x, _{n}y) + _{n}f(x + _{n}h, y +
_{n}h f(x, _{n}y)))
_{n}or, cleaning things up slightly: y_{n+1} = y_{n} + (h/2) (f(x, _{n}y) + _{n}f(x + _{n}h, y +
_{n}h f(x, _{n}y)))
_{n}
y) + _{n}f(x + _{n}h, y +
_{n}h f(x, _{n}y)))
_{n}
It's now time to implement these newly minted formulas in |
|||

ODE Laboratories: A Sabbatical Project by Christopher A. Barker©2009 San Joaquin Delta College, 5151 Pacific Ave., Stockton, CA 95207, USA e-mail: cbarker@deltacollege.edu |