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

**TODO**
- traduire en anglais ce qu'il reste à traduire...

## 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$


```python
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=(6, 6))
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


```python
import sympy as sp
sp.init_printing()

%matplotlib inline
```


```python
import numpy as np
import matplotlib.pyplot as plt
```

#### The function to integrate


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

#### Random points


```python
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


```python
# 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])])
```


```python
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)
```


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

#### The actual integral value


```python
x = sp.symbols("x")

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

#### The error ratio


```python
actual_s1 = float(integ.doit())

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

#### Graphical illustration


```python
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.

413f810111ef57b4160dce476d67fb7a8bd64e1b
ipynb
Jupyter Notebook
nb_sci_maths/maths_montecarlo_compute_integral_en.ipynb
jdhp-docs/python-notebooks
91a97ea5cf374337efa7409e4992ea3f26b99179
MIT
2017-05-03T12:23:36.000Z
2020-10-26T17:30:56.000Z
nb_sci_maths/maths_montecarlo_compute_integral_en.ipynb
jdhp-docs/python-notebooks
91a97ea5cf374337efa7409e4992ea3f26b99179
MIT
nb_sci_maths/maths_montecarlo_compute_integral_en.ipynb
jdhp-docs/python-notebooks
91a97ea5cf374337efa7409e4992ea3f26b99179
MIT
2020-10-26T17:30:57.000Z
2020-10-26T17:30:57.000Z
Qwen/Qwen-72B
1. YES
2. YES