Page - 86 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python

86 3 Computing Integrals # Draw n**2 random points in the rectangle x = np.random.uniform(x0, x1, n) y = np.random.uniform(y0, y1, n) # Compute sum of f values inside the integration domain f_mean = 0 num_inside = 0 # number of x,y points inside domain (g>=0) for i in range(len(x)): for j in range(len(y)): if g(x[i], y[j]) >= 0: num_inside += 1 f_mean += f(x[i], y[j]) f_mean = f_mean/float(num_inside) area = num_inside/float(n**2)*(x1 - x0)*(y1 - y0) return area*f_mean (See thefileMC_double.py.) Verification Asimple test case is to check the area of a rectangle Œ0;2 Œ3;4:5 embedded in a rectangle Œ0;3 Œ2;5 . The right answer is 3, butMonte Carlo integration is, unfortunately,never exact so it is impossible topredict theoutputof the algorithm. Allweknowis that the estimated integral should approach3 as the number of randompoints goes to infinity. Also, for a fixed number of points, we can run thealgorithmseveral timesandgetdifferentnumbers thatfluctuate around theexactvalue,sincedifferentsamplepointsareusedindifferentcalls to theMonte Carlo integrationalgorithm. The area of the rectangle can be computed by the integral R2 0 R4:5 3 dydx, so in this casewe identifyf.x;y/ D 1, and theg function can be specified as (e.g.) 1 if .x;y/ is inside Œ0;2 Œ3;4:5 and 1 otherwise. Here is an example on how wecanutilize theMonteCarlo_doublefunction to compute the area for different numberof samples: >>> from MC_double import MonteCarlo_double >>> def g(x, y): ... return (1 if (0 <= x <= 2 and 3 <= y <= 4.5) else -1) ... >>> MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 100) 2.9484 >>> MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 1000) 2.947032 >>> MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 1000) 3.0234600000000005 >>> MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 2000) 2.9984580000000003 >>> MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 2000) 3.1903469999999996 >>> MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 5000) 2.986515 Wesee that the valuesfluctuate around3, a fact that supports a correct implemen- tation,but inprinciple,bugscouldbehiddenbehind the inaccurateanswers. It ismathematically known that the standard deviation of theMonte Carlo es- timate of an integral converges asn 1=2, wheren is the number of samples. This