Deturck, Wilf - Lecture Notes on Numerical Analysis (523142), страница 9
Текст из файла (страница 9)
The name arises from the fact that the first derivative yn0 isbeing approximated by the slope of the chord that joins the two points (xn−1 , yn−1 ) and(xn+1 , yn+1 ), instead of the chord joining (xn , yn ) and (xn+1 , yn+1 ) as in Euler’s method.At first sight it seems that (2.5.5) can be used just like Euler’s method, because it is arecurrence formula in which we compute the next value yn+1 from the two previous valuesyn and yn−1 . Indeed the rules are quite similar, except for the fact that we can’t get startedwith the midpoint rule until we know two consecutive values y0 , y1 of the unknown functionat two consecutive points x0 , x1 .
Normally a differential equation is given together withjust one value of the unknown function, so if we are to use the midpoint rule we’ll need tomanufacture one more value of y(x) by some other means.40The Numerical Solution of Differential EquationsThis kind of situation will come up again and again as we look at more accurate methods,because to obtain greater precision without computing higher derivatives we will get thenext approximate value of y from a recurrence formula that may involve not just one or two,but several of its predecessors. To get such a formula started we will have to find severalstarting values in addition to the one that is given in the statement of the initial-valueproblem.To get back to the midpoint rule, we can get it started most easily by calculating y1 ,the approximation to y(x0 + h), from Euler’s method, and then switching to the midpointrule to carry out the rest of the calculation.Let’s do this, for example with the same differential equation (2.1.7) that we used toillustrate Euler’s rule, so we can compare the two methods.
The problem consists of theequation y 0 = 0.5y and the initial value y(0) = 1. We’ll use the same step size h = 0.05 asbefore.Now to start the midpoint rule we need two consecutive values of y, in this case at x = 0and x = 0.05. At 0.05 we use the value that Euler’s method gives us, namely y1 = 1.025(see Table 1).
It’s easy to continue the calculation now from (2.5.5).For instancey2 = y0 + 2h(0.5y1 )= 1 + 0.1(0.5 ∗ 1.025)(2.5.6)= 1.05125andy3 = y1 + 2h(0.5y2 )= 1.025 + 0.1(0.5 ∗ 1.05125)(2.5.7)= 1.0775625 .In the table below we show for each x the value computed from the midpoint rule, fromEuler’s method, and from the exact solution y(x) = ex/2 . The superior accuracy of themidpoint rule is apparent.x0.000.050.100.150.200.25...Midpoint(x)1.000001.025001.051251.077561.105131.13282...Euler(x)1.000001.025001.050631.076891.103811.13141...Exact(x)1.000001.025321.051271.077881.105171.13315...1.002.003.00...1.648472.717634.48032...1.638622.685064.39979...1.648722.718284.48169...5.0010.0012.17743148.3127411.81372139.5638912.18249148.413162.5 The midpoint and trapezoidal rules41y = f (x)abFigure 2.1: The trapezoidal ruletable 2Next, we introduce a third method of numerical integration, the trapezoidal rule.
Thebest way to obtain it is to convert the differential equation that we’re trying to solve intoan integral equation, and then use the trapezoidal approximation for the integral.We begin with the differential equation y 0 = f (x, y(x)), and we integrate both sidesfrom x to x + h, gettingZx+hy(x + h) = y(x) +f (t, y(t)) dt.(2.5.8)xNow if we approximate the right-hand side in any way by a weighted sum of values ofthe integrand at various points we will have found an approximate method for solving ourdifferential equation.The trapezoidal rule states that for an approximate value of an integralZbf (t) dt(2.5.9)awe can use, instead of the area under the curve between x = a and x = b, the area of thetrapezoid whose sides are the x axis, the lines x = a and x = b, and the line through thepoints (a, f (a)) and (b, f (b)), as shown in Figure 2.1.
That area is 12 (f (a) + f (b))(b − a).If we apply the trapezoidal rule to the integral that appears in (2.5.8), we obtainy(xn + h) ≈ y(xn ) +h(f (xn , y(xn )) + f (xn + h, y(xn + h)))2(2.5.10)in which we have used the “≈” sign rather than the “=” because the right hand side is notexactly equal to the integral that really belongs there, but is only approximately so.If we use our usual abbreviation yn for the computed approximate value of y(xn ), then(2.5.10) becomeshyn+1 = yn + (f (xn , yn ) + f (xn+1 , yn+1 )).(2.5.11)242The Numerical Solution of Differential EquationsThis is the trapezoidal rule in the form that is useful for differential equations.At first sight, (2.5.11) looks like a recurrence formula from which the next approximatevalue, yn+1 , of the unknown function, can immediately be computed from the previousvalue, yn .
However, this is not the case.Upon closer examination one observes that the next value yn+1 appears not only on theleft-hand side, but also on the right (it’s hiding in the second f on the right side).In order to find the value yn+1 it appears that we need to carry out an iterative process.First we would guess yn+1 (guessing yn+1 to be equal to yn wouldn’t be all that bad, butwe can do better). If we use this guess value on the right side of (2.5.11) then we wouldbe able to calculate the entire right-hand side, and then we could use that value as a new“improved” value of yn+1 .Now if the new value agrees with the old sufficiently well the iteration would halt, andwe would have found the desired value of yn+1 . Otherwise we would use the improved valueon the right side just as we previously used the first guess.
Then we would have a “moreimproved” guess, etc.Fortunately, in actual use, it turns out that one does not actually have to iterate toconvergence. If a good enough guess is available for the unknown value, then just onerefinement by a single application of the trapezoidal formula is sufficient. This is not thecase if a high quality guess is unavailable. We will discuss this point in more detail insection 2.9. The pair of formulas, one of which supplies a very good guess to the nextvalue of y, and the other of which refines it to a better guess, is called a predictor-correctorpair, and such pairs form the basis of many of the highly accurate schemes that are used inpractice.As a numerical example, take the differential equationy 0 = 2xey + 1(2.5.12)with the initial value y(0) = 1.
If we use h = 0.05, then our first task is to calculate y1 , theapproximate value of y(0.05). The trapezoidal rule asserts thaty1 = 1 + 0.025(2 + 0.1ey1 )(2.5.13)and sure enough, the unknown number y1 appears on both sides.Let’s guess y1 = 1. Since this is not a very inspired guess, we will have to iterate thetrapezoidal rule to convergence. Hence, we use this guess on the right side of (2.5.13),compute the right side, and obtain y1 = 1.056796. If we use this new guess the same way,the result is y1 = 1.057193. Then we get 1.057196, and since this is in sufficiently closeagreement with the previous result, we declare that the iteration has converged.
Then wetake y1 = 1.057196 for the computed value of the unknown function at x = 0.05, and we gonext to x = 0.1 to repeat the same sort of thing to get y2 , the computed approximation toy(0.1).In Table 3 we show the results of using the trapezoidal rule (where we have iterateduntil two successive guesses are within 10−4 ) on our test equation y 0 = 0.5y, y(0) = 1 asthe column Trap(x). For comparison, we show Midpoint(x) and Exact(x).2.6 Comparison of the methods43x0.000.050.100.150.200.25...Trap(x)1.000001.025321.051271.077891.105181.13316...Midpoint(x)1.000001.025001.051251.077561.105131.13282...Exact(x)1.000001.025321.051271.077881.105171.13315...1.002.003.00...1.648762.718424.48203...1.648472.717634.48032...1.648722.718284.48169...5.0010.0012.18402148.4508912.17743148.3127412.18249148.41316table 32.6Comparison of the methodsWe are now in possession of three methods for the numerical solution of differential equations.
They are Euler’s methodyn+1 = yn + hyn0 ,(2.6.1)the trapezoidal ruleyn+1 = yn +h 00),(y + yn+12 n(2.6.2)and the midpoint ruleyn+1 = yn−1 + 2hyn0 .(2.6.3)In order to compare the performance of the three techniques it will be helpful to havea standard differential equation on which to test them. The most natural candidate forsuch an equation is y 0 = Ay, where A is constant. The reasons for this choice are firstthat the equation is easy to solve exactly, second that the difference approximations arealso relatively easy to solve exactly, so comparison is readily done, third that by varying thesign of A we can study behavior of either growing or shrinking (stable or unstable) solutions,and finally that many problems in nature have solutions that are indeed exponential, atleast over the short term, so this is an important class of differential equations.We will, however, write the test equation in a slightly different form for expositoryreasons, namely asyy0 = − ;y(0) = 1 ;L>0(2.6.4)Lwhere L is a constant.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
