Calcule numérique avec la méthode Monte-Carlo (deuxième partie)

Jérémie Decock (www.jdhp.org)

Attention:

Veuillez noter que ce document est en cours de rédaction. Dans son état actuel, il n'est pas destiné à être lu par d'autres personnes que ses auteurs.

Ce document fait suite à l'article http://www.jdhp.org/docs/notebook/maths_montecarlo_compute_integral_fr.html où le calcul numérique à l'aide de la méthode Monte-Carlo à été introduit.

Application au calcul de la constante $\pi$

...

Exemple written in Python

Random points

In [1]:
N = 100000           # The number of random points

lower_bound = 0.0
upper_bound = 1.0

random_points = np.random.uniform(low=lower_bound, high=upper_bound, size=(N, 2))

Numerical computation of $\pi$ with Monte-Carlo

Points in the circle, i.e. $x^2 + y^2 \leq 1$

...

Compute $\pi$ (with $r=1$):

square_area = $r^2$ circle_area = $\pi * r^2$

quarter_circle_area = $\frac{\pi * r^2}{4}$ <=> quarter_circle_area = $\pi$ square_area / 4 <=> $\pi$ = quarter_circle_area / square_area 4

In [2]:
random_points_in = np.array([p for p in random_points if p[0]**2 + p[1]**2 <= 1.0])
In [3]:
ratio = float(len(random_points_in)) / float(N)
print('Percentage of points in the circle: {:.2f}%'.format(ratio * 100))

s2 = (upper_bound - lower_bound)**2
print("Box surface:", s2)
Percentage of points in the circle: 78.58%
Box surface: 1.0
In [4]:
s1 = ratio * s2
pi = 4. * s1
print("Pi (numerical computation using Monte-Carlo):", pi)
Pi (numerical computation using Monte-Carlo): 3.14316

The error ratio

In [5]:
error = math.pi - pi
print("Error ratio = {:.6f}%".format(abs(error / math.pi) * 100.))
Error ratio = 0.049890%

Graphical illustration

In [6]:
x_array = np.arange(lower_bound, upper_bound, 0.01)
y_array = np.sqrt(1. - x_array**2)

fig = plt.figure(figsize=(6, 6))
plt.axis('equal')
plt.plot(random_points[:,0], random_points[:,1], ',k')
plt.plot(random_points_in[:,0], random_points_in[:,1], ',r')
plt.hlines(y=0, xmin=lower_bound, xmax=upper_bound)
plt.plot(x_array, y_array, '-r', linewidth=2)

plt.show()