Calcule numérique avec la méthode Monte-Carlo (première partie)

TODO

traduire en français certaines phrases restées en anglais

Approximation numérique d'une surface avec la méthode Monte-Carlo

$\newcommand{\ounk}{{\color{red}{O_1}}}$ $\newcommand{\aunk}{{\color{red}{\mathcal{A}_1}}}$ $\newcommand{\nunk}{{\color{red}{\mathcal{N}_1}}}$ $\newcommand{\okn}{{\color{blue}{O_2}}}$ $\newcommand{\akn}{{\color{blue}{\mathcal{A}_2}}}$ $\newcommand{\nkn}{{\color{blue}{\mathcal{N}_2}}}$ Conceptuellement, pour calculer la surface $\aunk$ d'un objet $\ounk$ avec la methode Monte-Carlo, il suffit:

de placer cet objet $\ounk$ entièrement dans une figure géométrique $\okn$ dont on connait la surface $\mathcal \akn$ (par exemple un carré ou un rectangle)
tirer aléatoirement un grand nombre de points dans cette figure $\okn$ (tirage uniforme)
compter le nombre de points $\nunk$ tombés dans l'objet $\ounk$ dont on veut calculer la surface
calculer le rapport $\frac{\nunk}{\nkn}$ où $\nkn$ est le nombre total de points tiré aléatoirement (en multipliant par 100 ce rapport on obtient le pourcentage de points tombés dans l'objet $\ounk$)
appliquer ce rapport à la surface $\mathcal \akn$ de la figure englobante $\okn$ (le carré, rectangle, ... dont on connait la surface) pour obtenir la surface $\aunk$ recherchée: $\aunk \simeq \frac{\nunk}{\nkn} \mathcal \akn$



In [ ]:

    
import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0., 2. * np.pi, 100)
x = np.cos(t) + np.cos(2. * t)
y = np.sin(t)

N = 100
rand = np.array([np.random.uniform(low=-3, high=3, size=N), np.random.uniform(low=-3, high=3, size=N)]).T

fig, ax = plt.subplots(1, 1, figsize=(7, 7))
ax.plot(rand[:,0], rand[:,1], '.k')
ax.plot(x, y, "-r", linewidth=2)
ax.plot([-3, -3, 3, 3, -3], [-3, 3, 3, -3, -3], "-b", linewidth=2)
ax.set_axis_off()
ax.set_xlim([-4, 4])
ax.set_ylim([-4, 4])

plt.show()

Le même principe peut être appliqué pour calculer un volume.

Cette methode très simple est parfois très utile pour calculer la surface (ou le volume) de figures géometriques complexes. En revanche, elle suppose l'existance d'une procedure ou d'une fonction permettant de dire si un point est tombé dans l'objet $O_2$ ou non.

Application au calcul d'intégrales

Calculer l'intégrale d'une fonction, revient à calculer la surface entre la courbe décrivant cette fonction et l'axe des abscisses (les surfaces au dessus de l'axe des abscisses sont ajoutées et les surfaces au dessous sont retranchées).

Exemple written in Python



In [ ]:

    
import sympy as sp
sp.init_printing()

%matplotlib inline



In [ ]:

    
import numpy as np
import matplotlib.pyplot as plt

The function to integrate



In [ ]:

    
def f(x):
    return -x**2 + 3.

Random points



In [ ]:

    
N = 100000           # The number of random points

x_lower_bound = -4.0
x_upper_bound = 4.0
y_lower_bound = -16.0
y_upper_bound = 16.0

random_points = np.array([np.random.uniform(low=x_lower_bound, high=x_upper_bound, size=N),
                          np.random.uniform(low=y_lower_bound, high=y_upper_bound, size=N)]).T

Numerical computation of the integral with Monte-Carlo



In [ ]:

    
# Points between f and the abscissa
random_points_in_pos = np.array([p for p in random_points if 0 <= p[1] <= f(p[0])])
random_points_in_neg = np.array([p for p in random_points if 0 >  p[1] >= f(p[0])])



In [ ]:

    
ratio_pos = float(len(random_points_in_pos)) / float(N)
ratio_neg = float(len(random_points_in_neg)) / float(N)
print('Percentage of "positive" points between f and the abscissa: {:.2f}%'.format(ratio_pos * 100))
print('Percentage of "negative" points between f and the abscissa: {:.2f}%'.format(ratio_neg * 100))

s2 = (x_upper_bound - x_lower_bound) * (y_upper_bound - y_lower_bound)
print("Box surface:", s2)



In [ ]:

    
s1 = ratio_pos * s2 - ratio_neg * s2
print("Function integral (numerical computation using Monte-Carlo):", s1)

The actual integral value



In [ ]:

    
x = sp.symbols("x")

integ = sp.Integral(f(x), (x, x_lower_bound, x_upper_bound))
sp.Eq(integ, integ.doit())

The error ratio



In [ ]:

    
actual_s1 = float(integ.doit())

error = actual_s1 - s1
print("Error ratio = {:.6f}%".format(abs(error / actual_s1) * 100.))

Graphical illustration



In [ ]:

    
fig, ax = plt.subplots(1, 1, figsize=(12, 8))

x_array = np.arange(x_lower_bound, x_upper_bound, 0.01)
y_array = f(x_array)

plt.axis([x_lower_bound, x_upper_bound, y_lower_bound, y_upper_bound])
plt.plot(random_points[:,0], random_points[:,1], ',k')
plt.plot(random_points_in_pos[:,0], random_points_in_pos[:,1], ',r')
plt.plot(random_points_in_neg[:,0], random_points_in_neg[:,1], ',r')
plt.hlines(y=0, xmin=x_lower_bound, xmax=x_upper_bound)
plt.plot(x_array, y_array, '-r', linewidth=2)

plt.show()

Suite

Dans le document suivant, nous allons appliquer le même principe pour calculer la constante $\pi$ : http://www.jdhp.org/docs/notebook/maths_montecarlo_compute_pi_fr.html.