Download (direct link):
h t = 0.90 t= 1.0
0.1 14.02182 735.0991
0.05 14.27117 1.75863 x 105
0.01 14.30478 2.0913 x 102893
Chapter 8. Numerical Methods
While the values at t = 0.90 are reasonable and we might well believe that the solution has a value of about 14.305 at t = 0.90, it is clear that something strange is happening between t = 0.9 and t = 1.0. To help determine what is happening let us turn to some analytical approximations to the solution of the initial value problem (3). Note that on 0 < t < 1
Ó2 < t2 + y2 < 1 + y2. (4)
This suggests that the solution ó = ô() of
/= 1 + ó2, ó (0) = 1 (5)
and the solution ó = ô2 (t) of
Ó = ó2, y(0) = 1 (6)
are upper and lower bounds, respectively, for the solution ó = ô(t) of the original problem, since all these solutions pass through the same initial point. Indeed, it can be shown (for example, by the iteration method of Section 2.8) that ô2(³) < ô() < ô()
as long as these functions exist. The important thing to note is that we can solve Eqs. (5)
and (6) for ô1 and ô2 by separation of variables. We find that
/ n \ 1
ô1(Ã) = tan ^t + -j , ô2(0 = . (7)
Thus ô2(³) ^òî as t ^ 1, and ô() ^òî as t ^ n/4 = 0.785. These calculations show that the solution ofthe original initial value problem exists at least for 0 < t <n/4 and at most for 0 < t < 1. The solution of the problem (3) has a vertical asymptote for some t in n/4 < t < 1 and thus does not exist on the entire interval 0 < t < 1.
Our numerical calculations, however, suggest that we can go beyond t = n/4, and probably beyond t = 0.9. Assuming that the solution of the initial value problem exists at t = 0.9 and has the value 14.305, we can obtain a more accurate appraisal of what happens for larger t by considering the initial value problems (5) and (6) with y(0) = 1 replaced by y(0.9) = 14.305. Then we obtain
ô1 (t) = tan(t + 0.60100), ô2(t) = 1/(0.96991 - t), (8)
where only five decimal places have been kept. Thus ô() as t ^ n/2 -
0.60100 = 0.96980 and ô2(³) ^òî as t ^ 0.96991. We conclude that the asymptote of the solution of the initial value problem (3) lies between these two values. This example illustrates the sort of information that can be obtained by a judicious combination of analytical and numerical work.
Stability. The concept of stability is associated with the possibility that small errors that are introduced in the course of a mathematical procedure may die out as the procedure continues. Conversely, instability occurs if small errors tend to increase, perhaps without bound. For example, in Section 2.5 we identified equilibrium solutions of a differential equation as (asymptotically) stable or unstable, depending on whether solutions that were initially near the equilibrium solution tended to approach it or to depart from it as t increases. Somewhat more generally, the solution of an initial value problem is asymptotically stable if initially nearby solutions tend to approach the given solution, and unstable if they tend to depart from it. Visually, in an asymptotically stable problem the graphs of solutions will come together, while in an unstable problem they will separate.
8.5 More on Errors; Stability
If we are solving an initial value problem numerically, the best that we can hope for is that the numerical approximation will mimic the behavior of the actual solution. We cannot make an unstable problem into a stable one merely by solving it numerically. However, it may well happen that a numerical procedure will introduce instabilities that were not part of the original problem, and this can cause trouble in approximating the solution. Avoidance of such instabilities may require us to place restrictions on the step size h.
To illustrate what can happen in the simplest possible context consider the differential equation
dy/dt = ry, (9)
where r is a constant. Suppose that in solving this equation we have reached the point (tn, yn). Let us compare the exact solution of Eq. (9) that passes through this point, namely,
y = Óï exp[r(t — Ï)]> with numerical approximations obtained from the Euler formula
Óï+1 = Óï + hf(tn ’ Óï ) and from the backward Euler formula
Óï+1 = Óï + hf(tn+V Óï+1 Ó From the Euler formula (11) we obtain
Óï+1 = Óï + hryn = Óï (1 + rh)-Similarly, from the backward Euler formula (12) we have