Conte, de Boor - Elementary Numerical Analysis. An Algorithmic Approach (523140), страница 53
Текст из файла (страница 53)
Theerror of the corrected trapezoid rule isIf the above-mentioned rules for numerical integration do not give asatisfactory approximation to I(f), we could, of course, increase the degreek of the interpolating polynomial used. We discussed the dangers of suchan action in Sec. 6.7 and proposed there the use of piecewise-polynomialinterpolation as a more reasonable and certain means for achieving highaccuracy. Accordingly, we approximate I(f) by I(g k ), where g k (x) is apiecewise-polynomial function of “low” degree k which interpolates f(x).We discuss the resulting integration rules, usually called composite rules, inSec. 7.4.We have derived in this section five basic integration rules.
These arethe rectangle rule (7.24), the midpoint rule (7.25), the trapezoid rule (7.26),Simpson’s rule (7.28), and the corrected trapezoid rule (7.29). The corrected trapezoid rule is the only one of these requiring knowledge of thederivative of f(x), and this is an obvious disadvantage of this particularmethod. The error terms of these rules suggest that Simpson’s rule or thecorrected trapezoid rule should be preferred whenever the function f(x) issufficiently smooth.
There are, nevertheless, some functions for whichlower-order formulas yield better results than do higher-order formulas[see Exercise 7.2-2].Example 7.1 Apply each of the five rules given above to find estimates forwe set a = 0, b = 1, (a + b)/2 = ½, and from a table of values find thatf (0) = 1f(l) = e-1 = 0.36788f(½) = e-¼ = 0.77880We will also needf´(0) = 0f´(1) = -2e-1=-0.73576310DIFFERENTIATION AND INTEGRATIONWe can then calculate from the appropriate formulasR = l · e0 = 1M = 1 · e-1/4 = 0.77880T - ½[e0 + e-1 ] = 0.68394S = 1/6[e0 + 4e-1/4 + e-1 ] = 0.74718CT = ½[e0 + e-1 ] + 1/12[0 + 2e - 1 ] = 0.74525The value of the integral correct to five decimal places is I = 0.74682.
The correctedtrapezoid (CT) rule and Simpson’s (S) rule clearly give the best results, as might beexpected from a consideration of the error terms and the fact that the first fewderivatives of the functiondo not vary much in size.EXERCISES7.2-l Verify by direct integration that=(x - a) (x - (a + b) / 2 ) (x - b)7.2-2 Apply each of the five rules given in this section to find an approximation toI =Compare the results with the correct value I = sin 1 - cos 1 = 0.301169.7.2-3 The function f(x) is defined on the interval [0, 1] as follows:Calculate the results of applying the following rules to find(a) The trapezoid rule over the interval [0, 1](b) The trapezoid rule first over the interval [0, ½] and then over the interval [½ , 1](c) Simpson’s rule over the interval [0, 1](d) The corrected trapezoid rule over the interval [0, l]Account for the differences in the results.7.2-4 The corrected trapezoid rule can be derived more simply by observing that since p3(x)is a polynomial of degree 3, piv3(x) = 0, and hence that Simpson’s rule (7.28) can be used toevaluate I(p3) exactly.
HenceSince p3(x) interpolates f(x) at a, a, b, b we must have p3(a) = f(a), p3(b) = f(b). Show usingthe results of Sec. 2.7 on osculatory interpolation thatThen substitute into the expression for I(p3 ) above to derive the corrected trapezoid rule(7.296).7.2-5 Use Simpson’s rule to estimate the value of the integralObtain a7.2-6 use the trapezoid rule to estimate the value of the integralbound on the error of the trapezoid rule (7.26) and compare with the actual error.7.3NUMERICAL INTEGRATION: GAUSSIAN RULES3117.3 NUMERICAL INTEGRATION: GAUSSIAN RULESAll the rules derived in Sec.
7.2, except for the corrected trapezoid rule, canbe written in the form(7.30)where the weights A,, . . . , A, do not depend on the particular functiong(x). We have, so far, picked the nodes x0, . . . , xk somehow, for example,equispaced as in a table, and have then calculated the weights Ai as I(li ),all i. This guarantees that the rule is exact for polynomials of degree < k.But it is possible to make such a rule exact for polynomials of degree< 2k - 1, by choosing also the nodes appropriately. This is the basic ideaof gaussian rules.The resulting rules look more complicated than the rules derived inSec. 7.2.
Both nodes and weights for gaussian rules are, in general,irrational numbers. This fact may have deterred people from using theserules when calculations were done by hand. But, on a computer, it usuallymakes no difference whether one evaluates a function at x = 3 or atOnce the nodes and weights of such a ruleare stored in some form (for example, as in the subroutine LGNDREbelow), these rules are as easily used as the trapezoid rule or Simpson’srule. At the same time, these gaussian rules are usually much moreaccurate when compared with the rules of Sec.
7.2 on the basis of numberof function values used.We discuss these gaussian rules in the more general context of anintegralin which the integrand f(x) may not be often enoughdifferentiable to justify application of the rules of Sec. 7.2. For example,f(x) may behave like (x - a) α near a, for some α > -1, or a and/or bmay be infinite. In such situations, it is often possible to rewrite theintegral aswhere w(x) is a nonnegative integrable function, andis smooth. In the above example, this is the case with w(x) = (x - a) α .Other choices for w(x) are discussed below. The situation of a trouble-freeintegrand is also covered in this setup, by the simple choice w(x) = 1.Consider now the approximate evaluation of the weighted integral(7.3 1)by a rule of the form (7.30). We say that the rule (7.30) is exact for the312DIFFERENTIATION AND INTEGRATIONparticular function p(x) if substitution of p(x) for g(x) into (7.30) makes(7.30) an equality.
The trapezoid rulefor instance, is exact for all polynomials of degree < 1. To check this, weonly have to look at the error term for this rule,Since this error term involves the second derivative of g(x), and the secondderivative of any polynomial of degree < 1 is identically zero, it followsthat the error is zero whenever g(x) is a polynomial of degree < 1. Moregenerally, if the error term of (7.30) is of the formE = const g(r+1) (η) · (some function of x0, .
. . , xk )(7.32)then the rule (7.30) must be exact for all polynomials of degree < r.Hence, if we wish to construct a rule of the form (7.30) which, forfixed k, is exact for polynomials of as high a degree as possible, we shouldconstruct the rule in such a way that it has an error term of the form (7.32),with r as large an integer as possible.
This we can do, using a trick alreadyemployed in Sec. 7.2.As in Sec. 7.2, we use analytic substitution, picking points x0, . . . , xkin (a,b) and writingwhere p k(x) is the polynomial of degree < k which interpolates g(x) atx0, . . . , xk, andThis givesThe approximation I(pk) to I(g) is clearly of the form (7.30). For if wewrite pk(x) in Lagrange form (see Sec. 2.2),p k (x) = g(x 0 )l 0 (x) + g(x 1 )l 1 (x) + · · · + g(x k )l k (x)withtheni = 0, . . . ,k7.3NUMERICAL INTEGRATION: GAUSSIAN RULES313HenceI(pk) = A0g(x0) + A1g(x1) + ·where·· + Akg(xk)i = 0, . . . , k(7.33)(7.34)Next, consider the errorSuppose thatThen, as argued in Sec. 7.2,for any choice of xk+1. If now alsothen, by the same reasoning,Hence, in general, if for Certain x0, .
. . , xk+m,i = 0, . . . , m - 1 (7.35)then, for any choice of xk+m+1(7.36)Now recall from Sec. 6.6 that we can find, for many w(x), a polynomial Pk+1(x) such that(7.37)for all polynomials q(x) of degree < k (see Property 3 of orthogonalpolynomials in Sec. 6.6). Further, by Property 2 of orthogonal polynomials,we can writewhere ξ0, .
. . , ξk are the k + 1 distinct points in the interval (a,b) at314DIFFERENTIATION AND INTEGRATIONwhich Pk+1 , vanishes. Hence, if we setxj = ξj(7.38)j = 0, . . . ,kand let xk+j be arbitrary points in (a,b), j = 1, . . . , k + 1, then (7.35),and therefore (7.36), is satisfied for m = k. For then (7.35) is of the form(7.37), withi=0,...,m-lwhich, for m < k, are all polynomials of degree < k.
Therefore(7.39)To get this error into the form (7.32), we pick the xk+j’s asx k+j = ξj - 1j = l, . . . , k + 1Thenso thatis of one sign, i.e., nonnegative, on (a,b). Hence wecan apply the mean-value theorem for integrals (see Sec. 1.7) to getFinally, if g(x) is 2k + 2 times continuously differentiable, we can makeuse of Theorem 2.5 to express the error in the form(7.40)whereTo summarize, we have shown that if we choose the points x0, .