Projet 02 : Estimation de Pi par la méthode de Monte-Carlo#
La méthode de Monte-Carlo regroupe une famille d'algorithmes informatiques calculant une solution numérique approchée à l'aide de simulations aléatoires répétées. Conçue dans les années 1940 par Stanislaw Ulam, John von Neumann et Nicholas Metropolis dans le cadre du projet Manhattan, elle est aujourd'hui fondamentale en physique statistique, en finance quantitative et en science des données.
Dans ce projet, nous allons implémenter en langage C un estimateur probabiliste de la constante mathématique \(\pi\) et évaluer sa vitesse de convergence.
1. Principe mathématique#
Considérons le plan cartésien délimité par le carré unité \([0, 1] \times [0, 1]\) et le quart de disque centré à l'origine de rayon \(R = 1\) :
y ^
|
1 +---------+ (Carré de côté 1, Aire = 1)
|* * * * /|
|* * * / |
|* * * / | (Quart de disque de rayon 1,
|* * / | Aire = pi * 1^2 / 4 = pi / 4)
|* * / |
|* / |
0 +---+-----+---> x
0 1
-
L'aire du carré unité est : $\(\text{Aire}_{\text{carré}} = 1 \times 1 = 1\)$
-
L'aire du quart de disque est : $\(\text{Aire}_{\text{disque}} = \frac{\pi \times R^2}{4} = \frac{\pi}{4}\)$
Si l'on tire uniformément au hasard un point \(M(x, y)\) dans le carré \([0, 1] \times [0, 1]\), la probabilité \(P\) que ce point appartienne au quart de disque est le ratio des deux aires : $\(P(M \in \text{disque}) = \frac{\text{Aire}_{\text{disque}}}{\text{Aire}_{\text{carré}}} = \frac{\pi / 4}{1} = \frac{\pi}{4}\)$
Un point \(M(x, y)\) se situe à l'intérieur du disque si et seulement si sa distance euclidienne à l'origine est inférieure ou égale à \(1\) : $\(x^2 + y^2 \le 1\)$
En vertu de la Loi des Grands Nombres, si l'on génère un échantillon total de \(N\) points aléatoires et que \(N_{\text{dedans}}\) points vérifient \(x^2 + y^2 \le 1\), la fréquence empirique converge vers la probabilité théorique :
2. Génération aléatoire uniforme en C#
La fonction rand() de <stdlib.h> renvoie un entier pseudo-aléatoire compris entre 0 et la constante RAND_MAX (au moins 32767). Pour obtenir un réel en virgule flottante uniformément distribué sur l'intervalle \([0, 1]\), on utilise la conversion :
Avant toute simulation, le générateur doit être ensemencé avec une graine variable via srand(time(NULL)) de <time.h>.
3. Implémentation complète en langage C#
Le programme ci-dessous réalise une étude comparative de la précision en fonction du nombre de tirs, calcule l'erreur relative par rapport à la valeur théorique \(\pi \approx 3.141592653589793\) et mesure le temps d'exécution CPU avec clock().
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define PI_THEORIQUE 3.14159265358979323846
// Structure contenant les métriques d'une simulation
typedef struct
{
unsigned long total_tirs;
unsigned long points_interieurs;
double pi_estime;
double erreur_absolue;
double erreur_relative_pct;
double temps_secondes;
} ResultatMonteCarlo;
// Exécute une simulation Monte-Carlo avec N tirs
ResultatMonteCarlo EstimerPi(unsigned long nb_tirs)
{
ResultatMonteCarlo res;
res.total_tirs = nb_tirs;
res.points_interieurs = 0;
clock_t debut = clock();
for (unsigned long i = 0; i < nb_tirs; i++)
{
// Coordonnées uniformes dans [0, 1]
double x = (double)rand() / (double)RAND_MAX;
double y = (double)rand() / (double)RAND_MAX;
// Test d'appartenance au quart de disque : x^2 + y^2 <= 1
if (x * x + y * y <= 1.0)
{
res.points_interieurs++;
}
}
clock_t fin = clock();
// Calcul de l'estimation : pi ~ 4 * (N_interieur / N_total)
res.pi_estime = 4.0 * ((double)res.points_interieurs / (double)res.total_tirs);
res.erreur_absolue = fabs(res.pi_estime - PI_THEORIQUE);
res.erreur_relative_pct = (res.erreur_absolue / PI_THEORIQUE) * 100.0;
res.temps_secondes = (double)(fin - debut) / CLOCKS_PER_SEC;
return res;
}
int main(void)
{
// Initialisation de la graine pseudo-aléatoire
srand((unsigned int)time(NULL));
printf("================================================================================\n");
printf(" ESTIMATION DE PI PAR LA METHODE DE MONTE-CARLO (ALGOJUNGLE) \n");
printf(" Valeur cible pi = %.15f \n", PI_THEORIQUE);
printf("================================================================================\n");
printf("%-12s | %-12s | %-15s | %-15s | %-10s\n",
"Nb Tirs", "Pi Estime", "Erreur Absolue", "Erreur Relative", "Temps (s)");
printf("--------------------------------------------------------------------------------\n");
// Échelle de tests de 1 000 à 10 000 000 de tirs
unsigned long paliers[] = {
1000UL,
10000UL,
100000UL,
1000000UL,
10000000UL
};
int nb_paliers = sizeof(paliers) / sizeof(paliers[0]);
for (int i = 0; i < nb_paliers; i++)
{
ResultatMonteCarlo r = EstimerPi(paliers[i]);
printf("%-12lu | %-12.6f | %-15.6f | %9.4f %% | %7.4f s\n",
r.total_tirs,
r.pi_estime,
r.erreur_absolue,
r.erreur_relative_pct,
r.temps_secondes);
}
printf("================================================================================\n\n");
printf("Conclusion didactique :\n");
printf("L'erreur diminue proportionnellement a 1 / sqrt(N).\n");
printf("Pour gagner un chiffre decimal de precision, il faut multiplier le nombre de tirs par 100 !\n");
return EXIT_SUCCESS;
}
4. Analyse des performances et vitesse de convergence#
Vitesse de convergence en \(\mathcal{O}(1 / \sqrt{N})\)#
Contrairement aux méthodes d'intégration déterministes (méthode des rectangles ou de Simpson), la vitesse de convergence de Monte-Carlo obéit au Théorème Central Limite :
| Nombre de tirs \(N\) | Erreur typique observée | Gain de précision |
|---|---|---|
| \(1\,000\) | \(\approx 2.5\%\) | Point de départ |
| \(10\,000\) | \(\approx 0.8\%\) | Précision \(\times 3\) |
| \(100\,000\) | \(\approx 0.2\%\) | Précision \(\times 10\) |
| \(1\,000\,000\) | \(\approx 0.07\%\) | Précision \(\times 30\) |
| \(10\,000\,000\) | \(\approx 0.02\%\) | Précision \(\times 100\) |
Les limites du générateur standard rand()#
Le générateur linéaire congruentiel rand() de la bibliothèque standard C présente deux limites connues en calcul scientifique :
1. Sa période est relativement courte (\(2^{31}-1\) sur certains compilateurs, voire \(32767\) sur de vieux environnements 16 bits).
2. Pour des simulations à grande échelle (plusieurs milliards de tirages), des générateurs modernes comme le Mersenne Twister (MT19937) ou les algorithmes PCG (Permuted Congruential Generator) sont privilégiés pour éviter tout biais de corrélation statistique.
5. Défi d'application#
Modifiez le programme pour estimer par la même méthode le volume de la sphère unité dans l'espace à trois dimensions \([0, 1] \times [0, 1] \times [0, 1]\) avec la condition \(x^2 + y^2 + z^2 \le 1\). Sachant que le volume théorique du huitième de sphère est \(V = \frac{1}{8} \times \frac{4}{3}\pi = \frac{\pi}{6}\), comparez l'estimation obtenue avec la méthode 2D.
Commentaires