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

6.7 DoubleandTriple Integrals 167 It is mathematically known that the standard deviation of the Monte Carlo estimate of an integral convergesasn−1/2, wheren is the numberof samples. This kind of convergence rate estimate could be used to verify the implementation, but the topic isbeyondthescopeof thisbook. Test FunctionforFunctionwith RandomNumbers To makea test function,we need a unit test that has identical behavior each time we run the test. Thus, since the algorithm generates pseudo-randomnumbers, we apply the standard technique of fixing the seed (of the random number generator), so that the sequence of numbersgenerated is the same every time we run the algorithm.Assuming that the MonteCarlo_doublefunctionworks,wefix theseed,observea certain result, and take this resultas thecorrectresult.Providedthe test functionalwaysuses thisseed, we should get exactly this result every time theMonteCarlo_double function is called. Of course, this procedure does not test whether the MonteCarlo_double function works right now (but we hope our assumption of correctness is well founded!).Still, it isneverthelessusefulwhenfuturechangesaremade,sinceatany time we can confirm that MonteCarlo_doublegives the same answer as before. The test functioncanbewrittenas shownbelow. def test_MonteCarlo_double_rectangle_area(): """Check the area of a rectangle.""" def g(x, y): return (1 if (0 <= x <= 2 and 3 <= y <= 4.5) else -1) x0 = 0; x1 = 3; y0 = 2; y1 = 5 # embedded rectangle n = 1000 np.random.seed(8) # must fix the seed! I_expected = 3.121092 # computed with this seed I_computed = MonteCarlo_double( lambda x, y: 1, g, x0, x1, y0, y1, n) assert abs(I_expected - I_computed) < 1E-14 (See the fileMC_double.py.) IntegralOvera Circle The test above involvesa trivial functionf(x,y)=1.We should also test a non-constantf function and a more complicated domain. LetΩ be a circle at the originwith radius 2, and letf =√x2+y2. This choice makes it possibletocomputeanexactresult: inpolarcoordinates, ∫ Ωf(x,y)dxdy simplifies to 2π ∫2 0 r 2dr = 16π/3. We must be prepared for quite crude approximations that fluctuate around this exact result. As in the test case above, we experience better results with larger number of points. When we have such evidence for a working implementation,wecanturn the test intoaproper test function.Here isanexample: def test_MonteCarlo_double_circle_r(): """Check the integral of r over a circle with radius 2.""" def g(x, y): xc, yc = 0, 0 # center R = 2 # radius return R**2 - ((x-xc)**2 + (y-yc)**2) # Exact: integral of r*r*dr over circle with radius R becomes # 2*pi*1/3*R**3 import sympy