Méthode de Monte-Carlo pour un calcul d'aire - Python

Principe de la méthode de Monte-Carlo

Considérons pour commencer un cas évident : soit un rectangle $ABCD$, que nous divisons en deux parties égales :

Nous comprenons très facilement que, si nous prenons un point de manière aléatoire dans le rectangle $ABCD$, la probabilité pour qu’il appartienne à $AEFD$ est égale à :

$\dfrac {\text{aire}(AEFD)}{\text{aire}(ABCD)}=\dfrac 12$

Il en va de même pour une figure $\mathcal F$ quelconque inscrite dans le rectangle :

La probabilité $p$ qu’un point pris au hasard dans $ABCD$ appartienne à $\mathcal F$ est égale à :

$p=\dfrac {\text{aire}(\mathcal F)}{\text{aire}(ABCD)}$

• On en déduit :

$\text{aire}(\mathcal F) = p\times \text{aire}(ABCD)$

À retenir

Voilà en quoi consiste la méthode de Monte-Carlo :

• la figure dont on cherche l’aire $\mathcal A$ est incluse dans un rectangle dont on peut calculer facilement l’aire $\mathcal A_\text{rect}$ ;
• on prend au hasard un grand nombre $n$ de points appartenant au rectangle, et on compte ceux qui appartiennent à la figure, pour obtenir leur nombre $c$ ;
• on en déduit la fréquence $f=\frac cn$ d’appartenance à la figure ;
• plus on prend de points, plus la fréquence d’appartenance à la figure $f$ approche de la probabilité théorique (loi des grands nombres) ;
• nous obtenons ainsi l’approximation :

$\begin{aligned} \mathcal A&\approx f\times \mathcal A$\text{rect} \ &\approx \dfrac cn\times \mathcal A\text{rect} \end{aligned}

Programme Python

Prérequis

Soit $g$ la fonction définie et continue sur $[0\ ;\, +\infty[$ par :

$g(x)=x(x^2-3x+3)$

Nous nous intéressons à l’aire $\mathcal A$ du domaine $\mathscr D$ sous la courbe représentative de $g$, sur l’intervalle $[0\ ;\, 1]$.

• Nous avons, dans un repère $(O\ ;\, I,\,J)$, avec $K\,(1\ ;\, 1)$ :

$g$ est strictement croissante et positive sur $[0\ ;\, +\infty[$ (vous pouvez le montrer).
Et nous avons $g(0)=0$ et $g(1)=1$.

• Cela signifie donc que la partie d’aire $\mathcal A$ est incluse dans le carré $OIKJ$.

Remarquons que cette aire est égale à :

$\mathcal A = \int_{0}^{1} g(x) \text d x$

Attention

Dans toute cette fiche, nous travaillerons avec des fonctions positives sur l’intervalle considéré. Ainsi, le programme que nous écrirons ne permettra de travailler qu’avec de telles fonctions (mais vous pourrez aussi adapter pour que cela marche également avec des fonctions négatives sur l’intervalle).

Comme nous l’avons vu dans la première partie, nous prendrons de manière aléatoire un grand nombre de points appartenant au carré $OIKJ$, d’aire égale à $1\ \text{u.a.}$

Pour chaque point, nous regarderons s’il appartient à $\mathscr D$.
Pour savoir si un point $M\,(x$M\ ;\, yM) appartient à $\mathscr D$, il suffit de calculer $f(x$M)f(xM​) et de comparer le résultat à $y$M.

• En effet, logiquement :
• si $y$M \leq g(xM), alors $M$ est sous la courbe et appartient à $\mathscr D$ ;
• si $y$M > g(xM), alors $M$ est au-dessus de la courbe et n’appartient pas à $\mathscr D$.

Nous compterons ceux qui appartiennent à $\mathscr D$ et en déduirons la fréquence $f$ de tels points.

• Cela nous permettra de donner une approximation de $\mathcal A$ :

$\begin{aligned} \mathcal A&\approx f\times \text{aire}(OIKJ) \ &\approx f \footnotesize{\textcolor{#A9A9A9}{\text{ [car aire$(OIKJ)=1\ \text{u.a.}$]}}} \end{aligned}$

Programme Python

Pour choisir au hasard un point, nous utiliserons la fonction $\purple{\text{random()}}$ du module $\purple{\text{random}}$, pour générer de manière aléatoire son abscisse et son ordonnée, comprise entre $0$ et $1$.

• Nous commençons donc par importer le module :

• Définissons aussi notre fonction $\purple{\text{g}}$, qui prend en paramètre $\purple{\text{x}}$ :

• $g:x\mapsto x(x^2-3x+3)$ : nous demandons à la fonction de faire le calcul et de retourner le résultat.

• Notre fonction Python $\purple{\text{aire}}$ prend en paramètres :
• la fonction $\purple{\text{f}}$ à laquelle on s’intéresse (dans notre cas, il s’agira donc de $\purple{\text{g}}$) ;
• le nombre $\purple{\text{n}}$ de points que nous voulons générer aléatoirement.

• Nous initialisons le compteur $\purple{\text{c}}$, qui comptabilisera le nombre de points qui sont sous la courbe :

• Nous utilisons une boucle $\purple{\text{for}}$ pour générer les $\purple{\text{n}}$ points :

Pour générer un point, nous utilisons la fonction $\purple{\text{random()}}$ pour choisir aléatoirement deux nombres $\purple{\text{x}}$ et $\purple{\text{y}}$ compris entre $0$ et $1$, qui représenteront respectivement :

• son abscisse,
• et son ordonnée.

• Comme nous l’avons dit, nous regardons alors si $\purple{\text{y}}$ est inférieur ou égal à l’image de $\purple{\text{x}}$ par la fonction $\purple{\text{f}}$ :

• Dans ce cas, nous comptabilisons le point comme étant sous la courbe :

• Nous allons maintenant pouvoir donner l’approximation de l’aire.
• Comme nous l’avons dit en prérequis, elle vaut environ la fréquence, que l’on calcule en divisant le nombre de points sous la courbe par le nombre total de points :

À retenir

Le programme est terminé :

\begin{aligned} \textcolor{#A9A9A9}{\text{1}}&\quad\text{from random import $\ast$} \ \textcolor{#A9A9A9}{\text{2}}& \ \textcolor{#A9A9A9}{\text{3}}&\quad\text{def g(x):} \ \textcolor{#A9A9A9}{\text{4}}&\quad\qquad\text{return x $\ast$ (x$\ast\ast\,$2 - 3 $\ast$ x + 3)} \ \textcolor{#A9A9A9}{\text{5}}& \ \textcolor{#A9A9A9}{\text{6}}&\quad\text{def aire(f, n):} \ \textcolor{#A9A9A9}{\text{7}}&\quad\qquad\text{c = 0} \ \textcolor{#A9A9A9}{\text{8}}&\quad\qquad\text{for i in range(n):} \ \textcolor{#A9A9A9}{\text{9}}&\quad\qquad\qquad\text{x = random()} \ \textcolor{#A9A9A9}{\text{10}}&\quad\qquad\qquad\text{y = random()} \ \textcolor{#A9A9A9}{\text{11}}&\quad\qquad\qquad\text{if y <= f(x):} \ \textcolor{#A9A9A9}{\text{12}}&\quad\qquad\qquad\qquad\text{c = c + 1} \ \textcolor{#A9A9A9}{\text{13}}&\quad\qquad\text{print($^{\backprime\backprime}$L'aire sous la courbe vaut environ : $^{\prime\prime}$)} \ \textcolor{#A9A9A9}{\text{14}}&\quad\qquad\text{print(c / n, $^{\backprime\backprime}$u.a.$^{\prime\prime}$)} \end{aligned}

Exemple

Exécutez ce programme pour la fonction $\purple{\text{g}}$, avec $10^6$ points à choisir aléatoirement :

$\purple{\text{aire(g, 10$\ast\ast\,$6)}}$

• Faites-le plusieurs fois et donnez une approximation de $\mathcal A$.

Maintenant, calculez, comme vous l’avez appris en cours :

$\int_{0}^{1} x(x^2-3x+3) \text d x$

• Confrontez les résultats.

Astuce

$\purple{\text{aire}}$ reste valable pour toute fonction continue et positive sur $[0\ ;\, 1]$, si son maximum sur cet intervalle est inférieur ou égal à $1$. (Sinon, le carré $OIKJ$ ne contiendrait pas tout le domaine dont on cherche l’aire…)

• Vous pouvez essayer avec la fonction carrée, dont, pour confrontation, on calcule facilement l’intégrale sur $[0\ ;\, 1]$ :

$\begin{aligned} \int${0}^{1} x^2 \text d x &= \left[\dfrac 13x^3\right]0^1 \ &=\dfrac 13 \end{aligned}

Une fonction Python améliorée

Nous l’avons dit, la fonction $\purple{\text{aire}}$ que nous avons programmée ne fonctionne que si on travaille sur l’intervalle $[0\ ;\, 1]$ avec une fonction positive dont le maximum est inférieur ou égal à $1$.

Nous allons donner ci-dessous une fonction $\purple{\text{aire2}}$, qui sera plus générale.
Cette fonction prendra pour arguments :

• une fonction $\purple{\text{f}}$, comme dans le précédent programme ;
• les bornes $\purple{\text{a}}$ et $\purple{\text{b}}$ de l’intervalle sur lequel on souhaite travailler ($\purple{\text{a}} < \purple{\text{b}}$) ;
• la « hauteur » $\purple{\text{h}}$ du rectangle, qui doit donc inclure l’ensemble du domaine dont on cherche l’aire :
• on pourra déterminer $\purple{\text{h}}$ graphiquement, si on a une représentation graphique de $\purple{\text{f}}$, on choisira alors $\purple{\text{h}}$ supérieur ou égal au maximum de la fonction sur $[\purple{\text{a}}\ ;\, \purple{\text{b}}]$ (choisir $\purple{\text{h}}$ le plus proche possible du maximum) ;
• si on connaît le maximum, on pourra l’entrer directement ;
• si on sait juste que $\purple{\text{f}}$ admet un maximum sur $[\purple{\text{a}}\ ;\, \purple{\text{b}}]$ en $\purple{\text{m}}$, on peut mettre directement en argument, pour $\purple{\text{h}}$ : $\purple{\text{f(m)}}$ ;
• si $\purple{\text{f}}$ admet un maximum en $\purple{\text{a}}$ ou $\purple{\text{b}}$ (pour une fonction monotone sur $[\purple{\text{a}}\ ;\, \purple{\text{b}}]$, par exemple), on peut tout simplement mettre $\purple{\text{0}}$, le programme prendra la plus grande valeur entre $\purple{\text{f(a)}}$ et $\purple{\text{f(b)}}$ ;
• le nombre $\purple{\text{n}}$ de points que nous souhaitons.

Astuce

Pour générer un nombre aléatoire non plus compris entre $0$ et $1$, mais compris entre $p$ et $q$, on utilise la fonction présente dans le module $\purple{\text{random}}$ :

$\purple{\text{uniform(p, q)}}$

À retenir

\begin{aligned} \textcolor{#A9A9A9}{\text{1}}&\quad\text{from random import $\ast$} \\ \textcolor{#A9A9A9}{\text{2}}& \\ \textcolor{#A9A9A9}{\text{3}}&\quad\text{def aire2(f, a, b, h, n):} \\ \textcolor{#A9A9A9}{\text{4}}&\quad\qquad\text{l = b - a} \\ \textcolor{#A9A9A9}{\text{5}}&\quad\qquad\text{if h == 0:} \\ \textcolor{#A9A9A9}{\text{6}}&\quad\qquad\qquad\text{h = max(f(a), f(b))} \\ \textcolor{#A9A9A9}{\text{7}}&\quad\qquad\text{r = l $\ast$ h} \\ \textcolor{#A9A9A9}{\text{8}}&\quad\qquad\text{c = 0} \\ \textcolor{#A9A9A9}{\text{9}}&\quad\qquad\text{for i in range(n):} \\ \textcolor{#A9A9A9}{\text{10}}&\quad\qquad\qquad\text{x = uniform(a, b)} \\ \textcolor{#A9A9A9}{\text{11}}&\quad\qquad\qquad\text{y = uniform(0, h)} \\ \textcolor{#A9A9A9}{\text{12}}&\quad\qquad\qquad\text{if y <= f(x):} \\ \textcolor{#A9A9A9}{\text{13}}&\quad\qquad\qquad\qquad\text{c = c + 1} \\ \textcolor{#A9A9A9}{\text{14}}&\quad\qquad\text{print($^{\backprime\backprime}$L'aire sous la courbe vaut environ :$^{\prime\prime}$)} \\ \textcolor{#A9A9A9}{\text{15}}&\quad\qquad\text{print(r $\ast$ c / n, $^{\backprime\backprime}$u.a.$^{\prime\prime}$)} \end{aligned}

Attention

Ce programme Python ne fonctionne toujours que pour les fonctions positives sur $[a\ ;\, b]$.

Exemple

La fonction $x\mapsto 4x^2+3x$ est continue, positive et strictement croissante sur $[0\ ;\, +\infty[$.

• En vous servant du dernier programme, donnez une valeur approchée de :

$I = \int_{1}^{\frac 32} (4x^2+3x)\text d x$

• Calculez maintenant la valeur exacte de $I$.