Conte, de Boor - Elementary Numerical Analysis. An Algorithmic Approach (523140), страница 54
Текст из файла (страница 54)
. . , xkin (7.33) as the zeros of the polynomial pk+1(x) of degree k + 1 which isorthogonal with respect to the weight function w(x) over the interval (a,b)to any polynomial of degree < k, and if the coefficients Ai ( i = 0, . . . , k)in (7.33) are chosen according to (7.34), the resulting gaussian formula(7.33) will then be exact for all polynomials of degree < 2k + 1. Quadrature rules of this type are said to be “best possible” in the sense defined,and under the conditions given above.We now give some examples. First, let w(x) = 1.
If (a,b) is a finiteinterval, then the linear change of variables x = [(b - a) t + (b + a) ] / 27.3NUMERICAL INTEGRATION: GAUSSIAN RULES315can be used to change the limits of integration from (a,b) to (-1, 1). Withthis,(7.41)Assuming that this transformation has already been made, we consider theintegral (7.31) to be in the formSince w(x) = 1, the appropriate orthogonal polynomials are the Legendrepolynomials (see Example 6.6).
In this caseP1 (x) = xξ0 = 0etc. If we choose k = 1, thensubstituting into (7.33) and (7.40), we obtainand(7.42)wheresinceSubstituting these constants into (7.42), we obtain the gaussian two-pointquadrature formula(7.43)with the error(7.44)For k > 1, both the points ξi and the weights Ai become irrational.Their calculation, however, is straightforward. We record these nodes and316DIFFERENTIATION AND INTEGRATIONweights for k = 0, .
. . , 5 in the following FORTRAN subroutineLGNDRE. Note that (former) FORTRAN restrictions have forced us tonumber nodes and weights from 1 through NP = k + 1 rather than from 0through k. Thus, the input parameter NP specifies the number of pointsrather than the degree of the underlying polynomial.SUBROUTINE LGNDRE ( NP , POINT, WEIGHT )C SUPPLIES POINTS AND WEIGHTS FOR GAUSS-LEGENDRE QUADRATIUREC INTEGRAL(F(X), -1 .LE.
X .LE. 1) IS APPROXIMATELY EQUAL TOCSUM(F(POINT(I))*WEIGHT(I), I=l,...,NP) .INTEGER NP, IREAL POINT(NP),WEIGHT(NP)IF (NP .GT. 6) THENPRINT 600,NPFORMAT(' THE GIVEN NUMBER NP =',12,' IS GREATER THAN 6.'600*/' EXECUTION STOPPED IN SUBROUTINE L G N D R E . ' )STOPEND IFGO TO (1,2,3,4,5,6),NP1 POINT(l) = 0.WEIGHT(l) = 2.GO TO 992 POINT(2) = .57735 02691 89626 D0WEIGHT(2) = 1.GO TO 953 POINT(2) = 0.POINT(3) = .77459 66692 41483 D0WEIGHT(2) = .88888 88888 88888 9 D0WEIGHT(3) = .55555 55555 55555 6 D0GO TO 954 POINT(3) = .33998 10435 84856 D0POINT(4) = .86113 63115 54053 D0WEIGHT(3) = .65214 51548 62546 D0WEIGHT(4) = .34785 48451 37454 D0GO TO 955 POINT(3) = 0.POINT(4) = .53846 93101 85683 D0POINT(5) = .90617 98459 38664 D0WEIGHT(3) = .56888 88088 88888 9 D0WEIGHT(4) = .47862 86704 99366 D0WEIGHT(5) = .23692 68850 56189 D0GO TO 956 POINT(4) = .23861 91860 83197 D0POINT(5) = .66120 93864 66265 D0POINT(6) = .93246 95142 03152 D0WEIGHT(4) = .46791 39345 72691 D0WEIGHT(5) = .36076 15730 48139 D0WEIGHT(6) = .l7132 44923 79170 D0C95 DO 96 I=1,NP/2POINT(I) = -POINT(NP+l-I)WEIGHT(I) = WEIGHT(NP+l-I)96RETURN99ENDExample 7.2 For comparison purposes, we again wish to evaluate I =but using the gaussian five-point formula (k = 4).The required change of variables (7.41) here is x = ( t + 1)/2, soNaturally, we use a program to carry out the calculation.7.3CNUMERICAL INTEGRATION: GAUSSIAN RULES317EXAMPLE 7.2 GAUSSIAN INTEGRATIONREAL INTGRL,P(5),WEIGHT(5)F(T) = EXP(-(l.+T)**2/4.)/2.CALL LGNDRE ( 5, P, WEIGHT )INTGRL = WEIGHT(l)*(F(P(l))+F(P(5))) + WEIGHT(2)*(F(P(2))+F(P(4)))*+ WEIGHT(3)*F(P(3))PRINT 600,INTGRL600 FORMAT(' EXAMPLE 7.2.
GAUSS QUADRATURE'/' INTEGRAL = ',1PE14.7)STOPENDThis gives the outputINTEGRAL = 7.4682413-001To achieve comparable accuracy with the trapezoidal rule would require some 2,800subdivisions, whereas Simpson’s rule would require about 20 subdivisions.Example 7.3 Find an approximation tousing gaussian quadrature with k = 3. (The correct value is I = 0.79482518 . . .
.)We again transform to the interval [-1, 1], this time by the change of variablex = t + 2. This yieldsAfter changing the body of the program for Example 7.2 appropriately toF(T) = SIN(T + 2.)**2/(T + 2.)CALL LGNDRE (4, P, WEIGHT )INTGRL = WEIGHT (1)*(F(P(1)) + F(P(4))) + WEIGHT(2)*(F(P(2)) + F(P(3)))we obtain the outputINTEGRAL = 7.9482833-001Gaussian-type formulas are especially useful in dealing with singularis to be calculated, where f(x) has anintegrals.
If, for example,algebraic singularity at a and/or b, then one transforms the integral intowherew(x) = (1 - x ) α (1 + x) βfor appropriate exponents α and β. In this case, the ξi ‘s are the zeros of theappropriate Jacobi polynomial.
In the special case α = β = -½, these arejust the Chebyshev polynomials introduced in Example 6.7 and discussedin Sec. 6.1. For this special case, one gets the very attractive rule(7.45)for which all the weights Ai coincide, and for which the ξ i ’s are the318DIFFERENTIATION AND INTEGRATIONChebyshev points [see (6.18)](7.46)If the interval of integration is semi-infinite, it is at times of help totransform the integral intowithw(x) = x α e - xIn this case, the ξi ’s are the zeros of the appropriate Laguerre polynomial;see Example 6.9.
Finally, integrals of the formcan often be successfully estimated using the zeros of the appropriateHermite polynomial (see Example 6.8).For all these examples (and others), tables are available both for theξ i ’s and the weights Ai , the most recent, and probably most extensive, beingStroud and Secrest’s “Gaussian Quadrature Formulas” [20]. See also [27].EXERCISES7.3-l For which polynomials is Simpson’s rule exact?7.3-2 Construct a rule of the formwhich is exact for all polynomials of degree < 2.7.3-3 Calculatecorrect to four significant digits. [Hint: Transform the integral appropriately and use (7.45)and (7.46).]7.3-4 Find an estimate for7.3-5 Derive the weights Ai for the gaussian formula with k = 3, using the zeros ξ i given inLGNDRE.7.3-6 Use the gaussian five-point formula to obtain an estimate for the integrals given inExercises 7.4-3 and 7.4-4.7.3-7 Use Exercise 6.3-7 to show that (7.34) can also be written Ai =i = 0, .
. . , k. Conclude that gaussian weights are always positive.except that it7.3-8 Lobatto’s rule is a gaussian formula for integrating I =includes ± 1 as two fixed abscissas. It has the form [see (7.31)]Derive the Lobatto rule for the case k = 2 and show that it is exact for all polynomials ofdegree < 3.7.4NUMERICAL INTEGRATION: COMPOSITE RULES3197.3-9 Check out the subroutine LGNDRE by using it to calculatefor n = 0, 1, 2, .
. .For what values of n should the Gauss-Legendre rule on NP points give the integral exactly?7.4 NUMERICAL INTEGRATION: COMPOSITE RULESThe simple quadrature rules developed in the preceding sections to estimatewill usually not produce sufficiently accurate estimates, particularly whenthe interval [a, b] is reasonably large. It is customary in practice to dividethe given interval [a,b] into N smaller intervals and to apply the simplequadrature rules to each of these subintervals. We therefore subdivide theinterval [a, b] in such a way thata = x0 < x1 < x2 < · · · < xN = band we denote by gk(x) a piecewise-polynomial function (see Sec.
6.7) withbreakpoints {xi } (i = 1, . . . , N - 1). Furthermore, let Pi,k(x) ( i =1, . . . , N) denote the polynomial of degree < k which agrees with gk(x)on (xi-1, xi ). By the rules of integration we know thatand thatHence, approximating I(f) by I(gk) amounts to approximatingbyi=1,...,Nand summing the results. Evidently, on each subinterval (xi-1, xi ), we areproceeding just as in Secs. 7.2 and 7.3. In particular, we can apply any ofthe rules derived in Secs. 7.2 and 7.3 by substituting some polynomial forthe integrand, on each subinterval, and then summing the results.In the absence of any reason to do otherwise, we choose the xi ’s to beequally spaced,xi = a + ihWe also use, as in Sec.
2.6, the abbreviationf s = f(a + sh)so that fi = f(xi ), i = 0, . . . , N.320DIFFERENTIATION AND INTEGRATIONWe now consider specific examples. If we apply the rectangle rule(7.23) on each subinterval, we getfor the subinterval (xi-1,xi ). Summing, we obtain(7.47a)the composite rectangle rule (on N intervals). Its error is just the sum ofthe errors committed in each subinterval,whereIf f´(x) is continuous (as we assume), this can besimplified, using Theorem 1.2 in Sec. 1.7, as follows:so that, with Nh = b - a,some η(7.47b)We derive next the composite Simpson rule. Letting a = xi-1, b = xi ,and xi - xi-l = h in (7.28), we obtain for a single subintervalSumming for i = 1, .
. . , N, we obtainThe composite Simpson approximation SN can be simplified to yield(7.48a)7.4NUMERICAL INTEGRATION: COMPOSITE RULES321while the error term can be simplified, again using Theorem 1.2 of Sec. 1.7,toa < ξ < b(7.48 b)Note that in Simpson’s rule we must be able to evaluate the function at themidpoints x i-½ (i = 1, . . . , N) as well as at the breakpoints xi (i =0, 1, . .
. , N). This implies in particular that we always need an oddnumber of equally spaced points at which we know the value of theintegrand.In the same manner, one gets the composite midpoint rulefrom the midpoint rule (7.25), and the composite trapezoid rule(7.49b)from (7.26).From the corrected trapezoid rule (7.29), one obtains(7.50)Note that all the interior derivatives f´(xi ), i = 1, . . . , N - 1, cancel eachother when the results of applying the corrected trapezoid rule on eachsubinterval are summed. Hence the composite corrected trapezoid rule is,in fact, a corrected composite trapezoid rule, i.e.,(7.51)The corrected trapezoid rule has, of course, the disadvantage that thederivative of f(x) must be known or calculable [except when f(x) happenst o b e (b - a)-periodic].322DIFFERENTIATION AND INTEGRATIONIf any of these composite rules are to be applied, one has to determinefirst an appropriate N, or equivalently, an appropriate h = (b - a) /N.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
