Heath - Scientific Computing (523150), страница 69
Текст из файла (страница 69)
To illustrate a change of interval, we use a two-pointGaussian quadrature rule derived for the interval [−1, 1] in Example 8.4 to approximate theintegralZ1I(g) =2e−t dt0from Example 8.1. Using the linear transformation of variables just given, we gett=x+1,2so that the integral is approximated by!2 √1(−1/3)+1 + exp −I(g) ≈exp −22!2 √(1/ 3 ) + 1 ≈ 0.746595,2which is slightly more accurate than the result given by Simpson’s rule for this integral (seeExample 8.1) despite using only two points instead of three.8.3.3Gauss-Kronrod Quadrature RulesAs we have seen, one convenient way to obtain an error estimate is by using two differentquadrature rules.
Since Newton-Cotes quadrature rules use equally spaced nodes, rulesof different orders often have nodes in common. For example, the three nodes used inSimpson’s rule are the same as those used in the midpoint and trapezoid rules. We can takeadvantage of this fact to minimize the number of times that the integrand function mustbe evaluated in using multiple rules of different orders to estimate the error.We have seen that Gaussian quadrature rules are more accurate than Newton-Cotesrules for the same number of nodes, but unfortunately, Gaussian rules of different ordersdo not have any nodes in common (except that Gaussian rules of odd order always havethe midpoint as one node).
Thus, if we seek to estimate the error by using Gaussian rulesof different orders, we must evaluate the integrand function at the full set of nodes of bothrules.Avoiding this additional work is the motivation for Gauss-Kronrod quadrature rules.Such rules come in pairs: an n-point Gaussian rule Gn and a (2n + 1)-point Kronrod ruleK2n+1 whose nodes are optimally chosen subject to the constraint that all of the nodes ofGn are reused in K2n+1 . The (2n + 1)-point Kronrod rule is of polynomial degree 3n + 1,whereas a true (2n + 1)-point Gaussian rule would be of polynomial degree 4n + 1.8.4. COMPOSITE AND ADAPTIVE QUADRATURE255In using such a Gauss-Kronrod pair, the value of K2n+1 is taken as the approximation tothe integral, and a realistic but conservative estimate for the error, based partly on theoryand partly on experience, is given by(200|Gn − K2n+1 |)1.5 .Because they efficiently provide both high accuracy and a reliable error estimate, GaussKronrod rules are among the most effective methods for numerical quadrature, and theyform the basis for many of the quadrature routines available in major software libraries.The pair of rules (G7 , K15 ), in particular, has become a commonly used standard.8.48.4.1Composite and Adaptive QuadratureComposite Quadrature RulesIt is not feasible to use arbitrarily high-order quadrature rules in an attempt to attainarbitrarily high accuracy in evaluating an integral over a given interval.
A much betteralternative is to subdivide the original interval into subintervals, often called panels in thiscontext, then apply a lower-order quadrature rule in each panel. Summing all of thesepartial results then yields an approximation to the overall integral.This approach is equivalent to using piecewise interpolation to derive a composite, orcompound, quadrature rule over the given interval.
For example, if the interval [a, b] ispartitioned into n panels, [xi−1 , xi ], i = 1, . . . , n, with a = x0 < x1 < x2 < · · · < xn−1 <xn = b, then the composite midpoint rule is given byI(f ) ≈ M (f ) =nX(xi − xi−1 )fi=1xi−1 + xi2,and the composite trapezoid rule byI(f ) ≈ T (f ) =nXi=1(xi − xi−1 )f (xi−1 ) + f (xi ).2Composite quadrature rules offer a particularly simple means of estimating error by usingtwo rules of different order.
For example, we observed in Section 8.2.3 that halving theinterval length reduces the error in the midpoint or trapezoid rules by a factor of about 18 .For a given interval [a, b], however, halving the width of each panel doubles the number ofpanels, so the overall reduction in the error is by a factor of about 41 .If the number of panels is n, and hence the average panel width is h = (b − a)/n,then the dominant term in the remainder for the composite midpoint or trapezoid rules isO(nh3 ) = O(h2 ), so the accuracy of these rules is said to be of second order. Similarly, thecomposite Simpson’s rule is of fourth-order accuracy, meaning that the dominant term inits remainder is O(h4 ), and hence halving the panel width reduces the error by a factor of1about 16.256CHAPTER 8.
NUMERICAL INTEGRATION AND DIFFERENTIATION8.4.2Automatic and Adaptive QuadratureA composite quadrature rule with an error estimate can be used to produce an automaticquadrature procedure: simply continue to subdivide all of the panels, say, by half, until theoverall error estimate falls below the required tolerance. This approach usually works, butit may require substantially more work than methods tailored for the particular problem.A more intelligent approach is adaptive quadrature, in which the domain of integration isselectively refined to reflect the behavior of the particular integrand function.For example, one might apply a quadrature rule over the entire original interval.
If theerror tolerance is not met, subdivide the interval into two halves and apply the quadraturerule in each. From this point on, if the sum of the error estimates for the individual panelsstill exceeds the required tolerance, then the panel with the largest error is further halved,and so on until the error tolerance is eventually met, if possible.
In this way the integrandfunction tends to be sampled most densely in regions where it is most active, as shown byexample in Fig. 8.3. Such an adaptive strategy forms the basis for most library subroutinesfor one-dimensional integration.................... ................................. ........................................... .
.|......|.....|......|.....|......|......|.....|......|......|.....|......|......|.....|......|......|.....|...|...|...|..|...|...|...|...|...|..|...|...|...|...|..|...|...|...|...|..|...|...|...|...|...|..|...|...|...|...|...|..|...|...|...|...|...|..|...|...|...|...|..|...|...|...|...|...|..|...|...|...|...|...|..|...|...|...|...|..|...|...|...|...|.....|......|......|.....|......|......|.....|......|.....|......|......|.....|......|......|.....|......|.Figure 8.3: Typical placement of evaluation points by an adaptive quadrature routine.It may not be possible, however, to meet a given error tolerance in computing a givenintegral. The accuracy attainable is limited both by the precision of the arithmetic usedand by the accuracy with which the integrand function can be evaluated.
If the integrandis noisy, or if the error tolerance is unrealistically tight relative to the machine precision,then an adaptive quadrature routine may be unable to meet the error tolerance and willlikely expend a large number of function evaluations only to return a warning messagethat its subdivision limit was exceeded. Such a result should not be regarded as a faultof the adaptive routine but as a reflection of the difficulty of the problem or unrealisticexpectations on the part of the user, or both.Although adaptive quadrature procedures tend to be very effective in practice, they canbe fooled: both the approximate integral and the error estimate can be completely wrong.The reason is that the integrand function is sampled at only a finite number of points,so it is possible that significant features of the integrand may be missed.
For example,it may happen that the interval of integration is very wide, but all of the “interesting”behavior of the integrand is confined to a very narrow range. In this case, sampling bythe automatic routine may completely miss the interesting part of the integrand’s behavior,and the resulting value for the integral may be completely wrong. This situation may seemunlikely, but it can happen, for example, if we are trying to evaluate an integral over aninfinite interval and have truncated it unwisely (see Section 8.5.2).Another potential difficulty with adaptive quadrature routines is that they may be veryinefficient in handling discontinuities (finite jumps in the integrand) and integrable singular-8.5. OTHER INTEGRATION PROBLEMS257ities (points where the integrand becomes infinite but the integral still exists).
For example,an adaptive routine may expend a great many function evaluations in refining the regionaround a discontinuity of the integrand because it assumes that the integrand is smooth(but very steep). A good way to prevent this behavior is to call the quadrature routineseparately to compute the integral on either side of the discontinuity, thereby obviatingthe need for the routine to resolve the discontinuity.
A good strategy for dealing with asingularity is to obtain an analytic formula for the integral in a neighborhood around thesingularity and use the adaptive routine to compute the integral elsewhere.8.58.5.1Other Integration ProblemsIntegrating Tabular DataThus far we have assumed that the integrand function can be evaluated at any desiredpoint within the interval of integration. This assumption may not be valid if the integrand isdefined only by a table of its values at selected points.
A reasonable approach to integratingsuch tabular data is by piecewise interpolation. For example, integrating the piecewise linearinterpolant to tabular data gives a composite trapezoid rule.An excellent method for integrating tabular data is provided by Hermite cubic or cubicspline interpolation. In effect, the overall integral is computed by integrating each of thecubic pieces that make up the interpolant.
This facility is provided by some of the splineinterpolation packages mentioned in Section 7.4.8.5.2Infinite IntervalsAlthough some quadrature routines are capable of handling integrals over infinite or semiinfinite intervals, one may also be able to deal adequately with such problems using standardquadrature routines for finite intervals. A number of approaches are possible:• Replace the infinite limits of integration by finite values. Such finite limits should bechosen carefully so that any omitted tail is negligible or its contribution is estimated,if possible. But the remaining finite interval should not be so wide that an automaticquadrature routine will be fooled into sampling badly.• Transform the variable of integration so that the new interval is finite. Typical transformations include x = − log t or x = t/(1 − t). Care must be taken not to introducesingularities or other difficulties by such a transformation.• Apply a quadrature rule, such as Gauss-Laguerre or Gauss-Hermite, designed for aninfinite interval.8.5.3Double IntegralsThus far we have considered only one-dimensional integrals, where we wish to determinethe area under a curve over an interval.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
