c4-1 (779478), страница 2
Текст из файла (страница 2)
There are elegant ways, but the most straightforward is to write down thebasic form of the formula, replacing the numerical coefficients with unknowns, sayp, q, r, s. Without loss of generality take x0 = 0 and x1 = 1, so h = 1. Substitute inturn for f(x) (and for f1 , f2 , f3 , f4 ) the functions f(x) = 1, f(x) = x, f(x) = x2 ,and f(x) = x3 . Doing the integral in each case reduces the left-hand side to anumber, and the right-hand side to a linear equation for the unknowns p, q, r, s.Solving the four equations produced in this way gives the coefficients.Extended Formulas (Closed)If we use equation (4.1.3) N − 1 times, to do the integration in the intervals(x1 , x2 ), (x2 , x3 ), . . . , (xN−1 , xN ), and then add the results, we obtain an “extended”or “composite” formula for the integral from x1 to xN .Extended trapezoidal rule:Z xN1f(x)dx = h f1 + f2 + f3 +2x1(4.1.11)1(b − a)3 f 00+O· · · + fN−1 + fN2N2Here we have written the error estimate in terms of the interval b − a and the numberof points N instead of in terms of h.
This is clearer, since one is usually holdinga and b fixed and wanting to know (e.g.) how much the error will be decreasedSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Notice that the integral from a = x0 to b = x5 is estimated, using only the interiorpoints x1 , x2 , x3 , x4. In our opinion, formulas of this type are not useful for thereasons that (i) they cannot usefully be strung together to get “extended” rules, as weare about to do with the closed formulas, and (ii) for all other possible uses they aredominated by the Gaussian integration formulas which we will introduce in §4.5.Instead of the Newton-Cotes open formulas, let us set out the formulas forestimating the integral in the single interval from x0 to x1 , using values of thefunction f at x1 , x2 , .
. . . These will be useful building blocks for the “extended”open formulas.134Chapter 4.Integration of FunctionsZxNf(x)dx = hx1135f1 + f2 + f3 + f4 +1212135· · · + fN−2 + fN−1 + fN1212+O1N3(4.1.12)(We will see in a moment where this comes from.)If we apply equation (4.1.4) to successive, nonoverlapping pairs of intervals,we get the extended Simpson’s rule:ZxNf(x)dx = hx14241f1 + f2 + f3 + f4 +3333412· · · + fN−2 + fN−1 + fN333+O1N4(4.1.13)Notice that the 2/3, 4/3 alternation continues throughout the interior of the evaluation.
Many people believe that the wobbling alternation somehow contains deepinformation about the integral of their function that is not apparent to mortal eyes.In fact, the alternation is an artifact of using the building block (4.1.4). Anotherextended formula with the same order as Simpson’s rule isZxNf(x)dx = hx17233f1 + f2 + f3 + f4 + f5 +86242373(4.1.14)· · · + fN−4 + fN−3 + fN−2 + fN−1 + fN24681+ON4This equation is constructed by fitting cubic polynomials through successive groupsof four points; we defer details to §18.3, where a similar technique is used in thesolution of integral equations. We can, however, tell you where equation (4.1.12)came from.
It is Simpson’s extended rule, averaged with a modified version ofitself in which the first and last step are done with the trapezoidal rule (4.1.3). Thetrapezoidal step is two orders lower than Simpson’s rule; however, its contributionto the integral goes down as an additional power of N (since it is used only twice,not N times). This makes the resulting formula of degree one less than Simpson.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).by taking twice as many steps (in this case, it is by a factor of 4). In subsequentequations we will show only the scaling of the error term with the number of steps.For reasons that will not become clear until §4.2, equation (4.1.11) is in factthe most important equation in this section, the basis for most practical quadratureschemes.The extended formula of order 1/N 3 is:4.1 Classical Formulas for Equally Spaced Abscissas135Extended Formulas (Open and Semi-open)Equations (4.1.8) and (4.1.12) giveZ xN723f2 + f3 + f4 + f5 +f(x)dx = h1212x1723· · · + fN−3 + fN−2 + fN−112121+ON3Equations (4.1.9) and (4.1.13) giveZ xN13427f2 + 0 + f4 + f5 +f(x)dx = h12123x113274· · · + fN−4 + fN−3 + 0 + fN−1312121+ON4(4.1.16)(4.1.17)The interior points alternate 4/3 and 2/3.
If we want to avoid this alternation,we can combine equations (4.1.9) and (4.1.14), givingZ xN11155f2 − f3 + f4 + f5 + f6 + f7 +f(x)dx = h2468x111155· · · + fN−5 + fN−4 + fN−3 − fN−2 + fN−186241+ON4(4.1.18)We should mention in passing another extended open formula, for use wherethe limits of integration are located halfway between tabulated abscissas.
This one isknown as the extended midpoint rule, and is accurate to the same order as (4.1.15):Z xNf(x)dx = h[f3/2 + f5/2 + f7/2 +x1(4.1.19)1· · · + fN−3/2 + fN−1/2 ] + ON2Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).We can construct open and semi-open extended formulas by adding the closedformulas (4.1.11)–(4.1.14), evaluated for the second and subsequent steps, to theextrapolative open formulas for the first step, (4.1.7)–(4.1.10). As discussedimmediately above, it is consistent to use an end step that is of one order lowerthan the (repeated) interior step.
The resulting formulas for an interval open atboth ends are as follows:Equations (4.1.7) and (4.1.11) giveZ xN313+O(4.1.15)f(x)dx = h f2 + f3 + f4 + · · ·+ fN−2 + fN−122N2x1136Chapter 4.Integration of FunctionsN=1234Figure 4.2.1. Sequential calls to the routine trapzd incorporate the information from previous calls andevaluate the integrand only at those new points necessary to refine the grid. The bottom line shows thetotality of function evaluations after the fourth call. The routine qsimp, by weighting the intermediateresults, transforms the trapezoid rule into Simpson’s rule with essentially no additional overhead.There are also formulas of higher order for this situation, but we will refrain fromgiving them.The semi-open formulas are just the obvious combinations of equations (4.1.11)–(4.1.14) with (4.1.15)–(4.1.18), respectively. At the closed end of the integration,use the weights from the former equations; at the open end use the weights fromthe latter equations.
One example should give the idea, the formula with error termdecreasing as 1/N 3 which is closed on the right and open on the left:Z xN723f2 + f3 + f4 + f5 +f(x)dx = h1212x1 (4.1.20)1351+O· · · + fN−2 + fN−1 + fN1212N3CITED REFERENCES AND FURTHER READING:Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), §25.4. [1]Isaacson, E., and Keller, H.B. 1966, Analysis of Numerical Methods (New York: Wiley), §7.1.4.2 Elementary AlgorithmsOur starting point is equation (4.1.11), the extended trapezoidal rule.
There aretwo facts about the trapezoidal rule which make it the starting point for a variety ofalgorithms. One fact is rather obvious, while the second is rather “deep.”The obvious fact is that, for a fixed function f(x) to be integrated between fixedlimits a and b, one can double the number of intervals in the extended trapezoidalrule without losing the benefit of previous work. The coarsest implementation ofthe trapezoidal rule is to average the function at its endpoints a and b. The firststage of refinement is to add to this average the value of the function at the halfwaypoint. The second stage of refinement is to add the values at the 1/4 and 3/4 points.And so on (see Figure 4.2.1).Without further ado we can write a routine with this kind of logic to it:Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).(total after N = 4).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
