Heath - Scientific Computing (523150), страница 80
Текст из файла (страница 80)
INITIAL VALUE PROBLEMS FOR ODESthe second-order explicit method with stepsize h = 0.25 to take a step from t1 = 0.25 tot2 = 0.5, obtaining the predicted value1ŷ2 = y1 + (3y10 − y00 )h = 0.9375 + 0.5(−1.3184 + 0)0.25 = 0.7727.2We evaluate f at this predicted value ŷ2 to obtain the corresponding derivative value ŷ20 =−0.5971. We can now use these predicted values in the corresponding implicit method (inthis case the trapezoid rule) to obtain the corrected solution value1y2 = y1 + (y20 + y10 )h = 0.9375 + 0.5(−0.5971 − 0.4395)0.25 = 0.8079.2We evaluate f again using this new value y2 to obtain the improved value y20 = −0.6528,which would be needed in taking further steps.
At this point we have completed the PECEprocedure for this step. The corrector could be repeated, if desired, until convergence isobtained.For comparison, the exact solution for this problem is y(t) = 1/(1 + t2 ), and hence thetrue solution at the integration points is y(0.25) = 0.9412 and y(0.5) = 0.8.One of the most popular pairs of multistep methods is the explicit fourth-order AdamsBashforth predictoryk+1 = yk +1000(55yk0 − 59yk−1+ 37yk−2− 9yk−3)h24and the implicit fourth-order Adams-Moulton correctoryk+1 = yk +100(9y 0 + 19yk0 − 5yk−1+ yk−2)h.24 k+1Backward differentiation formulas (BDF), due to Gear, form another important family ofimplicit multistep methods.
BDF methods, typified by the formulayk+1 =16 0(18yk − 9yk−1 + 2yk−2 ) + yk+1h,1111have stability properties that make them particularly effective for solving stiff equations.The general properties of multistep methods can be summarized as follows:• They are not self-starting, because several previous solution values are required. Thus,a special starting procedure must be used initially, such as a single-step method, untilenough values have been generated to begin using a multistep method of the desiredorder.• Changing stepsize is complicated, since the interpolation formulas are most convenientlybased on equally spaced intervals for several consecutive points.• A good local error estimate can be determined from the difference between the predictorand the corrector.• They are relatively complicated to program.• Being based on interpolation, they can efficiently provide solution values at output pointsother than the integration points.9.6.
SURVEY OF NUMERICAL METHODS FOR ODES297• Implicit methods have a much greater region of stability than explicit methods but mustbe iterated to convergence to realize this benefit fully (e.g., a PECE scheme is actuallyexplicit, albeit in a somewhat complicated way).• Although implicit methods are more stable than explicit methods, they are still notnecessarily unconditionally stable. Indeed, no multistep method of greater than secondorder is unconditionally stable, even if it is implicit.• A properly designed implicit multistep method can be very effective for solving stiffequations.The stability and accuracy of some of the most popular multistep methods are summarized in Table 9.1, where “stability threshold” indicates the left endpoint of the stabilityinterval for a scalar equation, and “error constant” indicates the coefficient of the hp+1term in the local truncation error, where p is the order of the method.
All of these Adamsmethods have α1 = 1, and αi = 0 for i > 1, so we list only the βi . We observe that theimplicit methods are both more stable and more accurate than the corresponding explicitmethods of the same order.Order1234Table 9.1: Properties of multistep methodsExplicit MethodsStabilityβ1β2β3β4threshold1−23/2−1/2−123/12 −16/12 5/12−6/1155/24 −59/24 37/24 −9/24−3/10Errorconstant1/25/123/8251/720Implicit MethodsOrder12349.6.5β011/25/129/24β1β2β31/28/1219/24−1/12−5/241/24Stabilitythreshold−∞−∞−6−3Errorconstant−1/2−1/12−1/24−19/720Multivalue MethodsAs we have seen, changing stepsize is difficult with multistep methods because the pasthistory of the solution is most easily maintained at equally spaced intervals.
Like multistepmethods, multivalue methods are based on polynomial interpolation, but they avoid manyof the implementation difficulties associated with multistep methods.One of the key ideas motivating multivalue methods is the observation that the interpolating polynomial itself can be evaluated at any point, not just at equally spaced intervals.The equal spacing associated with multistep methods is simply an artifact of the way themethods are represented as a linear combination of successive solution and derivative values298CHAPTER 9.
INITIAL VALUE PROBLEMS FOR ODESwith fixed weights.Another key idea in implementing multivalue methods is choosing the representation ofthe interpolating polynomial so that its parameters are essentially the values of the solutionand one or more of its derivatives at a single point tk . This approach is analogous to usinga Taylor, rather than Lagrange, representation of the polynomial. The solution is advancedin time by a simple transformation of this representation from one point to the next, andchanging the stepsize in doing so is easy. Multivalue methods turn out to be mathematicallyequivalent to multistep methods, but multivalue methods are more convenient and flexibleto implement, so most modern software implementations are based on them.Example 9.17 Multivalue Method. To make these ideas a bit more concrete, weconsider a four-value method for solving an ODEy 0 = f (t, y).Instead of representing the interpolating polynomial by its value at four different points, werepresent it by its value and the values of its first three derivatives at a single point tk ,yk hy 0kyk = (h2 /2)y 00 ,k(h3 /6)yk000where the solution value and derivative values indicated are approximations to those of thetrue solution.
For convenience, the derivatives are scaled to match the coefficients in aTaylor series expansion. By differentiating the Taylor seriesy(tk + h) = y(tk ) + hy 0 +h2 00 h3 000y + yk + · · ·26three times, we see that the corresponding values at the next point tk+1 = tk + h are givenapproximately by the transformationŷk+1 = Byk ,where the matrix B is given by10B=001100121013.31We have not yet used the differential equation, however, so we add a correction term to theforegoing prediction to obtain the final valueyk+1 = ŷk+1 + αr,where r is a fixed 4-vector and0α = h(f (tk+1 , yk+1 ) − ŷk+1).9.7.
SOFTWARE FOR ODE INITIAL VALUE PROBLEMS299For consistency, i.e., for the ODE to be satisfied, we must have r2 = 1; but the three remaining components of r can be chosen in various ways, resulting in different methods, analogousto the different choices of parameters in multistep methods. For example, the four-valuemethod with r = [ 38 1 34 16 ]T is equivalent to the implicit fourth-order Adams-Moultonmethod given in Section 9.6.4.We can now see why it is easy to change stepsize with a multivalue method: we needmerely rescale the components of yk to reflect the new stepsize.
Moreover, it is also easyto change the order of the method simply by changing the components of r. These twocapabilities, combined with sophisticated tests and strategies for deciding when to changeorder and stepsize, have led to the development of very powerful and efficient softwarepackages for solving ODEs based on variable-order/variable-step methods. Such routinesare analogous to adaptive quadrature routines (see Section 8.4.2) in that they automaticallyadapt to a given problem, varying the order and stepsize of the integration method asnecessary to meet the user-supplied error tolerance in an efficient manner.
Such routinesoften have options for solving either stiff or nonstiff problems, and some even detect stiffnessautomatically and select an appropriate method accordingly.The ability to change order easily also obviates the need for special starting procedures.With a variable-order/variable-step method, one can simply start with a first-order method,which requires no additional starting values, and let the automatic order/stepsize selectionprocedure increase the order as needed for the required accuracy.9.7Software for ODE Initial Value ProblemsTable 9.2 is a list of some of the software available for numerical solution of initial valueproblems for ordinary differential equations. Many of these routines have additional variants for special situations, such as root finding or sparse Jacobians.
Another importantcategory that we have not discussed is differential-algebraic systems, in which the solutionmust satisfy a system containing both differential and algebraic equations. The best-knownroutine for solving such problems is dassl, which is available from netlib.Software for solving an ODE y 0 = f (t, y) typically requires the user to supply the nameof a routine that computes the value of the function f for any given values of t and y.Additional input includes the number of equations in the system; the initial values of theindependent variable t and the vector y of dependent variables at the start of the integration;the value tout of the independent variable at which the integration is to stop; and absoluteor relative error tolerances, or both. Additional input, especially for a stiff ODE solver, mayinclude the name of a routine for computing the Jacobian of f and the name of an arrayto be used as workspace for storing such matrices.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
