Dans une autre page, nous avons expérimenté plusieurs méthodes d'intégration numérique, des plus simples aux plus précises. Certaines, comme la méthode des trapèzes apparaissent maintenant dans les cours sur l'intégration en terminale S. Les méthodes de Simpson et de Romberg sont utilisées classiquement en physique numérique. Alors pourquoi chercher une autre méthode ?
Toutes les méthodes citées ci-dessus partent de l'hypothèse que l'aire à calculer (la fonction à intégrer) n'est pas trop tourmentée et que la surface peut être facilement découpée en rectangles et autres trapèzes. Si l'aire est tarabiscotée, ces méthodes sont difficilement applicables, et encore plus lorsqu'il s'agit de volumes ou de dimensions supérieures ! Il fallu trouver autre chose. La solution vint des probabilités et de l'exploitation de l'aléatoire. Elle repose sur la méthode de Monte-Carlo, universellement exploitée dans le monde de la simulation pour sa puissance et son efficacité. Nous aborderons ici son utilisation pour réaliser une intégration, sachant que l'on peut faire des tas d'autres choses avec, comme simuler la désintégration d'une population de noyaux radioactifs, le refroidissement d'un verre, le comportement d'un gaz d'électrons, etc...
Notons toutefois un point important: la méthode de Monte-Carlo est une méthode moins précise que les méthodes d'intégration comme Simpson et Romberg (pour ne citer qu'elles), et même que la simple méthode des trapèzes... Ce manque de précision est du à son principe, mais aussi à l'usage obligatoire d'un générateur de nombres aléatoires, dont la qualité influe sur la précision du calcul. Or vous savez quoi en penser après la lecture de la page qui leur consacrée...
Cette méthode trouve toute sa puissance dans le calcul des intégrales multiples. En physique, on peut être amené à calculer des intégrales sur un nombre important de dimensions, par exemple dans les systèmes n-corps. On démontre que lorsque n dépasse 10, il n'existe aucune autre méthode connue qui soit aussi efficace que la méthode de Monte-Carlo.
Enfin, disons qu'il serait dommage de cantonner la méthode de Monte-Carlo au calcul d'intégrales ! C'est une méthode inspirée de la physique statistique, qui offre un moyen très puissant pour comprendre beaucoup de phénomènes dans tous les domaines de la physique et même au delà.
Pour la petite histoire, le nom de la méthode "Monte-Carlo" fut inventé par Nicholas C. Metropolis, également inventeur d'un algorithme éponyme célèbre d'amélioration de la méthode de Monte-Carlo. Il paraît que le terme lui fut soufflé par Stanslaw Ulam, célèbre topologiste, qui utilisait cette méthode pour ses travaux de physique appliqués à la bombe H.
Pour aborder le principe de l'intégration par la méthode de Monte-Carlo, j'emprunterai un problème que j'aime bien, à Rubin Landau (bien connu des physiciens numériciens) dans "Computational Physics".
Imaginez que vous soyez fermier et que vous vouliez mesurer la surface d'un étang dans un de vos champs. Les contours de cet étang sont loin d'être réguliers ! Et vous ne disposez que d'un tas de cailloux , d'une chaine d'arpenteur et de votre intelligence... Voici le moyen proposé par Landau (qui, semble-t-il, a emprunté cet exemple à H.Gould):
Bon, raisonnons maintenant: si vos jets ont été suffisamment aléatoires, il semble évident que le nombre de cailloux Se tombés dans l'eau sera proportionnel à la surface de l'étang et que l'on peut écrire Se / (Se + Ss) = Ae/Ac. Ce qui nous donne naturellement Ae = Ac*Se/(Se + Ss). Nous avons trouvé un moyen d'estimer de manière assez précise l'aire de l'étang !
La manip décrite ci dessus correspond presque parfaitement à l'algorithme de la méthode de Monte-Carlo, que je détaillerai ci-dessous: un lancé aléatoire, une condition par rapport à la surface à intégrer (ici le splash du caillou qui tombe dans l'eau), des sommations et un rapport.
On peut aussi aborder ce principe par une autre voie. Vous connaisez sans doute, c'est dans le cours de TS sur les intégrales, la définition de la valeur moyenne d'une fonction <f>. Elle s'écrit:
\begin{align} \int_a^bf(x)dx = (b - a) <f> \end{align}
Cela revient à dire que la valeur de l'intégrale de f(x) entre a et b est la surface du rectangle de longeur ab et de largeur <f>, la valeur moyenne de f entre a et b. L'algorithme de Monte-Carlo permet d'estimer <f> et ainsi de calculer la valeur de l'intégrale.
Nous allons nous limiter à l'intégration d'une fonction de dimension 1, c'est à dire à l'évaluation de l'aire comprise entre cette courbe et l'axe des abscisses pour un intervalle donné sur Ox. Il vous sera facile de généraliser le programme pour une fonction à deux variables ou plus, Maple étant particulièrement approprié à ce genre d'exercice.
L'algorithme proposé est fort simple:
j'ai choisis, par simplicité, d'utiliser une boucle à nombre constant d'itérations (N). Ce nombre détermine la précision du calcul, comme la qualité du générateur de nombres aléatoires. Pour obtenir de bons résultats, il faut donc tâtonner un peu ! Il existe une autre méthode : utiliser une boucle à condition déterminée par le calcul. Ici par exemple, on pourrait décider de sortir de la boucle lorsque le résultat du calcul d'aire ne changerait plus de façon significative. C'est généralement ce qu'on fait, comme par exemple dans le schéma de Romberg, exposé dans la page sur l'intégration numérique, où l'on fixe la précision souhaitée.
Le code source est disponible ici. Il reprend l'algorithme décrit ci-dessus, avec toutefois quelques ajouts:
J'ai utilisé le générateur de nombres aléatoires de Maple, qui génère des nombres aléatoires compris entre 0 et 1012 selon une distribution normale. Il existe sans doute mieux, mais pour ce que j'en fait ici, cela ira bien... En cas d'utilisation sérieuse, il faudrait revoir ça.
J'ai retrouvé dans le bestiaire des fonctions un peu bizarres cette fonction qui correspond bien au besoin de l'exemple:
\begin{align} f(x) = x(1 - x)sin^2(200x(1 - x)) \end{align}
dont voici la courbe représentative sur l'intervalle [0,1] :
Pas facile à découper en trapèzes ou en rectangles ! Et Maple dit qu'il existe une solution analytique (une somme de fonctions de Fresnel), ce qui nous permettra de comparer l'aire calculée par Maple avec celle calculée par Monte-Carlo.
Lançons le programme MonteCarloMaple.mpl, que vous aurez téléchargé, et voyons ce que cela donne pour N = 2 000.
Lors du premier essai, j'obtiens una aire calculée de 0.080500. Maple me fournit la solution analytique et une valeur numérique de 0.080498.
J'ai effectué plusieurs essais et et obtenu des résultats qui fluctuent : 0.078250, 0.078625, 0.075500, 0.079375, etc.. La précision pour 2 000 points n'est pas mauvaise, et le temps de calcul court ! En utilisant 5 000 points (modifiez la valeur de N dans le code), on constate que la précision n'est pas vraiment plus grande, alors que le temps de calcul est plus long. Idem pour 10 000 points. Il existe d'ailleurs un calcul qui permet d'estimer un majorant de la précision en fonction du nombre de points.
Je vous laisse faire vos expériences ! Essayez de trouver une fonction plus "tourmentée", ou encore une fonction à deux ou trois variables...
Un calcul hyperclassique consiste à évaluer la valeur de Pi en calculant l'aire sous la fonction qui va bien. Vous savez bien sur qu'il s'agit de :
\begin{align} f(x) = \sqrt(1 - x^2) \end{align}
que l'on va intégrer sur l'intervalle [0,1]. Le résultat est bien connu analytiquement et vaut \(\frac{\pi}{4}\).
Dans le programme MonteCarloMaple, il suffit de modifier la définition de f, en écrivant f := x -> sqrt(1 - x^2): et aussi les valeurs des bornes en posant xmin = ymin = 0 et xmax = ymax = 1.
Lançons le programme ainsi modifié et voyons ce que cela donne. Lors du premier essais, j'obtiens une aire de 0.784800, soit une valeur de pi = 3.1392. Au cours des essais suivants, j'obtiens des valeurs comprises entre 3.1212 et 3.1648, avec une valeur moyenne de 3.1484. Pas trop mauvais, mais pas vraiment génial... Je travaille avec 2 000 points. Avec 3 000 points c'est un peu mieux, mais guère plus !
La méthode de Monte-Carlo utilisée pour calculer la valeur d'une intégrale est une méthode simple mais relativement peu précise. Elle ne devra pas être utilisée si d'autres schémas numériques sont possibles (Simpson, Romberg ou leurs variantes). Mais elle est irremplacable pour calculer les intégrales multiples ou dans le cas de fonctions analytiques ou non, un peu tourmentées.
La méthode de Monte-Carlo est un pilier de la simulation numérique en physique statistique, nucléaire et dans bien d'autres domaines. Nous aurons sans doute l'occasion d'y revenir sur TangenteX.com !
Contenu et design par Dominique Lefebvre - tangenteX.com mars 2013 -- Vous pouvez me joindre par ou sur PhysiqueX
Cette œuvre est mise à disposition selon les termes de la Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Pas de Modification 3.0 France .