f4-6

4.6MultidimensionalIntegrals 155Sample 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).Carnahan, B., Luther, H.A., and Wilkes, J.O. 1969, Applied Numerical Methods (New York: Wiley), §§2.9–2.10. Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis , 2nd ed. (New York: McGraw-Hill), §§4.4–4.8. 4.6 Multidimensional Integrals Integrals of functions of several variables, overregions with dimensiongreater than one, are noteasy. Thereare two reasons for this. First, the numberof function evaluations needed to sample an N-dimensional space increases as the Nth power of the number needed to do a one-dimensional integral. If you need 30 function evaluations to do a one-dimensional integral crudely, then you will likely need on theorderof30000evaluationsto reachthesame crudelevelfora three-dimensional integral. Second, the region of integration in N-dimensional space is defined by anN−1dimensional boundary which can itself be terribly complicated: It need not be convex or simply connected, for example. By contrast, the boundary of a one-dimensionalintegral consists of two numbers, its upperand lower limits. The first question to be asked, when faced with a multidimensional integral, is, “can it be reduced analytically to a lower dimensionality?” For example, so-called iterated integrals of a function of one variable f(t)can be reduced to one-dimensional integrals by the formula /integraldisplayx 0dtn/integraldisplaytn 0dtn−1···/integraldisplayt3 0dt2/integraldisplayt2 0f(t1)dt1 =1 (n−1)!/integraldisplayx 0(x−t)n−1f(t)dt(4.6.1 ) Alternatively, the function may have some special symmetry in the way it depends on its independent variables. If the boundary also has this symmetry, then the dimension can be reduced. In three dimensions, for example, the integration of asphericallysymmetricfunctionoverasphericalregionreduces,inpolarcoordinates, to a one-dimensional integral. The next questions to be asked will guide your choice between two entirely different approaches to doing the problem. The questions are: Is the shape of the boundary of the region of integration simple or complicated? Inside the region, is the integrandsmooth and simple, or complicated, or locally strongly peaked? Does the problem require high accuracy, or does it require an answer accurate only to a percent, or a few percent? If your answers are that the boundary is complicated, the integrand is not strongly peaked in very small regions, and relatively low accuracyis tolerable, then yourproblemis a good candidatefor MonteCarlo integration . This methodis very straightforward to program, in its cruder forms. One needs only to know a region with simple boundaries that includesthe complicated region of integration, plus a method of determining whether a random point is inside or outside the region of integration. Monte Carlo integration evaluates the function at a random sample of 156 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).points,and estimates its integralbased onthat randomsample. We will discuss it in more detail, and with more sophistication, in Chapter 7. If the boundaryis simple, and the function is very smooth, then the remaining approaches, breaking up the problem into repeated one-dimensional integrals, or multidimensional Gaussian quadratures, will be effective and relatively fast [1].I f you require high accuracy, these approaches are in any case the onlyones available to you,since MonteCarlo methodsarebynatureasymptoticallyslow to converge. Forlowaccuracy,userepeatedone-dimensionalintegrationormultidimensional Gaussianquadratureswhentheintegrandisslowlyvaryingandsmoothintheregion of integration, Monte Carlo when the integrand is oscillatory or discontinuous, butnot strongly peaked in small regions. If the integrand isstronglypeakedin small regions,andyouknowwherethose regionsare,breaktheintegralupintoseveralregionssothattheintegrandis smoothineach,anddoeachseparately. Ifyoudon’tknowwherethestronglypeakedregions are,youmightaswell (atthelevelofsophisticationofthisbook)quit: Itis hopeless to expectan integrationroutineto search outunknownpocketsoflargecontributionin a huge N-dimensional space. (But see §7.8.) If, on the basis of the aboveguidelines,you decide to pursue the repeatedone- dimensional integration approach, here is how it works. For definiteness, we willconsider the case of a three-dimensional integral in x, y, z-space. Two dimensions, or more than three dimensions, are entirely analogous. The first step is to specify the region of integration by (i) its lower and upper limits in x, which we will denote x 1and x2; (ii) its lower and upper limits in yat a specified value of x, denoted y1(x)and y2(x); and (iii) its lower and upper limits inzat specified xand y, denoted z1(x, y )and z2(x, y ). In other words, find the numbers x1and x2, andthefunctions y1(x),y2(x),z1(x, y ), and z2(x, y )such that I≡/integraldisplay/integraldisplay/integraldisplay dx dy dzf (x, y, z ) =/integraldisplayx2 x1dx/integraldisplayy2(x) y1(x)dy/integraldisplayz2(x,y ) z1(x,y )dz f (x, y, z )(4.6.2 ) For example, a two-dimensional integral over a circle of radius one centered on the origin becomes /integraldisplay1 −1dx/integraldisplay√ 1−x2 −√ 1−x2dy f (x, y )( 4.6.3 ) Now we can define a function G(x, y )that does the innermost integral, G(x, y )≡/integraldisplayz2(x,y ) z1(x,y )f(x, y, z )dz (4.6.4 ) and a function H(x)that does the integral of G(x, y ), H(x)≡/integraldisplayy2(x) y1(x)G(x, y )dy (4.6.5 ) 4.6MultidimensionalIntegrals 157Sample 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).inner integration y x outer integration Figure 4.6.1. Function evaluations for a two-dimensional integral over an irregular region, shown schematically. The outer integration routine, in y, requests values of the inner, x, integral at locations along the yaxis of its own choosing. The inner integration routine then evaluates the function at xlocations suitable to it. This is more accurate in general than, e.g., evaluating the function on a Cartesian mesh of points. andfinally our answer as an integral over H(x) I=/integraldisplayx2 x1H(x)dx (4.6.6 ) To implement equations (4.6.4) –(4.6.6)in a program, one needs three separate copies ofa basic one-dimensionalintegrationroutine(andof anysubroutinescalled by it), one each for the x,y, and zintegrations. If you try to make do with only one copy, then it will call itself recursively, since (e.g.) the function evaluationsofHfor the xintegration will themselves call the integration routine to do the y integration(see Figure4.6.1). Inourexample,let us supposethatwe planto usethe one-dimensionalintegrator qgausof§4.5. Thenwemakethreeidenticalcopiesand call them qgausx,qgausy, and qgausz. The basic program for three-dimensional integration then is as follows: SUBROUTINE quad3d(x1,x2,ss) REAL ss,x1,x2,h EXTERNAL h C USES h,qgausx Returns as ssthe integral of a user-supplied function func over a three-dimensional region specified by the limits x1,x2, and by the user-supplied functions y1,y2,z1,a n d z2,a s defined in (4.6.2). call qgausx(h,x1,x2,ss)return END FUNCTION f(zz) REAL f,zz,func,x,y,z 158 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).COMMON /xyz/ x,y,z C USES func Called by qgausz . Calls func . z=zzf=func(x,y,z) return END FUNCTION g(yy) REAL g,yy,f,z1,z2,x,y,z EXTERNAL fCOMMON /xyz/ x,y,z C USES f,qgausz,z1,z2 Called by qgausy . Calls qgausz . REAL ssy=yy call qgausz(f,z1(x,y),z2(x,y),ss) g=ssreturnEND FUNCTION h(xx) REAL h,xx,g,y1,y2,x,y,z EXTERNAL g COMMON /xyz/ x,y,z C USES g,qgausy,y1,y2 Called by qgausx . Calls qgausy . REAL ss x=xxcall qgausy(g,y1(x),y2(x),ss) h=ss returnEND The necessaryuser-suppliedfunctionshavethe followingcalling sequences: FUNCTION func(x,y,z) The 3-dimensional function to be integrated FUNCTION y1(x) FUNCTION y2(x)FUNCTION z1(x,y)FUNCTION z2(x,y) CITED REFERENCES AND FURTHER READING: Stroud,A.H.1971, ApproximateCalculationofMultipleIntegrals (EnglewoodCliffs,NJ:Prentice- Hall). [1] Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), §7.7, p. 318. Johnson, L.W., and Riess, R.D. 1982, Numerical Analysis , 2nd ed. (Reading, MA: Addison- Wesley), §6.2.5, p. 307. 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), equations 25.4.58ff.