Heath - Scientific Computing (523150), страница 77

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 77)

This result meansthat Euler’s method is first-order accurate.y....................0 ............... ......... ............................. .... ............... ..... ............ ..... ........0.... .... ................. .... ....... .. ....... .........

....... ......... ...... ..................... ............ ................... .... ................. ............ ............ ................. ..................................................... ............ ............ ............ ................. .......... .......... .......... .............. .................. ............................................. ......... ......... ............... ............................ .........

.......... ................. ...................... .......... ......... ...... ........... ........................ ........... ................. .. .......... ...................... ........... ............ ............ .......................... ............ ........... ............ .............................. .............. ..............

.............. ................................ ................ ...................... ................ .................................. ................................................. .................. ............................................................................................................. ................................................................................. .................................................................... ....................................................................................................................................................................................................................................................................................yy = −yt0t1t2t3t4global errortFigure 9.7: Local and global errors in Euler’s method for the ODE y 0 = −y.9.3.2Stability of a Numerical MethodThe concept of stability for numerical solutions of ODEs is analogous to, but distinct from,the concept of stability of the ODE itself.

Recall that an ODE is stable if its solutioncurves do not diverge from each other with time. Similarly, a numerical method is said tobe stable if small perturbations do not cause the resulting numerical solutions to divergefrom each other without bound (recall the general notion of stability in Section 1.2.7). Suchdivergence of numerical solutions could be due to instability of the ODE being solved, butit can also be due to the numerical method itself, even when solving a stable ODE.Example 9.9 Stability of Euler’s Method. From the derivation in Example 9.8 we seethat the global error is the sum of the local error and what might be termed the propagatederror .

To characterize the latter, note that by the Mean Value Theorem we can writef (tk , yk ) − f (tk , y(tk )) = J(ξ)(yk − y(tk ))for some (unknown) value ξ, so that we can express the global error at step k + 1 asEk+1 = (1 + hk J)Ek + Lk+1 .9.3. ACCURACY AND STABILITY285Thus, the global error is multiplied at each step by the factor (1 + hk J), which is called theamplification factor or growth factor . If |1 + hJ| < 1, then the errors do not grow, and themethod is stable. This condition is equivalent to requiring hJ to lie in the interval (−2, 0).If this is not the case, then the errors grow and the method is unstable.

Note that suchinstability could be due to instability of the ODE (i.e., J > 0), but it can also occur for astable equation (J < 0) if h > −2/J. We will see a dramatic example of such numericalinstability for a stable equation in Example 9.11.For a system of equations, the amplification factor for Euler’s method is the matrix(I + hJ ), and the condition for stability of the method is ρ(I + hJ ) < 1, which is satisfiedif the eigenvalues of hJ lie inside a circle in the complex plane of radius 1 and centered at−1 [notice that this includes the interval (−2, 0) of the single-equation case].In general, the amplification factor depends on the particular ODE being solved (whichdetermines the Jacobian J), the particular numerical method used (which determines theform of the amplification factor), and the stepsize h.An alternative approach to assessing the accuracy and stability of a numerical methodis to apply the method to the linear ODE y 0 = λy with initial condition y(0) = y0 , whoseexact solution is given by y(t) = y0 eλt .

This will enable us to determine the accuracy ofthe method by comparing the computed and exact solutions and to determine stability bycharacterizing the growth factor of the numerical solution.For example, applying Euler’s method to this equation using a fixed stepsize h, we haveyk+1 = yk + λyk h = (1 + λh)yk ,which means thatyk = (1 + λh)k y0 .Provided λ < 0, the exact solution decays to zero as t increases, as will the computedsolution if |1 + λh| < 1. This result agrees with our earlier stability analysis because J = λfor this ODE. We also note that the growth factor 1 + λh agrees with the series expansioneλh = 1 + λh +(λh)2 (λh)3++ ···26through terms of first order in h, and hence Euler’s method is first-order accurate.

Especiallyfor more complicated numerical methods, a linear ODE is easier to work with than a generalODE, and it produces essentially the same stability result if we equate λ with the Jacobian Jat a given point. An important caveat, however, is that λ is constant, whereas the JacobianJ varies for a nonlinear equation, and hence the stability can potentially change.9.3.3Stepsize ControlIn choosing a stepsize h for advancing the numerical solution of an ODE we would like totake as large a step as possible to minimize computational cost, but we must also take intoaccount both stability and accuracy.

Obviously, to yield a meaningful solution, the stepsizemust obey any stability restrictions imposed by the method being used. In addition, alocal error estimate is needed to ensure that the desired accuracy is attained. With Euler’s286CHAPTER 9. INITIAL VALUE PROBLEMS FOR ODESmethod, for example, we know that the local error is approximately (y 00 /2)h2 , and hencewe should choose the stepsize so thath ≤ (2tol/|y 00 |)1/2 ,where tol is the specified local error tolerance.

Of course, we do not know the value of y 00 ,but we can estimate it by a difference quotient of the formy 00 ≈0yk0 − yk−1.tk − tk−1Other methods of obtaining local error estimates are based on the difference between resultsobtained using methods of different orders or different stepsizes.9.4Implicit MethodsEuler’s method is an explicit method in that it uses only information at time tk to advancethe solution to time tk+1 .

This may appear to be a virtue, but we saw that Euler’s methodhas a rather limited stability interval of (−2, 0). A larger stability region can be obtainedby using information at time tk+1 , which makes the method implicit. The simplest exampleis the backward Euler method ,yk+1 = yk + f (tk+1 , yk+1 )hk .This method is implicit because we must evaluate f with the argument yk+1 before we knowits value.

This statement simply means that a value for yk+1 that satisfies the precedingequation must be determined, and if f is a nonlinear function of y, as is often the case, thenan iterative solution method, such as fixed-point iteration or Newton’s method, must beused. A good starting guess for the iteration can be obtained from an explicit integrationmethod, such as Euler’s method, or from the solution at the previous time step.Example 9.10 Backward Euler Method.

Consider the nonlinear ODEy 0 = −y 3with initial condition y(0) = 1. Using the backward Euler method with a stepsize of h = 0.5,we obtain the equationy1 = y0 + f (t1 , y1 )h = 1 − 0.5y13for the solution value at the next step. This nonlinear equation for y1 is already set up tosolve by fixed-point iteration, repeatedly substituting successive values for y1 on the righthand side, or we could use any other method from Chapter 5, such as Newton’s method. Inany case, we need a starting guess for y1 , for which we could simply use the previous solutionvalue, y0 = 1, or we could use an explicit method to produce a starting guess for the implicitmethod. Using Euler’s method, for example, we would obtain y1 = y0 − 0.5y03 = 0.5 as astarting guess for the iterative solution of the implicit equation. The iterations eventuallyconverge to the final value y1 ≈ 0.7709.9.4. IMPLICIT METHODS287Given the extra trouble and computation in using an implicit method, one might wonderwhy we would bother.

The answer is that implicit methods generally have a significantlylarger stability region than comparable explicit methods. To determine the stability of thebackward Euler method, we apply it to the linear ODE y 0 = λy, obtainingyk+1 = yk + λyk+1 h,or(1 − λh)yk+1 = yk ,so thatyk =11 − λhky0 .Thus, to mimic the exponential decay of the exact solution when λ < 0, we must have|1/(1 − λh)| < 1. Moreover, the growth factor1= 1 + λh + (λh)2 + · · ·1 − λhagrees with the expansion for eλh through terms of order h, so the backward Euler methodis first-order accurate.More generally, the amplification factor for the backward Euler method for a scalarequation is 1/(1 − hJ), which is less than 1 in magnitude for any positive h provided thatJ < 0.

Характеристики

Список файлов книги