Calcule numérique avec la méthode Monte-Carlo (première 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.

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

Conceptuellement, pour calculer la surface $\mathcal A_1$ d'un objet $O_1$ avec la methode Monte-Carlo, il suffit:

  1. de placer cet objet $O_1$ entièrement dans une figure géométrique $O_2$ dont on connait la surface $\mathcal A_2$ (par exemple un carré ou un rectangle)
  2. tirer aléatoirement un grand nombre de points dans cette figure $O_2$ (tirage uniforme)
  3. compter le nombre de points $N_1$ tombés dans l'objet $O_1$ dont on veut calculer la surface
  4. calculer le rapport $\frac{N_1}{N}$ où $N$ 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 $O_1$)
  5. appliquer ce rapport à la surface $\mathcal A_2$ de figure englobante $O_2$ (le carré, rectangle, ... dont on connait la surface) pour obtenir la surface $\mathcal A_1$ recherchée: $\mathcal A_1 \simeq \frac{N_1}{N} \mathcal A_2$

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 [3]:
import numpy as np
import matplotlib.pyplot as plt

The function to integrate

In [4]:
def f(x):
    return -x**2 + 3.

Random points

In [5]:
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 [6]:
# 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 [7]:
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)
Percentage of "positive" points between f and the abscissa: 2.72%
Percentage of "negative" points between f and the abscissa: 10.11%
Box surface: 256.0
In [8]:
s1 = ratio_pos * s2 - ratio_neg * s2
print("Function integral (numerical computation using Monte-Carlo):", s1)
Function integral (numerical computation using Monte-Carlo): -18.915840000000003

The actual integral value

Out[9]:
$$\int_{-4.0}^{4.0} - x^{2} + 3.0\, dx = -18.6666666666667$$

The error ratio

Error ratio = 1.334857%

Graphical illustration

In [11]:
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_en.html.