Page - 166 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python 3.6, Volume Second Edition

166 6 ComputingIntegralsandTestingCode f_mean = f_mean + f(x[i], y[j]) f_mean = f_mean/num_inside area = num_inside/(n**2)*(x1 - x0)*(y1 - y0) return area*f_mean (See thefileMC_double.py.) Verification A simple 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, but Monte Carlo integration is, unfortunately,never exact so it is impossible to predict the output of the algorithm. All we know is that the estimated integral should approach 3 as the numberofrandompointsgoes to infinity.Also, forafixednumberofpoints,wecan run the algorithm several times and get different numbers that fluctuate around the exact value, since different sample points are used in different calls to the Monte Carlo integrationalgorithm. The area of the rectangle can be computed by the integral ∫2 0 ∫4.5 3 dydx, so in this case we identify f(x,y) = 1, and the g 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, using samplesofdifferentsizes,onhowwecanutilize theMonteCarlo_doublefunction tocompute thearea: In [1]: from MC_double import MonteCarlo_double In [2]: def g(x, y): ...: return (1 if (0 <= x <= 2 and 3 <= y <= 4.5) else -1) ...: In [3]: MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 100) Out[3]: 2.9484 In [4]: MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 1000) Out[4]: 2.947032 In [5]: MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 1000) Out[5]: 3.0234600000000005 In [6]: MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 2000) Out[6]: 2.9984580000000003 In [7]: MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 2000) Out[7]: 3.1903469999999996 In [8]: MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 5000) Out[8]: 2.986515 Note the compact if-else construction in the definition of g. It is a one-line alternative to, forexample, if (0 <= x <= 2) and (3 <= y <= 4.5): return 1 else: return -1 From the output, we see that the values fluctuate around 3, a fact that supports a correct implementation,but inprinciple,bugscouldbehiddenbehindtheinaccurate answers.