T’es sûr du résultat ?

Piste rouge 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

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

    le 5 mai 2015 à 09:39, par B !gre

    Une autre conclusion qui peut être tirée de cet article est l’importance du calcul formel, à côté du calcul numérique. Par définition, le calcul formel consiste à effectuer des calculs de manière exacte, contrairement au calcul numérique qui travaille avec des valeurs approchées.

    Il n’est pas question de dire que l’un est mieux que l’autre. Le calcul numérique est souvent (pas toujours !) plus rapide que le calcul formel, mais le calcul formel permet de donner des réponses pour des problèmes numériquement instables.

    Répondre à ce message
    • T’es sûr du résultat ?

      le 5 mai 2015 à 11:22, par Aurélien Alvarez

      Je partage pleinement ce point de vue, comme je l’ai d’ailleurs rapidement écrit juste avant la note 6. Théorie, calcul formel, calcul numérique, etc. tout se complète et rien ne s’oppose !

      Répondre à ce message
    • T’es sûr du résultat ?

      le 7 mai 2015 à 11:10, par Jérôme Buzzi

      Le bug est partout où la main de l’homme a posé le pied, y compris dans les logiciels de calcul formel.

      Répondre à ce message
      • T’es sûr du résultat ?

        le 7 mai 2015 à 11:23, par B !gre

        En effet, écrire un logiciel sans bug est très difficile, en particulier pour des logiciels qui font des calculs sophistiqués comme les logiciels de calcul formel. Deux remarques cependant :

        1. On remarque dans la liste des bugs proposée dans le lien précédent un très grand nombre de bugs affectant des logiciels propriétaires. Il y en a plus en proportion que dans les logiciels libres dans cette liste (je ne dis pas que c’est le cas dans l’absolu, je pense que personne n’en sait rien). Ces bugs ont tendance à rester des années dans ces logiciels propriétaires, alors qu’ils sont corrigés beaucoup plus rapidement dans les logiciels libres comme SageMath déjà cité dans l’article. Par exemple, le premier bug de la liste affectant SageMath a été corrigé depuis longtemps. D’où l’importance d’utiliser des logiciels libres ! On peut aussi rêver du jour où on saura utiliser des techniques telles que la preuve assistée par ordinateur pour éliminer ces bugs. Ce n’est sans doute pas pour tout de suite malheureusement.

        2. Le type d’erreurs dont il est question dans l’article me semble (très) différent du type d’erreurs reportées dans la liste en lien. En effet, du côté des logiciels de calcul formel, il s’agit d’une liste de bugs qui ont donc vocation à être corrigés (et qui existent aussi bien sûr dans les logiciels numériques, et dans tout logiciel comme je l’écrivais au dessus). Les erreurs signalées par l’article sont beaucoup plus intrinsèques. Le calcul numérique ne peut par nature pas représenter l’ensemble des réels, et les approximations sont une nécessité. L’article montre le besoin de bien gérer ces approximations, au risque de trouver des résultats aberrants. Mais ce ne sont pas à proprement parler des bugs. Il s’agit plutôt d’une mécompréhension possible par un utilisateur du principe même du calcul flottant.

        Répondre à ce message
  • T’es sûr du résultat ?

    le 5 mai 2015 à 11:11, par le_cheveulu

    Petite question. Dans le programme Python que vous faites fonctionner, ne faudrait-il pas mettre la variable « pas » au lieu de « h » ? C’est peut-être une explication de cette incrémentation bizarre du pas.

    Répondre à ce message
    • T’es sûr du résultat ?

      le 5 mai 2015 à 11:17, par Aurélien Alvarez

      Merci c’est corrigé. J’avais en effet initialement appelé la variable « pas » dans mon programme mais un relecteur m’a suggéré de la nommer « h » comme on fait traditionnellement. Désolé pour cet oubli de correction.
      Bien cordialement, Aurélien Alvarez.

      Répondre à ce message
      • T’es sûr du résultat ?

        le 5 mai 2015 à 13:50, par le_cheveulu

        J’ai lancé votre programme et effectivement l’incrémentation en h crée des trucs bizarres. Vous avez une explication de ce qui se passe dans la machine à ce sujet ?

        Répondre à ce message
        • T’es sûr du résultat ?

          le 5 mai 2015 à 14:26, par Aurélien Alvarez

          Je ne comprends pas votre question. Qu’entendez-vous par « l’incrémentation en h crée des trucs bizarres » ?

          Répondre à ce message
  • T’es sûr du résultat ?

    le 5 mai 2015 à 16:10, par le_cheveulu

    Par incrémentation, j’entends le ligne de code t = t + h qui devrait donner les valeurs successives 0 ; 0.1 ; 0.2 ; 0.3 ; 0.4 etc... Je ne comprend pas ce qui se passe dans la machine pour donner 0.30000000000000004 ou 0.7999999999999999.

    Répondre à ce message
    • T’es sûr du résultat ?

      le 5 mai 2015 à 18:17, par Aurélien Alvarez

      Ah d’accord. Eh bien la réponse tient dans un nouvel exercice : écrire 0.1 en base 2...

      Répondre à ce message
      • T’es sûr du résultat ?

        le 6 mai 2015 à 09:58, par le_cheveulu

        Donc la machine passe en binaire pour faire simple addition. On peut donc aussi faire la remarque que le choix d’un langage de programmation a aussi sont importance. Si on avait programmé en C par exemple, il aurait fallut typer les variables (en « real » par exemple pour h) et nous n’aurions pas eu ce problème.

        Répondre à ce message
  • T’es sûr du résultat ?

    le 5 mai 2015 à 22:18, par Bruno Duchesne

    Merci Aurélien pour cet exemple simple. Dans le même genre, on pourra consulter une conférence grand public et un texte de Sylvie Bodo : « Pourquoi mon ordinateur calcule faux ».

    Répondre à ce message
  • T’es sûr du résultat ? Expérience vécue

    le 6 mai 2015 à 01:32, par Rémi Peyre

    En guise de commentaire, une petite expérience personnelle vécue pas plus tard que cette année : à l’examen du cours de statistique dont j’étais enseignant assistant, il s’agissait d’évaluer la variance due à un certain facteur. Dans son corrigé, l’enseignant principal effectue un calcul avec 6 chiffres de précision. (Dans le détail, il manipulait les quantités suivantes, correspondant à des tailles d’élèves : $181,\!100$, $179,\!647$, $178,\!641$, $177,\!867$ et $179,\!012$ — on m’objectera que les deux premiers chiffres sont toujours essentiellement les mêmes, de sorte que la précision réelle est plutôt de 4 chiffres ; mais cela ne change pas grand-chose à la suite). (Je précise encore que, notant respectivement $y_1, \ldots, y_5$ ces valeurs, le calcul consistait à évaluer $10 y_1^2 + 17 y_2^2 + 39 y_3^2 + 15 y_4^2 - 81 y_5^2$). Puis il écrit sur son corrigé : « on trouve $86,\!58$, mais si on n’avait pas arrondi les données, on aurait trouvé $75,\!49$ ». Je saute alors aussitôt sur mon clavier pour lui faire remarquer qu’il s’est forcément trompé quelque part, car pour un calcul aussi simple, on ne peut pas passer comme cela d’une précision de quatre chiffres à une précision de même pas un chiffre... Eh bien si !! Piégé par l’habitude que les problèmes d’arrondis soient rares en pratique (et presque jamais très supérieurs à la précision des données initiales), j’avais cru à tort qu’il n’y en aurait pas ici ; mais c’est bien l’enseignant principal qui avait raison : la différence s’expliquait uniquement par les arrondis...!

    Addendum : Sachant que $y_5$ était en fait défini comme $(10 y_1 + 17 y_2 + 39 y_3 + 15 y_4) \mathbin{/} 81$, on pouvait montrer que la formule calculée correspondait aussi à $10 (y_1 - y_5)^2 + 17 (y_2 - y_5)^2 + 39 (y_3 - y_5)^2 + 15 (y_4 - y_5)^2$. Ce qui est intéressant alors, c’est qu’avec cette formule-là, les problèmes d’arrondis auraient eu beaucoup moins d’impact...! Moralité : la rapidité n’est pas le seul critère de choix entre deux formules de calcul numérique ; la résistance à la propagation des erreurs compte aussi...! (L’auteur lui-même évoquait d’ailleurs déjà cette idée quand il parlait de la méthode de Runge-Kutta).

    .

    Cordialement,

    Répondre à ce message
    • T’es sûr du résultat ? Expérience vécue

      le 9 mai 2015 à 13:17, par Aurélien Alvarez

      Merci Rémi pour cet exemple tout simple ; je viens d’essayer à mon tour et effectivement la première formule donne 86,58 alors que la deuxième 75,49. Un peu inquiétant tout ça s’il s’avérait que de pareilles erreurs ne soient pas contrôlées dans certains systèmes critiques (en aviation par exemple pour ne citer qu’un exemple)...

      Répondre à ce message
  • 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
  • T’es sûr du résultat ?

    le 9 mai 2015 à 11:31, par Jean-François Colonna

    Je me permets de rappeler mon article « Un ordinateur est-il une parfaite machine à calculer ? » publié en 2010 qui montrait déjà ces « anomalies », en particulier sur des exemples beaucoup plus simples. Sur mon site, on peut voir de nombreuses expériences montrant l’influence des méthodes d’intégration utilisées, mais aussi celles des options de compilation (je rappelle qu’à cause des erreurs d’arrondi, l’addition et la multiplication ne sont plus associatives).

    Répondre à ce message
    • T’es sûr du résultat ?

      le 9 mai 2015 à 13:03, par Aurélien Alvarez

      Bonjour Jean-François, ton article m’avait échappé et je viens de rajouter une note dans mon texte pour le mentionner. Bien cordialement, Aurélien.

      Répondre à ce message
  • T’es sûr du résultat ?

    le 10 mai 2015 à 23:01, par Pierre Lecomte

    J’ai lu avec beaucoup d’intérêt ce bel article !

    Ce qui est amusant, à propos de l’équation x’=x^2 dont les solutions (autres que celle pour laquelle x(0)=0) présentent une singularité (i.e. ne sont pas définies sur tout R), c’est que cette singularité est artificielle.

    D’abord, avec le changement de variable y=1/x, on transforme l’équation en y’=-1 dont les solutions sont définies sur tout R --- mais, d’accord, le changement de variable est singulier en l’origine, ce qui n’est peut être pas une bonne chose.

    Il existe cependant un cadre, plus abstrait sans doute, où toute singularité disparaît et dans lequel on comprend pourquoi elles se créent lorsqu’on revient à l’équation initiale x’=x^2.

    En deux mots : sur la droite projective réelle PR qui est compacte, tout champ de vecteur X est complet, i.e. ses courbes intégrales maximales sont définies sur R tout entier. La droite à deux cartes canoniques et, dans celles-ci, l’expression locale d’une courbe intégrale de X est solution une équation différentielle x’=f(x), où f est l’expression locale de X. Une courbe intégrale maximale de X est définie sur R mais ses expressions locales ne le sont peut-être pas car la courbe en question n’est pas nécessairement entièrement contenue dans le domaine de la carte où on l’étudie. Si c’est le cas, les solutions de l’équation x’=f(x) présentent alors des singularités en des valeurs de t qui correspondent au moment où la courbe quitte le domaine de carte.

    C’est ce qui se passe avec l’équation x’=x^2. Il existe un champ de vecteur X de PR dont les expressions locales dans les deux cartes canoniques sont x^2 et -1 ! Ses courbes intégrales maximales sortent toutes, sauf une, du domaine de la première carte, générant ainsi les singularités des solutions de x’=x^2, mais toutes sont contenues en entier dans le domaine de la seconde carte ! Naturellement, vous l’aurez soupçonné, le changement de coordonnées locales, lorsqu’on passe d’une carte à l’autre, est précisément la fonction y=1/x...

    En fait, toutes les équations de la forme x’=a+bx+cx^2 où a,b,c sont constants donnent lieu au même phénomène : elles sont des expressions locales de champs de vecteurs sur PR et les singularités éventuelles de leurs solutions ne sont dues qu’au passage en coordonnées locales. Tout étant explicite, on peut même trouver les courbes intégrales maximales des champs de vecteurs en question et en déduire les solutions de ces équations !

    Je ne vois pas encore clair à propos de l’équation x’=t^2+x^2 qui n’est pas autonome. Mais il se pourrait bien qu’elle donne lieu au même genre de considérations.

    Bien cordialement,

    Pierre Lecomte

    Répondre à ce message
    • T’es sûr du résultat ?

      le 11 mai 2015 à 11:14, par Nils Berglund

      Juste une petite remarque technique : les équations du type $x’=a(t)+b(t)x+c(t)x^2$ sont appelées équations de Riccati, et on sait les transformer en équations linéaires du second ordre, voir ici. Les singularités des solutions de l’équation de Riccati correspondent alors aux zéros des solutions de l’équation du second ordre. C’est un phénomène bien connu des physiciens, car il apparaît notamment dans l’étude des points tournants de l’équation de Schrödinger stationnaire. L’exemple le plus simple est celui de l’équation d’Airy, qui correspond à l’équation de Riccati $x'(t)=x(t)^2+t$.

      Répondre à ce message
      • T’es sûr du résultat ?

        le 11 mai 2015 à 11:28, par Pierre Lecomte

        Merci beaucoup pour ces informations ! Cordialement,

        Pierre Lecomte

        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é ?

Suivre IDM