Conte, de Boor - Elementary Numerical Analysis. An Algorithmic Approach (523140), страница 63
Текст из файла (страница 63)
In this case the solution will be output atthe M equally spaced points XBEGIN + k ∆ X where ∆X = (XMXBEGIN)/ M. Internally, DVERK will automatically select the properstep sizes so as to achieve the required accuracy. The step size normallywill vary as the integration proceeds. The subroutine also keeps track ofthe number of function evaluations required to find the solution at XEND.DVERK is a high-order-accuracy routine which requires a minimum ofeight function evaluations per integration step.
The number of functionevaluations actually used is stored in C(24) and can on option be outputtedas we have done above.As applied to the differential equationy (1) = -1which we considered in Example 8.2; the complete program and the resultsare given below.CUSE OF DVERK TO SOLVE EXAMPLE 8.2 .INTEGER IER,IND,K,N,NWREAL C(24),TOL,W(1,9),X,XEND,Y(l)DATA N , X ,Y(l), TOL ,IND,NW*/ 1 , 1.,-l. ,l.E-7, 1 , 1 /EXTERNAL FCNlDO 10 K=1,4XEND = 1. + FLOAT(K)/4.CALL DVERK ( N, FCNl, X, Y, XEND, TOL, INC, C, NW, W, IER )PRINT 600, XEND,Y(l),C(24)FORMAT(llX,F8.6,5X,El6.8,5X,Fll.0)60010 CONTINUE372THE SOLUTION OF DIFFERENTIAL EQUATIONSSTOPENDSUBROUTINE FCNl ( N, X, Y, YPRIME )REAL X,Y(l),YPRIME(l)YPRIME(1) = (1./X - Y(l))/X - Y(l)*Y(l)RETURNENDOUTPUTXY(l)FCN EVALS1.251.501.752.00-0.79999999-0.66666664-0.57142854-0.4999999616.24.32.40.The results are comparable in accuracy with those obtained using theclassical fourth-order method with a fixed step size of h = 1/32.
Since theclassical fourth-order method requires four function evaluations per step, atotal of 128 function evaluations was required to achieve about seven-decimal-place accuracy. By contrast, DVERK requires only 40 functionevaluations for the same accuracy. Note that 16 function evaluations wererequired for the output at x = 1.25, indicating that the step h = 1/4 was toolarge and apparently had to be halved to achieve 1.10-7 accuracy.In Sec. 8.12 we will illustrate the use of DVERK to solve a system offirst-order differential equations.EXERCISES8.6-l Suppose we are using a Runge-Kutta method of order 2 and step size control based oninterval halving to solve a differential equation. If we are using a step h = 0.1 and the error-4criterion ε = 10 -6 , and we find that Dn = 10 at a point x = xn what step should be usedfor the next integration step?8.6-2 Write a program for the Runge-Kutta method of order 2 with step size controlrestricted to doubling or halving.
Apply this program to solve the equation of Example 8.2with ε = 10 -6 .8.6-3 Check with your computing center to see whether they carry the IMSL collection ofsubroutines. Use subroutine DVERK to solve the following differential equations. In eachcase set TOL = 10-7 and request output at the XEND valuesXEND = XO + K(XM-XO)/10(a) y´ = x - l + y/xXO = 1, XM = 2, y(X0) = 2(b) y´ = xy2XO = l, XM = 4, y(XO) = 1K = 1, 2, . . .,108.7MULTISTEP FORMULAS3738.7 MULTISTEP FORMULASThe Taylor algorithm of order k and the Runge-Kutta methods are bothexamples of one-step methods.
They require information about the solution at a single point x = xn, from which the methods proceed to obtain yat the next point x = xn+1 . Multistep methods make use of informationabout the solution at more than one point. Let us assume that we havealready obtained approximations to y´ and y at a number of equally spacedpoints, say x0, x1, . . . , xn. One class of multistep methods is based on theprinciple of numerical integration. If we integrate the differential equationy´ = f(x,y) from xn to xn+1 , we will haveor(8.44)To carry out the integration in (8.44) we now approximate f(x,y(x)) by apolynomial which interpolates f(x,y(x)) at the (m + 1) pointsxn, xn-1, xn-2, .
. . , xn-m. If we use the notationf(xk,y(xk)) = fkwe can use the Newton backward formula (see Exercise 2.6-8) of degree mfor this purpose:Inserting this into (8.44) and noting that dx = h ds, we obtain(8.45)where(8.45a)From the definition of the binomial function given in Chap. 2 we caneasily compute the γk, the first few of which areγ0 = 1γ 1 = 1/2γ 2 = 5/12γ 3 = 3/8γ 4 = 251/720Formula (8.45) is known as the Adams-Bashford method. The simplestcase, obtained by setting m = 0 in (8.45), again leads to Euler’s method. Ingeneral, the use of (8.45) requires the value of y´ = f at the m + 1 pointsxn, xn-1, . .
. , xn-m. From these we can form the differencesfrom (8.45) we can compute yn+1 ; from thedifferential equation we can compute fn+1 = f(xn+1, yn+1). We now relabel374 THE SOLUTION OF DIFFERENTIAL EQUATIONSthe point xn+1 as xn, form a new line of differences, and repeat the process.For m = 3, which is commonly used in practice, the difference table isxn-3yn-3 fn-3∆f n - 3xn-2yn-2∆2fn-2fn-2∆ 3f n - 3∆f n - 2xn-1yn-1∆2 f n-2fn-1∆f n - 1xnynfnand (8.45) specializes to(8.46)In practice, it is more convenient computationally to work withordinates instead of differences. From the definition of the forwarddifference operator ∆ we find thatSubstituting in (8.46) and regrouping, we obtain(8.47)The local error of (8.46) may be derived as follows: From Exercise2.6-8 we know that the error of Newton’s backward formula with n = 3and k = 0 isThe error of (8.46) is then given bySincedoes not change sign on the interval [0, 1], there exists apoint ξ between xn-3 and xn+1 , such that(8.48)8.7MULTISTEP FORMULAS375To use (8.47) we must have four starting values.
These starting valuesmust be obtained from some independent source. To illustrate how (8.47)is used, we carry out a few steps of the integration of the equationy´ = -y2y(1) = 1with h = 0.1. The exact solution of this problem is y = 1/ x. In the tablebelow, the first four starting values are obtained from the exact solution,and the remaining entries by (8.47).xnynf n = -y n 2y(x n ) = 1/x n1.01.11.21.31.41.51.61.000000000.909090910.833333330.769230770.714436320.666860300.62524613-1.00000000-0.82644628-0.69444444-0.59171598-0.51041926-0.44470266-0.390932720.909090910.833333330.769230770.714285710.666666670.62500000The values yn computed by formula (8.47) are seen to be in error byabout two units in the fourth decimal place. Using the local error estimate(6.43) and the fact thatl < x < 2we obtain the error boundThis bound is about twice as large as the errors encountered in going fromone step to the next.A number of other formulas of the multistep type can be derivedsimilarly, using numerical integration.
Instead of integrating f(x,y) in(8.43) from xn to xn+1, we could, for example, integrate from xn-p to xn+1for some integer p > 0. If we again interpolate at the m + 1 pointsxn, xn-1, · · · , xn-m with Newton’s backward formula, we obtain(8.49)The case p = 0 yields the Adams-Bashforth formula (8.44). Some especially interesting formulas of this type are those corresponding to m = 1,p = 1 and to m = 3, p = 3. These formulas together with their local-error376THE SOLUTION OF DIFFERENTIAL EQUATIONSterms are(8.50)(8.51)Formula (8.50), which is comparable in simplicity to Euler’s method, has amore favorable discretization error.
Similarly (8.51), which requires knowledge of f(x,y) at only three points, has a discretization error comparablewith that of the Adams-Bashforth method (8.47). It can be shown that allformulas of the type (8.49) with m odd and m = p have the property thatthe coefficient of the mth difference vanishes, thus yielding a formula ofhigher order than might be expected. On the other hand, these formulasare subject to greater instability, a concept which will be developed later.A major disadvantage of multistep formulas is that they are notself-starting.
Thus, in the Adams-Bashforth method (8.47), we must havefour successive values of f(x,y) at equally spaced points before thisformula can be used. These starting values must be obtained by someindependent method. We might, for example, use Taylor’s algorithm orone of the Runge-Kutta methods to obtain these starting values. We mustalso be assured that these starting values are as accurate as necessary forthe overall required accuracy. A second disadvantage of the Adams-Bashforth method is that, although the local discretization error isthecoefficient in the error term is somewhat larger than for formulas of theRunge-Kutta type of the same order. Runge-Kutta methods are generally,although not always, more accurate for this reason.
On the other hand, themultistep formulas require only one derivative evaluation per step, compared with four evaluations per step with Runge-Kutta methods, and aretherefore considerably faster and require less computational work.Example 8.5 Solve the equationy´ = x + yy(0) = 0from x = 0 to x = 1, using the Adams-Bashforth method.A FORTRAN program and the results for this problem are given below. The exactsolution of this problem is y = ex - 1 - x. The first four starting values are computed,using this solution. The first column of the results gives the values of xn with h = 1/32, thesecond column gives y n as computed by formula (8.47), the third column gives thevalue y(xn) as computed from the solution, and the fourth column gives the error en =yn - y(xn)The results are correct to about six significant figures, which is approximately whatwould be expected from the error formula (8.48).
Since the accumulated discretizationerror iswe would expect to reduce the error by 1/16 if the step size h were halved.8.7MULTISTEP FORMULASFORTRAN PROGRAM FOR EXAMPLE 8.5C ADAMS-BASHFORTH METHODINTEGER I,N,NSTEPSREALERROR,F(4),H,XBEGIN,XN,YBEGIN,YNCSOLN(X) = EXP(X) - 1. - XCC ** INITIALIZEPRINT 600608 FORMAT('lADAMS-BASHFORTH METHOD'/* '0',4X,'N',l3X,'XN',l5X,'YN',l3X,'Y(XN)',12X,'ERROR'/)NSTEPS = 32H = 1./NSTEPSYBEGIN = 0.XBEGIN = 0.CC ** COMPUTE FIRST FOUR POINTS USING EXACT SOLUTIONF(1) = XBEGIN + YBEGINN=0ERROR = 0.PRINT 601, N,XBEGIN,YBEGIN,YBEGIN,ERROR601 FORMAT(' ',I3,4X,4El7.8)DO 20 N=l,3XN = XBEGIN + N*HYN = SOLN(XN)F(N+l) = XN + YNPRINT 601, N,XN,YN,YN,ERROR20 CONTINUECC ** BEGIN ITERATIONDO 50 N=4,NSTEPSYN = YN + (H/24.)*(55.'F(4)-59.*F(3)+37.*F(2)-9.*F(l))XN = XBEGIN + N*HF(1) = F(2)F(2) = F(3)F(3) = F(4)F(4) = XN + YNYOFXN = SOLN(XN)ERROR = YN - YOFXNPRINT 601, N,XN,YN,YOFXN,ERROR50 CONTINUESTOPENDCOMPUTER RESULTS FOR EXAMPLE 8.5NXN0123456789101112130.0.31250000E-010.6250000E-010.93750000E-010.12500000E-000.1562500E-000.18750000E-000.21875000E-000.25000000E-000.28125000E-000.31250000E-000.34375000E-000.37500000E-000.40625000E-00YN0.0.49340725E-030.19944459E-020.45351386E-020.81484411E-020.12868421E-010.18730211E-010.25770056E-010.34025350E-010.43534677E-010.54337843E-010.66475919E-010.79991280E-010.94927646E-01Y(XN)0.0.49340725E-030.19944459E-020.45351386E-020.81484467E-020.12868434E-010.18730238E-010.25770098E-010.34025416E-010.43534756E-010.54337934E-010.66476032E-010.79991400E-01094927788E-01ERROR0.0.0.0.-0.55879354E-08-0.12922101E-07-0.26309863E-07-0.41676685E-07-0.65192580E-07-0.78696758E-07-0.90803951E-07-0.11269003E-06-0.12014061E-06-0.14156103E-06377378THE SOLUTION OF DIFFERENTIAL EQUATIONSCOMPUTER RESULTS FOR EXAMPLE 8.5 (continued)NXN141516171819202122232425262728293031320.43750000E-00046875000E-000.50000000E 000.53125000E 000.56250000E 000.59375000E 000.62500000E 000.65625000E 000.68750000E 000.71875000E 000.75000000E 000.78125000E 000.81250000E 000.84375000E 000.87500000E 00090625000E 000.93750000E 00096875000E 000.09999999E 01YN0.11133012E-000.12924525E-000.14872105E-000.16980705E-000.19255438E-000.21701577E-000.24324562E-000.27130008E-000.30123707E-000.33311634E-000.36699954E-00040295030E-000.44103424E-000.48131907E-000.52387466E 000.56877308E 000.61608872E 00066589829E 000.71828098E 00Y(XN)ERROR0.11133029E-000.12924545E-000.14872126E-000.16980730E-000.19255446E-000.21701607E-000.24324594E-000.27130044E-000.30123746E-000.33311677E-000.36700001E-000.40295079E-000.44103476E-00048131964E-000.52387527E 000.56877375E 000.61608934E 00066589907E 000.71828181E 00-0.16111881E-06-0.19185245E-06-0.21234155E-06-0.24400651E-06-0.26822090E-06-0.29988587E-06-0.31664968E-06-0.34645200E-06-0.39115548E-06-0.42840838E-06-0.46566129E-06-0.49173832E-06-0.52526593E-06-0.56624413E-06-0.61094761E-06-066310167E-06-0.71525574E-06-0.77486038E-06-0.82701445E-06EXERCISES8.7-l Using (8.45a), derive the coefficients γ k (k = 1, .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
