Les problèmes du calcul numérique avec Python : exemples et quelques explications
Cet article peut être librement diffusé et son contenu réutilisé suivant la licence CC-by-sa
Si vous tapez 0.1 + 0.1 + 0.1 - 0.3 dans la console Python, qu’attendez-vous en retour ? Sur un ordinateur standard, la réponse est 5.551115123125783e-17, autrement dit $5,551115123125783\,\times\, 10^{-17}$, ce qui n’est pas bien gros mais pas 0 quand même. Voilà déjà de quoi instiller un doute sur la fiabilité de l’utilité de Python pour les calculs, d’autant que n’importe quelle calculatrice, à la question précédente, donne la bonne réponse - nous y reviendrons.
Mais voici un premier exemple de ce qui peut se produire à cause de ces approximations.
Fixons deux entiers strictement positifs $p$ et $q$, et considérons la suite $(u_n)_n$ définie par son premier terme $u_0$ et la relation de récurrence $u_{n+1} = (q+1)\,u_n - p$. Si $u_0$ vaut $\dfrac{p}{q}$, en réduisant au même dénominateur on obtient immédiatement $u_1 = u_0$ et, en itérant, on voit que la suite est constante. Qu’en est-il avec Python ? Demandons-lui de calculer les 100 premiers termes, et la différence avec la valeur théorique ; nous n’afficherons qu’un terme sur 10, pour mieux voir. Voici le programme [1] :
- def suite(p,q):
- u = p/q
- for i in range(100): # i varie de 0 à 99
- u = (q+1) * u - p
- if (i+1) % 10 == 0 : # l'indice mathématique est décalé de 1
- print('Le terme de rang '+str(i+1)+' vaut : '+str(u))
- print(' écart avec la valeur théorique : ',u-p/q)
Et voici son exécution pour $p=2, q=3$ :
>>> suite(2,3)
Le terme de rang 10 vaut : 0.6666666666278616
écart avec la valeur théorique : -3.8805070268210784e-11
Le terme de rang 20 vaut : 0.6666259765625
écart avec la valeur théorique : -4.069010416662966e-05
Le terme de rang 30 vaut : -42.0
écart avec la valeur théorique : -42.666666666666664
Le terme de rang 40 vaut : -44739242.0
écart avec la valeur théorique : -44739242.666666664
Le terme de rang 50 vaut : -46912496118442.0
écart avec la valeur théorique : -46912496118442.664
Le terme de rang 60 vaut : -4.9191317529892135e+19
écart avec la valeur théorique : -4.9191317529892135e+19
Le terme de rang 70 vaut : -5.1580834970224175e+25
écart avec la valeur théorique : -5.1580834970224175e+25
Le terme de rang 80 vaut : -5.4086425609737785e+31
écart avec la valeur théorique : -5.4086425609737785e+31
Le terme de rang 90 vaut : -5.671372782015641e+37
écart avec la valeur théorique : -5.671372782015641e+37
Le terme de rang 100 vaut : -5.946865386274833e+43
écart avec la valeur théorique : -5.946865386274833e+43
C’est beaucoup plus gênant que le $10^{-17}$ précédent ! Si on creuse un peu, on voit que le problème n’est pas systématique : pour certaines valeurs de $p$ et $q$, la suite est bien constante. Voici quelques tests, où nous n’affichons que le terme de rang 100 :
>>> suite(1,2)
Le terme de rang 100 vaut : 0.5
écart avec la valeur théorique : 0.0
Cette fois tout est normal.
>>> suite(2,7)
Le terme de rang 100 vaut : -3.2308060613090455e+73
écart avec la valeur théorique : -3.2308060613090455e+73
On retrouve le problème précédent. Serait-ce dû au fait que notre premier terme n’est pas décimal ?
>>> suite(2,5)
Le terme de rang 100 vaut : 3.868423350539693e+61
écart avec la valeur théorique : 3.868423350539693e+61
Ah ben non, c’est plus compliqué que ça. Il va falloir se plonger dans la représentation des nombres en machine.
Le codage des nombres non entiers
Ce qui suit s’applique à Python sur un ordinateur ; les calculatrices n’ont en général pas le même problème, nous y revenons dans un autre onglet.
Commençons par un exercice qu’on peut donner en première (programme actuel ou futur) à condition sans doute de le décomposer : calculer $S= \sum_{n=1}^{+\infty}\,\dfrac{1}{2^{4n}}$ . En mettant $\dfrac{1}{2}$ en facteur (ou en ajoutant 1 qu’on retranchera ensuite), on fait démarrer la somme à $n=0$ et on peut utiliser la formule bien connue $1+q+q^2+\ldots +q^n = \dfrac{1-q^{n+1}}{1-q}$, ici pour $q=\dfrac{1}{2^4}$ ; puis en passant à la limite, on trouve $S = \dfrac{1}{15}$. Cette fraction peut donc s’écrire en binaire : $S = (0,000100010001001\ldots)_2$ , le bloc 0001 se répétant à l’infini, ce que nous conviendrons de noter $(0,\overline{0001})_2$.
Pour $S’= \sum_{n=1}^{+\infty}\,\dfrac{1}{2^{4n+1}}$ on peut utiliser un calcul analogue au précédent, ou plus économiquement remarquer que $S’$ est la moitié de $S$, donc $\dfrac{1}{30}$. Or en binaire, diviser par 2 revient à décaler la virgule vers la gauche, donc $S’$ s’écrit $(0,0000100010001001\ldots)_2$ c’est -à-dire $(0,0\overline{0001})_2$.
Ainsi $\dfrac{1}{10} = \dfrac{3}{30} = \dfrac{2}{30} + \dfrac{1}{30} = \dfrac{1}{15} + \dfrac{1}{30}$ s’écrit en binaire $(0,00011001100110011\ldots)_2$ c’est-à-dire $(0,0\overline{0011})_2$.
Voyons maintenant comment un ordinateur peut s’en accommoder. Ne pouvant stocker une infinité de chiffres, il doit les tronquer. Pour des raisons évidentes de compatibilité matérielle et logicielle, la façon de tronquer est standardisée, régie par une norme internationale, IEEE 754, bien décrite sur Wikipédia. Regardons ce que donne cette norme sur notre $0,1 = (0,00011001100110011\ldots)_2$ et un ordinateur qui travaille avec des « mots » de 64 bits, comme presque tous les ordinateurs personnels actuels. La virgule (flottante, c’est bien pour cela que les nombres - non entiers - des ordinateurs sont appelés « flottants ») est placée après le premier chiffre non nul, donc après le 1 qui est en quatrième position (quatrième après la « vraie » virgule, celle de la représentation ci-dessus). Ce chiffre 1-là n’est pas enregistré dans l’ordinateur : un chiffre binaire non nul, c’est forcément 1, donc inutile de perdre de la place pour ça. Les 52 chiffres suivants sont enregistrés, cela forme la « mantisse » de notre nombre ; tous les suivants sont ignorés. L’exposant, ici -4, est codé sur les 11 bits restants, sous une forme un peu compliquée qui n’a pas d’importance ici. En résumé, la décomposition binaire ci-dessus de 0,1 est tronquée juste après le 56-ième chiffre.
Vérifions-le... avec Python, mais prudemment cette fois !
- def undixieme(N):
- '''retourne (0,00011001100110011...)_2 tronqué au N-ième chiffre
- après la virgule'''
- puiss = 1/2
- i = 1
- S = 0
- while i < N :
- puiss = puiss / 2
- i += 1
- if i % 4 == 0 or i % 4 == 1 : # le chiffre de rang i est un 1
- S += puiss
- return S
Testons :
>>> 0.1 - undixieme(56)
1.3877787807814457e-17
>>> 0.1 - undixieme(57)
0.0
Ainsi, Python remplace tous les nombres (non entiers) par des valeurs approchées, sauf ceux qui ont en binaire un développement fini : à partir d’un certain rang il n’y a plus que des 0. Ces nombres sont appelés « dyadiques », ce sont les équivalents des décimaux quand on remplace la base 10 par la base 2. Ce sont les nombres qui peuvent s’écrire comme quotient d’un entier par une puissance de 2. Pour tous les autres, avec Python - comme avec tout système de calcul numérique - il y a approximation, donc risque d’erreur.
Une première conséquence : à partir du moment où on travaille avec des flottants, il faut absolument éviter de tester une égalité. Par exemple [4], selon la réciproque du théorème de Pythagore, le triangle de côtés $\sqrt{3}, \sqrt{5}$ et $2\sqrt{2}$ est rectangle. Mais si l’on teste :
>>> from math import sqrt
>>> a, b, c = sqrt(3), sqrt(5), 2*sqrt(2)
>>> a**2 + b**2 == c**2
False
Il faut donc remplacer un tel test par une estimation $|a^2+b^2-c^2| < \varepsilon$ pour un $\varepsilon > 0$ très petit.
Si on pouvait rester entier...
Il y a un domaine - parmi bien d’autres - où Python fait beaucoup mieux que les calculatrices, mais aussi que le tableur : le calcul sur les nombres entiers. Là, l’unique limitation de Python [5] est la taille de la mémoire de l’ordinateur utilisé.
Une première application vaut d’être signalée, même si elle n’est pas liée à notre sujet : le célèbre problème des grains de riz sur l’échiquier, dit aussi problème de Sissa. On peut le donner en collège et montrer ainsi les limites du tableur : avec une cellule par case de l’échiquier, les onze dernières contiennent des multiples de 10, alors que, par construction, on ne peut avoir que des puissances de 2. Alors que Python donne instantanément le bon résultat, même si on double la taille de l’échiquier, et au-delà.
Nous allons utiliser cette capacité de Python, pour creuser la question.
Revenons à notre problème de conversion du décimal au binaire, en reprenant la représentation de 0,1 en flottant : $0,1 = (0,00011001100110011...)_2$ tronqué à 56 chiffres après la virgule.
On va pouvoir déterminer l’approximation faite, en n’utilisant que des entiers grâce au module fractions
. Celui-ci définit un type « Fraction » sur lequel sont définies les quatre opérations usuelles, entre fractions ou avec des entiers ; la fonction float
permet d’obtenir une valeur décimale approchée d’un tel objet.
- from fractions import *
- def undixieme_dec(N):
- '''retourne (0,00011001100110011...)_2 tronqué au N-ième chiffre
- après la virgule, calculé sous forme de fraction'''
- puiss = Fraction(1,2)
- i = 1
- S = 0
- while i < N :
- puiss = puiss / 2
- i += 1
- if i % 4 == 0 or i % 4 == 1 : # le chiffre de rang i est un 1
- S += puiss
- return(S,float(S))
Testons :
>>> undixieme_dec(56)
(Fraction(7205759403792793, 72057594037927936), 0.09999999999999999)
>>> undixieme_dec(57)
(Fraction(14411518807585587, 144115188075855872), 0.1)
Pourtant, on l’a dit, Python tronque au 56ème chiffre, donc pour lui $0,1\simeq \dfrac{7205759403792793}{72057594037927936}$ !
On peut aussi laisser moins de travail à Python, et faire à la main, avec la formule déjà rappelée pour la somme des termes d’une suite géométrique, le calcul du nombre tronqué : $S= \sum_{n=1}^{14}\,\dfrac{1}{2^{4n}} + \sum_{n=1}^{13}\,\dfrac{1}{2^{4n+1}}$ puisque le dernier rang est $56 = 4 \times 14$ . On trouve $\dfrac{2^{56}+2^{55}-9}{15\times 2^{56}}$ qui après simplification par 15 redonne la fraction ci-dessus.
Nous savons donc maintenant quel est le nombre dyadique qui, pour Python, remplace 0.1, et nous connaissons son expression fractionnaire. Mais cela n’éclaire en rien le phénomène d’« explosion » de la suite récurrente. Ce sera l’objet du prochain onglet, où nous utiliserons à nouveau les capacités de Python à calculer avec les entiers.
Pourquoi cette « explosion » ?
Revenons à notre suite $u_{n+1} = (q+1)\,u_n - p$. Pour $p=2, q=3$ par exemple, nous savons maintenant que Python va devoir tronquer l’écriture binaire de $\dfrac{p}{q}$, donc dès le départ introduire une erreur. Que se passe-t-il ensuite ? Continuons d’expérimenter, et toujours avec Python, mais avec le module fractions
déjà utilisé, qui ne calcule qu’avec des entiers, donc sans risque d’erreur.
Reprenons le calcul des valeurs de la suite :
- from fractions import *
- def suite_frac(p,q):
- u = Fraction(p,q)
- for i in range(100): # i varie de 0 à 99
- u = (q+1) * u - p
- if (i+1) % 10 == 0 : # l'indice mathématique est décalé de 1
- approx = float(u)
- print('Le terme de rang '+str(i+1)+' vaut : '+str(approx))
- print(' écart avec la valeur théorique : ',approx-p/q)
Exécutons ce programme pour $p=2, q=3$ , en gardant seulement le dernier affichage :
>>> suite_frac(2,3)
Le terme de rang 100 vaut : 0.6666666666666666
écart avec la valeur théorique : 0.0
Python a calculé avec les rationnels donc n’a pas fait d’erreur. Mais avec ce nouvel outil, si nous modifions un peu notre programme pour prendre pour $u_0$ un rationnel très proche mais différent de $\dfrac{p}{q}$, que se passe-t-il ?
- def suite_frac_pt_initial(p,q,p0,q0):
- '''comme suite_frac mais en partant d'un autre point initial p0/q0 '''
- u = Fraction(p0,q0)
- print('Valeur initiale : '+str(float(u)))
- for i in range(100): # i varie de 0 à 99
- u = (q+1) * u - p
- if (i+1) % 10 == 0 : # l'indice mathématique est décalé de 1
- approx = float(u)
- print('Le terme de rang '+str(i+1)+' vaut : '+str(approx))
- print(' écart avec la valeur théorique : ',approx-p/q)
Exécutons-le pour $p=2, q=3$ encore, et à partir d’une fraction très proche :
>>> suite_frac_pt_initial(2,3,2000000001,3000000001)
Valeur initiale : 0.6666666667777777
Le terme de rang 10 vaut : 0.6667831751110723
écart avec la valeur théorique : 0.00011650844440569408
Le terme de rang 20 vaut : 122.83462526772179
écart avec la valeur théorique : 122.16795860105512
Le terme de rang 30 vaut : 128102390.02472664
écart avec la valeur théorique : 128102389.35805997
Le terme de rang 40 vaut : 134325091023517.77
écart avec la valeur théorique : 134325091023517.1
Le terme de rang 50 vaut : 1.4085006664507546e+20
écart avec la valeur théorique : 1.4085006664507546e+20
Le terme de rang 60 vaut : 1.4769199948242665e+26
écart avec la valeur théorique : 1.4769199948242665e+26
Le terme de rang 70 vaut : 1.54866286049285e+32
écart avec la valeur théorique : 1.54866286049285e+32
Le terme de rang 80 vaut : 1.6238907076041507e+38
écart avec la valeur théorique : 1.6238907076041507e+38
Le terme de rang 90 vaut : 1.70277282261673e+44
écart avec la valeur théorique : 1.70277282261673e+44
Le terme de rang 100 vaut : 1.7854867152481602e+50
écart avec la valeur théorique : 1.7854867152481602e+50
Ainsi l’explosion se produit même en calcul exact, dès lors que la valeur initiale s’écarte un peu de celle attendue.
En fait, c’est un phénomène classique pour ces suites récurrentes du type $u_{n+1} = f(u_n)$ et qui n’a rien à voir avec Python. Voyons ce qui peut se passer pour une telle suite [6]. Nous supposerons la fonction $f$ dérivable. Si notre suite récurrente a une limite, disons $l$, en passant à la limite dans la relation qui définit $(u_n)$ on a par continuité de $f$ l’égalité $l=f(l)$. Ainsi $l$ doit être un point fixe de $f$.
Pour mieux voir ce qui peut se passer, nous allons étudier le cas de la fonction $f,\quad f(x) = \dfrac{x^2}{2}-3\dfrac{x}{2}+2$ : elle a deux points fixes (et deux seulement) en $x=1$ et $x=4$. Le programme qui suit va nous permettre de visualiser le phénomène.
- import matplotlib.pyplot as plt
- def f(x):
- '''la fonction qu'on veut representer'''
- return(x**2/2 - 3*x/2 + 2)
- def graphe(f,a,b,N):
- '''trace le graphe de la fonction f entre a et b avec N segments'''
- lx = [a+i*(b-a)/N for i in range(N+1)]
- ly = [f(x) for x in lx]
- return(plt.plot(lx,ly,'m',label ="$f$"))
- def suite_rec(f,u0,a,b,N):
- '''dessine sur [a,b] la suite des segments reliant les points qui
- representent la suite u_{n+1}=f(u_n) partant de u0 avec N termes'''
- u,v = u0,0
- # on crée les listes des abscisses et ordonnées des segments
- lx = [u] ; ly = [v]
- v = f(u) ; lx.append(u) ; ly.append(v) # premier segment vertical
- for i in range(N):
- u = v ; lx.append(u) ; ly.append(u) # segment horizontal
- v = f(u) ; lx.append(u) ; ly.append(v) # segment vertical
- plt.plot(lx,ly,'r',label = "$u_{n+1} =f(u_n)$")
- # on trace la droite d'équation y = x
- plt.plot([a,b],[a,b],'b',label = "$y=x$")
- #on trace le graphe de f sur le même intervalle
- graphe(f,a,b,30)
- plt.legend(loc='upper left') # ces legendes doivent avoir ete creees (option 'label' de plot)
- plt.title("$u_{0} = $ "+str(u0)+", $N = $"+str(N))
- # on ajoute les axes
- plt.axhline(color = 'k')
- plt.axvline(color='k')
- # et une grille
- plt.grid()
- # on affiche !
- plt.show()
Un premier test :
Un second :
Nous invitons le lecteur à télécharger le programme et à poursuivre les tests s’il le souhaite. Pour résumer, le point fixe 1 est « attractif » ($|f’(1)| < 1$), l’autre point fixe 4 est « répulsif » ($|f’(4)| > 1$) : pour une valeur initiale $u_0$ même proche de 4 (mais inférieure à 4), c’est vers 1 que la suite va converger. Pour $u_0$ plus grand que 4, elle diverge, et très rapidement :
C’est exactement ce qui se passe pour le phénomène dont nous étions partis, de divergence massive de la suite $(u_n)_n$ définie par la relation de récurrence $u_{n+1} = (q+1)\,u_n - p$ pour certaines valeurs du premier terme $u_0=\dfrac{p}{q}$ : ici la fonction $f$ est affine : $f(x)=(q+1)\,x - p$ donc sa dérivée est partout $q+1$ qui, pour les valeurs que nous avions prises, est plus grand que 1.
Pour nous résumer, à ce stade : ce problème semble en fait dû à la combinaison de deux phénomènes : l’un, informatique, est la conversion de la fraction en binaire ; l’autre, mathématique, est l’instabilité du « système dynamique » $u_{n+1} = f(u_n)$ au voisinage d’un point $a$ où $|f’(a)|\geq 1$ : Python prend une valeur approchée de $u_0$ et cela suffit pour faire diverger la suite.
Pour une dernière vérification, nous allons reprendre le même problème et notre premier programme, mais en modifiant $f$ de façon à trouver un système dynamique stable. Pour le cas $p=2, q=3$ par exemple, il suffit de prendre $p’=-\dfrac{1}{3}$ et $q’=-\dfrac{1}{2}$ de sorte que le premier terme $u_0=\dfrac{p’}{q’}$ vaut toujours $\dfrac{2}{3}$, mais cette fois la dérivée de $f$ est $q’+1=\dfrac{1}{2}$. On obtient :
>>> suite(-1/3,-1/2)
Le terme de rang 100 vaut : 0.6666666666666666
écart avec la valeur théorique : 0.0
C’est bien ça : Python fait toujours une approximation en convertissant en binaire, mais comme le point fixe $\dfrac{2}{3}$ est maintenant attractif, on ne s’écarte pas de cette valeur et la suite apparaît constante.
Et les calculatrices ?
La calculatrice NumWorks a un module Python intéressant et régulièrement amélioré. Et c’est un « vrai » Python, la preuve : à nos deux questions de départ, 0.1 + 0.1 + 0.1 - 0.3 et la suite récurrente, il a exactement les mêmes réponses (fausses) qu’un ordinateur ! De même lorsqu’on lui demande si le triangle de côtés $\sqrt{3}, \sqrt{5}$ et $2\sqrt{2}$ est rectangle.
Mais si l’on demande 0.1+0.1+0.1-0.3 au module « calculs » de cette même machine, la réponse est bien 0, et notre triangle de test est bien rectangle.
Sur une (assez ancienne) TI-Nspire CAS, en mode « calcul exact », toutes les réponses sont mathématiquement correctes, même pour notre suite récurrente : il n’y a pas d’approximation numérique, donc en partant de $\dfrac{p}{q}$ la suite reste constante. La calculatrice fait ici du calcul formel et non numérique. En passant en mode « approché », la suite récurrente diverge comme avec l’ordinateur, par contre les réponses restent correctes pour 0.1+0.1+0.1-0.3 et pour le triangle rectangle test.
En fait, les calculatrices utilisent le système « décimal codé binaire » (DCB) : pour un nombre entré en système décimal, chacun de ses chiffres est converti en binaire, et les calculs sont menés en traitant séparément chaque bloc de « bits » représentant un chiffre décimal. Ainsi par exemple 0.1 est représenté exactement, sans troncature ni arrondi. Et les calculs sont menés, de fait, en système décimal. Ce codage DCB est peu performant en termes d’occupation de la mémoire : il faut quatre bits pour coder chaque chiffre décimal, mais avec ces quatre bits on a $2^4= 16$ valeurs possibles, et là on n’en utilise [7] que 10. Mais pour les systèmes qui, comme les calculatrices, ne sont utilisés que sur des nombres et avec un format décimal pour l’entrée et la sortie, il a l’avantage de la simplicité.
La calculatrice reste cependant limitée par le nombre de chiffres avec lesquels elle travaille, 14 pour la TI. Donc pour notre suite récurrente, comme elle remplace (en mode « approché », toujours) $\dfrac{2}{3}$ par 0.66666666666667, cela suffit pour faire diverger le calcul, à cause de l’instabilité du système. On peut le vérifier avec notre petit programme de l’onglet précédent : pour suite_frac_pt_initial(2,3,p0,q0)
avec $p_0=6666666666667$ et $q_0=10^{13}$, le centième terme vaut $5.4 \times 10^{46}$ pour Python et $1.3 \times 10^{46}$ pour la TI : l’ordre de grandeur est le même mais Python a calculé avec plus de chiffres.
Pour aller un peu plus loin
Intéressons-nous à un cas plus compliqué, la suite récurrente d’ordre 2 et non linéaire définie par :
$u_0 = \dfrac{3}{2},\quad u_1 = \dfrac{5}{3}$
$u_{n+2} = 2003 - \dfrac{6002}{u_{n+1}} + \dfrac{4000}{u_n\,u_{n+1}}$
(source : Jean-Michel Muller, Vincent Lefèvre, « Erreurs en arithmétique des ordinateurs » — Images des Mathématiques, CNRS, 2009 [8])
On l’étudie expérimentalement avec le programme :
- x, y = 3/2, 5/3
- for i in range(20):
- z = y
- y = 2003 - 6002/y + 4000/(x*y)
- x = z
- print(y)
On obtient :
>>>
1.800000000000182
1.8888888890910494
1.941176684634911
1.9699174994632358
2.2085111769288233
204.74869262199354
1982.5318590274205
1999.9824122480152
1999.999982429244
1999.9999999824379
1999.9999999999825
2000.0
2000.0
2000.0
2000.0
2000.0
2000.0
2000.0
2000.0
2000.0
Il semble donc que la suite tende vers 2000. Mais s’agissant de calcul avec des « flottants », nous avons appris à être prudents.
Solution avec le calcul fractionnaire
Voyons d’abord ce que donne le calcul « formel » avec le module fractions
- from fractions import *
- def suite_frac(p,q,r,s,N):
- u = Fraction(p,q)
- v = Fraction(r,s)
- for i in range(N):
- w = v
- v = 2003 - 6002/v + 4000/(u*v)
- u = w
- if (i+1) % 10 == 0 :
- approx = float(v)
- print('Le terme de rang '+str(i+1)+' vaut : '+str(v))
- print(' valeur approchée : '+str(approx))
Testons :
>>> suite_frac(3,2,5,3,100)
Le terme de rang 10 vaut : 4097/2049
valeur approchée : 1.9995119570522206
Le terme de rang 20 vaut : 4194305/2097153
valeur approchée : 1.9999995231630692
Le terme de rang 30 vaut : 4294967297/2147483649
valeur approchée : 1.9999999995343387
Le terme de rang 40 vaut : 4398046511105/2199023255553
valeur approchée : 1.9999999999995453
Le terme de rang 50 vaut : 4503599627370497/2251799813685249
valeur approchée : 1.9999999999999996
Le terme de rang 60 vaut : 4611686018427387905/2305843009213693953
valeur approchée : 2.0
Le terme de rang 70 vaut : 4722366482869645213697/2361183241434822606849
valeur approchée : 2.0
Le terme de rang 80 vaut : 4835703278458516698824705/2417851639229258349412353
valeur approchée : 2.0
Le terme de rang 90 vaut : 4951760157141521099596496897/2475880078570760549798248449
valeur approchée : 2.0
Le terme de rang 100 vaut : 5070602400912917605986812821505/2535301200456458802993406410753
valeur approchée : 2.0
On constate déjà que le résultat est différent de ce qu’il était quand Python calculait avec les flottants.
Essayons avec d’autres valeurs initiales :
>>> suite_frac(1,2,5,4,20)
Le terme de rang 10 vaut : 9229840145297874162306378414416051/4614920072648937081153189199027
valeur approchée : 2000.0
Le terme de rang 20 vaut : 9451356308785023142201731496378820040650956108684973117189208447155/4725678154392511571100865748189410020325478054342486558586214579
valeur approchée : 2000.0
>>> suite_frac(1,1,1,1,20)
Le terme de rang 10 vaut : 1
valeur approchée : 1.0
Le terme de rang 20 vaut : 1
valeur approchée : 1.0
Partant de valeurs de $u_0$ et $u_1$ proches de 1, on obtient donc trois résultats différents. Où est le problème ?
Brève étude mathématique
La limite de $(u_n)$, si elle existe, est un point fixe de $f$,
$f(x) = 2003 - \dfrac{6002}{x} + \dfrac{4000}{x^2}$, autrement dit une solution de
$x^3 - 2003\, x^2 + 6002\, x - 4000 = 0$
Les racines plus ou moins évidentes de ce polynôme sont 1, 2 et 2000. Vérifions avec sympy
:
- >>> import sympy
- >>> import sympy.solvers as sol
- >>> x = sympy.symbols('x')
- >>> sol.solve(x**3-2003*x**2+6002*x-4000,x)
- [1, 2, 2000]
On commence à voir apparaître un phénomène analogue à celui décrit dans l’article : l’instabilité numérique amène la convergence vers un autre point fixe que celui qu’on aurait dû obtenir. Étudions de plus près ces points fixes.
La fonction $(u_{n-1},u_n)\longmapsto (u_n,u_{n+1})$ est ici $f : (x,y)\longmapsto (y,2003 - \dfrac{6002}{y} + \dfrac{4000}{x\,y})$.
L’analogue de la dérivée pour les fonctions d’une seule variable est ici la matrice jacobienne $J(x,y)$ : matrice $2\times 2$ des dérivées partielles de $f$. Dériver $y$ par rapport à $x$ ou à $y$, c’est facile, mais pour la deuxième composante de $f$, on peut la demander à sympy
:
- >>> x, y, t = sympy.symbols('x, y, t')
- >>> from sympy import Derivative
- >>> f2 = 2003 - 6002/y + 4000/(x*y)
- >>> d1 = Derivative(f2,x).doit()
- >>> d2 = Derivative(f2,y).doit() ; d1, d2
- (-4000/(x**2*y), 6002/y**2 - 4000/(x*y**2))
Ce qui va nous intéresser, ce sont les valeurs propres de la jacobienne, autrement dit les racines de son polynôme caractéristique, en $(x,y)$ :
>>> carac = t**2 - d2 * t -d1
Demandons maintenant à sympy
les racines de ce polynôme en $t$ :
>>> vp = sol.solve(carac,t) ; vp
[(3001*x - sqrt(9006001*x**2 - 12004000*x - 4000*y**3 + 4000000) - 2000)/(x*y**2),
(3001*x + sqrt(9006001*x**2 - 12004000*x - 4000*y**3 + 4000000) - 2000)/(x*y**2)]
On n’est guère avancé. Mais ces valeurs propres nous intéressent en fait seulement aux points $(x,y) = (u,u)$ où $u$ est l’un des points fixes 1, 2 et 2000. Il suffit donc de demander à sympy
de remplacer $x$ et $y$ par chacune de ces valeurs, et de donner les racines du polynôme en $t$ qui en résulte :
- >>> for u in [1,2,2000]:
- print('Les valeurs propres de la Jacobienne en ('+str(u)+','+str(u)+') sont : ')
- print(sol.solve(carac.subs({x:u,y:u})))
- Les valeurs propres de la Jacobienne en (1,1) sont :
- [2, 2000]
- Les valeurs propres de la Jacobienne en (2,2) sont :
- [1/2, 1000]
- Les valeurs propres de la Jacobienne en (2000,2000) sont :
- [1/2000, 1/1000]
Ainsi en $(1,1)$ les deux valeurs propres sont plus grandes que 1 : ce point est répulsif.
En $(2000,2000)$, elles sont toutes deux inférieures (en valeur absolue) à 1. Donc 2000 est très attractif.
Le cas de $(2,2)$ est, avec ce critère, indéterminé : une seule valeur propre au-delà de 1.
Étude expérimentale de la convergence
On va voir de plus près, en fonction des valeurs initiales $u_0 = \dfrac{p}{q}$ et $u_1 = \dfrac{r}{s}$ vers quelle valeur converge vraiment (avec le calcul exact sur les rationnels, donc) la suite. Modifions un peu notre programme, pour ne garder que le résultat au bout d’un nombre déterminé d’itérations :
- def sfrac(p,q,r,s,Niter):
- '''calcule avec le module Fractions les Niter premiers termes
- de la suite et retourne le dernier, en valeur approchée'''
- u = Fraction(p,q)
- v = Fraction(r,s)
- for i in range(Niter):
- w = v
- v = 2003 - 6002/v + 4000/(u*v)
- u = w
- return(float(v))
Étude graphique
Pour fixer les idées, on va considérer l’intervalle $I = [\frac{1}{2},\frac{9}{2}]$ et, pour des valeurs de $(u_0,u_1)$ dans $I^2$ régulièrement espacées, calculer avec sfrac
la limite de $(u_n)$. On pourra alors représenter ce maillage du carré, avec des couleurs différentes selon la limite trouvée pour $(u_n)$.
- def cvgce(Niter,Npts,eps):
- '''divise l'intervalle en Npts régulièrement espacés, et en chaque point du maillage
- de I^2 ainsi obtenu, calcule lim de u_n avec Niter itérations ;
- le point du maillage est rangé dans une liste qui est fonction de la limite obtenue,
- identifiée à eps près à 1, 2, 2000 ou autre'''
- lx1, ly1, lx2, ly2, lx2000, ly2000, lxa, lya = [], [], [], [], [], [], [], []
- for i in range(Npts) :
- pi, qi = Npts - 1 + 8*i, 2*(Npts - 1)
- for j in range(Npts) :
- rj, sj = Npts - 1 + 8*j, 2*(Npts - 1)
- zij = sfrac(pi,qi,rj,sj,Niter)
- if abs(zij - 1) < eps :
- lx1.append(pi/qi)
- ly1.append(rj/sj)
- elif abs(zij - 2) < eps :
- lx2.append(pi/qi)
- ly2.append(rj/sj)
- elif abs(zij - 2000) < eps :
- lx2000.append(pi/qi)
- ly2000.append(rj/sj)
- else :
- lxa.append(pi/qi)
- lya.append(rj/sj)
- return(lx1, ly1, lx2, ly2, lx2000, ly2000, lxa, lya)
Et pour la représentation :
- import matplotlib.pyplot as plt
- def figure_cvgce(Niter,Npts,eps):
- '''dessine le résultat de cvgce, en mettant en évidence les points
- où la limite est autre que 2000'''
- lx1, ly1, lx2, ly2, lx2000, ly2000, lxa, lya = cvgce(Niter,Npts,eps)
- plt.plot(lx1,ly1,'o',color = 'r',markersize = 8)
- plt.plot(lx2,ly2,'o',color = 'y',markersize = 8)
- plt.plot(lx2000,ly2000,'o',markersize = 1)
- plt.plot(lxa,lya,'o',color = 'k',markersize = 8)
- plt.show()
La commande figure_cvgce(20,121,0.01)
donne l’image
Là on voit le point fixe $(1,1)$, et quelques points initiaux pour lesquels la suite tend vers 2, parmi un océan de points pour lesquels elle tend vers 2000. Mais c’est parce que $Npts$ est bien choisi : avec
figure_cvgce(20,120,0.01)
on obtient
où tous les points sont « attirés » par 2000...
Étude plus fine avec les fractions
Déjà, pour $u_0 = u_1 = a \pm 10^{-k}, \quad a = 1$ ou 2, pour quelles valeurs de $k$ a-t-on convergence vers $a$ ?
- def seuil(a,eps,k):
- '''attention ça n'a d'intérêt que si a = 1 ou 2 et eps = 1 ou -1'''
- u = Fraction(a+eps*10**(-k))
- p, q = u.numerator, u.denominator
- return sfrac(p,q,p,q,20)
- def test_seuil(a,eps,tol=0.1):
- k = 1
- while abs(seuil(a,eps,k) - a) > tol :
- k += 1
- return k
Testons quatre cas ensemble :
>>> test_seuil(1,1), test_seuil(2,1), test_seuil(1,-1), test_seuil(2,-1)
(16, 16, 17, 16)
Il faut donc partir très près d’un des points fixes 1 et 2, pour espérer y rester. Et si l’on essaye avec $u_0 = a, \quad u_1 = a + 10^{-k}$ :
- def seuil_var(a,eps,k):
- '''attention ça n'a d'intérêt que si a = 1 ou 2 et eps = 1 ou -1'''
- u = Fraction(a+eps*10**(-k))
- p, q = u.numerator, u.denominator
- return sfrac(a,a,p,q,20)
c’est bien pire ! Même avec k = 500 la réponse est encore 2000.
Si enfin on regarde de plus près ce qui se passe pour les premiers termes, avec cette variante
- def suite_frac_detail(p,q,r,s,N):
- u = Fraction(p,q)
- v = Fraction(r,s)
- for i in range(N):
- w = v
- v = 2003 - 6002/v + 4000/(u*v)
- u = w
- approx = float(v)
- print('Le terme de rang '+str(i+1)+' vaut : '+str(v))
- print(' valeur approchée : '+str(approx))
avec nos valeurs initiales, on obtient :
>>> suite_frac_detail(3,2,5,3,10)
Le terme de rang 1 vaut : 9/5
valeur approchée : 1.8
Le terme de rang 2 vaut : 17/9
valeur approchée : 1.8888888888888888
Le terme de rang 3 vaut : 33/17
valeur approchée : 1.9411764705882353
Le terme de rang 4 vaut : 65/33
valeur approchée : 1.9696969696969697
Le terme de rang 5 vaut : 129/65
valeur approchée : 1.9846153846153847
Le terme de rang 6 vaut : 257/129
valeur approchée : 1.9922480620155039
Le terme de rang 7 vaut : 513/257
valeur approchée : 1.9961089494163424
Le terme de rang 8 vaut : 1025/513
valeur approchée : 1.9980506822612085
Le terme de rang 9 vaut : 2049/1025
valeur approchée : 1.9990243902439024
Le terme de rang 10 vaut : 4097/2049
valeur approchée : 1.9995119570522206
>>>suite_frac_detail(5,3,3,2,10)
Le terme de rang 1 vaut : -1195/3
valeur approchée : -398.3333333333333
Le terme de rang 2 vaut : 2403591/1195
valeur approchée : 2011.3732217573222
Le terme de rang 3 vaut : 4807208383/2403591
valeur approchée : 2000.010976493089
Le terme de rang 4 vaut : 9614416817967/4807208383
valeur approchée : 2000.000010810224
Le terme de rang 5 vaut : 19228833636037135/9614416817967
valeur approchée : 2000.000000010727
Le terme de rang 6 vaut : 38457667272074475471/19228833636037135
valeur approchée : 2000.0000000000107
Le terme de rang 7 vaut : 76915334544148951352143/38457667272074475471
valeur approchée : 2000.0
Le terme de rang 8 vaut : 153830669088297902705105487/76915334544148951352143
valeur approchée : 2000.0
Le terme de rang 9 vaut : 307661338176595805410212612175/153830669088297902705105487
valeur approchée : 2000.0
Le terme de rang 10 vaut : 615322676353191610820425227625551/307661338176595805410212612175
valeur approchée : 2000.0
Donc permuter $u_0$ et $u_1$ suffit à ce que 2000 l’emporte à nouveau.
En résumé, le couple de valeurs initiales $(\frac{3}{2},\frac{5}{3})$ est un cas très particulier : des valeurs relativement éloignées de 2 qui donnent une limite de 2. Sinon, 2000 est tellement attractif qu’il faut prendre $u_0=u_1$ et très proches de 1 ou 2 (de l’ordre de $10^{-16}$) pour que la limite soit 1 ou 2. Mais prédire le problème, ou plutôt l’éventualité d’un problème, suppose une étude mathématique déjà assez poussée.
Il est important de souligner ici que le problème présenté n’est en rien particulier au langage Python : comme on l’a montré, la cause est la traduction en langage binaire, commune à tous les langages utilisés pour le calcul numérique. L’utilisateur doit être conscient de cette particularité, d’abord pour respecter des précautions élémentaires : déjà, comme on l’a dit, ne jamais tester l’égalité de deux flottants. Mais aussi parce que ses conséquences peuvent être bien plus dramatiques que la divergence d’une suite qui devrait converger, voir par exemple ici. Les mathématiques, plus précisément l’analyse numérique, peuvent souvent dire dans quelles conditions un système risque ou non de diverger. Mais cela ne dispense pas, pour tout calcul confié à une machine, de contrôler à l’arrivée la vraisemblance du résultat...