Deturck, Wilf - Lecture Notes on Numerical Analysis (523142), страница 16
Текст из файла (страница 16)
A total of p + 2 points are involved in theformula, and so we can call (2.12.1) a p + 2-point formula.The trapezoidal rule, for example, arises by taking p = 0, a0 = 1, b0 = 1/2, b−1 = 1/2.Euler’s method has p = 0, a0 = 1, b−1 = 0, b0 = 1, whereas for the midpoint rule we havep = 1, a0 = 0, a−1 = 1, b1 = 0, b0 = 2, b−1 = 0.We can recognize an explicit, or noniterative, formula by looking at b1 . If b1 is nonzero,then yn+1 appears implicitly on the right side of (2.12.1) as well as on the left, and theformula does not explicitly tell us the value of yn+1 .
Otherwise, if b1 = 0, the formula isexplicit. In either case, if p > 0 the formula will need help getting started or restarted,whereas if p = 0 it is self-starting.70The Numerical Solution of Differential EquationsThe general linear multistep formula contains 2p + 3 constantsa0 , .
. . , a−pand b1 , b0 , b−1 , . . . , b−p .These constants are chosen to give the method various properties that we may deem to bedesirable in a particular application. For instance, if we value explicit formulas, then wemay set b1 = 0 immediately, thereby using one of the 2p + 3 free parameters, and leaving2p + 2 others.One might think that the remaining parameters should be chosen so as to give thehighest possible accuracy, in some sense. However, for a fixed p, the more accuracy wedemand, the more we come into conflict with stability, another highly desirable feature.Indeed, if we demand “too much accuracy” for a fixed p, it will turn out that no stableformulas exist. An important theorem of the subject, due to Dahlquist, states roughly thatwe cannot use more than about half of the 2p + 3 “degrees of freedom” to achieve highaccuracy if we want to have a stable formula (and we do!).First let’s discuss the conditions of accuracy.
These are usually handled by asking thatthe equation (2.12.1) should be exactly true if the unknown function y happens to be apolynomial of low enough degree. For instance, suppose y(x) = 1 for all x. Substituteyk = 1 and yk0 = 0 for all k into (2.12.1), and there follows the conditionpXa−i = 1.(2.12.2)i=0Notice that in all three of the methods we have been studying, the sum of the a’s isindeed equal to 1.Now suppose we want our multistep formula to be exact not only when y(x) = 1, butalso when y(x) = x. We substitute yk = kh, yk0 = 1 for all k into (2.12.1) and obtain, aftersome simplification, and use of (2.12.2),−pXi=0ia−i +pXb−i = 1.(2.12.3)i=−1The reader should check that this condition is also satisfied by all three of the methods wehave studied.In general let’s find the condition for the formula to integrate the function y(x) = xrand all lower powers of x exactly for some fixed value of r.
Hence, in (2.12.1), we substituteyk = (kh)r and yk0 = r(kh)r−1 . After cancelling the factor hr , we get(n + 1)r =pXi=0a−i (n − i)r + rpXb−i (n − i)r−1 .(2.12.4)i=−1Now clearly we do xr and all lower powers of x exactly if and only if we do (x + c)r andall lower powers of x exactly. The conditions are therefore translation invariant, so we canchoose one special value of n in (2.12.4) if we want.2.12 The big leagues71Let’s choose n = 0, because the result simplifies then tor(−1) =pXi=0pXi a−i − rrir−1 b−ir = 0, 1, 2, .
. . .(2.12.5)i=−1A small technical remark here is that when r = 1 and i = 0 we see a 00 in the secondsum on the right. This should be interpreted as 1, in accordance with (2.12.3).By the order (of accuracy) of a linear multistep method we mean the highest power of xthat the formula integrates exactly, or equivalently, the largest number r for which (2.12.5)is true (together with its analogues for all numbers between 0 and r).The accuracy conditions enable us to take the role of designer, and to construct accurateformlas with desirable characteristics. The reader should verify that the trapezoidal rule isthe most accurate of all possible two-point formulas, and should seach for the most accurateof all possible three-point formulas (ignoring the stability question altogether).Now we must discuss the question of stability.
Just as in section 2.6, stability will bejudged with respect to the performance of our multistep formula on a particular differentialequation, namelyyy0 = −(L > 0)(2.12.6)Lwith y(0) = 1.To see how well the general formula (2.12.1) does with this differential equation, whosesolution is a falling exponential, substitute yk0 = −yk /L for all k in (2.12.1), to obtain(1 + τ b1 )yn+1 −pX(a−i − τ b−i )yn−i = 0(2.12.7)i=0where as in section 2.6, we have written τ = h/L, the ratio of the step size to the relaxationlength of the problem.Equation (2.12.7) is a linear difference equation with constant coefficients of order p + 1.If as usual with such equations, we look for solutions of the form αn , then after substitutionand cancellation we obtain the characteristic equation(1 + τ b1 )αp+1 −pX(a−i − τ b−i )αp−i = 0.(2.12.8)i=0This is a polynomial equation of degree p + 1, so it has p + 1 roots somewhere in thecomplex plane.
These roots depend on the value of τ , that is, on the step size that we useto do the integration.We can’t expect that these roots will have absolute value less than 1 however large wechoose the step size h. All we can hope for is that it should be possible, by choosing h smallenough, to get all of these roots to lie inside the unit disk in the complex plane.Just for openers, let’s see where the roots are when h = 0. In fact, if when h = 0 someroot has absolute value strictly greater than 1, then there is no hope at all for the formula,72The Numerical Solution of Differential Equationsbecause for all sufficiently small h there will be a root outside the unit disk, and the methodwill be unstable.Now when h = 0, τ = 0 also, so the polynomial equation (2.12.8) becomes simplyαp+1 −pXa−i αp−i = 0.(2.12.9)i=0This is also a polynomial equation of degree p + 1.
However, its coefficients don’t dependon the step size, or even on the values of the various b’s in the formula. The coefficients,and therefore the roots, depend only on the a’s that we use.Hence, our first necessary condition for the stability of the formula (2.12.1) is that thepolynomial equation (2.12.9) should have no roots of absolute value greater than 1.
Thereader should now check the three methods that we have been using to see if this conditionis satisfied.Next we have to study what happens to the roots when h is small and positive. Qualitatively, here’s what happens. One of the roots of (2.12.9) is always α = 1, because condition(2.12.2) is always satisfied in practice, and since all it asks is that we correctly integrate theequation y 0 = 0, which is surely not an excessive request.The one root of (2.12.9) that is = 1 is our good friend, and we’ll call it the principalroot. As h increases slightly away from 0 the principal root moves slightly away from 1.Let α1 (h) denote the value of the principal root at some small value of h.
Then α1 (0) = 1.Now for small positive h, it turns out that the principal root tries as hard as it can to be agood approximation to e−τ . This means that the portion of the solution of the differenceequation (2.12.7) that comes from the one root α1 (h) isα1 (h)n = (“nearly” e−τ )n= “nearly” e−nτ(2.12.10)= “nearly” e−nh/L .But e−nh/L is the exact solution of the differential equation (2.12.6) that we’re trying tosolve.
Therefore the principal root is trying to give us a very accurate approximation tothe exact solution.Well then, what are all of the other p roots of the difference equation (2.12.7) doing forus? The answer is that they are trying as hard as they can to make our lives difficult. Thehigh quality of our approximate solution derives from the nearness of the principal root toe−τ . This high quality is bought at the price of using a p + 2-point multistep formula, andthis forces the characteristic equation to be of degree p + 1, and hence the remining rootshave to be there, too.People have invented various names for these other, non-principal, roots of the differenceequations. One that is printable is “parasitic,” so we’ll call them the parasitic roots.
Wewould be happiest if they would all be zero, but failing that, we would like them to be assmall as possible in absolute value.2.12 The big leagues73A picture of a good multistep method in action, then, with some small value of h, showsone root of the characteristic equation near 1, more precisely, very near e−τ , and all of theother roots near the origin.Now let’s try to make this picture a bit more quantitative.
We return to the polynomialequation (2.12.8), and attempt to find a power series expansion for the principal root α1 (τ )in powers of τ . The expansion will be in the formα1 (τ ) = 1 + q1 τ + q2 τ 2 + · · ·(2.12.11)where the q’s are to be determined. If we substitue (2.12.11) into (2.12.8), we get(1 − τ b1 )(1 + q1 τ + q2 τ 2 + · · ·)p+1 −pX(a−i − τ b−i )(1 + q1 τ + q2 τ 2 + · · ·)p−i = 0. (2.12.12)i=0Now we equate the coefficient of each power of τ to zero. First, the coefficient of thezeroth power of τ is1−pXa−i(2.12.13)i=0and, according to (2.12.2), this is indeed zero if our multistep formula correctly integratesa constant function.Second, the coefficient of τ isb1 + (p + 1)q1 −pX[a−i (p − i)q1 + b−i ] = 0.(2.12.14)i=0However, if we use (2.12.2) and (2.12.3), this simplifies instantly to the simple statementthat q1 = −1.So far, we have shown by direct calculation that if our multistep formula is of orderat least 1 (i.e., if (2.12.4) holds for r = 0 and r = 1), then the expansion (2.12.11) of theprincipal root agrees with the expansion of e−τ through terms of first degree in τ .Much more is true, although we will not prove it here: the expansion of the principal rootagrees with the expansion of e−τ through terms of degree equal to the order of the multistepmethod under consideration.
The proof can be done just by continuing the argument thatwe began above, but we omit the details.Thus the careful determination of the a’s and the b’s so as to make the formula asaccurate as possible all result in the principal root being “nearly” e−τ . Equal care must betaken to assure that the parasitic roots stay small.We illustrate these ideas with a new multistep formula,yn+1 = yn−1 +h 00+ 4yn0 + yn−1).(y3 n+1(2.12.15)This, like the midpoint rule, is a three-point method (p = 1). It is iterative, becauseb1 = 1/3 =6 0, and it happens to be very accurate.
Indeed, we can quickly check that74The Numerical Solution of Differential Equationsthe accuracy conditions (2.12.5) are satisfied for r = 0, 1, 2, 3, 4. The method is offourth order, and in fact its error term can be found by the method of section 2.8 to be−h5 y (v) (X)/90, where X lies between xn−1 and xn+1 .Now we examine the stability of this method. First, when h = 0 the equation (2.12.9)that determines the roots is just α2 − 1 = 0, so the roots are +1 and −1.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
