f4-1

124 Chapter4. IntegrationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).of various orders, with higher order sometimes, but not always, giving higher accuracy. “Rombergintegration,”which is discussed in §4.3,is a generalformalism for making use of integration methods of a variety of different orders, and we recommend it highly. Apart from the methods of this chapter and of Chapter 16, there are yet other methods for obtaining integrals. One important class is based on function approximation. We discuss explicitly the integration of functions by Chebyshev approximation (“Clenshaw-Curtis” quadrature) in §5.9. Although not explicitly discussedhere,yououghttobeable tofigureouthowto do cubicsplinequadrature using the output of the routine splinein§3.3. (Hint: Integrate equation 3.3.3 over xanalytically. See [1].) Some integrals related to Fourier transforms can be calculated using the fast Fourier transform (FFT) algorithm. This is discussed in §13.9. Multidimensionalintegrals are anotherwhole multidimensionalbag of worms. Section 4.6 is an introductorydiscussion in this chapter; the importanttechnique of Monte-Carlo integration is treated in Chapter 7. CITED REFERENCES AND FURTHER READING: Carnahan, B., Luther, H.A., and Wilkes, J.O. 1969, Applied Numerical Methods (New York: Wiley), Chapter 2. Isaacson,E.,andKeller,H.B.1966, AnalysisofNumericalMethods (NewYork:Wiley),Chapter7. Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), Chapter 4. Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), Chapter 3. Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis , 2nd ed. (New York: McGraw-Hill), Chapter 4. Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), §7.4. Kahaner,D.,Moler,C.,andNash,S.1989, NumericalMethods andSoftware (EnglewoodCliffs, NJ: Prentice Hall), Chapter 5. Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for Mathematical Computations (Englewood Cliffs, NJ: Prentice-Hall), §5.2, p. 89. [1] Davis, P., and Rabinowitz, P. 1984, Methods of Numerical Integration , 2nd ed. (Orlando, FL: Academic Press). 4.1 Classical Formulas for Equally Spaced Abscissas Where wouldany bookon numericalanalysis be withoutMr. Simpsonand his “rule”? The classical formulas for integrating a function whose value is known at equally spaced steps have a certain elegance about them, and they are redolent withhistoricalassociation. Throughthem,themodernnumericalanalystcommuneswith the spirits of his or her predecessors back across the centuries, as far as the time of Newton, if not farther. Alas, times dochange; with the exception of two of the most modest formulas (“extended trapezoidal rule,” equation 4.1.11, and “extended 4.1ClassicalFormulasforEquallySpacedAbscissas 125Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).x0 xNxN + 1 open formulas use these points closed formulas use these pointsx1x2h Figure 4.1.1. Quadrature formulas with equally spaced abscissas compute the integral of a function between x0and xN+1. Closed formulas evaluate the function on the boundary points, while open formulas refrain from doing so (useful if the evaluation algorithm breaks down on the boundary points). midpointrule, ”equation4.1.19,see §4.2),the classical formulas are almost entirely useless. They are museum pieces, but beautiful ones. Some notation: We have a sequence of abscissas, denoted x0,x1,...,x N, xN+1which are spaced apart by a constant step h, xi=x0+ih i =0 ,1,...,N +1 ( 4.1.1 ) A function f(x)has known values at the xi’s, f(xi)≡fi (4.1.2 ) We want to integrate the function f(x)between a lower limit aand an upper limit b, where aand bare each equal to one or the other of the xi’s. An integration formula that uses the value of the function at the endpoints, f(a)orf(b), is called aclosedformula. Occasionally, we want to integrate a function whose value at one or both endpoints is dif ficult to compute (e.g., the computation of fgoes to a limit of zero over zero there, or worse yet has an integrable singularity there). In this case we want an openformula, which estimates the integral using only xi’s strictly between aand b(see Figure 4.1.1). The basic building blocks of the classical formulas are rules for integrating a function over a small number of intervals. As that number increases, we can find rules that are exact for polynomials of increasingly high order. (Keep in mind that higher order does not always imply higher accuracy in real cases.) A sequence ofsuch closed formulas is now given. Closed Newton-CotesFormulas Trapezoidal rule: /integraldisplayx2 x1f(x)dx=h/bracketleftbigg1 2f1+1 2f2/bracketrightbigg +O(h3f/prime/prime)( 4.1.3 ) Here the error term O()signifies that the true answer differs from the estimate by anamountthatistheproductofsomenumericalcoef ficienttimes h3timesthevalue 126 Chapter4. IntegrationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).of the function ’s second derivative somewhere in the interval of integration. The coefficient is knowable, and it can be found in all the standard references on this subject. The point at which the second derivative is to be evaluated is, however, unknowable. If we knew it, we couldevaluate the functionthere and have a higher- order method! Since the product of a knowable and an unknowableis unknowable,we will streamline our formulas and write only O(), instead of the coef ficient. Equation(4.1.3)isatwo-pointformula( x 1and x2). Itis exactforpolynomials up to and including degree 1, i.e., f(x)= x. One anticipates that there is a three-pointformulaexactuptopolynomialsofdegree2. Thisistrue;moreover,bya cancellationofcoef ficientsduetoleft-rightsymmetryoftheformula,thethree-point formulais exact for polynomialsup to and includingdegree3, i.e., f(x)= x3: Simpson’s rule: /integraldisplayx3 x1f(x)dx=h/bracketleftbigg1 3f1+4 3f2+1 3f3/bracketrightbigg +O(h5f(4))( 4.1.4 ) Here f(4)means the fourth derivative of the function fevaluated at an unknown place in the interval. Note also that the formula gives the integral over an interval of size 2h, so the coef ficients add up to 2. There is no lucky cancellation in the four-point formula, so it is also exact for polynomials up to and including degree 3. Simpson’s3 8rule: /integraldisplayx4 x1f(x)dx=h/bracketleftbigg3 8f1+9 8f2+9 8f3+3 8f4/bracketrightbigg +O(h5f(4))(4.1.5 ) Thefive-point formula again bene fits from a cancellation: Bode’s rule: /integraldisplayx5 x1f(x)dx=h/bracketleftbigg14 45f1+64 45f2+24 45f3+64 45f4+14 45f5/bracketrightbigg +O(h7f(6))(4.1.6 ) This is exact for polynomials up to and including degree 5. At this point the formulas stop being named after famous personages, so we will not go any further. Consult [1]for additional formulas in the sequence. Extrapolative Formulasfor a SingleInterval We are going to depart from historical practice for a moment. Many texts would give, at this point, a sequence of “Newton-Cotes Formulas of Open Type. ” Here is an example: /integraldisplayx5 x0f(x)dx=h/bracketleftbigg55 24f1+5 24f2+5 24f3+55 24f4/bracketrightbigg +O(h5f(4)) 4.1ClassicalFormulasforEquallySpacedAbscissas 127Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Notice that the integral from a=x0tob=x5is estimated, using only the interior points x1,x2,x3,x4. In our opinion, formulas of this type are not useful for the reasonsthat(i)theycannotusefullybestrungtogethertoget “extended”rules,aswe are aboutto do with the closed formulas,and (ii) for all otherpossible uses theyare dominatedby the Gaussian integrationformulaswhich we will introducein §4.5. Instead of the Newton-Cotes open formulas, let us set out the formulas for estimating the integral in the single interval from x0tox1, using values of the function fatx1,x2,.... These will be useful building blocks for the “extended” open formulas. /integraldisplayx1 x0f(x)dx=h[f1]+ O(h2f/prime)( 4.1.7 ) /integraldisplayx1 x0f(x)dx=h/bracketleftbigg3 2f1−1 2f2/bracketrightbigg +O(h3f/prime/prime)( 4.1.8 ) /integraldisplayx1 x0f(x)dx=h/bracketleftbigg23 12f1−16 12f2+5 12f3/bracketrightbigg +O(h4f(3))( 4.1.9 ) /integraldisplayx1 x0f(x)dx=h/bracketleftbigg55 24f1−59 24f2+37 24f3−9 24f4/bracketrightbigg +O(h5f(4))(4.1.10 ) Perhaps a word here would be in order about how formulas like the above can bederived. Thereareelegantways,butthemoststraightforwardistowritedownthe basic form of the formula, replacing the numerical coef ficients with unknowns, say p, q, r, s. Withoutlossofgeneralitytake x0=0and x1=1,so h=1. Substitutein turn for f(x)(and for f1,f2,f3,f4) the functions f(x)=1,f(x)= x,f(x)= x2, and f(x)= x3. Doing the integral in each case reduces the left-hand side to a number, and the right-hand side to a linear equation for the unknowns p, q, r, s. Solving the four equations produced in this way gives the coef ficients. Extended Formulas(Closed) If we use equation (4.1.3) N−1times, to do the integration in the intervals (x1,x2),(x2,x3),..., (xN−1,x N),andthenaddtheresults,weobtainan “extended” or“composite ”formula for the integral from x1toxN. Extended trapezoidal rule: /integraldisplayxN x1f(x)dx=h/bracketleftbigg1 2f1+f2+f3+ ···+fN−1+1 2fN/bracketrightbigg +O/parenleftbigg(b−a)3f/prime/prime N2/parenrightbigg (4.1.11 ) Herewehavewrittentheerrorestimateintermsoftheinterval b−aandthenumber of points Ninstead of in terms of h. This is clearer, since one is usually holding aand bfixed and wanting to know (e.g.) how much the error will be decreased 128 Chapter4. IntegrationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).by taking twice as many steps (in this case, it is by a factor of 4). In subsequent equationswe will show onlythescalingof theerrortermwiththe numberofsteps. For reasons that will not become clear until §4.2, equation (4.1.11) is in fact the most important equation in this section, the basis for most practical quadrature schemes. Theextended formula of order 1/N3is: /integraldisplayxN x1f(x)dx=h/bracketleftbigg5 12f1+13 12f2+f3+f4+ ··· +fN−2+13 12fN−1+5 12fN/bracketrightbigg +O/parenleftbigg1 N3/parenrightbigg (4.1.12 ) (We will see in a moment where this comes from.) If we apply equation (4.1.4) to successive, nonoverlapping pairsof intervals, we get the extended Simpson’s rule: /integraldisplayxN x1f(x)dx=h/bracketleftbigg1 3f1+4 3f2+2 3f3+4 3f4+ ···+2 3fN−2+4 3fN−1+1 3fN/bracketrightbigg +O/parenleftbigg1 N4/parenrightbigg (4.1.13 ) Notice that the 2/3, 4/3 alternation continues throughout the interior of the evalu- ation. Many people believe that the wobbling alternation somehow contains deepinformation about the integral of their function that is not apparent to mortal eyes. In fact, the alternation is an artifact of using the building block (4.1.4). Another extended formula with the same order as Simpson ’s rule is /integraldisplay xN x1f(x)dx=h/bracketleftbigg3 8f1+7 6f2+23 24f3+f4+f5+ ···+fN−4+fN−3+23 24fN−2+7 6fN−1+3 8fN/bracketrightbigg +O/parenleftbigg1 N4/parenrightbigg(4.1.14 ) This equationis constructedby fitting cubicpolynomialsthroughsuccessive groups of four points; we defer details to §18.3, where a similar technique is used in the solution of integral equations. We can, however, tell you where equation (4.1.12) came from. It is Simpson ’s extended rule, averaged with a modi fied version of itself in which the first and last step are done with the trapezoidal rule (4.1.3). The trapezoidalstepis twoorderslowerthanSimpson ’srule;however,itscontributionto the integral goes down as an additionalpower of N(since it is used only twice, not Ntimes). This makes the resulting formulaof degree oneless than Simpson. 4.1ClassicalFormulasforEquallySpacedAbscissas 129Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Extended Formulas(Openand Semi-open) We can construct open and semi-open extendedformulas by adding the closed formulas (4.1.11) –(4.1.14), evaluated for the second and subsequent steps, to the extrapolative open formulas for the first step, (4.1.7) –(4.1.10). As discussed immediately above, it is consistent to use an end step that is of one order lower than the (repeated) interior step. The resulting formulas for an interval open at both ends are as follows: Equations (4.1.7) and (4.1.11) give /integraldisplayxN x1f(x)dx=h/bracketleftbigg3 2f2+f3+f4+···+fN−2+3 2fN−1/bracketrightbigg +O/parenleftbigg1 N2/parenrightbigg (4.1.15 ) Equations (4.1.8) and (4.1.12) give /integraldisplayxN x1f(x)dx=h/bracketleftbigg23 12f2+7 12f3+f4+f5+ ··· +fN−3+7 12fN−2+23 12fN−1/bracketrightbigg +O/parenleftbigg1 N3/parenrightbigg(4.1.16 ) Equations (4.1.9) and (4.1.13) give /integraldisplayxN x1f(x)dx=h/bracketleftbigg27 12f2+0+13 12f4+4 3f5+ ···+4 3fN−4+13 12fN−3+0+27 12fN−1/bracketrightbigg +O/parenleftbigg1 N4/parenrightbigg(4.1.17 ) The interior points alternate 4/3 and 2/3. If we want to avoid this alternation, we can combine equations (4.1.9) and (4.1.14), giving /integraldisplayxN x1f(x)dx=h/bracketleftbigg55 24f2−1 6f3+11 8f4+f5+f6+f7+ ···+fN−5+fN−4+11 8fN−3−1 6fN−2+55 24fN−1/bracketrightbigg +O/parenleftbigg1 N4/parenrightbigg (4.1.18 ) We should mention in passing another extended open formula, for use where thelimitsofintegrationarelocatedhalfwaybetweentabulatedabscissas. Thisoneis knownas the extendedmidpointrule , andis accurateto thesame orderas (4.1.15): /integraldisplayxN x1f(x)dx=h[f3/2+f5/2+f7/2+ ···+fN−3/2+fN−1/2]+ O/parenleftbigg1 N2/parenrightbigg (4.1.19 ) 130 Chapter4. Integrationof FunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).N = 1 234 (total after N = 4) Figure 4.2.1. Sequential calls to the routine trapzdincorporate the information fromprevious calls and evaluate the integrand only at those new points necessary to re fine the grid. The bottom line shows the totality of function evaluations after the fourth call. The routine qsimp, by weighting the intermediate results, transforms the trapezoid rule into Simpson ’s rule with essentially no additional overhead. There are also formulas of higher order for this situation, but we will refrain from giving them. Thesemi-openformulas arejusttheobviouscombinationsofequations(4.1.11) – (4.1.14) with (4.1.15) –(4.1.18), respectively. At the closed end of the integration, use the weights from the former equations; at the open end use the weights from the latter equations. One example should give the idea, the formulawith error term decreasing as 1/N3which is closed on the right and open on the left: /integraldisplayxN x1f(x)dx=h/bracketleftbigg23 12f2+7 12f3+f4+f5+ ··· +fN−2+13 12fN−1+5 12fN/bracketrightbigg +O/parenleftbigg1 N3/parenrightbigg (4.1.20 ) CITED REFERENCES AND FURTHER READING: Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions , Applied Mathe- matics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), §25.4. [1] Isaacson, E., and Keller, H.B. 1966, Analysis of Numerical Methods (New York: Wiley), §7.1. 4.2 Elementary Algorithms Ourstartingpointis equation(4.1.11),theextendedtrapezoidalrule. Thereare two facts about the trapezoidal rule which make it the starting point for a variety of algorithms. One fact is rather obvious, while the second is rather “deep.” Theobviousfactisthat,fora fixedfunction f(x)tobeintegratedbetween fixed limits aand b, one can double the number of intervals in the extended trapezoidal rule without losing the bene fit of previous work. The coarsest implementation of the trapezoidal rule is to average the function at its endpoints aand b. Thefirst stage of re finementis to add to this averagethe value of the functionat the halfway point. The second stage of re finement is to add the values at the 1/4 and 3/4 points. And so on (see Figure 4.2.1). Without further ado we can write a routine with this kind of logic to it: