Approximation d'une surface par la méthode de Monté-Carlo¶
Enoncé
Soit les fonctions définies sur $[0;2]$ par :
$f(x)=\dfrac{4x}{x²+1} $ et $ g(x)=x(x-1,2)$
Les courbes $\mathcal{C}_f$ et $\mathcal{C}_g$ délimitent une surface représentée ci-dessus, hachurée en orange.
Une méthode (la méthode de Monte-Carlo) pour calculer une valeur approchée de l’aire d’une surface consiste en un tirage aléatoire successif de points et de compter le nombre de points situés dans la surface cherchée.
Ainsi le rapport : $\dfrac{\text{nombre de points dans la surface}}{\text{nombre de points total}}\times \text{surface du rectangle}$
donnera une approximation de cette aire pour un grand nombre de tirage.
(1) Ecrire une fonction test_position prenant en argument x,y,f,g qui teste si le point de cordonnées (x,y) est dans la surface ou pas,
elle renvoie True ou False.
(2) Ecrire le programme principal déterminant la valeur approximative de cette surface.
Il est courant d'utiliser des nombres aléatoires pour estimer des quantités numériques. Nous allons ici utiliser des variables aléatoires de loi uniforme afin d'estimer l'aire d'une surface. Pour cela nous allons:
- générer des points aléatoirement et uniformément dans un rectangle de côté 2 $\times$ 2,5;
- calculer la proportion de points situés entre les deux courbes représentatives des fonctions $f$ et $g$. Ce qui nous permettra de donner une estimation de son aire;
Estimation¶
La fonction random de la librairie random permet de générer un nombre aléatoirement et uniformément entre 0 et 1.
from random import*
def test_position(x,y,f,g):
if y<f(x) and y>g(x):
return True
else:
return False
### Les fonctions ###
f=lambda x:4*x/(x**2+1)
g=lambda x:x*(x-1.2)
### Les variables ###
compteur=0
N=100000
### Le corps du programme ###
for i in range(N):
x=uniform(0,2)
y=uniform(-0.5,2)
if test_position(x,y,f,g):
compteur+=1
print('le compteur:',compteur)
print('aire approximative :',compteur/n*5)
Convergence¶
Nous allons maintenant utiliser les listes afin de stocker au fur et à mesure l'estimation de cette surface.
%matplotlib inline
import matplotlib.pyplot as plt
from random import *
def test_position(x,y,f,g):
if y<f(x) and y>g(x):
return True
else:
return False
### Les variables ###
N=1000
compteur=0
estimation = []
### Les fonctions ###
f=lambda x:4*x/(x**2+1)
g=lambda x:x*(x-1.2)
### Le corps du programme ###
for i in range(N):
x=uniform(0,2)
y=uniform(-0.5,2)
if test_position(x,y,f,g):
compteur+=1
estimation.append(5*compteur/(i+1))
plt.plot([0,N],[2.95,2.95],color='m')
plt.text(N+2,2.95,"Surface",color='m',fontsize=14)
plt.title("Estimation de la surface avec {} points: S=$\simeq${}".format(N,estimation[-1]),color="#1e7fcb",fontsize=14)
plt.plot(estimation)
plt.show()
Calcul de surface avec le calcul formel¶
from sympy import*
x=Symbol("x")
integrate(4*x/(x**2+1)-x*(x-1.2),(x,0,2))
Visualisation de l'erreur¶
On ajoute une visualisation de l'erreur commise au fur et à mesure.
%matplotlib inline
import matplotlib.pyplot as plt
from random import *
from sympy import*
### figure ###
fig, (ax1, ax2) = plt.subplots(1, 2,figsize=(15, 6))
def test_position(x,y,f,g):
if y<f(x) and y>g(x):
return True
else:
return False
def calcul_formel():
x=Symbol("x")
return integrate(4*x/(x**2+1)-x*(x-1.2),(x,0,2))
### Les variables ###
N=500
compteur=0
estimations = []
erreurs = []
### Les fonctions ###
f=lambda x:4*x/(x**2+1)
g=lambda x:x*(x-1.2)
### Le corps du programme ###
surface=float(calcul_formel())
for i in range(N):
x=uniform(0,2)
y=uniform(-0.5,2)
if test_position(x,y,f,g):
compteur+=1
estimation=5*compteur/(i+1)
estimations.append(estimation)
erreurs.append(abs(estimation-calcul_formel()))
ax1.plot([0,N],[surface,surface],color='m')
ax1.text(N+2,2.95,"Surface",color='m',fontsize=14)
ax1.set_title("Estimation de la surface avec {} points: S$\simeq${}".format(N,estimations[-1]),color="#1e7fcb",fontsize=14)
ax2.set_title("Evolution de l'erreur avec {} points".format(N),color="#1e7fcb",fontsize=14)
ax1.plot(estimations)
ax2.plot(erreurs,'r')
ax2.plot([0,N],[0,0],color='orange')
plt.show()
Animation¶
Simulation de points générés aléatoirement et approximation de la surface
import matplotlib.pyplot as plt
import matplotlib.animation
from random import random as rd
from IPython.display import HTML
import numpy as np
### Les constantes ###
N = 200
### Les fonctions ###
f=lambda x:4*x/(x**2+1)
g=lambda x:x*(x-1.2)
### paramètres de la figure ###
fig, (ax1, ax2) = plt.subplots(1, 2,figsize=(15, 6))
ax1.set_xlim(( 0, 2))
ax1.set_ylim((-0.5,2))
ax2.set_xlim(( 0, N))
ax2.set_ylim((1, 4))
ax1.set_aspect('equal')
ax2.set_aspect(N//4)
Xin,Yin,Xout,Yout = [],[],[],[]
compteur = 0
surf = []
pointsIn, = ax1.plot([],[],'go')
pointsOut, = ax1.plot([],[],'ro')
estimation = ax2.text(N//2-N//4,1.2,'Estimation',color="#1e7fcb",fontsize=14)
courbe, = ax2.plot([],[])
def test_position(x,y,f,g):
if y<f(x) and y>g(x):
return True
else:
return False
def init():
pointsIn.set_data([], [])
pointsOut.set_data([], [])
estimation.set_text("0")
courbe.set_data([], [])
ax1.grid(True)
x = np.linspace(0, 2, 30)
y1 = f(x)
y2=g(x)
ax1.plot(x, y1)
ax1.plot(x,y2)
ax2.plot([0,N],[2.95,2.95],'m-')
ax2.text(N+0.5,2.95,'Surface',color="m",fontsize=14)
ax2.set_title('5 x Proportion de points dans la surface')
return (pointsIn,)
def animate(i):
global compteur
x=uniform(0,2)
y=uniform(-0.5,2)
if test_position(x,y,f,g):
compteur+=1
Xin.append(x)
Yin.append(y)
pointsIn.set_data([Xin,Yin])
else:
Xout.append(x)
Yout.append(y)
pointsOut.set_data([Xout,Yout])
surf.append(5*compteur/(i+1))
estimation.set_text("Surface $\simeq${} avec {} points".format(int(1000*surf[i])/1000,i+1))
courbe.set_data(range(i+1),surf)
return (pointsIn,)
ani = matplotlib.animation.FuncAnimation(fig, animate, frames=N,init_func=init,blit=True)
# l'un ou l'autre
HTML(ani.to_jshtml())
#HTML(ani.to_html5_video())