Download (direct link):
Óï+1 = Óï + hryn+1’
Óï = Óï [1 + rh + (rh)2 + •••].
1 — rh
Finally, evaluating the solution (10) at t + h, we find that
Óï+1 = Óï exp(rh) = Óï
1 + rh +
Comparing Eqs. (13), (14), and (15) we see that the errors in both the Euler formula and in the backward Euler formula are of order h2, as the theory predicts.
Now suppose that we change the value yn to yn + S. Think, if you wish, of S as the error that has accumulated by the time we reach t = tn. In going one more step to tn+1, the question is whether this error increases or decreases.
For the exact solution (15) the change in yn+1 due to the change S in yn is just
S exp(rh). This quantity is less than S if exp(rh) < 1, that is, if r < 0. This confirms our conclusion in Chapter 2 that Eq. (9) is asymptotically stable if r < 0 and unstable if r > 0.
For the backward Euler method the change in yn+1 in Eq. (14) due to S is S/(1 — rh). For r < 0 this quantity is always nonnegative and never greater than 1. Thus, if the differential equation is stable, then so is the backward Euler method for an arbitrary step size h.
Chapter 8. Numerical Methods
On the other hand, for the Euler method, the change in yn+1 in Eq. (13) due to S is 5(1 + rh). If we recall that r < 0 and require that |1 + rhi < 1, then we find that h must satisfy h < 2/|r |. Thus the Euler method is not stable for this problem unless h is sufficiently small.
The restriction on the step size h in using the Euler method in the preceding example is rather mild unless ir | is quite large. Nonetheless, the example illustrates that it may be necessary to restrict h in order to achieve stability in the numerical method, even though the initial value problem itself is stable for all values of h. Problems for which a much smaller step size is needed for stability than for accuracy are called stiff. The backward differentiation formulas described in Section 8.4 (of which the backward Euler formula is the lowest order example) are the most popular formulas for solving stiff problems. The following example illustrates the kind of instabilities that can occur in trying to approximate the solution of a stiff problem.
A Stiff Problem
Consider the initial value problem
Ó = — 100y + 100? + 1, y(0) = 1.
Find numerical approximations to the solution for 0 < t < 1 using the Euler, backward Euler, and Runge-Kutta methods. Compare the numerical results with the exact solution.
Since the differential equation is linear, it is easy to solve, and the solution of the initial value problem (16) is
Ó = ô(ã) = e~100t + t.
Some values of the solution ô(´), correct to six decimal places, are given in the second column of Table 8.5.3, and a graph of the solution is shown in Figure 8.5.2. There is a thin layer (sometimes called a boundary layer) to the right of t = 0 in which the exponential term is significant and the solution varies rapidly. Once past this layer, however, ô(´) = t and the graph of the solution is essentially a straight line. The width of the boundary layer is somewhat arbitrary, but it is certainly small. At t = 0.1, for example, exp(—100t) = 0.000045.
FIGURE 8.5.2 The solution of the initial value problem (16).
8.5 More on Errors; Stability
If we plan to approximate the solution (17) numerically, we might intuitively expect that a small step size will be needed only in the boundary layer. To make this expectation a bit more precise recall from Section 8.1 that the local truncation errors for the Euler and backward Euler methods are proportional to ô" (t). For this problem ô"(²) = 104 e-100t, which varies from a value of 104 at t = 0 to nearly zero for t > 0.2. Thus a very small step size is needed for accuracy near t = 0, but a much larger step size is adequate once t is a little larger.
On the other hand, the stability analysis in Eqs. (9) through (15) also applies to this problem. Since r = —100 for Eq. (16), it follows that for stability we need h < 0.02 for the Euler method, but there is no corresponding restriction for the backward Euler method.
Some results using the Euler method are shown in columns 3 and 4 of Table 8.5.3. The values obtained for h = 0.025 are worthless due to instability, while those for h = 0.01666... are reasonably accurate for t > 0.2. However, comparable accuracy for this range of t is obtained for h = 0.1 using the backward Euler method, as shown by the results in column 7 of the table.
The situation is not improved by using instead of the Euler method a more accurate one, such as Runge-Kutta. For this problem the Runge-Kutta method is unstable for h = 0.033 ... but stable for h = 0.025, as shown by the results in columns 5 and 6 in Table 8.5.3.
The results given in the table for t = 0.05 and for t = 0.1 show that in the boundary layer a smaller step size is needed to obtain an accurate approximation. You are invited to explore this matter further in Problem 3.