Stabilité de Von Neumann

La stabilité d'un schéma numérique est intimement liée à l'erreur numérique. Un schéma de différences finies est dit stable si les erreurs commises en un pas de temps ne font pas augmenter les erreurs au fil des itérations. Si les erreurs diminuent et finissent par s'estomper, le schéma numérique est dit stable. Si au contraire, l'erreur croît à chaque itération, le schéma est dit instable. La stabilité d'un schéma peut être déterminée grâce à l'analyse de von Neumann. Pour des problèmes dépendants du temps, la stabilité garantit qu'une méthode numérique produise une solution bornée lorsque la solution de l'équation différentielle exacte est bornée. La stabilité d'un schéma peut être ardue, surtout lorsque l'équation considérée est non linéaire.

Dans certains cas, la stabilité de von Neumann est suffisante et nécessaire pour la stabilité au sens de Lax–Richtmyer (comme utilisée dans le théorème de Lax): l'EDP et le schéma aux différences finies sont linéaires ; l'EDP est à coefficients constants avec des conditions de bord périodiques et dépend de deux variables indépendantes ; et le schéma n'utilise pas plus de deux niveaux de temps[4]. La stabilité de von Neumann est nécessaire dans une bien plus vaste variété de cas. La relative simplicité de cette méthode implique qu'elle est souvent utilisée à la place d'une analyse de stabilité plus détaillée afin de donner une bonne idée des restrictions sur les tailles des pas.

La méthode de von Neumann est basée sur la décomposition de l'erreur en séries de Fourier. Afin d'illustrer cette procédure, considérons l'équation de la chaleur uni-dimensionnelle

${\displaystyle {\frac {\partial u}{\partial t}}=\alpha {\frac {\partial ^{2}u}{\partial x^{2}}}}$

définie sur ${\displaystyle L}$, qui peut dans ce cas être discretisée de la manière suivante

${\displaystyle \quad (1)\qquad u_{j}^{n+1}=u_{j}^{n}+r\left(u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}\right)}$

où

${\displaystyle r={\frac {\alpha \,\Delta t}{\Delta x^{2}}}}$

et la solution ${\displaystyle u_{j}^{n}}$ de l'équation discrète approxime la solution analytique ${\displaystyle u(x,t)}$ de l'EDP sur les points de grille.

Définissons l'erreur d'arrondi ${\displaystyle \epsilon _{j}^{n}}$ par

${\displaystyle \epsilon _{j}^{n}=N_{j}^{n}-u_{j}^{n}}$

où ${\displaystyle u_{j}^{n}}$ est la solution de l'équation discrétisée (1) qui serait implémentée en l'absence d'erreurs d'arrondi, et ${\displaystyle N_{j}^{n}}$ est la solution numérique obtenue par précision arithmétique finie. Comme la solution exacte ${\displaystyle u_{j}^{n}}$ doit vérifier la solution discrétisée, l'erreur ${\displaystyle \epsilon _{j}^{n}}$ doit elle aussi vérifier l'équation discrétisée[5]. Ainsi

${\displaystyle \quad (2)\qquad \epsilon _{j}^{n+1}=\epsilon _{j}^{n}+r\left(\epsilon _{j+1}^{n}-2\epsilon _{j}^{n}+\epsilon _{j-1}^{n}\right)}$

est une relation de récurrence pour l'erreur. Les équations (1) et (2) montrent que l'erreur et la solution numérique ont le même comportement en fonction du temps. Pour des équations différentielles linéaires avec des conditions de bord périodiques, la variation spatiale de l'erreur peut être décomposée en une série de Fourier finie sur l'intervalle ${\displaystyle L}$, par

${\displaystyle \quad (3)\qquad \epsilon (x)=\sum _{m=1}^{M}A_{m}e^{ik_{m}x}}$

où le nombre d'onde ${\displaystyle k_{m}={\frac {\pi m}{L}}}$ avec ${\displaystyle m=1,2,\ldots ,M}$ et ${\displaystyle M=L/\Delta x}$. La dépendance du temps de l'erreur est incluse en supposant que l'amplitude de l'erreur ${\displaystyle A_{m}}$ est une fonction du temps. Sachant que l'erreur tend à croître ou décroître de manière exponentielle avec le temps, il est raisonnable de supposer que l'amplitude varie exponentiellement avec le temps ; d'où

${\displaystyle \quad (4)\qquad \epsilon (x,t)=\sum _{m=1}^{M}e^{at}e^{ik_{m}x}}$

où ${\displaystyle a}$ est une constante.

Comme l'équation des différences de l'erreur est linéaire, il est suffisant de considérer la croissance de l'erreur pour un terme choisi:

${\displaystyle \quad (5)\qquad \epsilon _{m}(x,t)=e^{at}e^{ik_{m}x}}$

Les caractéristiques de stabilité peuvent être étudiées en utilisant cette forme de l'erreur, sans perte de généralité. Pour trouver la variation de l'erreur en fonction du temps, on substitue l'équation (5) dans l'équation (2), après avoir noté que

• ${\displaystyle {\begin{aligned}\epsilon _{j}^{n}&=e^{at}e^{ik_{m}x}\\\epsilon _{j}^{n+1}&=e^{a(t+\Delta t)}e^{ik_{m}x}\\\epsilon _{j+1}^{n}&=e^{at}e^{ik_{m}(x+\Delta x)}\\\epsilon _{j-1}^{n}&=e^{at}e^{ik_{m}(x-\Delta x)},\end{aligned}}}$

pour arriver à (après simplification)

${\displaystyle \quad (6)\qquad e^{a\Delta t}=1+{\frac {\alpha \Delta t}{\Delta x^{2}}}\left(e^{ik_{m}\Delta x}+e^{-ik_{m}\Delta x}-2\right).}$

En utilisant les identités suivantes

${\displaystyle \qquad \cos(k_{m}\Delta x)={\frac {e^{ik_{m}\Delta x}+e^{-ik_{m}\Delta x}}{2}}\qquad {\text{et}}\qquad \sin ^{2}{\frac {k_{m}\Delta x}{2}}={\frac {1-\cos(k_{m}\Delta x)}{2}}}$

on peut réécrire (6) comme

${\displaystyle \quad (7)\qquad e^{a\Delta t}=1-{\frac {4\alpha \Delta t}{\Delta x^{2}}}\sin ^{2}(k_{m}\Delta x/2)}$

Définissons le facteur d'amplitude

${\displaystyle G\equiv {\frac {\epsilon _{j}^{n+1}}{\epsilon _{j}^{n}}}}$

La condition nécessaire et suffisante pour que l'erreur reste bornée est ${\displaystyle \vert G\vert \leq 1.}$ Cependant,

${\displaystyle \quad (8)\qquad G={\frac {e^{a(t+\Delta t)}e^{ik_{m}x}}{e^{at}e^{ik_{m}x}}}=e^{a\Delta t}}$

Ainsi, des équations (7) et (8) on tire que la condition de stabilité est donnée par

${\displaystyle \quad (9)\qquad \left\vert 1-{\frac {4\alpha \Delta t}{\Delta x^{2}}}\sin ^{2}(k_{m}\Delta x/2)\right\vert \leq 1}$

Afin que la condition ci-dessus soit vraie pour tout ${\displaystyle \sin ^{2}(k_{m}\Delta x/2)}$, on a

${\displaystyle \quad (10)\qquad {\frac {\alpha \Delta t}{\Delta x^{2}}}\leq {\frac {1}{2}}}$

L'équation (10) donne la condition de stabilité pour le schéma FTCS appliqué à l'équation de la chaleur uni-dimensionnelle. Elle dit que pour un ${\displaystyle \Delta x}$ donné, la valeur de ${\displaystyle \Delta t}$ doit être assez petite pour vérifier l'équation (10).