Lorsque l'on met en place une méthode numérique pour résoudre un problème scientifique il importe d'avoir une approche critique des résultats. En ce qui concerne la précision, il est utile de pouvoir apprécier l’erreur produite par la méthode numérique utilisée. Cette erreur est essentiellement de deux types : l’erreur d'arrondi et l'erreur d'approximation.
Dans une machine, un nombre réel $x$ est codé et stocké dans un format dit en « virgule flottante ». Il est représenté dans une base $b$ ($b = 2$ pour un ordinateur), par son signe ($+$ ou $-$), une mantisse $m$ et un exposant $e$ : $$\begin{array}{l} x = \pm m b^e\\ m = d. dd\dots d\\ e = dd\dots d\end{array}$$ où $d\in\{0, 1,\dots , b - 1\}$ représente un chiffre. Pour rendre cette représentation unique, on utilise une mantisse normalisée : le premier chiffre avant le point décimal de la mantisse est non nul (problème pour coder le $0$ !).
Les exemples qui suivent ne prennent pas en compte l’offset introduit pour le codage signé de l’exposant $e$ suivant la norme IEE754. Pour une description plus précise, vous devez vous référer au cours.
Considérons une mantisse à $n$ chiffres, deux situations sont à envisagées :
Ainsi, si $n = 4$, alors $m = 3.456$. Il y a une erreur d'arrondi.
Le nombre $\dfrac{1}{5}= 0.2$ a un développement infini en puissance de $2$, son codage en base $2$ sera aussi tronqué.
Ces arrondis et ces troncatures des nombres codés ont pour origine la mantisse $m$ qui est toujours codée avec une précision finie.
C’est particulièrement vrai pour les nombres irrationnels tels que $\pi$, $e$, $\sqrt{2}$, $\dots$ et plus généralement pour les nombres non dyadiques (qui ne sont pas de la forme $k/2^n$). Il en est de même pour l'exposant $e$ qui a une dynamique de codage nécessairement limitée, si bien que les nombres trop grands, ou trop proches de zéro, ne peuvent être représentés.
En Python, le codage des nombres est conforme à la norme IEEE754 en double précision c’est-à-dire sur $64$ bits, soit $8$ octets.
Objectif de la séance : plusieurs types et différentes causes caractérisent les erreurs numériques. Nous nous intéresserons d'abord à la précision de la machine, puis plus généralement aux erreurs générées par le codage des nombres en précision finie ainsi que la propagation de ces erreurs (smearing).
En début de votre script, importer la bibliothèque math. En mode préfixé, dir(math) retourne la liste des fonctions de la bibliothèque. De la même façon, la bibliothèque des fonctions de plot : import matplotlib.pyplot as plt. En fin de document, vous avez un rappel des fonctions nécessaires pour une représentation correcte de vos résultats.
Q1. Arrondi, troncature, ... bref, que vous retourne Python suite au traitement de l’opération $0.2+0.1-0.1-0.2$ ?
# Cellule à compléter
Python, mais aussi Matlab, $\dots$ codent les nombres réels (type float
) avec une précision finie. Par conséquent, il existe un entier $k$ tel
que:
$a = 1 + 2^{-k}$ est codé sans erreur et pour tout $q > k$, le codage du réel $a$ est tel que $a = 1 + 2^{-q} = 1$.
L’égalité signifiant que le codage de $a$ sera confondu avec $1$ car $2^{-q}$ est trop petit pour pouvoir être codé.
On appelle précision machine le nombre $\varepsilon = 2^{-k}$. Sa connaissance est d’une grande importante, ainsi lorsque qu’un résultat $X$ est effectué en $N$ multiplications, on montre que l’erreur relative $\dfrac{\delta X}{X}$ varie comme $\sqrt{N}\varepsilon$.
L'objectif de cet exercice est de calculer la précision machine $\varepsilon$.
Q2. Proposer un algorithme qui retourne la valeur de la précision machine.
Ainsi par exemple, on part de $\eta = 2^0 = 1$, et tant que
$a=\eta + 1$ est différent de $1$, on divise $\eta$ par $2$.
En déduire le nombre de décimales exactes.
# Cellule à compléter
Quelques règles pour limiter les erreurs d'arrondi :
- utiliser un système d'unités tel que les valeurs traitées sont de l'ordre de l'unité;
- éviter de soustraire des nombres comparables, minimiser les soustractions;
- chercher à sommer par paquet les quantités de même ordre de grandeur.
Considérons la suite récurrente d’ordre $2$ $u_n$ définie par $u_0$, $u_1$ donnés dans $\mathbb{R}$ et : $$\forall n\geq 2,\quad u_n = 2003 -\dfrac{6002}{u_{n-1}}+\dfrac{4000}{u_{n-1}u_{n-2}}$$
On montre par récurrence que le terme général de la suite $u_n$ peut être mis sous la forme suivante : $$u_n=\dfrac{a + 2^{n+1} b + 2000^{n+1}c}{a + 2^n b + 2000^n c}$$
Les coefficients $a$, $b$, $c$ se déterminent à partir des termes $u_0$ et $u_1$ de la suite $(u_n)$.
On peut vérifier que si $c\neq 0$, la suite tend vers $2000$ et que si $c = 0$ et $b\neq 0$ elle tend vers $2$.
Ainsi pour $a = b = 1$ et $c = 0$ on a $u_0 =\dfrac{3}{2}$, $u_1 =\dfrac{5}{3}$ et $\displaystyle\lim_{n\to+\infty} u_n = 2$. L'évolution de la convergence de la suite est représentée sur la figure suivante.
Q3. Proposer le code d’une fonction de nom cvu(N)
d’argument N
qui retourne un tableau contenant les N
premiers termes de la suite $(u_n)$ de premiers termes $u_0 =\dfrac{3}{2}$ et $u_1 =\dfrac{5}{3}$. On calculera ces termes avec la formule de récurrence.
Vous justifierez ce résultat en vous aidant de la formule explicite.
# Cellule à compléter
def cvu(N) :
# Cellule pour tester
Cet exemple est une question du DM1 pour lequel on vous propose de faire le code afin d’étudier la convergence d'une suite associée à l’estimation de $\pi$. L’énoncé de la question est en partie une reprise du DM1.
Le nombre $\pi$ est connu depuis l’antiquité, en tant que méthode de calcul du périmètre du cercle ou de l’aire du disque. Le problème de la quadrature du cercle étudié par les anciens Grecs consiste à construire un carré de même aire qu’un cercle donné à l’aide d’une règle et d’un compas. Ce problème resta insoluble jusqu’au 19ème siècle, où la démonstration de la transcendance de $\pi$ montra que le problème ne peut être résolu en utilisant une règle et un compas.
L’aire d’un cercle de rayon $r$ est aujourd’hui connue, ... Il fut un temps où il n’en était rien. Parmi les solutions alors proposées pour l’approcher, une méthode consiste à construire un polygone dont le nombre de côté augmenterait jusqu’à ce qu’il devienne équivalent au cercle circonscrit. C’est Archimède, vers 250 avant J-C, qui appliquera cette propriété au calcul des décimales du nombre $\pi$, en utilisant à la fois un polygone inscrit et circonscrit au cercle. Il utilise ainsi un algorithme pour le calcul et parvient à l’approximation de $\pi$ dans l’intervalle $\left[3 +\dfrac{10}{71}, 3 +\dfrac{1}{7} \right]$ en faisant tendre le nombre de côtés jusqu’à $96$.
Regardons l’algorithme de calcul par les polygones inscrits. On considère un cercle de rayon $r = 1$ et on note $A_n$ l’aire associée au polygône inscrit à $n$ côtés. En notant $\alpha_n =\dfrac{2\pi}{n}$, l'aire $A_n$ est égale à $n$ fois l'aire du triangle $ABC$ représentée sur la figure ci-dessous.
Nous avons : $$A_n=n\cos\left(\dfrac{\alpha_n}{2}\right)\sin\left(\dfrac{\alpha_n}{2}\right) =\dfrac{n}{2}\left(2\cos\left(\dfrac{\alpha_n}{2}\right)\sin\left(\dfrac{\alpha_n}{2}\right)\right)=\dfrac{n}{2}\sin(\alpha_n)=\dfrac{n}{2}\sin\left(\dfrac{2\pi}{n}\right)$$
L’estimation de $\pi$ repose sur la recherche d’une relation de récurrence sur $A_n$. Elle est construite simplement à partir de relations trigonométriques élémentaires, soit: $$\sin\left(\dfrac{\alpha_n}{2}\right)=\sqrt{\dfrac{1-\cos(\alpha_n)}{2}}==\sqrt{\dfrac{1-\sqrt{1-\sin^2(\alpha_n)}}{2}}\qquad\qquad(*)$$
Choisissons $n$ telle que $n = 2^k$, c’est une relation dyadique qui existe entre deux itérées successifs; l’angle $\alpha_n$ est divisé par deux alors que le nombre de triangles isocèles est multiplié par deux.
Comme $A_n=\dfrac{n}{2}\sin\left(\dfrac{2\pi}{n}\right)$ on a $A_{2^k}=2^{k-1}s_k$ en posant $s_k=\sin\left(\dfrac{\pi}{2^k}\right)$. De plus la relation $(*)$ nous donne la relation de récurrence pour calculer les $s_k$: $$s_k=\sqrt{\dfrac{1-\sqrt{1-s_{k-1}^2}}{2}}$$
On a $\displaystyle\lim_{k\to+\infty} A_k = \pi$, ce n’est toutefois pas du tout ce que l’on observe sur machine à cause du codage sur un nombre fini de bits, comme on va le constater.
La récurrence précédente est initialisée par $n=2$ donc $k=2$ et $s_2=\sin(\pi/2) =1$. On calcule les $s_k$ avec la formule de récurrence puis $A_{2^k}$.
Q4. Décliner le code python de la fonction EstDePi(k)
d’argument k
qui
retourne l'aire $A_{2^k}$.
Modifier ensuite le code pour que la fonction renvoie les $N$ premiers de la suite $A_{2^k}$ puis utiliser là pour les représenter graphiquement.
Commenter la courbe obtenue.
# Cellule à compléter
from math import sqrt
def EstDePi(k) :
# Cellule pour tester
# Cellule à compléter
def EstDePi2(N) :
# Cellule pour tester
# Cellule à compléter
import matplotlib.pyplot as plt
plt.figure() # Initialise une nouvelle figure
u =
plt.plot(u,'r')
plt.show()
On voit qu'à partir de $n=25$ (à peu près) les termes calculés ne sont plus du tout corrects: ils sont proches de $0$ et non de $\pi$. En effet lorsque $s_k^2$ devient plus petit que la précision machine, on a $1-s_k^2=1$ et donc $s_{k+1}=0$. Le suite devient stationnaire égale à $0$ et donc on a aussi $A_{2^k}=0$.
On s’intéresse à la convergence suivante : $$\lim_{n\to+\infty}\left(1+\dfrac{1}{n}\right)^n=e$$
Q5. Calculer la quantité $\left|\left(1 + \dfrac{1}{n}\right)^n - e\right|$ pour $n = 10^k$ et $k = 1$, $\dots$, $16$.
Une représentation semi-logarithmique de l’erreur relative
d’estimation de $e$ semble pertinente ...
# Cellule à compléter
from math import exp
u =
abscisses =
plt.figure() # Initialise une nouvelle figure
plt.semilogx(abscisses, u, 'o-r') # plot semilog suivant x
# on obtient la courbe de u en fonction de ln(abscisses)
Q6. Expliquer le phénomène observé.
Soit la suite définie par :
$$S_n = \int_0^\pi \left(\dfrac{x}{\pi}\right)^{2n}\sin(x)\hbox{d}x$$
Remarquons que la suite $(Sn )$ est positive puisque la fonction intégrée est positive sur l'intervalle.
On démontre la relation de récurrence suivante :
$$\forall n\geq1,\quad S_n=1-\dfrac{2n(2n - 1)}{\pi^2}S_{n-1}$$
On peut aussi montrer que la suite $(S_n )$ tend vers $0$ lorsque $n$ tend vers $+\infty$.
De plus $S_0=\displaystyle\int_0^\pi\sin(x)\hbox{d}x=\big[-\cos(x)\big]_0^\pi=2$.
Q7. Ecrire une fonction S(n)
qui renvoie la liste des valeurs de $S_0$, $S_1$, $\dots$, $S_n$ calculées avec la formule de récurrence.
Commenter le résultat obtenu pour $n=30$.
# Cellule à compléter
from math import pi
def S(n) :
# Cellule pour tester
La suite $(|S_n|)$ semble tendre vers $+\infty$ et non pas vers $0$. C'est embêtant.
Pour terminer ce TP, on souhaite analyser l’effet d’accumulation d’erreurs associées à l’absorption. Pour cela, intéressons-nous à la série harmonique.
Le terme général $u_n$ de la série harmonique est défini par $u_n = \dfrac{1}{n}$. La $n$-ième somme partielle de la série harmonique a pour expression : $$S_n=\displaystyle\sum_{k=1}^n\dfrac{1}{k}$$
Pour les deux questions suivantes, vous compléterez le tableau ci-dessous.
Q8. Ecrire le code python qui retourne le résultat de la somme pour des indices croissants $k = 1$, $\dots$, $N$.
# Cellule à compléter
def somme(N):
Q9. Modifier le code précédent de manière à réaliser la sommation pour des indices décroissants $k = N$, $\dots$, $1$.
# Cellule à compléter
def somme2(N):
# Cellule pour tester
Un algorithme « bien pensé » peut limiter voire compenser les erreurs. Comparativement à l’approche naïve précédente, l’algorithme de Kahan (ou sommation compensée), aussi connu sous le nom patronymique conjoint de Kahan-Babuška, réduit de manière significative l’erreur numérique sur la somme d’une séquence de réels représentés en virgule flottante, et d’améliorer ainsi la précision. Voici Le pseudocode de cet algorithme est donné dans la suite.
Q10. Analyser le pseudo-code de l’algorithme de Kahan.
Q11. Décliner son code Python puis recalculer la somme pour les valeurs précédentes de N.
Analyser les résultats et conclure.
# Cellule à compléter
def kahan(N) :
# Cellule pour tester