Approximation de Pi par la méthode de Monte Carlo

Il est impossible en informatique de générer, sans intervention d’un système de mesure physique, de suites de nombres aléatoires. Les seules suites que l’on peut générer sont pseudo-aléatoires. Elles ressemblent à des suites aléatoires mais sont déterministes.

Certains problèmes informatiques peuvent être résolus simplement en utilisant des suites aléatoires et parfois seulement pseudo-aléatoires. C’est le cas de l’approximation de \(\pi\). Cette dernière peut être approchée en utilisant la méthode de Monte Carlo avec une suite de nombre pseudo-aléatoire.

Cette méthode repose sur le principe suivant : générer un grand nombre n de points de coordonnées comprise entre \((0,0)\) et \((m,m)\)\(m\) est une valeur fixée. Le rapport entre le nombre de points dans le cercle de centre \((0,0)\) et de rayon \(m\) avec le nombre total \(n\) de points générés est égal à \(\pi/4\).

L’objectif de ce TD est de créer un programme permettant d’approcher la valeur de \(\pi\) par la méthode de Monte Carlo. Pour cela, nous aurons besoin 1) d’une fonction permetant de générer des nombres pseudo-aléatoires ; 2) d’une fonction permetant de tester si un point appartient au disque de centre \((0,0)\) et de rayon \(m\) ; 3) d’une fonction générant un grand nombre \(n\) de points, comptabilisant le nombre de points \(d\) appartenant au cercle de rayon \(m\) afin de renvoyer le rapport \(4d/n = \pi\).

Génération pseudo-aléatoire

Les méthodes existantes pour générer des suites de nombres pseudo-aléatoires reposent sur des propriétés d’arithmétique. Nous en présentons une. Elle permet de générer un nouveau pseudo-aléatoire à partir du précédent. Nous ne décrirons pas ici son fonctionnement, pour plus d’informations reportez-vous sur la page wikipedia suivante : http://en.wikipedia.org/wiki/Linear_congruential_generator.

La suite de nombres pseudo-aléatoires est définie de la manière suivante : \(U_{n+1} \equiv U_{n}*a+c \textrm{ mod } m\). Elle repose sur des valeurs de \(a\), \(c\) et \(m\). Pour que la suite soit réellement proche d’une suite aléatoire, nous ne devons pas choisir ces valeurs arbitrairement. Nous les fixerons au moment de l’utilisation du générateur. Notons que cette suite génère des valeurs pseudo-aléatoire entre \(1\) et \(m\) et quelle nécessite une valeur initiale pour fonctionner : la graine.

Nous souhaitons définir une fonction next_int qui prend en paramètre : p (la valeur précédente de la suite \(U_n\)) ainsi que a, c et m et renvoie la valeur suivante de la suite (\(U_{n+1}\)).




(python_next_int)




(python_next_int_test)

Disque de rayon m

Nous souhaitons maintenant écrire une fonction qui nous permettra de tester si un point de coordonnées \((x,y)\) appartient ou non au disque de centre \((0,0)\) et de rayon \(m\). Pour cela il suffit de vérifier l’égalité suivante \(x*x+y*y \leq m*m\).

Écrire une fonction is_in_disk prenant 3 paramètres x, y, m et renvoie True si le point de coordonnées x, y est dans le disque de rayon m.




(python_is_in_disk)




(python_is_in_disk_test)

Programme général

Nous souhaitons maintenant assembler les briques de base que nous avons réalisées et testées. Pour cela, nous allons écrire une fonction qui prend en paramètres n, a, c, m, s et renvoie une estimation de la valeur de \(\pi\). La valeur n correspondra au nombre de points que nous allons générer. Les valeurs a, c, m et s sont les paramètres de notre générateur pseudo-aléatoire et la graine.




(python_pi)

Il est intéressant de constater qu’en fonction des paramètres utilisés pour la génération des nombres pseudo-aléatoire, la convergence est plus ou moins bonne...

**   PI   **             3.141592
** CONNUS **
Générateur de Knuth :            3.1324  estimate_pi(10000, 6364136223846793005, 1442695040888963407, 2**64, 10)
Générateur "Numerical Recipes" : 3.118   estimate_pi(10000, 1664525, 1013904223, 2**32, 10)
Générateur Borland C/C++ :       3.1488  estimate_pi(10000, 22695477, 1, 2**32, 10)
Générateur RtlUniform :          3.1336  estimate_pi(10000, 2147483629, 2147483587, 2**32-1, 10)
** AU PIF **
Générateur au pif 1 :            4.0     estimate_pi(10000, 1, 2, 7565577, 10)
Générateur au pif 2 :            3.126   estimate_pi(10000, 2147483587, 22695477, 2**32, 10)
Générateur au pif 3 :            0.0056  estimate_pi(10000, 2, 4, 2**32, 10)
Générateur au pif 4 :            3.1592  estimate_pi(10000, 2147483629, 245353423542542, 2**32-1, 10)
Next Section - Le jeu de Nim