T’es sûr du résultat ?

Piste rouge Le 5 mai 2015  - Ecrit par  Aurélien Alvarez Voir les commentaires (22)

Les simulations numériques jouent un rôle très important dans certaines branches des mathématiques et seront l’objet de ce mois. Plus exactement, nous regarderons deux petits exemples de nature bien différente... même si bien sûr il y aurait infiniment plus à dire pour évoquer les nombreux aspects de cette branche des mathématiques au service de presque toutes les sciences aujourd’hui.

L’expérimentation numérique en mathématiques et plus généralement dans les sciences est un sujet passionnant et un outil fort utile, devenu indispensable pour certains scientifiques [1]. Malgré la puissance vertigineuse de calcul de nos ordinateurs aujourd’hui, et encore plus de certains centres de calculs, on aurait tort d’oublier complètement la théorie et de trop se moquer de comment fonctionne la machine, au risque d’avoir quelques surprises... [2]

Une équation différentielle pour commencer

Les équations différentielles sont omniprésentes en sciences tant elles sont utilisées pour modéliser des phénomènes physiques, biologiques, économiques, etc. Les mathématiciens les étudient aussi pour elles-mêmes et ce depuis l’invention du calcul différentiel-intégral au XVIIème siècle, jonglant entre un point de vue théorique, une approche qualitative du comportement des solutions ou des méthodes numériques permettant de calculer des valeurs approchées des solutions.

Voyons un petit exemple, l’équation différentielle suivante :
\[x' = t^2 + x^2.\]

Il s’agit d’une équation différentielle du 1er ordre dont l’inconnue est une fonction $x$ de la variable réelle $t$. La notation $x'$ désigne la dérivée de la fonction $x$ par rapport à $t$. On pourra par exemple penser à la fonction $x$ comme la position d’un mobile et à la variable $t$ comme le temps : $x'$ s’interprètera alors comme la vitesse du mobile. [3]

On qualifie cette équation de non-linéaire du fait de la présence du terme $x^2$ dans le second membre. Une telle équation est également dite non-autonome puisque la variable $t$ apparaît dans le second membre de l’équation : c’est donc une équation de la forme $x' = f(t,x)$. À part ça, la fonction $f$ qui entre en jeu est une fonction plutôt sympathique puisqu’il s’agit d’une fonction polynomiale de deux variables : $f(u,v) = u^2 + v^2$ [4].

D’après la théorie, si on se donne une condition initiale, il existe une unique solution qui part de cette condition initiale. C’est le fameux théorème de Cauchy-Lipschitz, théorème illustré en images dans le chapitre 2 de Chaos.

PNG - 319.6 ko

On se propose de calculer la valeur de la solution $x$ au temps $t = 1$, étant donné la condition initiale $x(0) = 1$ au temps $t = 0$. Voilà une question tout à fait raisonnable à laquelle n’importe quel ordinateur devrait pouvoir apporter une réponse en une fraction de seconde. Du moins une solution approchée avec une grande précision. C’est ce que nous allons essayer de faire.

GIF - 9.8 ko
Graphique illustrant la méthode d’Euler

La méthode d’Euler est probablement la plus simple des méthodes de résolution numérique des équations différentielles. Cette méthode fournit une suite d’approximations de la solution. On se donne un pas d’intégration $h$ [5] et, connaissant une valeur approchée $x(t)$ de la solution au temps $t$, on calcule une valeur approchée de la solution au temps $t + h$ par la formule suivante :
\[x(t + h) = x(t) + x'(t) \cdot h,\]
où $x'(t) = t^2 + x(t)^2$ dans notre cas.

Essayons d’implémenter tout ça dans un langage de programmation. Lequel ? Disons Python qui est un langage de programmation très utilisé, multiplateforme dont la syntaxe est particulièrement simple.

## Équation différentielle
def equaDiff(t,x):
   return t*t + x*x

## Schéma d'Euler
def euler(t,x,h):
   return x + equaDiff(t,x)*h

## Calcul d'une trajectoire
def trajectoire(tempsDepart,tempsArret,condInit,h):
   t = tempsDepart
   orbite = [condInit]
   while(t < tempsArret):
       condInit = euler(t,condInit,h)
       orbite.append(condInit)
       t = t + h
   return orbite

########################

tempsDepart = 0
tempsArret = 1
condInit = 1
h = 0.1

orbite = trajectoire(tempsDepart,tempsArret,condInit,h)

print(orbite)

Commentons un peu ce programme. La première fonction equaDiff, c’est tout simplement notre fonction $f$. La deuxième fonction, c’est notre schéma d’Euler qui calcule donc la valeur au temps $t + h$. Enfin la fonction trajectoire calcule la trajectoire entre $tempsDepart$ et $tempsArret$ (c’est-à-dire $t=0$ et $t=1$ dans notre cas), étant donné une condition initiale ($x(0)=1$ pour nous) et un pas $h$. Les quatre paramètres du programme sont initialisés sous la ligne des #### et on stocke le résultat du calcul de la fonction trajectoire dans une liste appelée orbite. Enfin on demande à notre ordinateur d’afficher cette liste. Et on obtient :

[1, 1.1, 1.2220000000000002, 1.3753284000000003, 1.5734812207846565, 1.8370655360008539, 2.1995465143570643, 2.7193470012390955, 3.507831812553902, 4.802320215070421, 7.189548159877822, 12.45850843419808]

On note que le premier nombre de cette liste, c’est bien 1, la valeur de notre solution $x$ en $t=0$. La valeur cherchée $x(1)$ est le dernier élément de cette liste, à savoir environ 12.45. Nous avons donc répondu à la question :-).

PNG - 79.1 ko

L’exemple qui jette le doute

Regardons à présent un exemple que j’aime bien, dû à Siegfried Rump en 1988 [6]. Il s’agit de la fonction de deux variables suivante :

\[f(u,v) = (333.75-u^2) \cdot v^6 + u^2 \cdot (11 \cdot u^2 \cdot v^2-121 \cdot v^4-2) + 5.5 \cdot v^8 + u/(2 \cdot v)\]

qu’on souhaite évaluer sur le couple $(77617,33096)$. Ci-dessous le petit programme Python qui va bien :

def f(u,v):
   return (333.75-u**2)*v**6 + u**2*(11*u**2*v**2-121*v**4-2) + 5.5*v**8 + u/(2*v)

print(f(77617,33096))

et qui trouve la valeur 1.1726039400531787 sur mon ordinateur. Jusque là, pas de quoi s’affoler n’est-ce pas ? Mais imaginons que vous soyiez très courageux et que vous fassiez le calcul à la main. Après tout, il n’y a que des nombres rationnels qui interviennent dans cette fonction, des opérations élémentaires, donc aucun obstacle théorique à un calcul exact. Mais il est probablement plus instructif d’apprendre à utiliser un logiciel de calcul formel [7] plutôt que d’essayer d’élever à la main des nombres à la puissance 8... Quelle que soit votre méthode, vous trouverez une valeur exacte de $-\frac{54767}{66192}$, soit environ... −0.8273960599 !!

Non seulement mon ordinateur a calculé une valeur très éloignée du résultat mais en plus, le signe n’est même pas le bon ! Autant dire que le résultat calculé par mon ordinateur est complètement faux !! J’encourage ceux qui pensent que mon ordinateur est bon pour la poubelle à essayer sur leur propre machine :-). Il se peut d’ailleurs qu’ils trouvent quelque chose d’encore différent... Ce qui est sûr, c’est qu’il y a un truc bizarre qui se passe...

L’énigme du jour, laissée à la sagacité des lecteurs, ce sera donc d’essayer de comprendre la raison pour laquelle le calcul sur un ordinateur standard aboutit à un résultat aussi erroné, malgré une fonction $f$ donnée par une formule tout à fait explicite et somme toute relativement simple.

Retour à notre équation différentielle

L’exemple de Rump que nous venons de voir a de quoi jeter le doute sur nos esprits... En particulier, est-on bien sûr de la valeur de 12.45 pour $x(1)$ que nous avons calculée précédemment ?

Première remarque : si on relit la liste des valeurs approchées pour $x$ calculées par l’ordinateur, on remarque que cette liste contient 12 termes ! Or nous sommes partis de $t=0$ pour aller à $t=1$ par pas de 0.1. Pourquoi n’obtient-on pas une liste de 11 nombres ??

Modifions légèrement notre programme pour faire afficher les couples $(t,x(t))$. Oh surprise, on obtient :

[0, 1], [0.1, 1.1], [0.2, 1.2220000000000002], [0.30000000000000004, 1.3753284000000003], [0.4, 1.5734812207846565], [0.5, 1.8370655360008539], [0.6, 2.1995465143570643], [0.7, 2.7193470012390955], [0.7999999999999999, 3.507831812553902], [0.8999999999999999, 4.802320215070421], [0.9999999999999999, 7.189548159877822], [1.0999999999999999, 12.45850843419808]

On remarque que la variable $t$ ne prend pas successivement les valeurs 0, 0.1, 0.2, etc. jusqu’à 1.0 mais qu’il y a des erreurs d’arrondis : on voit apparaître des nombres comme 0.30000000000000004 ou 0.7999999999999999 au lieu des nombres 0.3 et 0.8 attendus. Dès lors, on comprend que le programme a calculé une valeur de plus par rapport à ce qui était attendu du fait que 0.9999999999999999 est bien strictement inférieur à 1 !

Bon soit, ça ne semble pas bien méchant. Juste que la base 10 qui nous est si familière n’est probablement pas la base préférée de notre ordinateur... Mais du coup, on en déduit que la valeur recherchée $x(1)$ est plutôt de l’ordre de 7.18 et non de 12.45 comme nous l’avions conclu un peu rapidement lors de notre premier essai !

Puisqu’on a fait tout ce travail et, quitte à pratiquer un certain doute hyperbolique à la lumière de ces derniers rebondissements, on se propose de relancer notre programme avec un pas plus petit pour calculer une valeur approchée de notre solution avec davantage de précision. Prenons un pas de 0.01, soit 10 fois plus petit. On s’attend donc à trouver une centaine de valeurs dans notre liste. Une fraction de seconde plus tard, notre ordinateur nous répond :

[1, 1.01, 1.02020, 1.03061, 1.04124, 1.05210, 1.06319, 1.07453, 1.08613, 1.09799, 1.11012, 1.12255, 1.13527, 1.14830, 1.16166, 1.17535, 1.18939, 1.20379, 1.21857, 1.23374, 1.24933, 1.26534, 1.28179, 1.29870, 1.31610, 1.33399, 1.35241, 1.37138, 1.39092, 1.41105, 1.43180, 1.45320, 1.47528, 1.49807, 1.52160, 1.54591, 1.57103, 1.59701, 1.62388, 1.65170, 1.68050, 1.710344, 1.74127, 1.77336, 1.80665, 1.84123, 1.87716, 1.91451, 1.95337, 1.99383, 2.03599, 2.07994, 2.12580, 2.17370, 2.22376, 2.27613, 2.33096, 2.38843, 2.44872, 2.51205, 2.57864, 2.64873, 2.72261, 2.80058, 2.88298, 2.97019, 3.06264, 3.16079, 3.26519, 3.37643, 3.49519, 3.62225, 3.75850, 3.90495, 4.06276, 4.23330, 4.41814, 4.61911, 4.83840, 5.07859, 5.34275, 5.634606, 5.95865, 6.32043, 6.72680, 7.18635, 7.71001, 8.31185, 9.01029, 9.82989, 10.80408, 11.97947, 13.42282, 15.23301, 17.56211, 20.65522, 24.93063, 31.15521, 40.87110, 57.58517, 90.75550]

Saperlipopette ! On pensait améliorer d’un petit epsilon notre résultat de 7.18 et on se retrouve avec une valeur $x(1)$ qui semble de plus en plus grande ! Bizarre bizarre... Recommençons avec un pas de 0.005, donc cinq fois plus précis encore et n’affichons que les derniers termes de la liste :

[..., 17.10690, 18.57440, 20.30377, 22.36936, 24.87572, 27.97419, 31.89148, 36.98138, 43.82410, 53.43152, 67.710867, 90.63942, 131.72175, 218.47972, 457.15156, 1502.09426]

1502.09 maintenant !! Un dernier essai avec un pas de 0.001 et on trouve :

[..., 172.47233, 202.21998, 243.11384, 302.21912, 393.55647, 548.44412, 849.23603, 1570.43881, 4036.71785, 20331.80986, 433714.30329, 188541811.18441, 35548203106509.33, 1.2636747441371875e+24, 1.5968738589701863e+45, 2.5500061214623346e+87, 6.5025312194953795e+171, inf, inf, inf, inf, inf, inf, inf, inf, inf, inf, inf, inf, inf, inf, inf]

Là c’est le bouquet, on excède la capacité de calcul de ma machine après avoir atteint des valeurs stratosphériques de $10^{171}$... Gloups !

Et si on essayait un autre schéma numérique que la méthode d’Euler pour intégrer notre équation. Par exemple l’une des méthodes de Runge-Kutta, disons celle d’ordre quatre [8].

Notre programme est le même, sauf qu’on change la fonction euler par la fonction rungeKutta.

## Équation différentielle
def equaDiff(t,x):
   return t*t + x*x

## Schéma de Runge-Kutta
def rungeKutta(t,x,h):
   k1 = equaDiff(t,x)
   k2 = equaDiff(t,x+k1*h/2)
   k3 = equaDiff(t,x+k2*h/2)
   k4 = equaDiff(t,x+k3*h)
   return x+(k1+2*k2+2*k3+k4)*h/6}}

## Calcul d'une trajectoire
def trajectoire(tempsDepart,tempsArret,condInit,h):
   t = tempsDepart
   orbite = [condInit]
   while(t < tempsArret):
       condInit = rungeKutta(t,condInit,h)
       orbite.append(condInit)
       t = t + h
   return orbite

########################

tempsDepart = 0
tempsArret = 1
condInit = 1
h = 0.01

orbite = trajectoire(tempsDepart,tempsArret,condInit,h)

print(orbite)

Pas la peine d’aller bien loin, avec un pas de 0.01, on est déjà dans le décor !

[..., 9.02087, 9.92347, 11.02514, 12.40004, 14.16437, 16.51119, 19.78639, 24.67730, 32.77133, 48.73988, 94.63913, 653.84819, 410496238508.723, 2.645061132895814e+151, inf]

Que se passe-t-il donc avec notre résolution numérique ? D’ailleurs est-ce bien un problème numérique ? Il est peut-être temps de se rappeler ce que dit précisément le théorème de Cauchy-Lipschitz...

Le théorème affirme en effet qu’il existe une unique solution passant en $t=0$ par la valeur $x_0=1$. Et que cette solution est bien définie sur tout un intervalle ouvert contenant $t=0$. Mais le théorème ne dit rien sur la taille de cet intervalle ; en particulier rien ne permet d’affirmer qu’il s’étende jusqu’à $t=1$, encore moins que la solution soit définie sur toute la droite réelle [9]. Notre solution peut parfaitement exploser en temps fini et c’est d’ailleurs ce qui se passe ici. On peut en effet démontrer que la solution explose quelque part entre $\pi/4$ et 1, ce que l’on visualise facilement quand on cherche à tracer le graphe de la solution calculée [10].

JPEG - 15.3 ko
Graphe de la solution passant par 1 à l’instant initial
JPEG - 14.3 ko
Zoom montrant l’explosion de la solution

En fait, un œil averti avait sans doute remarqué dès le début qu’il risquait bien d’y avoir un problème avant $t=1$. Pourquoi ? Regardons l’équation différentielle plus simple
\[x' = x^2.\]

Cette équation a le mérite qu’on peut donner des formules explicites pour ces solutions. Il est d’ailleurs très facile de vérifier que la fonction $g$ définie par
\[g(t) = \frac{1}{1-t}\]
vérifie bien $x' = x^2$ et $x(0)=1$. Et il est tout aussi immédiat de remarquer que la droite verticale $t=1$ est une asymptote de cette fonction. Du coup, on n’est pas vraiment surpris de voir la solution de notre équation différentielle $x' = t^2 + x^2$ exploser avant $t=1$ puisqu’elle est supposée croître encore plus vite que la fonction $g$ (à cause du terme $t^2$ qui est strictement positif sur $]0,1[$), ce que, graphiquement parlant, on vérifie facilement.

JPEG - 23.3 ko
Les solutions des équations $x'=x^2$ (rouge) et $x'=x^2+t^2$ (vert) passant par 1 en $t=0$

En guise de conclusion

La morale de ce qui précède, c’est qu’il ne faut jamais se fier trop vite au résultat du calcul d’une solution approchée. Surtout quand celle-ci n’existe pas !

La morale de l’exemple de Rump, c’est que connaître un peu le fonctionnement de sa machine (en particulier la représentation des nombres et la gestion de la mémoire) n’est pas forcément un luxe ! En plus d’être un sujet tout à fait intéressant en soi.

Post-scriptum :

Un merci amical à Jean-Marc Dellac pour son dessin pythonesque. Et merci également aux relecteurs Laurent Dietrich, alainfa, Quentin Gendron, Pierre Monmarché et Clément Caubel pour leurs remarques et critiques constructives qui m’ont permis d’améliorer ce texte.

Article édité par Serge Cantat

Notes

[1À ce propos, mentionnons un article (p. 41-50) de Jean-René Chazottes et Marc Monticelli retraçant une brève histoire de ce domaine dans les colonnes de la Gazette des mathématiciens de la SMF.

[2Dans le même esprit que cet article, on pourra également lire Un ordinateur est-il une parfaite machine à calculer ? de Jean-François Colonna.

[3En mécanique, on tombe généralement sur des équations du 2ème ordre à cause du principe fondamental de la dynamique $F = m.a$, qui exprime que les forces s’exerçant sur un mobile sont proportionnelles à l’accélération, la dérivée de la vitesse.

[4Remarquons que $f(\lambda\cdot x,\lambda\cdot y)=\lambda^2\cdot f(x,y)$ pour tout réel $\lambda$ ; l’équation différentielle est pour cela dite homogène.

[5ce nombre $h$ a donc vocation à être petit...

[6S. Rump. Algorithms for verified inclusion. Reliability in Computing, Perspectives in Computing, pages 109–126. Academic Press, New York, 1988.

[7par exemple Sage

[8Ces méthodes de discrétisation numérique pour calculer des approximations de solutions d’équation différentielle ordinaire datent des années 1900 et sont appréciées pour leur convergence rapide et leur grande stabilité.

[9Une solution suffisante pour assurer que la solution fût définie pour tout réel aurait été que la fonction $f$ fût lipschitzienne par rapport à la seconde variable. Malheureusement la fonction $x \mapsto x^2$ définie sur $\mathbf{R}$ est loin de satisfaire une telle hypothèse...

[10Des calculs numériques plus précis semblent montrer que c’est autour de 0.97 que la solution explose en fait.

Partager cet article

Pour citer cet article :

Aurélien Alvarez — «T’es sûr du résultat ?» — Images des Mathématiques, CNRS, 2015

Crédits image :

Graphique illustrant la méthode d’Euler - Article de Wikipédia sur la méthode d’Euler
img_14106 - Jean-Marc Dellac

Commentaire sur l'article

Voir tous les messages - Retourner à l'article

  • T’es sûr du résultat ?

    le 6 mai 2015 à 22:08, par FDesnoyer

    Merci pour cet article.
    Enseignant dans le secondaire, je me suis amusé à planter joyeusement des calculatrices en leur demandant d’évaluer des nombres dérivés : pour 0.2 on obtient 0.200000001 alors que le calcul « à la main » est trivial.

    De même, la méthode d’Euler classique donne des résultats aberrants avec y’=y / y(0)=5 (exemple soumis par l’ancien directeur de l’IREM de Clermont-Ferrand).

    Dans la même veine : u_0=0.3 et u_n+1=17u_n-4.8, n’importe quel élève de 1S vous montre qu’elle est constante. Excel n’aime pas du tout !!!

    (On pourrait aussi remarquer que la méthode d’Euler pour les EDS n’est pas améliorée par l’utilisation de termes supplémentaires dans le DL, théorème de 2004 dû à un mathématicien allemand dont le nom m’échappe.)

    Bref, vaste sujet ! Merci encore pour cet article !

    Répondre à ce message

Laisser un commentaire

Forum sur abonnement

Pour participer à ce forum, vous devez vous enregistrer au préalable. Merci d’indiquer ci-dessous l’identifiant personnel qui vous a été fourni. Si vous n’êtes pas enregistré, vous devez vous inscrire.

Connexions’inscriremot de passe oublié ?