Conte, de Boor - Elementary Numerical Analysis. An Algorithmic Approach (523140), страница 55
Текст из файла (страница 55)
Ifsome information about the size of the derivative appearing in the errorterm is available, one simply determines h or N so as to guarantee an errorless than a prescribed tolerance.Example 7.4 Determine N so that the composite trapezoid rule (5.33) gives the value ofcorrect to six digits after the decimal point, assuming thatcan be calculatedaccurately, and compute the approximation.In this example, f(x) =a = 0, b = 1, h - l/N; hence the error in thecomposite trapezoid rule is -f´´( η) N - 2 / 1 2 , f o r s o m e η(a,b).
Since we do not knowη, the best statement we can make is that the error is in absolute value no bigger thanWe computeFurther, f´´´(x) =4x(3 - 2x 2 ), which vanishes at x = 0 and x = ±max |f´´(x)| on [0, 1] must occur at x = 0 or at the end points x = 0, 1: thusHenceWe are therefore guaranteed six-place accuracy (after the decimal point) if we choose Nsuch thatororThe computer output below shows this to be a slight overestimate for N.As computed on an IBM 7094 in both single precision (SP) and double precision(DP), the results for various values of N are:NI(SP)I(DP)ERROR( S P )ERROR( D P )501002004008007.4679947E-017.4681776E-017.4682212E-017.4682275E-017.4682207E-017.4670061D-017.4681800D-017.4682260D-017.4682375D-017.4682404D-012.466E-056.37 E-062.01 E-061.56 E-062.06 E-062.452D-056.13 D-061.53 D-063.8 D-079.D-08The value of I correct to eight significant figures is I = 0.74682413.
It thus appearsthat in single-precision arithmetic we cannot obtain six-place accuracy, no matter howmany subdivisions we take. Indeed, for N = 800, the results are worse than those forN = 400. This shows that round-off error has affected the last three figures. Thedouble-precision results show that for N = 400 we have six-decimal-place accuracy,somewhat earlier than predicted above.7.4NUMERICAL INTEGRATION: COMPOSITE RULES323The FORTRAN program is:FORTRAN PROGRAM FOR EXAMPLE 7.4 (SINGLE PRECISION)CEXAMPLE 7.4 .
TRAPEZIOD RULE .INTEGER I,NREAL A,B,H,TF(X) = EXP(-X*X)1 PRINT 601601 FORMAT(' EXAMPLE 7.4 TRAPEZOIDAL INTEGRATION')READ 501, A,B,N501 FORMAT(2E20.0,15)IF (N .LT. 2)STOPT = F(A)/2.H = (B- A)/FLOAr(N)DO 2 I=l,N-12T = F(A + FLOAT(I)*H) + TT = (F(B)/2. + T)*HPRINT 602, A,B,N,T602 FORMAT(' INTEGRAL FROM A = ',lPE14.7,' 'TO B = ',E14.7,' FOR N = ',I5,' IS ',E14.7)*GO TO 1ENDIf we use the corrected trapezoid rule (7.50) instead, the required N dropsdramatically. We now have the error bounded byOne calculateshenceFor six-place accuracy, it is therefore sufficient thatororso that only 14 subintervals are required as compared with 578 for the compositetrapezoid rule without the differential end correction.As this example illustrates, higher-order formulas can reduce thenecessary number of function evaluations tremendously over lower-orderrules if the higher-order derivatives of the integrand are approximately thesame size as the lower-order derivatives.
Gaussian rules, in particular, canbe very effective.In the absence of information about the size of the appropriatederivative of f(x), it is possible only to apply the composite rules forvarious values of N, thus producing a sequence I N of approximations toI(f) which, theoretically, converges to I(f) as Nif f(x) is sufficientlysmooth. One terminates this process when the difference between successive estimates becomes “sufficiently small.” The dangers of such a procedure have been discussed in Sec.
1.6. An added difficulty arises in this case324DIFFERENTIATION AND INTEGRATIONfrom round-off effects, which increase with increasing N. The computerresults in Example 7.4 show this very clearly.Example 7.5 Write a program for the corrected trapezoid rule and solve the problem ofExample 7.4 using this program.FORTRAN PROGRAMC600110610EXAMPLE 7.5 . CORRECTED TRAPEZOID RULEINTEGER I,NREAL A,B,CORTRP,H,THAPF(X) = EXP(-X*X)FPRIME(X) = -2.*X*F(X)DATA A,B /0., 1. /PRINT 660FORMAT(9X,'N',7X,'TRAPEZOID SUM',7X,'CORR.TRAP.SUM')DO 10 N = 10,15H = (B - A)/FLOAT(N)TRAP = (F(A) + F(B))/2.DO 1 I=l,N-1TRAP = TRAP + F(A + FLOAT(I)*H)TRAP = H*TRAPCORTRP = TRAP + H*H*(FPRTME(A) - FPRIME(B))/l2.PRINT 610, N,TRAP,CORTRPFORMAT(I10,2E20.7)STOPENDSingle precision outputN101112131415TRAPEZOID SUM0.7462108E 000.7463173E 000.7463983E 000.7464612E 000.7465112E 000.7465516E 00CORR.TRAP.SUM0.7468239E 000.7468240E 000.7468240E 000.7468240E 000.7468240E 000.7468241E 00Double precision outputN101112131415TRAPEZOID SUM7.4621080E-017.4631727E-017.4639825E-0174646126E-017.4651126E-017.4655159E-01CORR.TRAP.SUM7.4682393E-017.4682399E-017.4682403E-017.4682406E-017.4682408E-017.4682409E-01Example 7.6 Write a program for Simpson’s rule and solve the problem of Example 7.4using this program in both single precision and double precision.The FORTRAN program and the results obtained on an IBM 7094 are givenbelow for N = 25, 50, and 100 subdivisions.
Note that the results in single precision areagain worse for N = 50, 100 than for N = 25, indicating round-off-error effects. Thedouble-precision results are all correct to the number of figures given. On comparingthese results with those of Examples 7.4 and 7.5, we see that both Simpson’s rule and thecorrected trapezoid rule are much more efficient than the trapezoid rule.7.4CNUMERICAL INTEGRATION: COMPOSITE RULES325PROGRAM FOR EXAMPLE 7.6 . SIMPSON'S RULE .INTEGER I,NREALA,B,H,HALF,HOVER2,S,XF(X) = EXP(-X*X)PRINT 600600 FORMAT(' EXAMPLE 7.6 SIMPSON''S RULE'/)1 READ 501, A,B,N501 FORMAT(2E20.0,15)IF (N .LT. 2)STOPH = (B - A)/FLOAT(N)HOVER2 = H/2.S = 0.HALF = F(A + HOVER2)DO 2 I=l,N-1X = A + FLOAT(I)*HS = S + F(X)2HALF = HALF + F(X+HOVER2)S = (H/6.)*(F(A) + 4.*HALF + 2.*S + F(B))PRINT 602, A,B,N,S602 FORMAT(' INTEGRAL FROM A = ',lPE14.7,' TO B = ',E14.7,*' FOR N = ',I5,' IS ',E14.7)GO TO 14 FORMAT(2E20.0,15)ENDCOMPUTER RESULTS FOR EXAMPLE 7.6NI(SP)E R R O R( S P )I(DP)E R R O R( D P )25501007.4682406E-017.4682400E-017.4682392E-017.
E-071.3E-062.1E-067.4682413D-017.4682413D-017.4682413D-010.0.0.Finally, composite rules based on gaussian formulas can also bederived. To be consistent with the composite rules already discussed, werestrict ourselves to definite integrals of the formWe again subdivide the interval (a,b) into N equally spaced panels so thatxi = a + ihi = 0, 1, .
. . , N withh = (b - a) /NWe wish to apply gaussian quadrature to the integral over the ith interval,i.e., to(7.52)The gaussian weights and points based on Legendre polynomials given inSec. 7.3 assume that the limits of integration are from -1 to +1. Hencewe first make the linear change of variablesw i t h x i - ½ = (xi+ xi-1 )/2326DIFFERENTIATION AND INTEGRATIONand substitute into (7.52) to obtainwhereWe now approximate the integral Ii with the gaussian formula on k + 1points to obtainIiA0gi (ξ0) + A1gi (ξ1) + · · · + Akgi (ξk )(7.53)where the weights and abscissas are taken from LGNDRE in Sec.
7.3.Finally, on summing over the N subintervals we obtainwhich from (7.53) gives the approximation(7.54 a)Notice that the weights are independent of i.According to the error equation (7.40), the error over the single panel(xi-1, xi ) is expressible in the formfor some η i in [-1, 1]but this means thatxi-1 < η´i < xiHence the error over the interval (a,b) can be expressed as(7.54b)Example 7.7 Evaluate the integral I =using gaussian quadrature withk = 3 and N = 2 subdivisions of the interval [1, 3]. See Example 7.3.CPROGRAM FOR EXAMPLE 7.7. COMPOSITE FOUR-POINT GAUSS-LEGENDRE.INTEGER I,NREAL A,B,H,HOVER2,P1,P2,POINT(2),S,S1,S2,WEIGHT(2),XDATA POINT,WEIGHT / .33998 10436, .86113 63116,*.65214 51549, .34785 48451 /F(X) = SIN(X)**2/XPRINT 600600 FORMAT(' EXAMPLE 7.7 FOUR-POINT GAUSS-LEGENDRE'/)1 READ 501, A,B,N7.4NUMERICAL INTEGRATION: COMPOSITE RULES327501FORMAT(2E20.0,I5)STOPIF (N .LT.
1)H = (B - A)/FLOAT(N)HOVER2 = H/2.P1 = POINT(l)*HOVER2P2 = POINT(2)*HOVER2Sl = 0.S2 = 0.DO 2 I=l,NX = A + FLOAT(I)*H - HOVER2Sl = Sl + F(-Pl+X) + F(P1+X)2S2 = S2 + F(-P2+X) + F(P2+X)S = HOVER2+(WEIGHT(l)*Sl + WEIGHT(2)*S2)PRINT 602, A,B,N,S602 FORMAT(' INTEGRAL FROM A = ',1PE14.7,' TO B = ',E14.7,*' FOR N = ',I3,' IS ',E14.7)GO TO 1ENDThe answer, as obtained on a UNIVAC 1110 in single precision, is 0.794825 17,which is in error by less than 3 in the last place.EXERCISES7.4-l Derive the composite trapezoid rule TN (7.49) and the composite midpoint rule MN(7.48).7.4-2 Derive the composite corrected trapezoid rule CTN (7.50) and verify that the interiorderivatives f´(xi) (i = 1, .
. . , N - 1) cancel out in the sum.7.4-3 Write a program for the composite Simpson rule. Inputs to the program should be f(x),the interval [a,b] and the number of subdivisions N. Use this program to calculatewith N = 10 and N = 20 subdivisions.7.4-4 Use the program for Simpson’s rule to calculate an approximation to the integralswhich are correct to six decimal places.
Do this by starting with N = 10 and doubling N untilyou are satisfied that you have the required accuracy.7.4-5 Write a program for the corrected trapezoid rule. In this case input will consist off(x),f´(x), [a,b], and N. Apply this program to the integral in Exercise 7.4-3 and compare theresults with those given by Simpson’s rule.7.4-6 Write a program for the composite gaussian rule (7.54a) using k = 3. Use it to evaluatethe integral in Exercise 7.4-3 first with N = 2 and then with N = 4 subdivisions. Compare theamount of computational effort and the accuracy obtained with those required by Simpson’srule.7.4-7 The error function erf(x ) is defined byUse the gaussian composite rule for k = 3 to evaluate erf(0.5) again with N = 2 and N = 4subdivisions.
Estimate the accuracy of your result and compare with the correct valueerf(0.5) = 0.520499876.328DIFFERENTIATION AND INTEGRATION7.4-8 The determination of the condensation of a pure vapor on the outside of a cooledhorizontal tube requires that the mean heat-transfer coefficient Q be computed.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
