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

3.7 DoubleandTriple Integrals 87 kind of convergence rate estimate could be used to verify the implementation, but this topic isbeyondthe scopeof thisbook. Test function for functionwith randomnumbers Tomake a test function, we need a unit test that has identical behavior each timewe run the test. This seems difficultwhen randomnumbers are involved, because these numbers are different every timewerun thealgorithm,andeach runhenceproducesa (slightly)different result. A standardway to test algorithms involving randomnumbers is to fix the seed of the randomnumber generator. Then the sequence of numbers is the same every timewerunthealgorithm.Assuming that theMonteCarlo_doublefunction works, wefix the seed, observe a certain result, and take this result as the correct result. Provided the test functionalways uses this seed,we shouldget exactly this resultevery timetheMonteCarlo_doublefunctioniscalled.Our test functioncan thenbewrittenas 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 thefileMC_double.py.) Integral overacircle The test above involves a trivial functionf.x;y/D 1. We shouldalso test anon-constantf functionandamorecomplicateddomain. Let˝ beacircleat theoriginwith radius2, and letf Dpx2Cy2. Thischoicemakes it possible to compute anexact result: in polar coordinates, R ˝ f.x;y/dxdy simpli- fies to2 R2 0 r2dr D 16 =3.Wemustbeprepared forquite crudeapproximations thatfluctuatearound this exact result. As in the test case above,weexperiencebet- ter resultswith largernumberofpoints.Whenwehavesuchevidenceforaworking implementation,wecanturnthe test intoapropertest 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 r = sympy.symbols(’r’) I_exact = sympy.integrate(2*sympy.pi*r*r, (r, 0, 2))