Thompson - Computing for Scientists and Engineers (523188), страница 27
Текст из файла (страница 27)
The program uses a while loop controlled by n, the number of steps inthe integration interval, with a zero for n being used to terminate program execution. The analytical value for the working-function integral, IntAnalyt, is obtainedfrom (4.4) and Table 4.1, which are coded in function YWInt. The trapezoid ruleis coded in function TrapInt and its value in the main program is IntByTrap.The function TrapInt integrates the function y, which here is coded as theworking-function expression (4.1).
In the C function for y, the roots are given either as exact decimals or as ratios of (floating-point) integers. The roots are therebycomputed to machine accuracy. Similarly, maximum accuracy is gained in the Cfunction YWInt for the working-function integral.Exercise 4.21(a) Code and test the trapezoid parts of program Trapezoid and SimpsonIntegrals. You can test the program by temporarily substituting for the working function any linear function of x, and substituting in the analytic-integralfunction the corresponding integral.
The safest way to code these temporarychanges is to add them just before the return statements in the two C functions. By doing this you do not have to disturb the code for the working function. Integration of a linear function should produce results accurate to machineroundoff error.(b) Run the program for the integral from a = - 0.2 up to 0.2, which Figure 4.1 suggests has a very small value because of the zero of the function atx = 0. Thus, for h = 0.1 you need to input n = 4, and so on for three successively halved values of h. Check your results against mine in Table 4.9.140NUMERICAL DERIVATIVES AND INTEGRALS(c) Use (4.44) to estimate the error in the trapezoidal rule for each choice of nand h, using the second derivative given in Table 4.3, namely the value atx = 0, as an estimator of the appropriate derivative.
Show that this over-predictsthe actual error by a factor of about 6. Suggest why. nTABLE 4.9 Trapezoid integrals for the working function (4.1) from -0.2 to 0.2 with variousstepsizes h. The analytical value of the integral is 9.9109 x 10-3.Notice, from comparing Table 4.9 with Table 4.4, relatively how much moreaccurate the simplest formula for numerical integration (trapezoid rule) is comparedwith the simplest formula for numerical differentiation (central-difference formula).This justifies our remarks at the start of this section about integrals being consideredas antiderivatives.We return to the trapezoid rule and consider its application to functions that areless pathological than our working function at the end of the next subsection, whereit is also compared with the Simpson rule.Simpson formula and program for integralsIn the trapezoid integration formula we assumed straight-line segments betweenpairs of “knots” (points where the function is to be evaluated), and this collapsed thegeneral integration formula to the basic trapezoid formula.
If we instead assumeparabolic segments between triples of knots, as shown in Figure 4.8, we can derivethe Simpson formula and an estimate of its error.Simpson rule formula. The method of deriving this formula is the same as for thetrapezoid formula (Exercise 4.20), except that having three values of the functionallows more of the (usually) unknown derivatives to be eliminated. Try it and see.Exercise 4.22(a) To derive the Simpson integration rule, start from the general integral expression (4.40), expand about x0 = a + h, then solve for the second derivativein terms of the function values at x = a, x = a + h, and x = a + 2h.
Keep1414.6 NUMERICAL INTEGRATION METHODSthe term containing the fourth derivative. Thus derive the basic Simpson formula for integration(4.45)(b) Show that if y is a cubic or lower-order polynomial, then the Simpson formula is exact.(c) Use additivity of integration and the basic integration formula (4.45) to derive the composite Simpson rule for n intervals, where n must be even,(4.46)where n2 = n /2, and the error estimate for the composite Simpson rule is(4.47)in which the second derivative is bounded by its maximum value in the intervalx = a to x = a + nh.
nThe parabolic shape that the Simpson formula assumes for the function in the basicinterval of width 2h is shown as the dashed curve in Figure 4.8 for our workingfunction (4.1). A different parabola is assumed for each triple of points, so that therather poor fit of the dashed to the solid curve outside the range [-O.1,O.l] is not ofimmediate concern.
The discontinuities of the parabolic curves between successivepanels of width 2h that occur in the Simpson formula are removed if one uses thecubic-spline technique developed in Chapter 5.Simpson rule program. The overall program structure for Trapezoid and Simpson Integrals was described in the preceding subsection on the trapezoid formulaand program, so here we describe only parts specific to the Simpson formula. Oninput the value of n is checked; only if n is even is the Simpson algorithm used toproduce the value IntBySimp, which is compared with IntAnalyt to obtain thevalue error for output.The function SimpInt uses the composite Simpson rule (4.46) to estimate numerically the integral of the function y. In case this function is to be used in anotherprogram, there is a trap to ensure that the number of points, n, is even.
For thecomposite rule (4.46) n > 2 is also required by the program. (For n = 2 thebasic Simpson formula (4.45) may be used.) The check for even n in SimpInt142NUMERICAL DERIVATIVES AND INTEGRALScan be tested by temporarily disabling the test made in the main program.
Note thatthe function is coded to resemble closely the formula (4.46), rather than to be veryefficient. If you need a speedy version, you should recode most of SimpInt.Now that you know how the Simpson integration formula is supposed to work,it’s time to plug and play.Exercise 4.23(a) Test the Simpson rule parts of the program Trapezoid and Simpson Integrals. Temporarily substitute for the working function any function of x ashigh as third order, and substitute in the analytic-integral function the corresponding integral. Code these temporary changes just before the returnstatements in the two C functions, so as not to disturb the code for the workingfunction. Integration of a cubic or lower-order polynomial should produce results accurate to machine roundoff error, according to (4.47).(b) Run the program for the Simpson integral from a = -0 .
2 up to 0 . 2,which has a very small value because of the zero of the working function atx = 0, as seen in (4.1) and Figure 4.1. For h = 0.1 input n = 4, and so onfor three successively halved values of h. Check your results against those inTable 4.10.(c) Use (4.47) to estimate the error in the Simpson rule for each choice of nand h, using the fourth derivative in Table 4.3 at x = 0 as an estimator of theappropriate derivative.
Show that this overpredicts the actual error by a factorof about 2. Suggest why. nTABLE 4.10 Simpson-rule integrals for the working function (4.1) from -0.2 to 0.2 with various stepsizes h. The analytical integral has value 9.9109 x 10-3.By comparing Tables 4.9 and 4.10 for the trapezoid and Simpson formulas, respectively, notice how the latter improves much more rapidly than the former as hdecreases, since there is an h 5 dependence of the error in the Simpson rule, compared with an h 3 dependence in the trapezoid-rule error. This occurs in spite of themuch larger fourth derivative than second derivative (Table 4.3) for the very wigglypolynomial that is our working function.4.6 NUMERICAL INTEGRATION METHODS143Integrals with cosinesNow that we have a working program that has been torture tested against a badlybehaved function, it is interesting to try a more-usual function such as the cosine.Since the cosine has derivatives that are bounded in magnitude by unity, the errorestimates for the trapezoid and Simpson rules, (4.44) and (4.47), will be dominatedby their dependence on h.
How about trying your integration program on the cosinefunction?Exercise 4.24(a) Modify the program Trapezoid and Simpson Integrals by replacing inthe C function y ( x )yeach h. Make a table and a log-log plot of the absolute value of the error againsth, as in Table 4.11 and Figure 4.9.(b) Assume that the graph for each method (trapezoid and Simpson) on the loglog plot is a straight line. Thus show that the error has a power-law dependenceon stepsize h, and that the overall scale of the error estimates the factor precedingthe power of h.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
