7.6 Simple Monte Carlo Integration 307 #include "nrutil.h" ns... Set to the number of sample points desired. den=... Set to the constant value of the density. Sw=SwX=sWy=sWz=0.0; Zero the various sums to be accumulated. varv=varx=vary=varz=0.0; v01=3.0*7.0*2.0; Volume of the sampled region. for(j=1;j<=m;j++) x=1.0+3.0*ran2(&1dum); Pick a point randomly in the sampled re- y=(-3.0)+7.0*ran2(&idum); gion. z=(-1.0)+2.0*ran2(&idum); 1f(z*z+SQR(sqrt(x*x+y*y)-3.0)<1.0)[ Is it in the torus? sw +den; If so,add to the various cumulants. swx +x*den; swy +=y*den; SwZ +z*den: varw +SQR(den); varx SQR(x*den); 100 vary +SQR(y*den); varz +SQR(z*den); w-vol*sw/n; The values of the integrals (7.6.5). x-vol*swx/n; y=vol*swy/n; (Nort server z=vol*swz/n; dw=vol*sqrt((varw/n-SQR(sw/n))/n); and their corresponding error estimates. computer, Americ Press. dx=vol*sqrt((varx/n-SQR(swx/n))/n); ART dy=vol*sqrt((vary/n-SQR(swy/n))/n); dz=vol*sqrt((varz/n-SQR(swz/n))/n) Programs A change of variable can often be extremely worthwhile in Monte Carlo SCIENTIFIC integration.Suppose,for example,that we want to evaluate the same integrals, but for a piece-of-torus whose density is a strong function of z,in fact varying 6 according to p(z,y,z）=e52 (7.6.6) 1920 COMPUTING (ISBN One way to do this is to put the statement Numerica 10621 den=exp(5.0*z) uctio 43106 inside the if (...block,just before den is first used.This will work,but it is a poor way to proceed.Since(7.6.6)falls so rapidly to zero as z decreases(down (outside Recipes to its lower limit-1).most sampled points contribute almost nothing to the sum North Software. of the weight or moments.These points are effectively wasted,almost as badly as those that fall outside of the region W.A change of variable,exactly as in the transformation methods of $7.2,solves this problem.Let ds=e5:dz so that s=e5,=(5s (7.6.7) Then pdz ds,and the limits-1<z<1 become.00135 <s<29.682.The program fragment now looks like this7.6 Simple Monte Carlo Integration 307 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). #include "nrutil.h" ... n=... Set to the number of sample points desired. den=... Set to the constant value of the density. sw=swx=swy=swz=0.0; Zero the various sums to be accumulated. varw=varx=vary=varz=0.0; vol=3.0*7.0*2.0; Volume of the sampled region. for(j=1;j<=n;j++) { x=1.0+3.0*ran2(&idum); Pick a point randomly in the sampled re￾y=(-3.0)+7.0*ran2(&idum); gion. z=(-1.0)+2.0*ran2(&idum); if (z*z+SQR(sqrt(x*x+y*y)-3.0) < 1.0) { Is it in the torus? sw += den; If so, add to the various cumulants. swx += x*den; swy += y*den; swz += z*den; varw += SQR(den); varx += SQR(x*den); vary += SQR(y*den); varz += SQR(z*den); } } w=vol*sw/n; The values of the integrals (7.6.5), x=vol*swx/n; y=vol*swy/n; z=vol*swz/n; dw=vol*sqrt((varw/n-SQR(sw/n))/n); and their corresponding error estimates. dx=vol*sqrt((varx/n-SQR(swx/n))/n); dy=vol*sqrt((vary/n-SQR(swy/n))/n); dz=vol*sqrt((varz/n-SQR(swz/n))/n); A change of variable can often be extremely worthwhile in Monte Carlo integration. Suppose, for example, that we want to evaluate the same integrals, but for a piece-of-torus whose density is a strong function of z, in fact varying according to ρ(x, y, z) = e5z (7.6.6) One way to do this is to put the statement den=exp(5.0*z); inside the if (...) block, just before den is first used. This will work, but it is a poor way to proceed. Since (7.6.6) falls so rapidly to zero as z decreases (down to its lower limit −1), most sampled points contribute almost nothing to the sum of the weight or moments. These points are effectively wasted, almost as badly as those that fall outside of the region W. A change of variable, exactly as in the transformation methods of §7.2, solves this problem. Let ds = e5zdz so that s = 1 5 e5z, z = 1 5 ln(5s) (7.6.7) Then ρdz = ds, and the limits −1 <z< 1 become .00135 <s< 29.682. The program fragment now looks like this