Heath - Scientific Computing (523150), страница 79
Текст из файла (страница 79)
To illustrate the second-order Taylor seriesmethod, we use it to solve the ODEy 0 = f (t, y) = −2ty 2 ,with initial value y(0) = 1. We differentiate f to obtain for this problemy 00 = ft (t, y) + fy (t, y)f (t, y) = −2y 2 + (−4ty)(−2ty 2 ) = 2y 2 (4t2 y − 1).Taking a step from t0 = 0 to t1 = 0.25 using stepsize h = 0.25, we obtainy1 = y0 + y00 h +y000 2h = 1 + 0 − 0.0625 = 0.9375.2Continuing with another step from t1 = 0.25 to t2 = 0.5, we obtainy2 = y1 + y10 h +y100 2h = 0.9375 − 0.1099 − 0.0421 = 0.7856.2For comparison, the exact solution for this problem is y(t) = 1/(1 + t2 ), and hence the truesolution at the integration points is y(0.25) = 0.9412 and y(0.5) = 0.8.9.6.2Runge-Kutta MethodsRunge-Kutta methods are single-step methods that are similar in motivation to Taylor seriesmethods but do not require the computation of higher derivatives.
Instead, Runge-Kuttamethods simulate the effect of higher derivatives by evaluating f several times between tkand tk+1 .Example 9.13 Derivation of a Runge-Kutta Method. The basic idea of Runge-Kuttamethods is best illustrated by example, the simplest of which is Heun’s method . Recall fromSection 9.6.1 that the second derivative of y is given byy 00 = ft + fy f,where each function is evaluated at (t, y). We can approximate the term on the right byexpanding f in a Taylor series in two variablesf (t + h, y + hf ) = f + hft + hfy f + O(h2 ),from which we obtainft + fy f =f (t + h, y + hf ) − f (t, y)+ O(h2 ).h292CHAPTER 9. INITIAL VALUE PROBLEMS FOR ODESWith this approximation to the second derivative, the second-order Taylor series methodgiven in Section 9.6.1 becomesf (tk + hk , yk + hk f (tk , yk )) − f (tk , yk ) 2hk2hkf (tk , yk ) + f (tk + hk , yk + hk f (tk , yk ))= yk +hk ,2yk+1 = yk + f (tk , yk )hk +which can be implemented in the form1yk+1 = yk + (k1 + k2 ),2wherek1 = f (tk , yk )hk ,k2 = f (tk + hk , yk + k1 )hk .Heun’s method, which is of second-order accuracy, is analogous to the implicit trapezoidrule but remains explicit by using the Euler prediction yk + k1 instead of yk+1 in evaluatingf at tk+1 .Example 9.14 Heun’s Method.
To illustrate the use of Heun’s method, we use it tosolve the ODEy 0 = −2ty 2 ,with initial value y(0) = 1. Taking a step from t0 = 0 to t1 = 0.25 using stepsize h = 0.25,we obtaink1 = f (t0 , y0 )h = 0 and k2 = f (t0 + h, y0 + k1 )h = −0.125,so that1y1 = y0 + (k1 + k2 ) = 1 − 0.0625 = 0.9375.2Continuing with another step from t1 = 0.25 to t2 = 0.5, we obtaink1 = f (t1 , y1 )h = −0.1099so thatand k2 = f (t1 + h, y1 + k1 )h = −0.1712,1y2 = y1 + (k1 + k2 ) = 0.9375 − 0.1406 = 0.7969.2For comparison, the exact solution for this problem is y(t) = 1/(1 + t2 ), and hence the truesolution at the integration points is y(0.25) = 0.9412 and y(0.5) = 0.8.The best-known Runge-Kutta method is the classical fourth-order scheme1yk+1 = yk + (k1 + 2k2 + 2k3 + k4 ),69.6.
SURVEY OF NUMERICAL METHODS FOR ODES293wherek1 = f (tk , yk )hk ,k2 = f (tk + hk /2, yk + k1 /2)hk ,k3 = f (tk + hk /2, yk + k2 /2)hk ,k4 = f (tk + hk , yk + k3 )hk .This method is analogous to Simpson’s rule; indeed it is Simpson’s rule if f depends onlyon t. For an illustration of the use of the classical fourth-order Runge-Kutta method forsolving a system of ODEs, see Example 10.1.Runge-Kutta methods have a number of virtues. To proceed to time tk+1 , they requireno history of the solution prior to time tk , which makes them self-starting at the beginningof the integration, and also makes it easy to change stepsize during the integration. Thesefacts also make Runge-Kutta methods relatively easy to program, which accounts in partfor their popularity.Unfortunately, classical Runge-Kutta methods provide no error estimate on which tobase the choice of stepsize.
More recently, however, Fehlberg devised a Runge-Kutta methodthat uses six function evaluations per step to produce both fifth-order and fourth-orderestimates of the solution, whose difference provides an estimate for the local error. Thisapproach has led to automatic Runge-Kutta solvers that are effective for many problemsbut are relatively inefficient for stiff problems or when very high accuracy is required. It ispossible, however, to define implicit Runge-Kutta methods with superior stability propertiesthat are suitable for solving stiff equations.9.6.3Extrapolation MethodsExtrapolation methods are based on the use of a single-step method to integrate the ODEover a given interval, tk ≤ t ≤ tk+1 , using several different stepsizes hi and yielding resultsdenoted by Y (hi ).
This gives a discrete approximation to a function Y (h), where Y (0) =y(tk+1 ). An interpolating polynomial or rational function Ŷ (h) is fit to these data, andŶ (0) is then taken as the approximation to Y (0).We saw another example of this approach in Richardson extrapolation for numericaldifferentiation and integration (see Section 8.8). Extrapolation methods are capable ofachieving very high accuracy, but they tend to be much less efficient and less flexible thanother methods for ODEs, so they are used mainly when extremely high accuracy is requiredand cost is not a significant factor.9.6.4Multistep MethodsMultistep methods use information at more than one previous point to estimate the solutionat the next point.
For this reason, they are sometimes called methods with memory. Linearmultistep methods have the formyk+1 =nXi=1αi yk+1−i + hnXi=0βi f (tk+1−i , yk+1−i ).294CHAPTER 9. INITIAL VALUE PROBLEMS FOR ODESThe parameters αi and βi are determined by polynomial interpolation. If β0 = 0, themethod is explicit, but if β0 6= 0, the method is implicit.Example 9.15 Derivation of Multistep Methods. To illustrate the derivation ofmultistep methods, we derive an explicit two-step method of the form0yk+1 = α1 yk + (β1 yk0 + β2 yk−1)h,where the parameters α1 , β1 , and β2 are to be determined. Using the method of undetermined coefficients, we will force the formula to be exact for the first three monomials. Ify(t) = 1, then y 0 (t) = 0, so that we have the equation1 = α1 · 1 + (β1 · 0 + β2 · 0)h.If y(t) = t, then y 0 (t) = 1, so that we have the equationtk+1 = α1 tk + (β1 · 1 + β2 · 1)h.If y(t) = t2 , then y 0 (t) = 2t, so that we have the equationt2k+1 = α1 t2k + (β1 · 2tk + β2 · 2tk−1 )h.All three of these equations must hold for any values of the ti , so we make the convenientchoice tk−1 = 0, h = 1 (hence tk = 1 and tk+1 = 2) and solve the resulting 3 × 3 linearsystem to obtain the values α1 = 1, β1 = 23 , β2 = − 12 .
Thus, the resulting explicit two-stepmethod is10yk+1 = yk + (3yk0 − yk−1)h,2and by construction it is of order two.Similarly, we can derive an implicit two-step method of the form0yk+1 = α1 yk + (β0 yk+1+ β1 yk0 )h.Again using the method of undetermined coefficients, we force the formula to be exact forthe first three monomials, obtaining the three equations1 = α1 · 1 + (β0 · 0 + β1 · 0)h,tk+1 = α1 tk + (β0 · 1 + β1 · 1)h,t2k+1 = α1 t2k + (β0 · 2tk+1 + β1 · 2tk )h.Making the convenient choice tk = 0, h = 1 (hence, tk+1 = 1), we solve the resulting 3 × 3linear system to obtain the values α1 = 1, β1 = 12 , β2 = 12 .
Thus, the resulting implicittwo-step method is1 0yk+1 = yk + (yk+1+ yk0 )h,2which we recognize as the trapezoid rule, and by construction it is of order two. Higherorder multistep methods can be derived in this same manner, forcing the desired formula to9.6. SURVEY OF NUMERICAL METHODS FOR ODES295be exact for as many monomials as there are parameters to be determined and then solvingthe resulting system of equations for those parameters.Alternatively, multistep methods can also be derived by numerical quadrature.
Forexample, sincey(tk+1 ) = y(tk ) +Ztk+10y (t) dt = y(tk ) +Ztkf (t, y(t)) dt,tkwe can takeyk+1 = yk +tk+1Ztk+1p(t) dt,tkwhere p(t) is a polynomial interpolating f (t, y) at the points (tk+1−n , yk+1−n ), . . . , (tk , yk ) foran explicit method of order n, or (tk+2−n , yk+2−n ), . . . , (tk+1 , yk+1 ) for an implicit methodof order n.Since multistep methods require several previous solution values and derivative values,how do we get started initially, before we have any past history to use? One strategy is to usea single-step method, which requires no past history, to generate solution values at enoughpoints to begin using a multistep method. Another option is to use a low-order methodinitially and gradually increase the order as additional solution values become available.As we saw with single-step methods, implicit multistep methods are usually more accurate and stable than explicit multistep methods, but they require an initial guess tosolve the resulting (usually nonlinear) equation for yk+1 .
A good initial guess is conveniently supplied by an explicit method, so the explicit and implicit methods are used as apredictor-corrector pair. One could use the corrector repeatedly (i.e., fixed-point iteration)until some convergence tolerance is met, but doing so may not be worth the expense. So,a fixed number of corrector steps, often only one, may be used instead, giving a PECE(predict, evaluate, correct, evaluate) scheme. Although it has no effect on the value of yk+1 ,0for futurethe second evaluation of f in a PECE scheme yields an improved value of yk+1use.Alternatively, the nonlinear equation for yk+1 given by an implicit multistep method canbe solved by Newton’s method or other similar iterative method, again with a good startingguess supplied by the solution at the previous step or by an explicit multistep method.
Inparticular, Newton’s method or a close variant of it is essential when using an implicitmultistep method designed for stiff ODEs, as fixed-point iteration will fail to converge forreasonable stepsizes.Example 9.16 Predictor-Corrector Method. To illustrate the use of a predictorcorrector pair, we use the two multistep methods derived in Example 9.15 to solve theODEy 0 = −2ty 2 ,with initial value y(0) = 1. The second-order explicit method requires two starting values,so in addition to the initial value y0 = 1 at t0 = 0 we will also use the value y1 = 0.9375 att1 = 0.25 obtained using the single-step Heun method in Example 9.14. We can now use296CHAPTER 9.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
