Heath - Scientific Computing (523150), страница 78
Текст из файла (страница 78)
Thus, the stability interval for the backward Euler method is (−∞, 0), or theentire left half of the complex plane in the case of a system of equations, and hence for astable equation the method is stable for any positive stepsize.
Such a method is said to beunconditionally stable (other terms sometimes used for this concept are absolutely stable,A-stable, or A0 -stable). The great virtue of an unconditionally stable method is that thedesired local accuracy places the only constraint on our choice of stepsize.
Thus, we may beable to take much larger steps than for an explicit method of comparable order and attainmuch higher overall efficiency despite requiring more computation per step.Although the backward Euler method is unconditionally stable, its first-order accuracyseverely limits its usefulness. We can obtain a method of higher-order accuracy by combiningthe Euler and backward Euler methods. In particular, averaging these two methods yieldsthe implicit trapezoid ruleyk+1 = yk +f (tk , yk ) + f (tk+1 , yk+1 )hk .2To determine the stability and accuracy of this method, we apply it to the linear ODEy 0 = λy, obtainingλyk + λyk+1yk+1 = yk +h,2which implies that1 + λh/2 kyk =y0 .1 − λh/2288CHAPTER 9. INITIAL VALUE PROBLEMS FOR ODESThus, the method is stable if |(1 + λh/2)/(1 − λh/2)| < 1, which is true for any positivevalue of h provided λ < 0.
In addition, the growth factor! 2 31 + λh/2λhλhλhλh=1+1++++ ···1 − λh/22222= 1 + λh +(λh)2 (λh)3++ ···24agrees with the expansion of eλh through terms of order h2 , and hence the trapezoid methodis second-order accurate.More generally, the trapezoid rule has amplification factor (1 + hJ/2)/(1 − hJ/2), whichis less than 1 in magnitude for any positive stepsize provided that J < 0. The resultingstability regions are the interval (−∞, 0) for a scalar equation and the entire left half of thecomplex plane for a system of equations. Thus, the trapezoid rule is unconditionally stableas well as second-order accurate.We have now seen two examples of implicit methods that are unconditionally stable,but not all implicit methods have this property.
Implicit methods generally have largerstability regions than explicit methods, but the allowable stepsize is not always unlimited.Implicitness is not sufficient to guarantee stability, and stability is not sufficient to guaranteeaccuracy.9.5Stiff Differential EquationsThe solution curves for a stable equation converge with time. This convergence has thefavorable property of damping errors in a numerical solution, but if it is too rapid, asillustrated in Fig. 9.8, then difficulties of a different type may arise. Such an equation issaid to be stiff .y...........................................................
.................................. ...................................................................... ....................................................................... ......................................................... ......................................................................................................................................................................................................................................................................................................................................................................tFigure 9.8: The family of solution curves for a typical stiff ODE.Formally, a stable system of ODEs is stiff if the eigenvalues of its Jacobian matrix Jhave greatly differing magnitudes.
There may be an eigenvalue with a large negative realpart (corresponding to a strongly damped component of the solution) or a large imaginarypart (corresponding to a rapidly oscillating component of the solution). Such a differentialequation corresponds to a physical process whose components have disparate time scales ora process whose time scale is small compared to the interval over which it is being studied.9.5. STIFF DIFFERENTIAL EQUATIONS289Some numerical methods are very inefficient for stiff equations because the rapidly varying component of the solution forces very small stepsizes to be used to maintain stability.Since the stability restriction depends on the rapidly varying component of the solution,whereas the accuracy restriction depends on the slowly varying component, the stepsizemay be much more severely restricted by stability than by the required accuracy. For example, Euler’s method with a fixed stepsize is unstable for solving a stiff equation, whereasthe implicit backward Euler method is stable for stiff problems.
Stiff ODEs need not bedifficult to solve numerically provided a suitable method is chosen.Example 9.11 Stiff ODE. To illustrate the numerical solution of a stiff ODE, considerthe equationy 0 = −100y + 100t + 101with initial condition y(0) = 1. The general solution of this ODE is y(t) = 1 + t + ce−100t ,and the particular solution satisfying the initial condition is y(t) = 1 + t (i.e., c = 0). Sincethe solution is linear, Euler’s method is theoretically exact for this problem.
However, toillustrate the effect of truncation or rounding errors, let us perturb the initial value slightly.With a stepsize h = 0.1, the first few steps for the given initial values are:tExact solutionEuler solutionEuler solution0.01.000.991.010.11.101.191.010.21.200.392.010.31.308.59−5.990.41.40−64.267.0The computed solution is incredibly sensitive to the initial value, as each tiny perturbationresults in a wildly different solution. An explanation for this behavior is shown in Fig. 9.9.Any point deviating from the desired particular solution, even by only a small amount,lies on a different solution curve, for which c 6= 0, and therefore the rapid transient ofthe general solution is present.
Euler’s method bases its projection on the derivative atthe current point, and the resulting large value causes the numerical solution to divergeradically from the desired solution. This behavior should not surprise us. The Jacobian forthis equation is J = −100, so the stability condition for Euler’s method requires a stepsizeh < 0.02, which we are violating.2.0•....................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................Euler solution1.0 •desired solution•desired solutiontransient solution0.00.10.2Figure 9.9: Unstable solution of stiff ODE using Euler method.By contrast, the backward Euler method has no trouble solving this problem.
In fact,the backward Euler solution is extremely insensitive to the initial value, as shown in thefollowing table,290CHAPTER 9. INITIAL VALUE PROBLEMS FOR ODEStExact solutionBE solutionBE solution0.01.000.002.000.11.101.011.190.21.201.191.210.31.301.301.300.41.401.401.40and illustrated in Fig. 9.10. Even with a very large perturbation in the initial value, byusing the derivative at the next point rather than the current point, the transient is quicklydamped out and the backward Euler solution converges to the desired solution curve afteronly a few steps.
This behavior is consistent with the unconditional stability of the backwardEuler method for a stable equation.2.0desired solution•BE solution•....................................................................................................................................................................................................................................................................................................................................................................................................................................................................................1.0 ..............................................................•0.00.10.2Figure 9.10: Stable solution of stiff ODE using backward Euler method.9.6Survey of Numerical Methods for ODEsHaving covered the basic concepts of solving ordinary differential equations numerically, wenow briefly survey each of the major categories of methods for such problems.9.6.1Taylor Series MethodsWe have already seen that Euler’s method can be derived from a Taylor series expansion.
Byretaining more terms in the Taylor series, we can generate higher-order single-step methods.For example, retaining one additional term in the Taylor seriesy(t + h) = y(t) + y 0 (t)h +y 00 (t) 2 y 000 (t) 3h +h + ···26gives the second-order methodyk+1 = yk + yk0 hk +yk00 2h .2 kNote, however, that this approach requires the computation of higher derivatives of y. Thesecan be obtained by differentiating y 0 = f (t, y) using the chain rule, e.g.,y 00 = ft (t, y) + fy (t, y)y 0 = ft (t, y) + fy (t, y)f (t, y),where the subscripts indicate partial derivatives with respect to the given variable.
As theorder increases, such expressions for the derivatives rapidly become too complicated to be9.6. SURVEY OF NUMERICAL METHODS FOR ODES291practical to compute, so Taylor series methods of higher order have not often been usedin practice. Recently, however, the availability of symbolic manipulation and automaticdifferentiation systems has made these methods more feasible.Example 9.12 Taylor Series Method.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
