Estimation d’aires avec la méthode de Monte Carlo
Activités #7

Ce TP/projet doit vous servir de support afin de rendre un rapport. Il est fortement conseillé de lire TOUT le TP/projet avant de le commencer.

Il s’agit d’un RAPPORT: vous ne devez pas seulement écrire du code, vous devez également le COMMENTER (à l’intérieur du code et à l’extérieur), faire attention aux NOMS des fonctions et des variables que vous utilisez (utilisez obligatoirement les noms donnés dans l’énoncé, choisissez des noms EXPLICITES pour vos variables et fonctions supplémentaires). Pensez à rédiger une introduction, une conclusion, et à faire des paragraphes pour vos explications.

Il est également essentiel de TESTER chacun des programmes que vous faites, d’expliquer vos tests, ainsi que les résultats que vous avez obtenu (ces résultats correspondaient-ils aux résultats que vous attendiez ? Vous devez faire part d’un esprit critique, aussi n’hésitez pas à mentionner si vous vous attendiez à d’autres résultats que ceux obtenus…). Faire des tests à l’intérieur d’un programme peut également vous aider à trouver vous-même votre erreur lorsqu’un programme ne fonctionne pas: pensez-y !

Si vous êtes bloqué à une question, parlez-en à vos professeurs (il est donc important de commencer le TP/projet le plus rapidement possible). Si, malgré tout, vous restez bloqués, expliquez ce qui vous a bloqué, n’hésitez pas également à exposer les idées que vous aviez pour la suite mais n’avez pas pu mettre en oeuvre. Certaines questions faciles ne nécessitent pas d’avoir fait la totalité des questions précédentes pour être faites: lisez donc bien l’intégralité du TP/projet. Le plus important n’est pas de faire la totalité du TP/projet, mais de faire BIEN ce que vous faites.

La qualité de la rédaction et la lisibilité du code seront prises en compte.

Le but de ce TP/projet est d’utiliser la méthode de Monte Carlo pour calculer une valeur approchée de l’aire d’un polygone. Afin de se familiariser avec la méthode, on se propose de l’utiliser dans un premier temps pour calculer une valeur approchée de π: cette première partie est indépendante du reste du TP/projet.

Ce TP/projet est long. Il est indispensable que:

• vous fassiez la totalité de la section 2 (méthode de Monte Carlo pour approximer π: estimation, simulations, calcul de l’erreur, représentation des résultats),
• vous compreniez la façon dont sont représentés les polygones (section 3: si cette partie vous bloque, vous ne pourrez pas faire la suite),
• vous fassiez la section 4 (dessin d’un polygone), qui est très facile et permet de tester si vous avez compris la section 3.

Pour la suite, il est à noter que:

• la section 5 (génération de polygones) n’est pas indispensable pour faire la suite du TP/projet, mais il est conseillé de l’utiliser à la fin des sections suivantes pour tester vos résultats,
• la section 6 est le coeur du TP/projet: la section 6.2 (définition de la fonction appartient) est difficile, et cette fonction est utilisée dans la section 6.3. Sans cette fonction, vous ne pourrez pas calculer l’aire approchée d’un polygone. TOUTEFOIS, cela ne vous empêche pas de:
  • répondre aux questions de la sous-partie « Principe » de la section 6.3 (ces questions ne sont pas liées à la programmation, elles sont là pour vérifier que vous avez compris le but du TP/projet),
  • définir une fonction dessin_simu pour laquelle les points sont tous de la même couleur, qu’ils soient à l’intérieur ou à l’extérieur du polygone.
• la section 7 (sur l’aire exacte) peut être faite indépendamment de la section 6 (hormis la devinette initiale),
• pour la section 8, il est indispensable d’avoir fait au préalable les sections 1, 6 et 7.

Il est important de tester les fonctions que vous programmez. De ce fait, le but du TP/projet n’est pas d’arriver à la fin de la section 8 (même si, évidemment, votre note risque d’être bien meilleure dans ce cas…), mais AVANT TOUT de faire BIEN ce que vous faites (commentaires, rédaction, tests, etc.) et de savoir répondre correctement aux questions faciles qui sont disséminées tout au long du TP/projet.

Le projet sera rendu dans une boîte de dépôt JALON sous la forme d’une archive ZIP contenant :

• Un fichier texte README comprenant les noms du binôme (NOM1 et NOM2) et quelques informations essentielles sur le projet.
• Le rapport au format PDF (au plus 10 pages) nommé NOM1_NOM2.pdf. Veuillez lire ce document pour préparer et écrire le rapport.
• Un dosser src contenant les fichiers source du projet commentés. Le dossier src doit au moins contenir un fichier main.R.

Le fichier main.R proposera une petite demonstration du projet. Attention, il est essentiel que le fichier main.R s’exécute sans erreur.

Ces instructions doivent impérativement être respectées.

Instruisez vous avec la page wikipedia sur PI. Comparer la qualité des estimations obtenues avec la méthode de Monte Carlo avec :

• les représentations fractionnaires,
• les séries,
• les produits,
• les formules mathématiques,
• …

N’oubliez pas qu’un ordinateur effectue des calcul approchés sur les nombres réels.

La suite du TP/projet utilise cette méthode dite de Monte Carlo, mais dans le but d’approximer l’aire d’un polygone. Un polygone est une figure géométrique plane, formée d’une suite cyclique de segments consécutifs. Les polygones considérés sont des polygones simples (c’est-à-dire dont les arêtes ne se croisent pas), par opposition au cas des polygones complexes qui n’est pas traité ici. Notez bien que les polygones considérés ne sont pas nécessairement convexes.

Dans la suite, on décrira un polygone à n côtés à l’aide d’une matrice de taille (n+1)*2. Plus pécisément, le polygone constitué des segments [(x_1,y_1);(x_2,y_2)], [(x_2,y_2);(x_3,y_3)], …, [(x_n,y_n);(x_1,y_1)] sera représenté par une matrice M telle que:

• M[i,1] = x_i et M[i,2] = y_i pour tout i compris entre 1 et n,
• M[n+1,1] = M[1,1] = x_1 et M[n+1,2] = M[1,2] = y_1.

On notera que:

• chaque ligne de la matrice M représente un sommet du polygone;
• l’ordre dans lequel sont données les lignes est important, car un segment du polygone correspond à deux points consécutifs (autrement dit deux lignes consécutives);
• la dernière ligne de la matrice est toujours identique à la première ligne: cette redondance d’information va permettre d’alléger le code par la suite.

Par exemple, on pourra utiliser le code suivant pour définir des polygones.

## Definition d'une fonction très utile
creer_polygone <- function (x,y) {
  matrix(c(x, x[1], y, y[1]), ncol=2,dimnames=list(c(), c("x","y")))
}

carre <- creer_polygone(c(10,10,90,90), c(30, 70, 70, 30))
## Une permutation cyclique des points donne le même polygone
carre <- creer_polygone(c(10,90,90,10), c(70, 70, 30, 30))
## En revanche, le code suivant ne définit pas un rectangle,
## mais un polygone dont les arêtes se croisent.
papillon <- creer_polygone(c(10,90,10,90), c(30,70,70,30))
## pour finir, voici un losange.
losange <- creer_polygone(c(50,10,50,90),c(30,50,70,50))

print(carre)

x  y
[1,] 10 70
[2,] 90 70
[3,] 90 30
[4,] 10 30
[5,] 10 70

Vous pouvez dessiner un polygone avec la fonction plot, lorsque celui-ci est donné sous la forme décrite précédemment.

plot(carre, type='l')
lines(papillon -1, type='b', col='firebrick')
lines(losange, type='l', col='darkblue')

Pour rappel sur les propriétés de ces polygones, il est conseillé de jeter un oeil à l’adresse: https://fr.wikipedia.org/wiki/Polygone_r%C3%A9gulier.

Définir une fonction reg_poly <- function(n, r=1) { ... } qui prend en argument un entier n, un réel strictement positif r (de valeur 1 par défaut), et qui renvoie un polygone p vérifiant:

• le polygone p a n côtés,
• il est inscrit dans un cercle de centre (0,0) et de rayon r.

Indications:

• Tous les sommets d’un polygone régulier peuvent être engendrés à partir d’un seul sommet et les images successives d’une rotation: quel est l’angle de cette rotation ? Donnez-le en radians et en fonction du nombre n de côtés du polygone.
• Quelles sont les coordonnées cartésiennes (x,y) d’un point dont les coordonnées polaires sont (r,theta)?
• Tester votre fonction reg_poly en dessinant un exemple de polygone régulier (vous pouvez utiliser la fonction dessin_polygone définie précédemment).

Quelle figure obtient-on lorsqu’on définit le polygone suivant ?

x <- c(0,0,9,11,11,9,8,11,9,6,3,3,8,9,9,8,2,2)
y <- c(0,12,12,10,7,5,5,0,0,5,5,7,7,8,9,10,10,0)
surprise <- creer_polygone(x,y)
plot(surprise,col="black", type="l")

Définir une fonction appartient(points, polygone) qui prend en arguments des points et un polygone, et renvoie pour chaque point TRUE si le point est à l’intérieur du polygone, FALSE sinon. On pourra s’appuyer sur une fonction auxiliaire appartient_poly(point, polygone) qui prend en arguments un point et un polygone, et renvoie TRUE si le point est à l’intérieur du polygone, FALSE sinon. On ne se souciera pas du résultat renvoyé par la fonction dans le cas où le point appartient à un des côtés du polygone. En théorie, les points seront tirés uniformément aléatoirement et la probabilité qu’un tel cas se produise sera nulle.

Indication: Lorsqu’un point est à l’intérieur d’un polygone, toute demi-droite partant de ce point possède un nombre impair d’intersections avec les côtés du polygone. Lorsqu’il est à l’extérieur, elle en possède nécessairement un nombre pair.

Vous pourrez utiliser le code suivant pout tester la fonction appartient.

## Réaliser un test de la fonction
carre <- creer_polygone(c(0, 0, 1, 1), c(0, 1, 1, 0))
cc <- seq(from=-0.25,to=1.25,by=0.25)
points <- do.call(rbind,lapply(cc, FUN=cbind, cc,deparse.level = 0))
pin <- appartient(points,carre);
## Dessiner le résultat du test
par(mar=c(2,2,3,2)+0.1)
plot(carre, type='l', main="Test de la fonction appartient", xlim=range(carre[,1],points[,1]), ylim=range(carre[,2],points[,2]))
points(points[pin,1], points[pin,2], col='firebrick', pch=20)
points(points[!pin,1], points[!pin,2], col='darkblue', pch=20)

Deviner l’aire du carré :

print(mc.poly(10000, carre))

Vérifier par le calcul que votre résultat est satisfaisant pour le losange et la carré.

Définir une fonction resultats qui prend en argument un polygone p, et renvoie un objet de type data frame, contenant des simulations obtenues à partir du polygone p. Par exemple, chaque ligne du data frame peut contenir:

• un nombre n de points utilisés pour les simulations,
• plusieurs estimations de l’aire de p, chacune obtenue par une simulation avec n points,
• le temps moyen mis pour obtenir ces estimations,
• l’erreur relative pour chaque simulation.

Représenter vos résultats à partir du data frame.

Le but de cette section est de créer un fichier .csv, qui contiendra les résultats obtenus pour plusieurs polygones. Il est conseillé de définir et d’utiliser des fonctions, plutôt que de représenter les résultats un par un.

Choisir une dizaine de polygones. Créer un data frame contenant, sur chaque ligne:

• un polygone p (ou une référence à un polygone, tous les polygones étant par exemple stockés dans une liste externe),
• l’aire exacte de ce polygone,
• pour différentes valeurs de n (par exemple n=10**j, pour j=1:3):
  • la moyenne m de plusieurs estimations de l’aire de p, chaque estimation ayant été obtenue par une simulation avec n points,
  • l’écart-type entre ces estimations,
  • le temps moyen mis pour obtenir ces estimations,
  • l’erreur relative entre la moyenne m et l’aire exacte.

Exporter ce data frame dans un fichier .csv, afin de pouvoir le réutiliser par la suite.

Représenter certains de vos résultats (par exemple moyenne et écart-type) à partir de la lecture du fichier .csv.