Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 36
Текст из файла (страница 36)
As that number increases, we can findrules that are exact for polynomials of increasingly high order. (Keep in mind thathigher order does not always imply higher accuracy in real cases.) A sequence ofsuch closed formulas is now given.Closed Newton-Cotes FormulasTrapezoidal rule:& x2x111f(x)dx = h f1 + f222+ O(h3 f )(4.1.3)Here the error term O( ) signifies that the true answer differs from the estimate byan amount that is the product of some numerical coefficient times h3 times the value132Chapter 4.Integration of Functionsof the function’s second derivative somewhere in the interval of integration. Thecoefficient is knowable, and it can be found in all the standard references on thissubject.
The point at which the second derivative is to be evaluated is, however,unknowable. If we knew it, we could evaluate the function there and have a higherorder method! Since the product of a knowable and an unknowable is unknowable,we will streamline our formulas and write only O( ), instead of the coefficient.Equation (4.1.3) is a two-point formula (x1 and x2 ). It is exact for polynomialsup to and including degree 1, i.e., f(x) = x. One anticipates that there is athree-point formula exact up to polynomials of degree 2. This is true; moreover, by acancellation of coefficients due to left-right symmetry of the formula, the three-pointformula is exact for polynomials up to and including degree 3, i.e., f(x) = x3 :Simpson’s rule:&x3x1411f(x)dx = h f1 + f2 + f3333+ O(h5 f (4) )(4.1.4)Here f (4) means the fourth derivative of the function f evaluated at an unknownplace in the interval.
Note also that the formula gives the integral over an intervalof size 2h, so the coefficients add up to 2.There is no lucky cancellation in the four-point formula, so it is also exact forpolynomials up to and including degree 3.Simpson’s&x4x138rule:9933f(x)dx = h f1 + f2 + f3 + f48888+ O(h5 f (4) )(4.1.5)The five-point formula again benefits from a cancellation:Bode’s rule:&x5f(x)dx = hx16424641414f1 + f2 + f3 + f4 + f54545454545+ O(h7 f (6) ) (4.1.6)This is exact for polynomials up to and including degree 5.At this point the formulas stop being named after famous personages, so wewill not go any further. Consult [1] for additional formulas in the sequence.Extrapolative Formulas for a Single IntervalWe are going to depart from historical practice for a moment.
Many textswould give, at this point, a sequence of “Newton-Cotes Formulas of Open Type.”Here is an example:&x5f(x)dx = hx0555555f1 + f2 + f3 + f424242424+ O(h5 f (4) )4.1 Classical Formulas for Equally Spaced Abscissas133Notice 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.&&x1x0f(x)dx = h[f1 ]x1f(x)dx = h&x0x1+ O(h2 f )13f1 − f222(4.1.7)+ O(h3 f )1623f1 − f2 +f(x)dx = h1212x0& x15955f1 − f2 +f(x)dx = h2424x05f312(4.1.8)+ O(h4 f (3) )379f3 − f42424(4.1.9)+ O(h5 f (4) )(4.1.10)Perhaps a word here would be in order about how formulas like the above canbe derived. 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:& xN1f(x)dx = h f1 + f2 + f3 +2x1(4.1.11)1(b − a)3 f +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 decreased134Chapter 4.Integration of Functionsby 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:&xNf(x)dx = hx1135f1 + f2 + f3 + f4 +1212· · · + fN−2 +135fN−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:&xNf(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 is&xNf(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.4.1 Classical Formulas for Equally Spaced Abscissas135Extended Formulas (Open and Semi-open)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) give& xN313+O(4.1.15)f(x)dx = h f2 + f3 + f4 + · · ·+ fN−2 + fN−122N2x1Equations (4.1.8) and (4.1.12) give& xN723f2 + f3 + f4 + f5 +f(x)dx = h1212x1723· · · + fN−3 + fN−2 + fN−112121+ON3Equations (4.1.9) and (4.1.13) give& 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), giving& 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):& xNf(x)dx = h[f3/2 + f5/2 + f7/2 +x1(4.1.19)1· · · + fN−3/2 + fN−1/2 ] + ON2136Chapter 4.Integration of FunctionsN=1234(total after N = 4)Figure 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.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
