Heath - Scientific Computing (523150), страница 67
Текст из файла (страница 67)
As an example, we approximate the integralI(f ) =Z12e−x dx0using each of the first three Newton-Cotes quadrature rules:M (f ) = (1 − 0) exp(−0.25) ≈ 0.7788011[exp(0) + exp(−1)] ≈ 0.683940T (f ) =21S(f ) =[exp(0) + 4 exp(−0.25) + exp(−1)] ≈ 0.7471806The integrand and the interpolating polynomial for each rule are shown in Fig. 8.1. Thecorrectly rounded result for this problem is 0.746824. It is somewhat surprising to see thatthe magnitude of the error from the trapezoid rule (0.062884) is about twice that from themidpoint rule (0.031977), and that Simpson’s rule, with an error of only 0.000356, seemsremarkably accurate considering the size of the interval over which it is applied. We willsoon see explanations for these phenomena.8.2.2Method of Undetermined CoefficientsAs we have seen, a quadrature rule can be derived directly by interpolating the integrandfunction by a polynomial at a set of points and then integrating the interpolant.
Analternative derivation that yields some additional insight is the method of undeterminedcoefficients. In deriving a quadrature rule of Newton-Cotes type on an interval [a, b], wetake the nodes x1 , . . . , xn as given and consider the weights w1 , . . . , wn as coefficients tobe determined. If we force the quadrature rule to integrate each of the polynomial basisfunctions exactly, then, by linearity, it will integrate any polynomial of degree n − 1 exactly.248CHAPTER 8. NUMERICAL INTEGRATION AND DIFFERENTIATION1.0 .........................................................................................................f.................................
.. ............ ........... ...... . ... ... ......................................... .......................... ...... .............. .................. ........ ........ ........ ........ ........ ............... ............. ........ ..................................... ........ ........ ........ ........ ........ ........ ........ ........ .................... ....... ........................ ....................... ............ ...... ........................................... ... .............................. ...... ..........................
... ............................... ...... ............................ ............................. ............................................... . ...................... .........................................................................................................................................................................................................................................................................................................................................................................midpointSimpsontrapezoid|0.50.01.02Figure 8.1: Integration of f (x) = e−x by Newton-Cotes quadrature rules.In this manner we obtain a system of n linear equations in n unknowns that determines theappropriate set of weights for the quadrature rule.Example 8.2 Method of Undetermined Coefficients.
We illustrate the method ofundetermined coefficients by deriving a three-point ruleI(f ) ≈ w1 f (x1 ) + w2 f (x2 ) + w3 f (x3 )on the interval [a, b] using the monomial basis. The three equally spaced points are x1 = a,x2 = (a + b)/2, and x3 = b, and the three monomial basis functions are 1, x, and x2 . Theresulting system of equations isw1 · 1 + w2 · 1 + w3 · 1 =Zb1 dx = x|ba = b − a,aw1 · a + w2 · (a + b)/2 + w3 · b =Zbx dx = (x2 /2)|ba = (b2 − a2 )/2,aw1 · a2 + w2 · ((a + b)/2)2 + w3 · b2 =Zbx2 dx = (x3 /3)|ba = (b3 − a3 )/3.aWritten in matrix form, we recognize this system of equations1aa21(a + b)/2((a + b)/2)2 1w1b−ab w2 = (b2 − a2 )/2 b2w3(b3 − a3 )/3as a Vandermonde system (recall Section 3.2).
Solving it by Gaussian elimination, we obtainthe weightsw1 = (b − a)/6,w2 = 2(b − a)/3,which we recognize as Simpson’s rule.w3 = (b − a)/6,8.2. NEWTON-COTES QUADRATURE8.2.3249Error EstimationThe error in the midpoint quadrature rule can be estimated by means of a Taylor seriesexpansion about the midpoint m = (a + b)/2 of the interval [a, b]:f (x) = f (m) + f 0 (m)(x − m) +f 00 (m)f 000 (m)f iv (m)(x − m)2 +(x − m)3 +(x − m)4 + · · · .2624When we integrate this expression from a to b, the odd-order terms drop out, yieldingf 00 (m)f iv (m)(b − a)3 +(b − a)5 + · · ·241920= M (f ) + E + F + · · · ,I(f ) = f (m)(b − a) +where we have used E and F to represent the first two terms in the error expansion for themidpoint rule.To derive a comparable error expansion for the trapezoid quadrature rule, we substitutex = a and x = b into the Taylor series, add the two resulting series together, observe onceagain that the odd-order terms drop out, solve for f (m), and substitute into the midpointformula to obtainI(f ) = T (f ) − 2E − 4F − · · · .Note thatT (f ) − M (f ) = 3E + 5F + · · · ,and hence the difference between the two quadrature rules provides an estimate for thedominant term in their error expansions,E≈T −M,3provided that the length of the interval, h = b − a, is sufficiently small that h5 h3 , andthe integrand f is such that f iv is reasonably well-behaved.
Under these assumptions, wemay draw several conclusions from the previous derivations:• Halving the length h of the interval decreases the error in either rule by a factor of about18.• The midpoint rule is about twice as accurate as the trapezoid rule, despite being basedon polynomial interpolation of degree one less.• The difference between the midpoint rule and the trapezoid rule can be used to estimatethe error in either of them.An appropriately weighted combination of the midpoint and trapezoid rules eliminates theE (i.e., h3 ) term from the error expansion:212M (f ) + T (f ) − F + · · ·3332= S(f ) − F + · · · ,3I(f ) =250CHAPTER 8.
NUMERICAL INTEGRATION AND DIFFERENTIATIONwhich provides an alternative derivation for Simpson’s rule as well as an expression for itserror term.Example 8.3 Error Estimation.R 1 2 We illustrate error estimation by computing theapproximate value for the integral 0 x dx. Using the midpoint rule, we obtain 211M (f ) = (1 − 0)= ,24and using the trapezoid rule we obtainT (f ) =11−0 2(0 + 12 ) = .22Thus, we have the estimateE≈T −M1/41== .33121We conclude that the error in M is about 12, and the error in T is about − 61 .In addition, we can now compute the approximate value for the integral given by Simpson’s rule,212 1 1 11S(f ) = M + T = · + · = ,333 4 3 23which is the exact value for this integral (as is to be expected since, by design, Simpson’srule is exact for quadratic polynomials). Thus, the error estimates for M and T are in factexactly correct for this simple problem (though this would not be true in general).8.2.4Polynomial DegreeThe accuracy of a quadrature rule is conveniently characterized by the notion of polynomialdegree.
A quadrature rule is said to be of polynomial degree d if it is exact (i.e., its remainderis zero) for every polynomial of degree d but is not exact for some polynomial of degreed + 1. Since an n-point Newton-Cotes rule is based on an interpolating polynomial of degreen − 1, we would expect such a rule to have polynomial degree at least n − 1, and we enforcedthis requirement by construction in the method of undetermined coefficients.Thus, we would expect the midpoint rule to have polynomial degree zero, the trapezoidrule degree one, Simpson’s rule degree two, and so on.
We saw from a Taylor series expansion, however, that the error for the midpoint rule depends on the second and higherderivatives of the integrand, which vanish for linear as well as constant polynomials. Thisimplies that the midpoint rule in fact integrates linear polynomials exactly, and hence itspolynomial degree is one rather than zero.
Similarly, the error for Simpson’s rule depends onthe fourth and higher derivatives, which vanish for cubics as well as quadratic polynomials,so that Simpson’s rule is of polynomial degree three rather than two.In general, an odd-order Newton-Cotes rule gains an extra degree beyond that of thepolynomial interpolant on which it is based. This phenomenon is due to cancellation ofpositive and negative errors, as illustrated for the midpoint and Simpson rules in Fig. 8.2,which, on the left, shows a linear polynomial and the constant function interpolating it at the8.3. GAUSSIAN QUADRATURE251midpoint and, on the right, a cubic and the quadratic interpolating it at the midpoint andendpoints. Integration of the linear polynomial by the midpoint rule yields two congruenttriangles of equal area.