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.
...
In [ ]:
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))
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 [ ]:
random_points_in = np.array([p for p in random_points if p[0]**2 + p[1]**2 <= 1.0])
In [ ]:
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)
In [ ]:
s1 = ratio * s2
pi = 4. * s1
print("Pi (numerical computation using Monte-Carlo):", pi)
In [ ]:
error = math.pi - pi
print("Error ratio = {:.6f}%".format(abs(error / math.pi) * 100.))
In [ ]:
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()