consistance vs stabilité d'un schéma numérique

L'analyse numérique  fait aujourd'hui partie intégrante des mathématiques, on l'enseigne dans de nombreux cursus, mais elle a mis du temps à s'imposer comme une domaine des mathématiques à part entière. Beaucoup ont considéré au début que les problèmes étudiés étaient trop appliqués pour permettre de dégager des concepts fondamentaux et de développer une approche suffisamment rigoureuse  pour faire partie des mathématiques. Lorsque j'étais étudiant  malgré des cours très complets sur le sujet  je ne me suis pas vraiment intéressé à l'analyse numérique car les concepts de base  m'en avaient échappé. Je n'ai commencé à m'intéresser à ce domaine  qu'en commençant à enseigner en DUT informatique  puis en cursus d'ingénieur , que je me suis rendu compte à quel point j’étais passé à côté de choses très intéressantes. J'aimerai  partager ici avec vous les quelques exemples  qui m'ont fait changer d'avis.

l'instabilité  du schéma explicite pour l'équation de Transport 1D  l'emporte  sur l'erreur  de consistance (écart entre les courbes rouges et bleu) à partir de t>1

Problème bien posé et schéma d'approximation numérique

Pour comprendre  l'analyse numérique il faut commencer par définir ce qu'est un schéma d'approximation numérique. Au départ  il faut comprendre que tout problème d'analyse numérique peut se voir d'un point de vu abstrait comme trouver une méthode de calcul approché  de $f(x)$

$f(x)$ = " la valeur d'une fonction $f$ pour des données numériques $x\in I\subset{\mathbb R}$" 

Avant de chercher à résoudre numériquement ce problème il faut s'assurer  que le problème est bien posé : que  résultat $f(x)$ est unique bien sûr mais aussi qu'il  dépend continûment des données $x$ de sorte que si on perturbe un peu les données le résultat  reste proche

$$\tilde{x}=x+\delta x\text{ alors }\displaystyle f(x)-f(\tilde{x})\mathop{\longrightarrow}_{\tilde{x}\to x}0$$

Quand c'est le cas on dit que le problème est  "bien posée au sens de Hadamard"  ou "bien conditionné". Ensuite on va donc remplacer le calcul  de $f(x)$ par celui de $f_n(x_j)$ où :

• on se restreint à des ensembles finis $(x_j)_{j\in J}$ de nombre caractérisés par un "pas" $\displaystyle \epsilon=\sup_{j\in J}\vert x_j-x_{j+1}\vert$ qu'on peut idéalement faire tendre vers 0
• on construit une suite d'applications $(f_n)_{\mathbb N}$ définies sur les $(x_j)_{j\in J}$  qui idéalement doit tendre vers $f$  quand $n\to\infty$

le but est d'obtenir la convergence du schéma vers la solution du problème de départ   :

$f_n(x_j)\to f(x)$  quand $n\to\infty$ et $\epsilon\to 0$ avec $x_j$  le point le plus proche de $x$ pour $j\in J$

Concrètement  je vais prendre un  exemple très simple :  calculer   la valeur de $e^x$    qui est celui qui m'a vraiment fait comprendre  les concepts de base de l'analyse numérique. On a alors que :

• ce problème revient à calculer la série $e^x=\sum_{k=0}^\infty {x^k\over k!}$
• c'est un problème bien posé car le résultat est unique et dépend continûment des données car  $\vert e^x-e^y\vert\leq e^M\vert x-y\vert$ pour $x,y\in [-M,M]$ (utiliser le TAF)
• le calcul numérique s'effectue en tronquant la série $\sum_{k=0}^n {x^k\over k!}=f_n(x)$  ce qui définit la suite d'applications $(f_n)$
• lors de l'application numérique  $x=\pi$ sera remplacé par son approximation par un nombre en virgule flottante $x_j$ le plus proche possible , par exemple avec des réels  double précision on aura $x_j=3.14159265358979312$  et le pas sera $\epsilon\approx 10^{-16}$

Les différentes erreurs numériques et le théorème de Lax

Lorsqu'on fait ce type de calculs on est confronté à trois types d'erreurs :

• les erreurs de données  qui peuvent souvent être  des données expérimentales comportant des imprécisions et qu'on ne peut pas  corriger
• l'erreur de méthode ou de consistance est celle qu'on fait quand on remplace $f(x_j)$  par $f_n(x_j)$  qui est  intrinsèque à l'algorithme utilisé , qu'on doit pouvoir majorer en fonction de $n$, et lorsqu'elle tend vers 0  avec $n$ on dira que le schéma est consistant
• les erreurs de troncature ou d'arrondis  qui vont apparaître  à chaque calcul de la valeur $f_n(x_j)$, elles sont faibles mais il y en a beaucoup , si on arrive à montrer  que le schéma ne les amplifie pas c'est à dire qu'elles restent inférieures à $C\times \epsilon$ alors le schéma est stable

Dans notre exemple du calcul de $e^\pi$

•  erreur de donnée  celle qui est faite quand on remplace $\pi$ par 3.14159265358979312  elle est d'environ $ 10^{-16}$ 
• l'erreur de  consistance s'évalue en calculant le reste de la série
  $$\displaystyle e^x-\sum_{k=0}^n {x^k\over k!}=\sum_{k=n+1}^\infty {x^k\over k!}\leq {M^{n+1}\over n!\vert n+1-M\vert}\neq 0$$
• les erreurs de troncature   elles  sont de l'ordre de $\epsilon$ sur chaque calcul donc au total
  $$ \left\vert \sum_{k=0}^n {x^k\over k!}\epsilon \right\vert
  \leq  e^M \times\epsilon  \mathop{\longrightarrow}_{\epsilon\to0}0 $$

On peut maintenant énoncer le résultat fondamental de l'analyse numérique:

Théorème de Lax Pour un problème bien posé  un schéma consistant et stable  est  convergent

La démonstration du théorème est très simple à écrire :

$$\begin{align*}
\vert f(x)-f_n(x_j)\vert
& = \vert f(x)-f(x_i)+f(x_i)-f_n(x_i)+f_n(x_i)-f_n(x_j)\vert \\
& \leq  \underbrace{\vert f(x)-f(x_i)\vert}_\text{bien posé $(\epsilon\to 0)$}
 +\underbrace{\vert f(x_i)-f_n(x_i)\vert}_\text{consistant $(n\to\infty)$}
 +\underbrace{\vert f_n(x_i)-f_n(x_j)\vert }_\text{stable $(\epsilon\to 0)$}\\
& \mathop{\longrightarrow}_{n\to\infty,~\epsilon\to 0}0+0+0=0 \text{ le schéma  converge!}
\end{align*}$$

Dans notre exemple  du calcul de $e^\pi$  toutes les conditions sont réunies pour appliquer le théorème de Lax, c'est pourquoi  on peut espérer  obtenir une bonne valeur approchée de $e^\pi$ par un calcul de la somme partielle de la série comme le montre les calculs ci-dessous faits avec scilab :

-->k=[1:20];
 
-->epi=1+sum((%pi.^k)./cumprod(k))
 epi  =
 
    23.140693  
 
-->exp(%pi)
 ans  =
 
    23.140693  
 
-->erreur=abs(epi-exp(%pi))
 erreur  =
 
    6.284D-10  
 
-->%eps
 %eps  =
 
    2.220D-16  

On peut même préciser l'erreur globale commise  avec notre méthode d'approximation :

$${\mathcal E}_{n,\epsilon}= \underbrace{M^{n+1}\over n!\vert n+1-M\vert}_\text{consistance}+ \underbrace{e^M \times\epsilon}_\text{troncature}$$

On peut ainsi se rendre compte que suivant la valeur de $\epsilon$ il existe une valeur optimale du paramètre $n$  à partir de laquelle  on améliore plus la précision du résultat. Dans notre cas  lorsque l'erreur de consistance devient plus faible que l'erreur de troncature et on ne gagne rien à calculer de nouveaux termes dans la série :
$${\vert M\vert^{n+1}\over n!\vert n+1-M\vert}\leq e^M \times 10^{-16} $$

• pour le calcul de $e^1$  on prend $M=1$  et on trouve que $n=16$
• pour le calcul de $e^\pi$  on prend $M=3$  et on trouve que $n=25$

 on peut  vérifier en traçant le graphe du logarithme de l'erreur  en fonction de $n$  que  cela donne assez précisément le $n$ optimal :

Conclusion 

Dans la pratique  la recherche d'un schéma d'approximation numérique commence par la construction d'une méthode consistante, dont il faut ensuite prouver la stabilité. Lorsque la stabilité du schéma  n'est pas vérifiée  les erreurs  d'arrondi vont s'accumuler et s'amplifier, le schéma instable n'est pas utilisable même s'il est consistant ! Ça se voit très bien sur l'exemple du schéma  explicite décentré l'équation de transport 1D, on a une erreur de troncature de la forme $a^n\epsilon$, avec $a>1$ ,de sorte que quand $n\to\infty$ l'erreur va tendre vers l'infini. C'est ce qui arrive dans cette animation à partir de $t>1$. Avant ce temps l’écart entre la solution exacte (en rouge) et numérique (en bleu)  est donné par l'erreur de consistance, mais  ensuite l'erreur de troncature prend le dessus et augmente exponentiellement ...