6.7 Double and Triple Integrals
167
It is mathematically known that the standard deviation of the Monte Carlo
estimate of an integral converges as n −1/2 , where n is the number of samples. This
kind of convergence rate estimate could be used to verify the implementation, but
the topic is beyond the scope of this book.
Test Function for Function with Random Numbers To make a test function, we
need a unit test that has identical behavior each time we run the test. Thus, since
the algorithm generates pseudo-random numbers, we apply the standard technique
of fixing the seed (of the random number generator), so that the sequence of
numbers generated is the same every time we run the algorithm. Assuming that the
MonteCarlo_double function works, we fix the seed, observe a certain result, and
take this result as the correct result. Provided the test function always uses this seed,
we should get exactly this result every time the MonteCarlo_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 is nevertheless useful when future changes are made, since at any
time we can confirm that MonteCarlo_double gives the same answer as before.
The test function can be written as shown below.
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 file MC_double.py.)
Integral Over a Circle The test above involves a trivial function f (x, y) = 1. We
should also test a non-constant f function and a more complicated domain. Let Ω
be a circle at the origin with radius 2, and let f =
x 2 + y 2 . This choice makes it
possible to compute an exact result: in polar coordinates,
Ω f (x, y)dxdy simplifies
to 2π
2
0 r 2 dr = 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, we can turn the test into a proper test function. Here is an example:
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
167
It is mathematically known that the standard deviation of the Monte Carlo
estimate of an integral converges as n −1/2 , where n is the number of samples. This
kind of convergence rate estimate could be used to verify the implementation, but
the topic is beyond the scope of this book.
Test Function for Function with Random Numbers To make a test function, we
need a unit test that has identical behavior each time we run the test. Thus, since
the algorithm generates pseudo-random numbers, we apply the standard technique
of fixing the seed (of the random number generator), so that the sequence of
numbers generated is the same every time we run the algorithm. Assuming that the
MonteCarlo_double function works, we fix the seed, observe a certain result, and
take this result as the correct result. Provided the test function always uses this seed,
we should get exactly this result every time the MonteCarlo_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 is nevertheless useful when future changes are made, since at any
time we can confirm that MonteCarlo_double gives the same answer as before.
The test function can be written as shown below.
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 file MC_double.py.)
Integral Over a Circle The test above involves a trivial function f (x, y) = 1. We
should also test a non-constant f function and a more complicated domain. Let Ω
be a circle at the origin with radius 2, and let f =
x 2 + y 2 . This choice makes it
possible to compute an exact result: in polar coordinates,
Ω f (x, y)dxdy simplifies
to 2π
2
0 r 2 dr = 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, we can turn the test into a proper test function. Here is an example:
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