Aller au menu - Aller au contenu

Icône Calculons pi !

Avatar
Par Avatar sebsheep
Mise à jour : 25/05/2010
Difficulté : Facile Facile Creative Commons BY
583 visites depuis 7 jours, dont 27 sur ce chapitre classé 194/786
Bon, je t'arrête tout de suite, pi je connais depuis que j'ai 10 ans, ça vaut 3.14159... !


Je suis d'accord avec vous ... Mais comment connaissez-vous cette valeur ? Grâce à la maîtresse qui vous a demandé d'apprendre par cœur ce nombre ? Mais cette maîtresse, comment elle l'a su ? Grâce à sa propre maîtresse ? :euh:

On voit bien que si l'on continue le raisonnement, il faut bien que quelqu'un, un jour, ait trouvé cette valeur. Ce calcul a toujours été un grand problème des mathématiques et plusieurs méthodes ont été imaginées ; on peut notamment citer celle d'Archimède dans la Grèce antique. Aujourd'hui encore, certains fous furieux tentent de découvrir des méthodes encore plus efficaces afin de connaître encore plus de décimales.

Ici, nous n'allons rien faire de compliqué, la méthode sera très simple à comprendre (satisfait ou remboursé !).
Retour en haut
Sommaire du chapitre :
Icône du chapitre
Chapitre précédent Sommaire Chapitre suivant
Retour en haut

Titi à la rescousse

Un jeu de fléchettes



Je vais vous montrer ici une technique très simple, tellement simple que c'est Titi, votre petit frère de 3 ans et demi, qui va nous servir de "calculateur". On va demander à Titi de jouer à un jeu, lui aussi très simple, consistant à envoyer des fléchettes sur la cible représentée ci-dessous :

cible vide

On ne va pas lui demander de viser, il en est incapable ; la fléchette doit juste arriver dans le carré. On donne à Titi 100 fléchettes et on le laisse jouer une heure ou deux. A la fin, la cible ressemblera à :

cible criblée de fléchette


Refléchissons un peu



Voilà, le plus gros du travail est fait (pas très fatigant, et en plus, on a occupé Titi tout l'après-midi) ! Bon, laissons Titi tranquille pour le moment et commençons à réfléchir... Lorsque Titi lançait ses fléchettes, le jet était complètement aléatoire, donc la probabilité qu'elles atterrissent sur la partie rouge est proportionnelle à la surface de cette partie rouge. Cela se comprend intuitivement : plus une partie est grande, plus on a de chance de "tomber" dedans. Plus précisément, on a :
p=\text{"Probabilit\'{e} que la fl\'{e}chette soit dans le rouge"} = \frac{Surface\ rouge}{Surface\ totale}


On cherche à connaitre p (vous comprendrez pourquoi dans un instant). Pour cela, répondons à quelques questions : comment est constituée la cible ? Réponse : un quart de disque rouge dans un carré noir. Maintenant, un peu plus dur : Quelle est l'aire de la partie rouge (on considère que le carré est de côté 1) ? Vous séchez ? Alors prenons calmement : un disque de rayon R a une surface de \pi R^2. Ici nous avons un quart de disque, donc sa surface est \frac{\pi R^2}{4} ; et R=1, donc, la surface est de \frac{\pi}{4}.

La surface totale de la cible est de 1x1=1 (c'est un carré de côté 1).

Donc nous avons : p=\frac{Surface\ rouge}{Surface\ totale} = \frac{\pi/4}{1} = \frac{\pi}{4}.

Mais nous avons un autre moyen d'estimer p. Cette estimation est toute simple, elle consiste à considérer que p est à peu près égale au rapport \frac{nombre\ de\ points\ dans\ le\ rouge}{nombre\ de\ points\ total}. Évidemment, plus le nombre de points est important, meilleure sera cette estimation. En m'amusant à compter combien de points sont dans la partie rouge, je trouve : 82. En remplaçant dans les formules trouvées plus tôt, on a donc : p=\frac{82}{100}=\frac{\pi}{4}. Une petite manipulation très simple donne : \pi = \frac{82}{100}*4=3.28.

Mais ton résultat est faux ! Elle n'est pas juste ta méthode !


Je dirais plutôt que mon résultat comporte une erreur ; souvenez vous, on a estimé la valeur p comme étant le rapport "nombre de points dans le rouge"/"nombre de points total". Il est donc normal que l'on ait qu'une estimation de pi. Pour réduire l'erreur, il suffit d'augmenter le nombre de tirages ; mais malheureusement, Titi est fatigué et ne peut plus rien faire... Nous allons donc demander à nos ordinateurs de lancer des fléchettes à notre place.
Retour en haut

Un piètre tireur : votre ordinateur

Nous allons donc implémenter ici un algorithme permettant de reproduire l'expérience que l'on vient de faire avec Titi ; grossièrement, il ressemblera à ceci :
Code : Autre
1
2
3
4
5

initialiser un compteur c à 0
lancer N fléchettes dans un carré de côté 1 :
    si la fléchette est dans le quart de cercle, on augmente le compteur de 1
    sinon, on ne fait rien
on retourne c*4/N


Les seuls points un peu compliqués sont les lignes 2 et 3.

Lancer une fléchette



Pour "lancer une fléchette dans un carré de côté 1", il suffit de tirer au hasard deux nombres compris entre 0 et 1. Pour ça, il faut implémenter cette fonction à partir de rand (référez-vous au tutoriel de natim sur le sujet pour tout comprendre) :
Code : C
1
2
3
4
double rand_0_1(void)
{
   return rand()/(double)RAND_MAX;
}


Maintenant, pour tirer effectivement notre fléchette, c'est tout simple :
Code : C
1
2
double x = rand_0_1();
double y = rand_0_1();


Les variables x et y contiennent alors respectivement l'abcisse et l'ordonnée de la fléchette.

Sommes-nous dans le disque ?



Image utilisateur

Pour répondre à cette question, il suffit de définir précisément ce qu'est un disque de rayon 1 et de centre O. On peut prendre la définition suivante :

Un disque de rayon 1 et de centre O est l'ensemble des points M tels que la distance OM soit inférieure ou égale à 1.

Pour mieux comprendre cette définition un peu abstraite, un petit exemple : prenons le point M(0.5,0.2) ; M est-il dans le cercle ?

Pour répondre, utilisons la définition : il faut que MO \leq 1. Pour simplifier les calculs, prenons le carré de cette inégalité : MO^2 \leq 1^2=1.
Or d'après notre cher Pythagore : MO² = (0.5-0)²+(0.2-0)² = 0.5² + 0.2² = 0.25 + 0.04 = 0.29.

0.29 est bien inférieur à 1 donc M est dans le disque.

Donc dans notre problème, la condition "le point (x,y) est dans le disque" s'exprimera :
Code : C
1
( x*x + y*y ) <= 1;


À vous de jouer !



Vous avez maintenant tous les ingrédients pour implémenter l'algorithme. Je vous propose donc de faire un programme qui demande à l'utilisateur combien de lancers il veut faire, puis qui affiche une valeur approchée de pi.

Voici un programme possible :
Secret (cliquez pour afficher)
Code : C
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

double rand_0_1(void)
{
	return rand()/(double) RAND_MAX;
}

int main(int argc, char *argv[])
{
	unsigned int n=0; // le nombre de jets
	unsigned int compteur=0; // le compteur de flechettes dans la cible
	unsigned int i=0; // le compteur de boucle
	double x,y;
	double pi;
	printf("Combien de jets voulez vous effectuer ? ");
	scanf("%d",&n);
	
	printf("Simulation en cours, cela peut prendre quelques minutes ... \n");

	for(i=0 ; i< n; i++)
	{
		// Jet de la flechette
		x = rand_0_1();
		y = rand_0_1();
		
		// La flechette est-elle dans le disque ?
		if( x*x + y*y <= 1 )
		{
			compteur++;
		}
	}
	
	pi = compteur*4 / (double) n; // "(double)" permet de convertir n d'un entier vers un double ;
                               // si il n'y est pas, c'est la division entière qui est effectuée (le résultat est un entier)
	printf("La simulation donne pour valeur de pi : %lf.\n",pi);
	return 0;
}


Vous remarquerez que je n'ai pas initialisé rand avec srand ; la suite de nombres générés sera donc toujours la même pour les différentes exécutions du programme. Cela n'a aucune importance ici et cela permet même de comparer les programmes. Par exemple pour 200 000 lancers (plus que Titi ne pourra jamais en faire de toute sa vie), j'obtiens 3.139080, ce que vous devriez obtenir également.
Retour en haut
Sans le savoir, vous venez d'implémenter votre première méthode de Monte Carlo ! Vous avez pu vous en apercevoir par vous-mêmes, cette méthode est terriblement simple, le cœur même de l'algorithme tient en 4 lignes. En revanche, elle est relativement peu efficace pour calculer pi. Effectivement, la théorie des Séries de Fourier permet de montrer que :
\frac{1}{1}+\frac{1}{4}+\frac{1}{9}+\frac{1}{16}+ ... +\frac{1}{n^2}+ ... = \frac{\pi ^2}{6}


Grâce à cette formule, en s'arrêtant à n = 100, on a 3.1302 comme valeur de pi, puis avec n = 3000, on obtient 3.14127, et les temps de calcul restent imperceptibles sur un ordinateur actuel. Monte Carlo est donc largement dépassé, d'autant qu'il existe d'autres formules donnant pi encore plus rapidement. Ceci étant, personne n'aurait le courage (moi y compris) de justifier et d'expliquer ces méthodes sur le SdZ. Justifier Monte Carlo est relativement simple (du moins dans l'idée), et c'est cette simplicité qui va être un des atouts de cette méthode.

Cette simplicité de l'algorithme a son importance, car même si les méthodes classiques sont plus performantes, il vaut mieux passer 4 heures de moins à coder plutôt que gagner 3 minutes à l'exécution. Ce qui compte c'est la rapidité (temps de codage compris) avec laquelle on obtient le résultat (évidement cette remarque n'a aucun sens dans un contexte "temps réel").

Mais une méthode n'est en général pas choisie uniquement pour sa simplicité, mais aussi pour son efficacité. Ici, le problème était en réalité trop simple pour qu'il soit intéressant d'y appliquer cette méthode. Pour faire une analogie avec le bricolage : pour visser des grosses vis, vous prenez un gros tournevis ; si vous prenez le même tournevis pour des toutes petites vis, vous allez avoir beaucoup de mal et ne serez pas très efficace. C'est exactement le même problème ici : Monte Carlo est un très gros tournevis, et pi, une toute petit vis !

Nous allons donc voir dans quels cas les méthodes "classiques" (des PMT : Petits et Moyens Tournevis) ne peuvent pas résoudre efficacement un problème et donc dans quels cas nous devrons utiliser Monte Carlo. Mais pour cela, il faut comprendre comment fonctionnent ces méthodes "classiques" !
Chapitre précédent Sommaire Chapitre suivant
Retour en haut

Partager

Retour en haut

Lire aussi

Icône
A la découverte de l'aléatoire et des probabilités
Icône
Simuler une distribution de cartes géométrique
Icône
Nombres flottants et processeurs
Icône
Savez-vous intégrer ?
Retour en haut
9 commentaires pour "Calculons pi !"
Note moyenne : 3.91 / 4 (32 votes)
Pseudo Commentaire
Hors ligne Rokil # Posté le 28/10/2009 à 14:41:28
Avatar
close

Bonjour, je débute dans la programmation en C, et je voudrais savoir comment on fait pour afficher le nombre d'essais et la valeur de Pi au cours de la simulation.
Un truc du genre :
Code : C
1
2
3
//Code...//
printf("%d ème essai , pi ≈ %ls", compteur, pi);
//Code...

Et le conde console devrait ressembler à ça :
Code : Console
Combien de jets voulez vous effectuer ? 1000
1 ème essai, pi = 2
2 ème essai, pi = 3,12 
...


Donc si vous pouviez m'aider, ce serait gentil ^^
                    
                        
Viendez sur #ircduzero, le chan qu'est bô !

Un jour, Torgi m'a mis 10% pour avoir complimenté une news. (cc @sdzfail) 



Compaq Presario 2100 :

-mobile AMD Athlon XP2500+

-192Mo RAM

-ATI Mobility Radeon

-SAMSUNG HM121HC (120GB)

                        
 

                    
                		
                

              
                

                  
                    
                      Hors ligne 
                    
                      sebsheep
                    
                  
                  
                    #
                    
                         Posté le 28/10/2009 à 17:39:57
                    
                    
                


                
                    

                
                
                    

Avatar

    

    Avis : Très bon
        
    



     Études : Universite Paris Sud 11






                		
                
                    
                    Je vais pas te donner la solution car c'est vraiment simple. Mais quand même des indication

Tu as la boucle ici :

Code : Autre
1
2
3
4
5
6
7
8
9
10
11
12

        for(i=0 ; i< n; i++)
        {
                // Jet de la flechette
                x = rand_0_1();
                y = rand_0_1();
                
                // La flechette est-elle dans le disque ?
                if( x*x + y*y <= 1 )
                {
                        compteur++;
                }
        }


Et juste après un petit calcul pour trouver pi et l'affichage :

Code : Autre
1
2

        pi = compteur*4 / (double) n;
        printf("La simulation donne pour valeur de pi : %lf.\n",pi);




Alors maintenant pose toi les questions suivantes dans l'ordre : 

A quel moment (endroit) dois-je insérer du code pour faire ce que je veux ?

Comment adapter l'affichage de fin de programme pour faire ce que je veux ?



N'hésite pas à me dire si tu bloques, mais réfléchis-y, c'est un très bon exercice.



</code>
                    
                        
"Il y a deux choses infinies dans le monde : l'univers et la bêtise humaine ... mais pour l'univers j'ai un doute" 


Image utilisateur
 Tuto sur l'aléatoire et les probabilités.

Un aperçu de Monte Carlo



                        
 

                    
                		
                

              
                

                  
                    
                      Hors ligne 
                    
                      Rokil
                    
                  
                  
                    #
                    
                         Posté le 29/10/2009 à 10:17:09
                    
                    
                


                
                    

                
                
                    

Avatar
	
	

		 
	




	
	close







    





                		
                
                    
                    Ok, merci, ça marche !



Effectivement, quand je mettais : 

Code : C
1
2
3
4
5
6
7
8
9
compteur++;
			
			printf("%lf\n",pi);
			
		}
	}


	pi = compteur*4 / (double) n;







Ca ne marchait pas ^^ ...



En tout cas, merci de m'avoir aidé !
                    
                        
Viendez sur #ircduzero, le chan qu'est bô !

Un jour, Torgi m'a mis 10% pour avoir complimenté une news. (cc @sdzfail) 



Compaq Presario 2100 :

-mobile AMD Athlon XP2500+

-192Mo RAM

-ATI Mobility Radeon

-SAMSUNG HM121HC (120GB)

                        
 

                    
                		
                

              
                

                  
                    
                      Hors ligne 
                    
                      Arthurus
                    
                  
                  
                    #
                    
                         Posté le 29/05/2010 à 14:06:48
                    
                    
                


                
                    

                
                
                    

Everyday I'm shuffling
Avatar
Validateurs


     Études : Ensimag






                		
                
                    
                    Salut.

À un moment donné tu dis que la surface du quart du cercle c'est Pi*R²/4, puis après tu supposes que R=1 pour continuer tes calculs.

Il n y a pas besoin de faire cette supposition, car en effet, quelque soit R, p = Pi/4, car en fait la surface totale sera égale à R², donc quoi qu'il en soit, on divisera toujours la surface du quart du cercle (qui est Pi*R²/4) par R², ce qui donne Pi/4.
                    
                        

                        
 

                    
                		
                

              
                

                  
                    
                      Hors ligne 
                    
                      Arthurus
                    
                  
                  
                    #
                    
                         Posté le 29/05/2010 à 14:25:49
                    
                    
                


                
                    

                
                
                    

Everyday I'm shuffling
Avatar
Validateurs


     Études : Ensimag






                		
                
                    
                    Voila en fait mon prog (qui avec ma remarque, devient beaucoup plus optimisé que le tien) :

Code : C
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int main (void) 
{
	//srand(time(NULL));
	long int i;
	long int nb = 200000;
	long int compt = 0;
	long int R = 10000;
	long int R_carre = R*R; // Pour ne pas le recalculer à chaque tour de boucle.
	for (i = 0; i<nb; i++) {
		int aleax = (rand() % (R+1));
		int aleay = (rand() % (R+1));
		if (aleax*aleax + aleay*aleay <= R_carre)
			compt++;
	}
	double pi = (double)compt * 4 / (double)nb;
	printf("pi = %lf\n",pi);
	return 0;
}





Code : Console
arthurus@Arthurus-Laptop:~/sdz$ ./toto 
pi = 3.141940

                    
                        

                        
 

                    
                		
                

              
            
        

        
            
Voir tous les commentaires
        


    

    

Ce tutoriel a été corrigé par les zCorrecteurs.