6.7 Double and Triple Integrals
165
3. count the fraction q of points that are inside Ω
4. approximate A(Ω)/A(R) by q, i.e., set A(Ω) = qA(R)
5. evaluate the mean of f , ¯
f , at the points inside Ω
6. estimate the integral as A(Ω) ¯
f
Note that A(R) is trivial to compute since R is a rectangle, while A(Ω) is
unknown. However, if we assume that the fraction of A(R) occupied by A(Ω) is
the same as the fraction of random points inside Ω, we get a simple estimate for
A(Ω).
To get an idea of the method, consider a circular domain Ω embedded in a
rectangle as shown below. A collection of random points is illustrated by black
dots.
Implementation A Python function implementing
Ω f (x, y)dxdy can be written
like this:
import numpy as np
def MonteCarlo_double(f, g, x0, x1, y0, y1, n):
"""
Monte Carlo integration of f over a domain g>=0, embedded
in a rectangle [x0,x1]x[y0,y1]. n^2 is the number of
random points.
"""
# 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 = num_inside + 1