1Calcul approche et arithmetique flottante
De nombreux problemes scientifiques n'ont pas de solution exacte calculable : resoudre \(\cos x = x\), calculer \(\int_0^1 e^{-x^2}\,\mathrm dx\), integrer l'equation du pendule. L'ingenierie numerique consiste a en calculer des solutions approchees par des algorithmes, avec une precision maitrisee.
| Probleme | Methodes | Chapitre |
|---|---|---|
| Resoudre \(f(x) = 0\) | dichotomie, Newton, secante | 21 |
| Calculer \(\int_a^b f(t)\,\mathrm dt\) | rectangles, trapezes, Simpson | 22 |
| Resoudre \(y' = f(t,y)\) | Euler, Runge-Kutta | 23 |
Toutes ces methodes suivent la meme demarche en trois temps :
- Discretiser : remplacer un objet continu (une fonction, un intervalle) par un nombre fini de valeurs, indexees par un pas \(h\) ou un nombre d'iterations \(n\).
- Iterer un calcul elementaire.
- Controler l'erreur : savoir comment elle depend de \(h\), et jusqu'ou il est utile de diminuer \(h\).
Les nombres flottants
Un flottant sur 64 bits s'ecrit \[ x = (-1)^s \times 1{,}m \times 2^{e-1023}, \] avec 1 bit de signe \(s\), 11 bits d'exposant \(e\) et 52 bits de mantisse \(m\). Cela donne environ 16 chiffres decimaux significatifs et un domaine allant de \(10^{-308}\) a \(10^{308}\).
import sys
sys.float_info.epsilon # 2.220446049250313e-16 : epsilon machine
sys.float_info.max # 1.7976931348623157e+308
sys.float_info.min # 2.2250738585072014e-308
1 + 1e-16 == 1 # True : 1e-16 est absorbe
1 + 1e-15 == 1 # FalseL'epsilon machine \(\varepsilon\) est le plus petit reel positif tel que \(1 + \varepsilon \neq 1\) en arithmetique flottante. En double precision, \(\varepsilon = 2^{-52} \approx 2{,}22\cdot10^{-16}\).
Consequence fondamentale : pour tout reel \(x\) representable, le flottant \(\mathrm{fl}(x)\) le plus proche verifie \[ \left|\frac{\mathrm{fl}(x) - x}{x}\right| \le \frac{\varepsilon}{2}. \] L'erreur relative est bornee, pas l'erreur absolue : les grands nombres sont approches grossierement en valeur absolue.
Les trois pieges
0.1 + 0.2 == 0.3 # False
0.1 + 0.2 # 0.30000000000000004En base 2, \(0{,}1 = 1/10\) n'a pas d'ecriture finie (10 n'est pas une puissance de 2), comme \(1/3\) en base 10. On ne teste jamais l'egalite :
from math import isclose
isclose(0.1 + 0.2, 0.3) # True
abs(a - b) <= 1e-12 * max(abs(a), abs(b)) # tolerance relative, a la main1e16 + 1 - 1e16 # 0.0 : le 1 a ete absorbe
s = 0.0
for _ in range(10**6):
s += 0.1
print(s) # 100000.00000133288, et non 100000.0Quand on ajoute \(y\) a \(x\) avec \(|y| < \varepsilon|x|\), le resultat vaut exactement \(x\) : \(y\) est perdu. C'est pourquoi il vaut mieux sommer une serie du plus petit terme au plus grand (chapitre 1).
from math import sqrt
x = 1e8
sqrt(x*x + 1) - x # 0.0 ! la vraie valeur est ~ 5e-09Soustraire deux nombres presque egaux annule les chiffres significatifs communs : il ne reste que du bruit d'arrondi. La parade consiste a reecrire l'expression :
1 / (sqrt(x*x + 1) + x) # 4.999999999999999e-09 : correct
# on a utilise l'identite sqrt(A) - B = (A - B^2)/(sqrt(A) + B)C'est le meme mecanisme que pour les racines d'un trinome (chapitre 1) et pour la differentiation numerique (chapitre 18).
2Erreur de methode, erreur d'arrondi, conditionnement
Erreur de methode et erreur d'arrondi
L'erreur totale d'un calcul numerique se decompose en deux contributions antagonistes :
\[ E_{\text{totale}}(h) \;\approx\; \underbrace{C h^p}_{\text{erreur de methode}} \;+\; \underbrace{\frac{\varepsilon}{h^q}}_{\text{erreur d'arrondi}} \]
- l'erreur de methode (ou de troncature) vient de la discretisation ; elle diminue quand \(h \to 0\) ;
- l'erreur d'arrondi vient des flottants ; elle augmente quand \(h\) diminue, car le nombre d'operations croit et/ou l'on divise par de petites quantites.
Une methode est d'ordre \(p\) si son erreur de methode est en \(O(h^p)\). Diviser \(h\) par 10 divise alors l'erreur par \(10^p\). L'ordre se lit experimentalement comme la pente de l'erreur en fonction de \(h\), en echelle log-log : \[ \ln E \approx \ln C + p\,\ln h. \]
import numpy as np
from math import exp
def ordre_experimental(erreur_de_h, h_values):
"""Estime l'ordre p par regression lineaire sur les logarithmes."""
h = np.array(h_values, dtype=float)
E = np.array([erreur_de_h(x) for x in h])
valides = E > 0
return np.polyfit(np.log(h[valides]), np.log(E[valides]), 1)[0]
f, x0, exact = exp, 1.0, exp(1.0)
decentree = lambda h: abs((f(x0+h) - f(x0)) / h - exact)
centree = lambda h: abs((f(x0+h) - f(x0-h)) / (2*h) - exact)
hs = [10**-k for k in range(2, 7)]
print(f"decentree : ordre {ordre_experimental(decentree, hs):.3f}") # 1.000
print(f"centree : ordre {ordre_experimental(centree, hs):.3f}") # 2.000Minimiser \(E(h) = Ch^p + \varepsilon/h\) donne \[ h^* = \left(\frac{\varepsilon}{pC}\right)^{\frac{1}{p+1}}, \qquad E(h^*) = \Theta\!\left(\varepsilon^{\frac{p}{p+1}}\right). \] Pour \(p = 1\) : \(h^* \approx \sqrt\varepsilon \approx 10^{-8}\) et une precision maximale de \(10^{-8}\). Pour \(p = 2\) : \(h^* \approx \varepsilon^{1/3} \approx 6\cdot10^{-6}\) et une precision de \(10^{-11}\).
Conclusion contre-intuitive : reduire indefiniment le pas degrade le resultat. Il existe une precision maximale atteignable, qui depend de l'ordre de la methode.
Conditionnement d'un probleme
Le conditionnement mesure la sensibilite de la solution aux perturbations des donnees, independamment de l'algorithme. Pour \(y = f(x)\), le conditionnement relatif est \[ \kappa(x) = \left|\frac{x f'(x)}{f(x)}\right|, \] et une perturbation relative \(\delta\) sur \(x\) produit une perturbation relative \(\approx \kappa\delta\) sur \(y\).
import numpy as np
# systeme bien conditionne
A = np.array([[2., 1.], [1., 3.]])
print(f"cond = {np.linalg.cond(A):.2f}") # 2.62 : tres bien
# systeme mal conditionne
B = np.array([[1., 1.], [1., 1.0000001]])
print(f"cond = {np.linalg.cond(B):.2e}") # 4.0e+07
b1 = np.array([2., 2.0000001])
b2 = np.array([2., 2.0000002]) # perturbation de 5e-8
print(np.linalg.solve(B, b1)) # [1. 1.]
print(np.linalg.solve(B, b2)) # [0. 2.] : reponse totalement differente !Une perturbation relative de \(5\cdot10^{-8}\) sur le second membre change completement la solution. Aucun algorithme ne peut y remedier : le probleme lui-meme est instable. Il faut le reformuler.
| Conditionnement | Stabilite | |
|---|---|---|
| Porte sur | le probleme | l'algorithme |
| Question | la solution depend-elle fortement des donnees ? | l'algorithme amplifie-t-il les erreurs d'arrondi ? |
| Remede | reformuler le probleme | changer d'algorithme |
| Exemple | matrice de Hilbert | Gauss sans pivot partiel |
Un algorithme stable applique a un probleme bien conditionne donne un resultat precis. Toute autre combinaison est risquee.
3Vitesses de convergence et validation d'un code
Convergence : vocabulaire
Soit \((x_n)\) convergeant vers \(\ell\), et \(e_n = |x_n - \ell|\).
| Type | Relation | Gain par iteration | Exemple |
|---|---|---|---|
| Lineaire | \(e_{n+1} \approx K e_n\), \(0<K<1\) | facteur constant | dichotomie (\(K = 1/2\)), point fixe |
| Quadratique | \(e_{n+1} \approx K e_n^2\) | le nombre de decimales double | Newton |
| Ordre \(\varphi\) | \(e_{n+1} \approx K e_n^{1{,}618}\) | intermediaire | secante |
from math import sqrt
exact = sqrt(2)
x = 1.0
print(f"{'n':>3}{'x_n':>22}{'erreur':>14}{'decimales':>12}")
for n in range(6):
e = abs(x - exact)
d = -1 if e == 0 else int(-__import__("math").log10(e))
print(f"{n:>3}{x:>22.16f}{e:>14.2e}{d:>12}")
x = 0.5 * (x + 2/x) # iteration de Newton (methode de Heron) n x_n erreur decimales
0 1.0000000000000000 4.14e-01 0
1 1.5000000000000000 8.58e-02 1
2 1.4166666666666665 2.45e-03 2
3 1.4142156862745097 2.12e-06 5
4 1.4142135623746899 1.59e-12 11
5 1.4142135623730951 0.00e+00 16Le nombre de decimales exactes suit 0, 1, 2, 5, 11, 16 : il double a chaque iteration. Cinq iterations suffisent a atteindre la precision machine. Une dichotomie demanderait 52 iterations pour le meme resultat.
Methodologie de validation
- Cas exactement soluble : tester d'abord sur un cas dont on connait la solution analytique. Si l'erreur n'est pas a la hauteur attendue, le code est faux.
- Ordre de convergence : diviser \(h\) par 2 et verifier que l'erreur est divisee par \(2^p\). Un ordre mesure different de l'ordre theorique revele presque toujours un bug.
- Invariants physiques : energie conservee, masse conservee, symetries. Un schema qui viole une loi de conservation est suspect.
- Deux methodes independantes : comparer les resultats de deux algorithmes differents.
- Sensibilite : perturber legerement les donnees et verifier que la solution ne change pas de facon disproportionnee.
import numpy as np
def teste_ordre(methode, exacte, hs):
"""Verifie l'ordre en divisant le pas par 2 et en observant le rapport des erreurs."""
print(f"{'h':>12}{'erreur':>14}{'rapport':>10}{'ordre':>9}")
precedent = None
for h in hs:
e = abs(methode(h) - exacte)
if precedent is None:
print(f"{h:>12.2e}{e:>14.3e}{'':>10}{'':>9}")
else:
r = precedent / e
print(f"{h:>12.2e}{e:>14.3e}{r:>10.2f}{np.log2(r):>9.2f}")
precedent = e
# integration par la methode des trapezes de sin sur [0, pi] : valeur exacte 2
from math import sin, pi
def trapezes(h):
n = max(1, int(round(pi / h)))
x = np.linspace(0, pi, n + 1)
y = np.sin(x)
return float((pi / n) * (y[0]/2 + y[1:-1].sum() + y[-1]/2))
teste_ordre(trapezes, 2.0, [pi/2**k for k in range(2, 9)]) h erreur rapport ordre
7.85e-01 1.045e-01
3.93e-01 2.579e-02 4.05 2.02
1.96e-01 6.433e-03 4.01 2.00
9.82e-02 1.608e-03 4.00 2.00
4.91e-02 4.019e-04 4.00 2.00
2.45e-02 1.005e-04 4.00 2.00
1.23e-02 2.512e-05 4.00 2.00Le rapport est constamment egal a 4, donc l'ordre vaut \(\log_2 4 = 2\) : la methode des trapezes est bien d'ordre 2. C'est le test le plus discriminant qui soit pour valider un code numerique.
L'essentiel du chapitre
- Un flottant double precision offre \(\approx16\) chiffres significatifs ; \(\varepsilon \approx 2{,}2\cdot10^{-16}\).
- Trois pieges : egalite entre flottants, absorption, cancellation catastrophique.
- Erreur totale = erreur de methode (\(\propto h^p\), decroissante) + erreur d'arrondi (\(\propto \varepsilon/h\), croissante).
- Il existe un pas optimal : diminuer \(h\) au-dela degrade le resultat.
- L'ordre \(p\) se mesure comme la pente de \(\ln E\) en fonction de \(\ln h\).
- Le conditionnement porte sur le probleme, la stabilite sur l'algorithme.
- Convergence lineaire (facteur constant) contre quadratique (les decimales doublent).