Introduction Notes Theory HOWTO Examples
Engineering Error Questions Matlab Maple
Consider the IVP
y(0) = 0
The solution to this IVP is the function
A plot of this solution is shown in Figure 1.
Figure 1. The solution to the given initial value problem.
If we use Euler's method to approximate a solution to this IVP using a step size of h = 0.1 for a few steps, we get the following values:
y(0.2) ≈ y2 = -0.0195162582
y(0.3) ≈ y3 = 0.1033409593
y(0.4) ≈ y4 = -0.0395932331
We immediately note a mistake with y2 which is negative, and rather than converging on the solution, the subsequent approximations appear to diverge. Figure 2 shows the solution and the first twenty iterations using Euler's method.
Figure 2. The first 20 approximations using Euler's method.
To demonstrate this more clearly, we can draw in the slopes followed by Euler's method at each step (drawing in the slopes), as is demonstrated in Figure 3.
Figure 3. Euler's method approximations together with slopes.
To understand this better, Figure 4 focuses on the solution and the directional field in the close vicinity of the solution.
Figure 4. The direction field plot for the given differential equation with the solution for the initial value y(0) = 0.
We note that even close to the solution, the slope of the direction field are very large and positive below the solution and negative above the solution. Such a differential equation is termed stiff. One method of solving such problems is to use a smaller value of h, however this may be unnecessarily expensive.
A better technique (rather than making h smaller) we will look at is a modification of Euler's method which is called backward-Euler's method.
Note: this problem was specifically tailored to fail for Euler's method with a step size of h = 0.1. While better techniques like 4th-order Runge Kutta will solve this problem satisfactorily, the problem can be made more stiff by substituting -21 with a larger negative number.
Backward-Euler's Method
For Euler's method, we approximate y(t1) with
where h = t1 − t0.
Instead, evaluate the differential equation at t = t1:
Replace the left-hand side with the backward divided-difference formula:
Next, replace y(ti) ≈ yi:
This gives us an implicit equation for y1:
We can now solve for y1 using any of our standard root-finding techniques. For example, substitute y1 = υ and rewrite the equation as
and let υ0 = y0 and υ1 = y0 + h f(t0, y0), that is, y(t0) and Euler's approximation to y(t1), respectively. We can now use the secant method to find solve for υ and we let y1 be this approximation.
We can then repeat this process arbitrarily often to approximate y2, y3, etc.
To demonstrate this in our example, we will work out two steps:
Finding y1
Given t0 = 0, t1 = 0.1, and y0 = 0, if we substitute these into:
we get the (rather simple) equation
In this case, we can trivially solve for υ = 0.02918830381, and therefore y1 = 0.02918830381.
Finding y2
To find y2, we use the equation
and substitute t2 = 0.2 and y1 = 0.02918830381 to get:
Again, this may be trivially solved for υ = 0.03582625132, and therefore we set y2 = 0.03582625132.
Further Solutions
If we continue this process, we get the approximations shown in Figure 5.
Figure 5. Approximations using the backward-Euler's method.
Higher-order Backward-Euler Methods
We have used the first-order backward divided-difference formula. There is no reason we could not have substituted the second-order backward divided-difference such as:
which, given two initial values y0 and y1, we could approximate y2, and then using y1 and y2, we could approximate y3, and so on.
Why does this work?
The critical question which you are probably asking is: why does this technique work when Euler's method does not? The following is a heuristic: consider the plot in Figure 6 which is a zoom of the solution shown in Figure 4.
Figure 6. A zoom of Figure 4.
Around the solution, all the direction arrows point in approximately the same direction of the solution, while further away, they tend to be more vertical.
If we apply Euler's method, we are following the slope, for example, in Figure 7 we follow the slope at the point (0.1, y(0.1)). This very quickly leads off into regions where the slope is large and negative, and therefore the next approximation with Euler's method will not point anywhere near the actual solution.
Figure 7. The tangent line at the point (0.1, y(0.1)).
If we apply Euler's method at (0.1, y(0.1)) for various values of h, Figure 8 shows the next iteration. Very quickly, the second iteration begins to diverge very quickly from the actual solution.
Figure 8. Two steps of Euler's method for values of h between 0.001 and 0.15.
Alternatively, the backward-Euler's method works as follows: going forward a distance h, what is the eligible point which could have had (0.1, y(0.1)) as a point on the solution? Figure 9 chooses a number of points around 0.2 and the lines passing through those points.
Figure 9. 21 various (0.2, y1) values and the backward paths.
Notice how only one of the points could have had (0.1, y(0.1)) as its predecessor. That point is approximately (0.2, 0.39) (the actual value being (0.2, 0.03902971461)). Figure 10 shows both Euler's method and the backward-Euler's method with h = 0.1.
Figure 10. Euler's and the backward-Euler's method at (0.1, y(0.1)) with h = 0.1.
The second-order backward-Euler's method uses a second-order backward divided difference approximation of the derivative. Given (0, 0) and (1, y(1)), Figure 11 shows the approximation of y(0.2) using this technique. In this case, the 2nd-order technique is not as good as the 1st-order technique demonstrated in Figure 10, however, for smaller values of h it becomes more accurate.
Figure 11. The second-order backward-Euler's method using (0, 0) and (0.1, y(0.1)).
Copyright ©2005 by Douglas Wilhelm Harder. All rights reserved.
