f4-5

140 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).so that /integraldisplayx=∞ x=af(x)dx=/integraldisplayt=e−a t=0f(−logt)dt t(4.4.8 ) The user-transparent implementation would be SUBROUTINE midexp(funk,aa,bb,s,n) INTEGER n REAL aa,bb,s,funk EXTERNAL funk This routine is an exact replacement for midpnt, except that bbis assumed to be infinite (valuepassednotactually used). Itisassumed thatthefunction funkdecreases exponen- tially rapidly at infinity. INTEGER it,jREAL ddel,del,sum,tnm,x,func,a,b func(x)=funk(-log(x))/x b=exp(-aa)a=0.if (n.eq.1) then The rest of the routine is exactly like midpntand is omitted. CITED REFERENCES AND FURTHER READING: Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), Chapter 4. Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), §7.4.3, p. 294. Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §3.7, p. 152. 4.5 Gaussian Quadratures and Orthogonal Polynomials Intheformulasof §4.1,theintegralofafunctionwas approximatedbythesum of its functional values at a set of equally spaced points, multiplied by certain aptlychosen weighting coefficients. We saw that as we allowed ourselves more freedom in choosing the coefficients, we could achieve integration formulas of higher and higher order. The idea of Gaussian quadratures is to give ourselves the freedom to choose not only the weighting coefficients, but also the location of the abscissas at whichthe functionis to be evaluated: Theywill no longerbe equallyspaced. Thus,wewill have twicethenumberofdegreesoffreedomat ourdisposal;it will turnout that we can achieveGaussian quadratureformulaswhose orderis, essentially, twice that of the Newton-Cotesformulawith the same numberof functionevaluations. Does this sound too good to be true? Well, in a sense it is. The catch is a familiar one, which cannot be overemphasized: High order is not the same as high accuracy. High order translates to high accuracy only when the integrand is very smooth, in the sense of being “well-approximated by a polynomial.” 4.5GaussianQuadraturesandOrthogonalPolynomials 141Sample 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).There is, however,one additional feature of Gaussian quadratureformulas that addstotheirusefulness: Wecanarrangethechoiceofweightsandabscissastomaketheintegralexactforaclass ofintegrands“polynomialstimes someknownfunction W(x)” rather than for the usual class of integrands “polynomials.” The function W(x)canthenbechosentoremoveintegrablesingularitiesfromthedesiredintegral. Given W(x), in other words, and given an integer N, we can find a set of weights w jand abscissas xjsuch that the approximation /integraldisplayb aW(x)f(x)dx≈N/summationdisplay j=1wjf(xj)( 4.5.1 ) is exact if f(x)is a polynomial. For example, to do the integral /integraldisplay1 −1exp(−cos2x)√ 1−x2dx (4.5.2 ) (notaverynaturallookingintegral,itmustbeadmitted),wemightwellbeinterested in a Gaussian quadrature formula based on the choice W(x)=1√ 1−x2(4.5.3 ) intheinterval (−1,1). (Thisparticularchoiceiscalled Gauss-Chebyshevintegration , for reasons that will become clear shortly.) Notice that the integration formula (4.5.1) can also be written with the weight function W(x)notovertlyvisible: Define g(x)≡W(x)f(x)andvj≡wj/W (xj). Then (4.5.1) becomes /integraldisplayb ag(x)dx≈N/summationdisplay j=1vjg(xj)( 4.5.4 ) Where did the function W(x)go? It is lurking there, ready to give high-order accuracytointegrandsoftheformpolynomialstimes W(x),andreadyto denyhigh- order accuracy to integrands that are otherwise perfectly smooth and well-behaved. When youfind tabulations of the weights and abscissas for a given W(x), you have to determine carefully whether they are to be used with a formula in the form of (4.5.1), or like (4.5.4). Hereisanexampleofaquadratureroutinethatcontainsthetabulatedabscissas and weights for the case W(x)=1andN=1 0. Since the weights and abscissas are, in this case, symmetric around the midpoint of the range of integration, thereare actually only five distinct values of each: SUBROUTINE qgaus(func,a,b,ss) REAL a,b,ss,funcEXTERNAL func Returns as ssthe integral of the function funcbetween aandb, by ten-point Gauss- Legendre integration: the function is evaluated exactly ten times at interior points in therange of integration. INTEGER j 142 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).REAL dx,xm,xr,w(5),x(5) The abscissas and weights. SAVE w,x DATA w/.2955242247,.2692667193,.2190863625,.1494513491,.0666713443/ DATA x/.1488743389,.4333953941,.6794095682,.8650633666,.9739065285/xm=0.5*(b+a) xr=0.5*(b-a) ss=0 Willbetwicetheaveragevalueofthefunction,sincetheten weights(fivenumbersaboveeachusedtwice)sumto2. do 11j=1,5 dx=xr*x(j) ss=ss+w(j)*(func(xm+dx)+func(xm-dx)) enddo 11 ss=xr*ss Scale the answer to the range of integration. return END The above routine illustrates that one can use Gaussian quadratures without necessarilyunderstandingthetheorybehindthem: Onejustlocatestabulatedweights and abscissas in a book (e.g., [1]or[2]). However, the theory is very pretty, and it willcomeinhandyifyoueverneedtoconstructyourowntabulationofweightsand abscissasforanunusualchoiceof W(x). Wewillthereforegive,withoutanyproofs, someusefulresultsthatwill enableyoutodothis. Severaloftheresultsassumethat W(x)does not change sign inside (a, b ), which is usually the case in practice. The theory behind Gaussian quadratures goes back to Gauss in 1814, who used continued fractions to develop the subject. In 1826 Jacobi rederived Gauss’s results by means of orthogonal polynomials. The systematic treatment of arbitrary weightfunctions W(x)usingorthogonalpolynomialsislargelyduetoChristoffelin 1877. To introduce these orthogonal polynomials, let us fix the interval of interest to be (a, b ). We can define the “scalar product of two functions fandgover a weight function W”a s /angbracketleftf|g/angbracketright≡/integraldisplayb aW(x)f(x)g(x)dx (4.5.5 ) The scalar product is a number, not a function of x. Two functions are said to be orthogonal if their scalar product is zero. A function is said to be normalized if its scalarproductwithitselfisunity. A setoffunctionsthatareallmutuallyorthogonal and also all individually normalized is called an orthonormal set. We can find a set of polynomials (i) that includes exactly one polynomial of order j, called pj(x), for each j=0,1,2,..., and (ii) all of which are mutually orthogonal over the specified weight function W(x). A constructive procedure for finding such a set is the recurrence relation p−1(x)≡0 p0(x)≡1 pj+1(x)=(x−aj)pj(x)−bjpj−1(x) j=0,1,2,...(4.5.6 ) where aj=/angbracketleftxpj|pj/angbracketright /angbracketleftpj|pj/angbracketrightj=0,1,... bj=/angbracketleftpj|pj/angbracketright /angbracketleftpj−1|pj−1/angbracketrightj=1,2,...(4.5.7 ) 4.5GaussianQuadraturesandOrthogonalPolynomials 143Sample 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).The coefficient b0is arbitrary; we can take it to be zero. The polynomials defined by (4.5.6) are monic, i.e., the coefficient of their leading term [ xjforpj(x)] is unity. If we divide each pj(x)by the constant [/angbracketleftpj|pj/angbracketright]1/2we canrendertheset ofpolynomialsorthonormal. Onealso encounters orthogonal polynomials with various other normalizations. You can convert froma given normalization to monic polynomials if you know that the coefficient of x jinpjisλj, say; then the monic polynomials are obtained by dividing each pj byλj. Note that the coefficients in the recurrence relation (4.5.6) depend on the adopted normalization. The polynomial pj(x)can be shown to have exactly jdistinct roots in the interval (a, b ). Moreover, it can be shown that the roots of pj(x)“interleave” the j−1roots of pj−1(x), i.e., there is exactly one root of the former in between each two adjacent roots of the latter. This fact comes in handy if you need to find all theroots: You can start with the one root of p 1(x)and then, in turn, bracket the roots of each higher j, pinningthem down at each stage more precisely by Newton’s rule or some other root-finding scheme (see Chapter 9). Why would you ever want to find all the roots of an orthogonal polynomial pj(x)? Because the abscissas of the N-point Gaussian quadrature formulas (4.5.1) and(4.5.4)withweightingfunction W(x)intheinterval (a, b )arepreciselytheroots of the orthogonal polynomial pN(x)for the same interval and weighting function. This is the fundamental theorem of Gaussian quadratures, and lets you find theabscissas for any particular case. Once you know the abscissas x 1,...,x N, you need to find the weights wj, j=1,...,N. One way to do this (not the most efficient) is to solve the set of linear equations  p0(x1)... p 0(xN) p1(x1)... p 1(xN) ...... p N−1(x1)... p N−1(xN)  w1 w2 ... wN = /integraltextb aW(x)p0(x)dx 0 ... 0 (4.5.8 ) Equation (4.5.8) simply solves for those weights such that the quadrature (4.5.1) givesthecorrectanswerfortheintegralofthefirst Northogonalpolynomials. Note that the zeros on the right-hand side of (4.5.8) appear because p1(x),...,p N−1(x) are all orthogonal to p0(x), which is a constant. It can be shown that, with those weights,theintegralofthe nextN−1polynomialsisalsoexact,sothatthequadrature is exact for all polynomials of degree 2N−1or less. Another way to evaluate the weights (though one whose proof is beyond our scope) is by the formula wj=/angbracketleftpN−1|pN−1/angbracketright pN−1(xj)p/prime N(xj)(4.5.9 ) where p/prime N(xj)is the derivative of the orthogonalpolynomial at its zero xj. ThecomputationofGaussianquadraturerulesthusinvolvestwodistinctphases: (i)thegenerationoftheorthogonalpolynomials p0,...,p N,i.e.,thecomputationof the coefficients aj,bjin (4.5.6); (ii) the determination of the zeros of pN(x), and thecomputationoftheassociatedweights. Forthecaseofthe“classical”orthogonal polynomials, the coefficients ajandbjare explicitly known (equations 4.5.10 – 144 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).4.5.14 below) and phase (i) can be omitted. However, if you are confronted with a “nonclassical” weight function W(x), and you don’t know the coefficients ajand bj, the construction of the associated set of orthogonal polynomials is not trivial. We discuss it at the end of this section. Computationof the Abscissas and Weights Thistaskcanrangefromeasytodifficult,dependingonhowmuchyoualready know about your weight function and its associated polynomials. In the case of classical, well-studied, orthogonal polynomials, practically everything is known,includinggoodapproximationsfortheirzeros. Thesecanbeusedasstartingguesses, enabling Newton’s method (to be discussed in §9.4) to converge very rapidly. Newton’s method requires the derivative p /prime N(x), which is evaluated by standard relations in terms of pNandpN−1. The weights are then convenientlyevaluatedby equation (4.5.9). For the following named cases, this direct root-finding is faster,by a factor of 3 to 5, than any other method. Here are the weight functions, intervals, and recurrence relations that generate the most commonlyused orthogonalpolynomialsand their correspondingGaussianquadrature formulas. Gauss-Legendre: W(x)=1 −1<x< 1 (j+1 )P j+1=( 2j+1 )xPj−jPj−1 (4.5.10 ) Gauss-Chebyshev: W(x)=( 1 −x2)−1/2−1<x< 1 Tj+1=2xTj−Tj−1 (4.5.11 ) Gauss-Laguerre: W(x)=xαe−x0<x< ∞ (j+1 )Lα j+1=(−x+2j+α+1 )Lα j−(j+α)Lα j−1 (4.5.12 ) Gauss-Hermite: W(x)=e−x2−∞ <x< ∞ Hj+1=2xHj−2jHj−1 (4.5.13 ) Gauss-Jacobi: W(x)=( 1 −x)α(1 +x)β−1<x< 1 4.5GaussianQuadraturesandOrthogonalPolynomials 145Sample 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).cjP(α,β ) j+1 =(dj+ejx)P(α,β ) j −fjP(α,β ) j−1 (4.5.14 ) where the coefficients cj,dj,ej, and fjare given by cj=2 (j+1 ) ( j+α+β+ 1)(2 j+α+β) dj=( 2j+α+β+1 ) ( α2−β2) ej=( 2j+α+β)(2j+α+β+ 1)(2 j+α+β+2 ) fj=2 (j+α)(j+β)(2j+α+β+2 )(4.5.15 ) We now give individual routines that calculate the abscissas and weights for these cases. First comes the most common set of abscissas and weights, those of Gauss-Legendre. The routine, due to G.B. Rybicki, uses equation (4.5.9) in the special form for the Gauss-Legendre case, wj=2 (1−x2 j)[P/prime N(xj)]2(4.5.16 ) Theroutinealsoscalestherangeofintegrationfrom (x1,x2)to(−1,1),andprovides abscissas xjand weights wjfor the Gaussian formula /integraldisplayx2 x1f(x)dx=N/summationdisplay j=1wjf(xj)( 4.5.17 ) SUBROUTINE gauleg(x1,x2,x,w,n) INTEGER nREAL x1,x2,x(n),w(n) DOUBLE PRECISION EPS PARAMETER (EPS=3.d-14) EPS is the relative precision. Giventhelowerandupperlimitsofintegration x1andx2,andgiven n,thisroutinereturns arrays x(1:n)andw(1:n)oflength n,containingtheabscissasandweightsoftheGauss- Legendre n-point quadrature formula. INTEGER i,j,mDOUBLE PRECISION p1,p2,p3,pp,xl,xm,z,z1 High precision is a good idea for this routine. m=(n+1)/2 The roots are symmetric in the interval, so we o n l yh a v et ofi n dh a l fo ft h e m . xm=0.5d0*(x2+x1) xl=0.5d0*(x2-x1) do 12i=1,m Loop over the desired roots. z=cos(3.141592654d0*(i-.25d0)/(n+.5d0)) Starting with the above approximation to the ith root, we enter the main loop of re- finement by Newton’s method. 1 continue p1=1.d0p2=0.d0 do 11j=1,n Loop up the recurrence relation to get the Leg- endre polynomial evaluated at z. p3=p2 p2=p1p1=((2.d0*j-1.d0)*z*p2-(j-1.d0)*p3)/j enddo 11 p1is nowthe desired Legendre polynomial. We nextcompute pp, its derivative, by a standard relation involving also p2, the polynomial of one lower order. pp=n*(z*p1-p2)/(z*z-1.d0) 146 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).z1=z z=z1-p1/pp Newton’s method. if(abs(z-z1).gt.EPS)goto 1 x(i)=xm-xl*z Scale the root to the desired interval, x(n+1-i)=xm+xl*z and put in its symmetric counterpart. w(i)=2.d0*xl/((1.d0-z*z)*pp*pp) Compute the weight w(n+1-i)=w(i) and its symmetric counterpart. enddo 12 return END Next we give three routines that use initial approximations for the roots given by Stroud and Secrest [2]. The first is for Gauss-Laguerre abscissas and weights, to be used with the integration formula /integraldisplay∞ 0xαe−xf(x)dx=N/summationdisplay j=1wjf(xj)( 4.5.18 ) SUBROUTINE gaulag(x,w,n,alf) INTEGER n,MAXITREAL alf,w(n),x(n)DOUBLE PRECISION EPS PARAMETER (EPS=3.D-14,MAXIT=10) Increase EPSifyou don’t havethis precision. C USES gammln Given alf,theparameter αoftheLaguerrepolynomials,thisroutinereturnsarrays x(1:n) andw(1:n)containingtheabscissasandweightsofthe n-pointGauss-Laguerrequadrature formula. The smallest abscissa is returned in x(1),t h el a r g e s ti n x(n). INTEGER i,its,jREAL ai,gammln DOUBLE PRECISION p1,p2,p3,pp,z,z1 High precision is a good idea for this routine. do 13i=1,n Loop over the desired roots. if(i.eq.1)then Initial guess for the smallest root. z=(1.+alf)*(3.+.92*alf)/(1.+2.4*n+1.8*alf) else if(i.eq.2)then Initial guess for the second root. z=z+(15.+6.25*alf)/(1.+.9*alf+2.5*n) else Initial guess for the other roots. ai=i-2z=z+((1.+2.55*ai)/(1.9*ai)+1.26*ai*alf/ * (1.+3.5*ai))*(z-x(i-2))/(1.+.3*alf) endif do 12its=1,MAXIT Refinement by Newton’s method. p1=1.d0 p2=0.d0 do11j=1,n Loop up the recurrence relation to get the Laguerre polynomial evaluated at z. p3=p2 p2=p1 p1=((2*j-1+alf-z)*p2-(j-1+alf)*p3)/j enddo 11 p1is nowthe desired Laguerre polynomial. We next compute pp, its derivative, by a standard relation involving also p2, the polynomial of one lower order. pp=(n*p1-(n+alf)*p2)/zz1=zz=z1-p1/pp Newton’s formula. if(abs(z-z1).le.EPS)goto 1 enddo 12 pause ’too many iterations in gaulag’ 1 x(i)=z Store the root and the weight. 4.5GaussianQuadraturesandOrthogonalPolynomials 147Sample 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).w(i)=-exp(gammln(alf+n)-gammln(float(n)))/(pp*n*p2) enddo 13 returnEND Next is a routine for Gauss-Hermite abscissas and weights. If we use the “standard” normalization of these functions, as given in equation (4.5.13), we findthat the computations overflow for large Nbecause of various factorials that occur. We can avoid this by using instead the orthonormal set of polynomials /tildewideH j. They are generated by the recurrence /tildewideH−1=0,/tildewideH0=1 π1/4,/tildewideHj+1=x/radicalbigg2 j+1/tildewideHj−/radicalBigg j j+1/tildewideHj−1 (4.5.19 ) The formula for the weights becomes wj=2 [/tildewideH/prime N(xj)]2(4.5.20 ) while the formula for the derivative with this normalization is /tildewideH/prime j=/radicalbig 2j/tildewideHj−1 (4.5.21 ) Theabscissasandweightsreturnedby gauherareusedwiththeintegrationformula /integraldisplay∞ −∞e−x2f(x)dx=N/summationdisplay j=1wjf(xj)( 4.5.22 ) SUBROUTINE gauher(x,w,n) INTEGER n,MAXITREAL w(n),x(n) DOUBLE PRECISION EPS,PIM4 PARAMETER (EPS=3.D-14,PIM4=.7511255444649425D0,MAXIT=10) Given n, this routine returns arrays x(1:n)andw(1:n)containing the abscissas and weights ofthe n-pointGauss-Hermitequadrature formula. The largestabscissa isreturned inx(1), the most negative in x(n). Parameters: EPSisthe relativeprecision, PIM4 =1/π1/4,MAXIT =maximumiterations. INTEGER i,its,j,m DOUBLE PRECISION p1,p2,p3,pp,z,z1 High precision is a good idea for this routine. m=(n+1)/2 The roots are symmetric about the origin, so we have to find only half of them. do13i=1,m Loop over the desired roots. if(i.eq.1)then Initial guess for the largest root. z=sqrt(float(2*n+1))-1.85575*(2*n+1)**(-.16667) else if(i.eq.2)then Initial guess for the second largest root. z=z-1.14*n**.426/z else if (i.eq.3)then Initial guess for the third largest root. z=1.86*z-.86*x(1) else if (i.eq.4)then Initial guess for the fourth largest root. z=1.91*z-.91*x(2) else Initial guess for the other roots. z=2.*z-x(i-2) 148 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).endif do12its=1,MAXIT Refinement by Newton’s method. p1=PIM4 p2=0.d0do 11j=1,n Loopuptherecurrence relationtogettheHermitepoly- nomial evaluated at z. p3=p2 p2=p1p1=z*sqrt(2.d0/j)*p2-sqrt(dble(j-1)/dble(j))*p3 enddo 11 p1is nowthe desired Hermite polynomial. We next compute pp, its derivative, by the relation (4.5.21) using p2, the polynomial of one lower order. pp=sqrt(2.d0*n)*p2z1=z z=z1-p1/pp Newton’s formula. if(abs(z-z1).le.EPS)goto 1 enddo 12 pause ’too many iterations in gauher’ 1 x(i)=z Store the root x(n+1-i)=-z and its symmetric counterpart. w(i)=2.d0/(pp*pp) Compute the weight w(n+1-i)=w(i) and its symmetric counterpart. enddo 13 return END Finally, here is a routine for Gauss-Jacobi abscissas and weights, which implement the integration formula /integraldisplay1 −1(1−x)α(1 +x)βf(x)dx=N/summationdisplay j=1wjf(xj)( 4.5.23 ) SUBROUTINE gaujac(x,w,n,alf,bet) INTEGER n,MAXIT REAL alf,bet,w(n),x(n) DOUBLE PRECISION EPSPARAMETER (EPS=3.D-14,MAXIT=10) Increase EPSifyou don’t havethis precision. C USES gammln Given alfandbet,theparameters αandβoftheJacobipolynomials,thisroutinereturns arrays x(1:n)andw(1:n)containingtheabscissasandweightsofthe n-pointGauss-Jacobi quadrature formula. The largest abscissa is returned in x(1), the smallest in x(n). INTEGER i,its,j REAL alfbet,an,bn,r1,r2,r3,gammln DOUBLE PRECISION a,b,c,p1,p2,p3,pp,temp,z,z1 High precision is a good idea for this routine. do13i=1,n Loop over the desired roots. if(i.eq.1)then Initial guess for the largest root. an=alf/nbn=bet/n r1=(1.+alf)*(2.78/(4.+n*n)+.768*an/n) r2=1.+1.48*an+.96*bn+.452*an*an+.83*an*bnz=1.-r1/r2 else if(i.eq.2)then Initial guess for the second largest root. r1=(4.1+alf)/((1.+alf)*(1.+.156*alf))r2=1.+.06*(n-8.)*(1.+.12*alf)/nr3=1.+.012*bet*(1.+.25*abs(alf))/n z=z-(1.-z)*r1*r2*r3 else if(i.eq.3)then Initial guess for the third largest root. r1=(1.67+.28*alf)/(1.+.37*alf) r2=1.+.22*(n-8.)/n 4.5GaussianQuadraturesandOrthogonalPolynomials 149Sample 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).r3=1.+8.*bet/((6.28+bet)*n*n) z=z-(x(1)-z)*r1*r2*r3 else if(i.eq.n-1)then Initial guess forthe second smallest root. r1=(1.+.235*bet)/(.766+.119*bet)r2=1./(1.+.639*(n-4.)/(1.+.71*(n-4.))) r3=1./(1.+20.*alf/((7.5+alf)*n*n)) z=z+(z-x(n-3))*r1*r2*r3 else if(i.eq.n)then Initial guess for the smallest root. r1=(1.+.37*bet)/(1.67+.28*bet) r2=1./(1.+.22*(n-8.)/n) r3=1./(1.+8.*alf/((6.28+alf)*n*n))z=z+(z-x(n-2))*r1*r2*r3 else Initial guess for the other roots. z=3.*x(i-1)-3.*x(i-2)+x(i-3) endifalfbet=alf+bet do 12its=1,MAXIT Refinement by Newton’s method. temp=2.d0+alfbet Start therecurrence with P0andP1to avoidadivi- sion by zero when α+β=0or−1. p1=(alf-bet+temp*z)/2.d0 p2=1.d0 do11j=2,n Loop up the recurrence relation to get the Jacobi polynomial evaluated at z. p3=p2 p2=p1 temp=2*j+alfbet a=2*j*(j+alfbet)*(temp-2.d0)b=(temp-1.d0)*(alf*alf-bet*bet+temp* * (temp-2.d0)*z) c=2.d0*(j-1+alf)*(j-1+bet)*temp p1=(b*p2-c*p3)/a enddo 11 pp=(n*(alf-bet-temp*z)*p1+2.d0*(n+alf)* * (n+bet)*p2)/(temp*(1.d0-z*z)) p1is nowthe desired Jacobi polynomial. We next compute pp, its derivative, by a standard relation involving also p2, the polynomial of one lower order. z1=z z=z1-p1/pp Newton’s formula. if(abs(z-z1).le.EPS)goto 1 enddo 12 pause ’too many iterations in gaujac’ 1 x(i)=z Store the root and the weight. w(i)=exp(gammln(alf+n)+gammln(bet+n)-gammln(n+1.)- * gammln(n+alfbet+1.))*temp*2.**alfbet/(pp*p2) enddo 13 returnEND LegendrepolynomialsarespecialcasesofJacobipolynomialswith α=β=0, butitisworthhavingtheseparateroutineforthem, gauleg,givenabove. Chebyshev polynomialscorrespondto α=β=−1/2(see§5.8). Theyhave analytic abscissas and weights: xj=c o s/parenleftbiggπ(j−1 2) N/parenrightbigg wj=π N(4.5.24 ) 150 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).Case of KnownRecurrences Turn now to the case where you do not know good initial guesses for the zeros of your orthogonal polynomials, but you do have available the coefficients ajandbjthat generate them. As we have seen, the zeros of pN(x)are the abscissas for the N-point Gaussian quadrature formula. The most useful computational formula for the weights is equation(4.5.9) above, since the derivative p /prime Ncan be efficiently computed by the derivative of (4.5.6) in the general case, or by special relations for the classical polynomials. Note that (4.5.9) isvalid as written only for monic polynomials; for other normalizations, there is an extra factor ofλ N/λ N−1, where λNis the coefficient of xNinpN. Except in those special cases already discussed, the best way to find the abscissas is not to use a root-finding method like Newton’s method on pN(x). Rather, it is generally faster to use the Golub-Welsch [3]algorithm, which is based on a result of Wilf [4]. This algorithm notes that if you bring the term xp jto the left-hand side of (4.5.6) and the term pj+1to the right-hand side, the recurrence relation can be written in matrix form as x p 0 p1 ... pN−2 pN−1 = a 0 1 b1a1 1 ...... b N−2aN−2 1 bN−1aN−1 · p 0 p1 ... pN−2 pN−1 + 0 0 ... 0 p N  or xp=T·p+p NeN−1 (4.5.25 ) HereTis a tridiagonal matrix, pis a column vector of p0,p1,...,p N−1, andeN−1is a unit vector with a 1 in the (N−1)st (last) position and zeros elsewhere. The matrix Tcan be symmetrized by a diagonal similarity transformation Dto give J=DTD−1= a0√b1 √b1a1√b2 ......√bN−2aN−2√bN−1 √bN−1aN−1 (4.5.26 ) The matrix Jis called the Jacobi matrix (not to be confused with other matrices named after Jacobi that arise in completely different problems!). Now we see from (4.5.25) thatp N(xj)=0is equivalent to xjbeing an eigenvalue of T. Since eigenvalues are preserved by a similarity transformation, xjis an eigenvalue of the symmetric tridiagonal matrix J. Moreover, Wilf [4]shows that if vjis the eigenvector corresponding to the eigenvalue xj, normalized so that v·v=1, then wj=µ0v2 j,1 (4.5.27 ) where µ0=/integraldisplayb aW(x)dx (4.5.28 ) and where vj,1is the first component of v. As we shall see in Chapter 11, finding all eigenvalues and eigenvectors of a symmetric tridiagonal matrix is a relatively efficient andwell-conditioned procedure. Weaccordingly give aroutine, gaucof,forfinding the abscissas and weights, given the coefficients a jandbj. Remember that if you know the recurrence relationfororthogonalpolynomialsthatarenotnormalizedtobemonic,youcaneasilyconvertit to monic form by means of the quantities λ j. 4.5GaussianQuadraturesandOrthogonalPolynomials 151Sample 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).SUBROUTINE gaucof(n,a,b,amu0,x,w) INTEGER n,NMAX REAL amu0,a(n),b(n),w(n),x(n) PARAMETER (NMAX=64) C USES eigsrt,tqli Computes the abscissas and weights for a Gaussian quadrature formula from the Jacobimatrix. On input, a(1:n)andb(1:n)are the coefficients of the recurrence relation for thesetofmonicorthogonal polynomials. Thequantity µ0≡/integraltextb aW(x)dxisinputas amu0. The abscissas x(1:n)are returned in descending order, with the corresponding weights inw(1:n). The arrays aandbare modified. Execution can be speeded up by modifying tqliandeigsrtto compute only the first component of each eigenvector. INTEGER i,j REAL z(NMAX,NMAX) do12i=1,n if(i.ne.1)b(i)=sqrt(b(i)) Set up superdiagonal of Jacobi matrix. do11j=1,n Setupidentitymatrixfor tqlitocomputeeigenvectors. if(i.eq.j)then z(i,j)=1. else z(i,j)=0. endif enddo 11 enddo 12 call tqli(a,b,n,NMAX,z)call eigsrt(a,z,n,NMAX) Sort eigenvalues into descending order. do 13i=1,n x(i)=a(i) w(i)=amu0*z(1,i)**2 Equation (4.5.12). enddo 13 return END OrthogonalPolynomialswithNonclassical Weights This somewhat specialized subsection will tell you what to do if your weight function is not one of the classical ones dealt with above and you do not know the aj’s and bj’s of the recurrence relation (4.5.6) to use in gaucof. Then, a method of finding the aj’s andbj’s is needed. Theprocedure of Stieltjes is to compute a0from (4.5.7), then p1(x)from (4.5.6). Knowing p0andp1, we can compute a1andb1from (4.5.7), and so on. But how are we to compute the inner products in (4.5.7)? The textbook approach is to represent each pj(x)explicitly as a polynomial in xand to compute the inner products by multiplying out term by term. This will be feasible if weknow the first 2Nmoments of the weight function, µ j=/integraldisplayb axjW(x)dx j =0,1,..., 2N−1( 4.5.29 ) However, the solution of the resulting set of algebraic equations for the coefficients ajandbj in terms of the moments µjis in general extremely ill-conditioned. Even in double precision, it is not unusual to lose all accuracy by the time N=1 2. We thus reject any procedure based on the moments (4.5.29). Sack and Donovan [5]discovered that the numerical stability is greatly improved if, instead of using powers of xas a set of basis functions to represent the pj’s, one uses some other known set of orthogonal polynomials πj(x), say. Roughly speaking, the improved stability occurs because the polynomial basis “samples” the interval (a, b )better than the power basis when the inner product integrals are evaluated, especially if its weight functionresembles W(x). 152 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).So assume that we know the modified moments νj=/integraldisplayb aπj(x)W(x)dx j =0,1,..., 2N−1( 4.5.30 ) where the πj’s satisfy a recurrence relation analogous to (4.5.6), π−1(x)≡0 π0(x)≡1 πj+1(x)=( x−αj)πj(x)−βjπj−1(x) j=0,1,2,...(4.5.31 ) and the coefficients αj,βjare known explicitly. Then Wheeler [6]has given an efficient O(N2)algorithm equivalent to that of Sack and Donovan for finding ajandbjvia a set of intermediate quantities σk,l=/angbracketleftpk|πl/angbracketright k, l≥− 1( 4.5.32 ) Initialize σ−1,l=0 l=1,2,..., 2N−2 σ0,l=νl l=0,1,..., 2N−1 a0=α0+ν1 ν0 b0=0(4.5.33 ) Then, for k=1,2,...,N −1, compute σk,l=σk−1,l+1−(ak−1−αl)σk−1,l−bk−1σk−2,l+βlσk−1,l−1 l=k, k +1,..., 2N−k−1 ak=αk−σk−1,k σk−1,k−1+σk,k +1 σk,k bk=σk,k σk−1,k−1 (4.5.34 ) Note that the normalization factors can also easily be computed if needed: /angbracketleftp0|p0/angbracketright=ν0 /angbracketleftpj|pj/angbracketright=bj/angbracketleftpj−1|pj−1/angbracketright j=1,2,...(4.5.35 ) You can find a derivation of the above algorithm in Ref. [7]. Wheeler’salgorithmrequiresthatthemodifiedmoments(4.5.30)beaccuratelycomputed. In practical cases there is often a closed form, or else recurrence relations can be used. Thealgorithmisextremelysuccessfulfor finiteintervals (a, b ). Forinfiniteintervals,thealgorithm does not completely remove the ill-conditioning. In this case, Gautschi [8,9]recommends reducing the interval to a finite interval by a change of variable, and then using a suitablediscretization procedure to compute the inner products. You will have to consult thereferences for details. We give the routine orthogfor generating the coefficients a jandbjby Wheeler’s algorithm, given the coefficients αjandβj, and the modified moments νj. To conform to the usual FORTRAN convention for dimensioning subscripts, the indices of the σmatrix are increased by 2, i.e., sig(k,l) =σk−2,l−2, while the indices of the vectors α,β,aand bare increased by 1. 4.5GaussianQuadraturesandOrthogonalPolynomials 153Sample 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).SUBROUTINE orthog(n,anu,alpha,beta,a,b) INTEGER n,NMAX REAL a(n),alpha(2*n-1),anu(2*n),b(n),beta(2*n-1) PARAMETER (NMAX=64) Computes the coefficients ajandbj,j =0,...N −1, of the recurrence relation for monicorthogonalpolynomialswithweightfunction W(x)byWheeler’salgorithm. Oninput, alpha(1:2*n-1) andbeta(1:2*n-1) arethecoefficients αjandβj,j =0,... 2N−2, of the recurrence relation for the chosen basis of orthogonal polynomials. The modifiedmoments ν jareinputin anu(1:2*n) .T h efi r s t ncoefficients arereturnedin a(1:n)and b(1:n). INTEGER k,lREAL sig(2*NMAX+1,2*NMAX+1)do 11l=3,2*n Initialization, Equation (4.5.33). sig(1,l)=0. enddo 11 do12l=2,2*n+1 sig(2,l)=anu(l-1) enddo 12 a(1)=alpha(1)+anu(2)/anu(1)b(1)=0. do 14k=3,n+1 Equation (4.5.34). do13l=k,2*n-k+3 sig(k,l)=sig(k-1,l+1)+(alpha(l-1)-a(k-2))*sig(k-1,l)- * b(k-2)*sig(k-2,l)+beta(l-1)*sig(k-1,l-1) enddo 13 a(k-1)=alpha(k-1)+sig(k,k+1)/sig(k,k)-sig(k-1,k)/sig(k-1,k-1)b(k-1)=sig(k,k)/sig(k-1,k-1) enddo 14 return END As an example of the use of orthog, consider the problem [7]of generating orthogonal polynomials with the weight function W(x)=−logxon the interval (0,1). A suitable set ofπj’s is the shifted Legendre polynomials πj=(j!)2 (2j)!Pj(2x−1) ( 4.5.36 ) The factor in front of Pjmakes the polynomials monic. The coefficients in the recurrence relation (4.5.31) are αj=1 2j=0,1,... βj=1 4(4−j−2)j=1,2,...(4.5.37 ) while the modified moments are νj=  1 j=0 (−1)j(j!)2 j(j+ 1)(2 j)!j≥1(4.5.38 ) A call to orthogwith this input allows one to generate the required polynomials to machine accuracyforverylarge N,andhencedoGaussianquadraturewiththisweightfunction. Before Sack and Donovan’s observation, this seemingly simple problem was essentially intractable. Extensions of Gaussian Quadrature There are many differentways in which the ideas of Gaussian quadraturehave been extended. One important extension is the case of preassigned nodes : Some pointsarerequiredtobeincludedinthesetofabscissas,andtheproblemistochoose 154 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).the weights and the remainingabscissas to maximize the degree of exactness of the the quadrature rule. The most common cases are Gauss-Radau quadrature, where one of the nodes is an endpoint of the interval, either aorb, andGauss-Lobatto quadrature,whereboth aandbarenodes. Golub [10]hasgivenanalgorithmsimilar togaucoffor these cases. The second important extension is the Gauss-Kronrod formulas. For ordinary Gaussian quadrature formulas, as Nincreases the sets of abscissas have no points in common. This means that if you compare results with increasing Nas a way of estimating the quadratureerror, you cannot reuse the previous function evaluations. Kronrod [11]posed the problem of searching for optimal sequences of rules, each of which reuses all abscissas of its predecessor. If one starts with N=m, say, and then adds nnew points, one has 2n+mfree parameters: the nnew abscissas and weights, and mnew weights for the fixed previous abscissas. The maximum degree of exactness one would expect to achieve would therefore be 2n+m−1. Thequestionis whetherthismaximumdegreeofexactnesscanactuallybeachieved in practice, when the abscissas are required to all lie inside (a, b ). The answer to this question is not known in general. Kronrod showed that if you choose n=m+1, an optimal extension can be found for Gauss-Legendre quadrature. Patterson [12]showed how to compute continued extensions of this kind. Sequences such as N=1 0 ,21,43,87,...are popularinautomaticquadratureroutines [13]thatattempttointegrateafunctionuntil some specified accuracy has been achieved. 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] Stroud, A.H., and Secrest, D. 1966, Gaussian Quadrature Formulas (Englewood Cliffs, NJ: Prentice-Hall). [2] Golub, G.H., and Welsch, J.H. 1969, Mathematics of Computation , vol. 23, pp. 221–230 and A1–A10. [3] Wilf, H.S. 1962, Mathematics for thePhysical Sciences (NewYork: Wiley),Problem9, p. 80.[4] Sack, R.A., and Donovan, A.F. 1971/72, Numerische Mathematik , vol. 18, pp. 465–478. [5] Wheeler, J.C. 1974, Rocky Mountain Journal of Mathematics , vol. 4, pp. 287–296. [6] Gautschi,W.1978,in RecentAdvancesinNumericalAnalysis ,C.deBoorandG.H.Golub,eds. (New York: Academic Press), pp. 45–72. [7] Gautschi, W.1981,in E.B. Christoffel , P.L. Butzer andF.Feh´ er,eds. (Basel: BirkhauserVerlag), pp. 72–147. [8] Gautschi, W.1990,in OrthogonalPolynomials ,P. Nevai, ed.(Dordrecht: Kluwer AcademicPub- lishers), pp. 181–216. [9] Golub, G.H. 1973, SIAM Review , vol. 15, pp. 318–334. [10] Kronrod, A.S. 1964, Doklady AkademiiNauk SSSR , vol. 154, pp. 283–286 (inRussian). [11] Patterson, T.N.L. 1968, Mathematics of Computation , vol. 22, pp. 847–856 and C1–C11; 1969, op. cit., vol. 23, p. 892. [12] Piessens, R., de Doncker, E., Uberhuber, C.W., and Kahaner, D.K. 1983, QUADPACK: A Sub- routine Package for Automatic Integration (New York: Springer-Verlag). [13] Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §3.6. Johnson, L.W., and Riess, R.D. 1982, Numerical Analysis , 2nd ed. (Reading, MA: Addison- Wesley), §6.5. 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