Méthode de Monte Carlo (Solution)

Les méthodes de Monté Carlo sont des méthodes probabilistes basées sur l'observation d'un grand nombre d'évènements.

Pour illustrer la méthode, nous allons construire une fonction permettant de calculer l'aire d'une surface du plan définie par une équation implicite. La méthode de Monté Carlo employée permettra de trouver une approximation de l'aire de la surface.

Les outils suivants vous seront peut être utiles :

• Wikipedia
• Wolfram Alpha

Le détail des choses à rendre est disponible sur Updago.

Pour expliquer le principe, nous allons prendre comme exemple de surface un disque de rayon 1 (pour réfléchir au problème, il est utile de le faire avec un exemple pour lequel la réponse est vérifiable…). La surface de ce disque étant connue, et égale à $\pi$, la méthode nous permettra donc de calculer une approximation de $\pi$.

En premier lieu, il faut chercher une zone simple (de préférence rectangulaire) d'aire $A$ connue, qui contient entièrement la surface à évaluer. Dans notre exemple, la zone rectangulaire sera : $-1\le x\le1$ et $-1\le y\le1$. Elle est délimitée par le carré dans lequel est inscrit notre disque de rayon 1.

Puis, on choisira $N$ points au hasard dans cette zone (c'est très facile pour une zone rectangulaire) et nous vérifierons combien sont dans la surface à évaluer. Cette opération de vérification doit être facile à faire. Ce sera effectivement le cas si la surface à calculer est définie par une équation implicite du type de $f(x,y)\leq a$ (dans l'exemple du disque : $f(x,y)=x^2+y^2 \leq 1$).

Admettons que sur les $N$ points choisis dans la zone rectangulaire, $n$ soient dans la surface définie par l'équation implicite.

L'aire cherchée vaudra approximativement : $\frac{n\timesA}{N}$

Plus $N$ sera grand, plus l'approximation sera précise.

Essayons sur notre exemple :

Avec 10000 points nous pourrions par exemple obtenir 7784 points dans le disque, ce qui nous donnerait une aire de : 3.1136

Toutes les fonctions seront entrées dans un fichier, que vous nommerez login_montecarlo.py. Ce fichier devra contenir uniquement des importations de module et vos fonctions, mais pas de script (le module pourra être importé silencieusement). De même, vos fonctions doivent avoir le nom indiqué dans le sujet. Enfin, pensez à donner une docstring claire pour chaque fonction écrite

Commencez par programmer l'équation implicite du disque. Pour cela, écrivez une fonction disque qui prend en paramètre x et y et renvoie 1 si le point (x,y) est dans le disque centré en 0 de rayon 1 et 0 sinon ((le choix des valeurs renvoyées simplifiera le comptage)) :

def  disque(x,y) :
    ''' Renvoie 1 si (x,y) est dans le disque
    centré en 0 de rayon 1 et renvoie 0 sinon.
    '''
    # Complétez ici...

Dans le shell, testez ensuite votre fonction sur quelques points pour vérifier qu'elle est correcte.

disque(0,0)
# => 1
disque(1.1,0) 
# => 0
disque(0.6,0.6) 
# => 1
disque(0.8,0.75) 
# => 0

Avant d'aborder le calcul de l'aire, nous allons voir comment passer une fonction en paramètre d'une autre fonction. Voici un exemple (inutile de le mettre dans votre fichier, il s'agit juste d'un exemple pour vous aidez à comprendre. Si toutefois vous voulez l'essayer, faites le dans un autre fichier, que vous ne rendrez pas) :

def foisdeux(x) :
    return 2*x
 
def foistrois(x) :
    return 3*x
 
def evaluation(fonc,x) :
    return fonc(x)

Nous pouvons ensuite écrire :

evaluation(foisdeux,4) 
# => 8
evaluation(foistrois,4) 
# => 12

Le premier paramètre de evaluation est une fonction

Si le passage d'une fonction en paramètre n'est pas parfaitement compris avec cet exemple, n'hésitez pas demander plus d'explications à l'encadrant.

La fonction de calcul d'aire que nous allons écrire devra marcher pour n'importe quelle équation implicite, du moment que cette équation est écrite sous forme d'une fonction prenant en paramètres x et y et renvoyant 0 ou 1 (c'est le cas de la fonction disque).

La fonction de calcul d'aire, que nous appellerons calcul_aire, aura donc pour paramètres :

• la fonction contenant l'équation implicite
• le nombre de points à utiliser pour calculer l'aire
• les bornes du rectangle englobant

Le rectangle englobant sera donné sous la forme d'une séquence (tuple ou liste) contenant : (xmin,xmax,ymin,ymax)

Nous pouvons à présent nous concentrer sur l'écriture de la fonction calcul_aire :

def calcul_aire(f,N,bornes) :
    ''' Calcule l'aire de la surface définie
        par l'équation implicite f(x,y), et incluse
        dans le rectangle défini par 'bornes'.
        Le calcul est réalisé en utilisant N points
    '''
    # Calculer l'aire du rectangle englobant => A
    # Initialiser le compteur de points nb à 0
    # Répéter N fois :
        # prendre un couple (x,y) au hasard dans les bornes
        # ajouter f(x,y) à nb (f(x,y) vaut 0 ou 1)
    # retourner nb/N*A    

Avec le module Matplotlib, vous avez la possibilité de visualiser facilement à l'écran les points dans et hors surface. Vous pouvez pour ceci refaire une fonction, similaire à calcul_aire, mais qui tracera un graphique au lieu de calculer l'aire numériquement. Pour réaliser ce graphique, vous aurez besoin de comprendre le code suivant :

import matplotlib.pyplot as plt
x=(1,2,4,5)        # abcisses
y=(1.2,1.3,0.5,-1) # ordonnées
plt.plot(x,y,'.b') # Trace les points dans le style '.b'
plt.show()

Dans le paramètre à plot qui définit le style, «.» indique qu'on trace des points, et b donne la couleur (blue). La liste des styles est disponible ici : [[http://matplotlib.sourceforge.net/api/pyplot_api.html#matplotlib.pyplot.plot]]

Vous pouvez afficher deux séries de points dans deux styles différents ainsi :

import matplotlib.pyplot as plt
xb,yb=(1,2,4,5) , (1.2,1.3,0.5,-1) 
xa,ya=(1.3,2.2,3.8,5.2) , (1.2,1.3,0.5,-1) 
plt.plot(xb,yb,'.b',xa,ya,'.r')
plt.show()

Vous pouvez maintenant utiliser votre fonction pour d'autres surface :

1. ellipse centrée en 0, et de rayons 2 et 1 (la fonction implicite sera codée dans ellipse(x,y)). L'ellipse soit être dans le sens horizontal (plus étalée sur l'axe des abscisses que des ordonnées)
2. cardioïde obtenue avec des cercles de rayon 2 (la fonction s'appellera cardioide(x,y), orientée dans ce sens :

Vous pouvez utiliser des nombres complexes en Python. Des exemples sont donnés dans la documentation : [[stu:python:python_notes#nombres_complexes|Aide mémoire / Notes sur Python 3]].

Il vous faudra environ 50000 tirages pour avoir une aire assez précise.