Stabilité d'un schéma numérique

En analyse numérique, la stabilité d’un schéma numérique aux différences finies est une propriété globale de l’algorithme qui en découle. Elle concerne essentiellement le comportement numérique qui se manifeste lorsque les pas de discrétisation ( $Δ t$ , $Δ x$ , etc.) tendent tous vers 0.

Sous certaines hypothèses, le théorème de Lax montre que la stabilité est une condition nécessaire et suffisante pour assurer la convergence.

Bien qu’un schéma numérique soit conçu pour tenter de résoudre un problème décrit par des équations aux dérivées partielles, la stabilité du schéma n’a aucun lien avec la solution exacte du problème.

La stabilité d’un schéma ne doit pas être confondue avec la stabilité de la solution du problème d’origine (par exemple la stabilité de Lyapunov des systèmes dynamiques).

La stabilité dans le théorème de Lax

Considérons un problème supposé être bien posé qui modélise un système évolutif caractérisé par

une condition initiale précisant son état d’origine (variables spatiales en $t = 0$ ),
des équations aux dérivées partielles et des conditions de bord auxquelles est soumis l’état du système au cours de son évolution.

Dans ce contexte, un schéma numérique procède de la manière suivante :

Discrétisation des variables spatiales (pas $Δ x$ ) pour établir une approximation numérique de l’ état d’origine.
Discrétisation de la variable temporelle (sur $[0, \, T]$ avec un pas $Δ t$ ) pour entreprendre un processus se déroulant par étapes successives au cours desquelles l’état numérique se transforme.

Notons $C (Δ t)$ l’opérateur de modification de l’état discret au cours d’une étape, ceci en supposant une relation liant $Δ x$ à $Δ t$ qui contraint $Δ x$ à converger vers 0 lorsque $Δ t$ fait de même.

La stabilité est définie par la propriété suivante ^[1] :

Il existe $\displaystyle \tau > 0$ tel que l’ensemble des opérateurs $C m (Δ t)$ est uniformément borné,

pour tout $\textstyle 0 < \Delta t \leqslant \tau$ et tout $\textstyle m \leqslant T / \Delta x$ .

On remarque que la stabilité est une qualité intrinsèque du schéma numérique ( $T$ est le seul élément du problème intervenant dans cette définition).

Lorsque cette propriété n’est pas satisfaite, le schéma est instable.

Exemple concret

Un problème typique qui peut être résolu par une méthode des différences finies est l’équation de la chaleur (présentée ici sous sa forme épurée) :

$\partial_t U (x,\,t) = \partial_x^2 U(x,\,t)$ pour $(x,\,t) \, \in \, [0,\,1]\times[0,\,T]$ , avec

$U(x,\,0) = U_0(x)$ donné comme condition initiale et

$U(0,\,t) = U(1,\,t) = 0$ pour les conditions aux limites.

Cette formulation modélise l’évolution de la température $\displaystyle U$ d’un barreau de métal (une dimension) isolé, préalablement chauffé et dont les extrémités sont maintenues à température nulle.

Afin de résoudre numériquement ce problème, considérons successivement deux schémas pour lesquels la dérivée temporelle est traitée par le schéma d'Euler. A partir des pas $\displaystyle \Delta x = 1 / N$ et $\displaystyle \Delta t = T / M$ , notons $u_i^j$ l’approximation de $U(i \, \Delta x,\, j \, \Delta t)$ .

Schéma explicite

Définissons $\displaystyle k = \Delta t / \Delta x^2$ et considérons le schéma explicite (ou progressif) suivant :

$u_i^0 \, = \, U_0(i \, \Delta x)$ pour tout $0 \leqslant i \leqslant N,$

$u_i^{j+1} \, = \, u_i^j + k \, (u_{i-1}^j \, - \, 2 \, u_i^j \, + \, u_{i+1}^j)$ pour tout $\displaystyle 0 < i < N,$

et $u_0^{j+1} \, = \, u_N^{j+1} \, = \, 0.$

Lorsque $k \leqslant 1/2$ , le schéma est stable pour la norme spatiale $\|.\|_\infty$ définie par

$\| u_.^j \|_\infty = \max_{0 \leqslant i \leqslant N} u_i^j.$

Preuve de la stabilité du schéma si $\, k \leqslant 1/2$

Avec l’hypothèse sur $\displaystyle k$ , il suffit de remarquer que $u_i^{j+1}$ est une pondération (avec des poids positifs ou nuls et de somme égale à 1) des valeurs $u_.^j$ de l’étape précédente. Ainsi $| \, u_i^{j+1} \, | \leqslant \| u_.^j \|_\infty$ pour tout $i$ .

Cependant, lorsque $\displaystyle k > 1/2$ , le schéma est instable pour toute norme spatiale.

Preuve de l’instabilité du schéma si $\displaystyle \, k > 1/2$

Notons $\vec u_j$ le vecteur dont les composantes sont les $u_i^j$ pour $0 < i < N$ .

Le schéma peut s’écrire $\vec u_{j+1} = A \, \vec u_j$ où $\displaystyle A$ est une matrice carrée symétrique de taille $N - 1$ (dont les valeurs propres sont réelles).

Pour prouver l’instabilité du schéma, il suffit de trouver un vecteur propre $\displaystyle \vec v$ de $\displaystyle A$ dont la valeur propre $\displaystyle \lambda_N$ est de module supérieur à $1$ . Ainsi, puisque $A \, \vec v = \lambda_N \, \vec v$ , alors $A^M \, \vec v = (\lambda_N)^M \, \vec v$ et par conséquent $\|A^M\| \geqslant |\lambda_N|^M.$

En choisissant $\displaystyle v_i = (-1)^i sin(i \pi/N)$ , on vérifie $\displaystyle v_{i-1} + v_{i+1} = - 2 \, v_i \, cos(\pi / N)$ . Ainsi, $\displaystyle \vec v$ est vecteur propre de $\displaystyle A$ pour la valeur propre

$\displaystyle \lambda_N = 1 - 2 \, k \, (1 + cos(\pi / N)) = 1 - 4 \, k \, cos^2(\pi / 2 N).$

Avec l’hypothèse sur $\displaystyle k$ (une valeur fixée) :

$\lim_{N \to \infty} \lambda_N = 1 - 4 \, k < -1$ et la suite $\|A^M\| \geqslant |\lambda_N|^M$ n’est pas bornée.

Remarque concernant l’instabilité

L’instabilité d’un schéma n’implique pas nécessairement que son application dans un cas particulier conduise nécessairement à une divergence.

L’exemple de conditions initiales nulles ( $U 0 (x) = 0$ ) le prouve puisque les solutions numérique et analytique sont identiquement nulles.

Un exemple plus instructif montre que, pour une version instable du schéma précédent, il existe des conditions initiales $U 0 (x)$ régulières pour lesquelles il y a convergence. En réalité, cette convergence n’est que « théorique», car il faut traiter les calculs numériques avec une précision infinie pour l’obtenir.

En pratique, même à l’aide d’outils de calcul offrant une grande précision relative, une dégénérescence se manifeste tôt ou tard. La figure ci-dessous illustre ce phénomène pour le problème précédent :

Preuve de la convergence théorique du schéma

Choisissons la condition initiale $U_0(x) = sin(\pi \,x)$ pour laquelle on vérifie que la solution analytique est la suivante :

$U(x,\, t) = sin(\pi \,x) e^{- \lambda \, t}$

où $\displaystyle \lambda = \pi^2$ .

Quels que soient $\displaystyle \Delta x = 1/N$ et $\displaystyle \Delta t = T/M$ , on définit toujours $\displaystyle k = \Delta t / \Delta x^2$ qui est supposé constant.

Les conditions initiales conduisent à $\displaystyle u_i^0 = sin(i \pi/N)$ .

On vérifie ensuite par induction que $\displaystyle u_{i-1}^j + u_{i+1}^j = 2 \, u_i^j \, cos(\pi / N)$ , ce qui entraîne

$\displaystyle u_i^{j+1} = \mu(k, \, N) \, u_i^j$ où $\mu(k, \, N) = 1 - 2 \, k \, (1 - cos(\pi / N)).$

A ce stade, on a montré que le rapport entre les solutions analytique et numérique s’écrit pour tout $i$ :

$\frac{u_i^j}{U(i \Delta x,\, j \Delta t)} = \phi^j(\Delta t)$ où $\phi(\Delta t) = e^{\lambda \,\Delta t} \, \mu (k, \, N).$

Montrons finalement $\lim_{\Delta t \to 0} \phi^j(\Delta t) = 1$ (uniformément en j) en utilisant le développement limité de Taylor des deux termes dont $\displaystyle \phi(\Delta t)$ est le produit :

$e^{\lambda \,\Delta t} = 1 + \pi^2 \Delta t + O(\Delta t^2),$

$\displaystyle cos(\theta) = 1 - \theta^2 / 2 + O(\theta^4),$

$\displaystyle cos(\pi / N) = cos(\pi (\Delta t / k)^{1/2}) = 1 - (\pi^2 \Delta t) / (2 k) + O(\Delta t^2),$

$\mu (k, \, N) = 1 - 2 \, k \, (1 - cos(\pi (\Delta t / k)^{1/2})) = 1 - \pi^2 \Delta t + O(\Delta t^2),$

$\displaystyle \phi(\Delta t) = (1 + \pi^2 \Delta t + O(\Delta t^2)) (1 - \pi^2 \Delta t + O(\Delta t^2)) = 1 + O(\Delta t^2),$

Par conséquent :

$| \phi^j(\Delta t) - 1| \leqslant |\phi^M(\Delta t) - 1| = (1 + O(\Delta t^2))^{(T / \Delta t)} - 1$

dont on déduit (quel que soit le choix de $\displaystyle k > 0$ ) : $\lim_{\Delta t \to 0} \phi^j(\Delta t) = 1$ avec une convergence en $\displaystyle O(\Delta t).$

Instabilité numérique d’un schéma appliqué dans un cas particulier où il est « théoriquement » convergent

Ce résultat est obtenu en résolvant le problème de la chaleur avec la condition initiale $U_0(x) = sin(\pi \,x)$ pour laquelle le schéma précédent est théoriquement convergent (19 pas d’espaces, k = 0.7, calculs avec 16 digits significatifs).

Pour les 50 premiers pas temporels, les résultats sont proches de la solution analytique. Après 60 pas de temps, les premières irrégularités se manifestent, puis le chaos s’installe rapidement.

Ce comportement s’explique sans doute par le fait que des erreurs négligeables du calcul apportent une très légère contribution à des composantes typiquement instables du schéma. On reconnaît d’ailleurs à la dernière étape visualisée les oscillations typiques de la fonction utilisée dans la preuve de l’instabilité (menu déroulant ci-dessus).

Schéma implicite

Toujours avec $\displaystyle k = \Delta t / \Delta x^2$ , considérons le schéma implicite (ou rétrograde) suivant :

$u_i^0 \, = \, U_0(i \, \Delta x)$ pour tout $0 \leqslant i \leqslant N,$

$u_i^{j+1} \, = \, u_i^j + k \, (u_{i-1}^{j+1} \, - \, 2 \, u_i^{j+1} \, + \, u_{i+1}^{j+1})$ pour tout $\displaystyle 0 < i < N,$

et $u_0^{j+1} \, = \, u_N^{j+1} \, = \, 0.$

Pour sa mise en œuvre, ce schéma nécessite de résoudre à chaque étape temporelle un système linéaire dont la matrice est symétrique tridiagonale et à à diagonale dominante (donc régulière). La matrice étant la même à chaque pas, une seule décomposition (LU, QR, Cholesky, etc) est suffisante.

Pour la norme spatiale $\|.\|_\infty$ , ce schéma est stable pour toute valeur de $\displaystyle k$ .

Preuve de la stabilité du schéma implicite

Notons $a$ (respectivement $b$ ) la valeur de l’indice $i$ pour laquelle $u_i^{j+1}$ est maximal (resp. minimal). On vérifie que le terme multiplié par $\displaystyle k$ est positif pour $i = a$ et négatif pour $i = b$ . Ainsi

$u_a^{j+1} \leqslant u_a^j,$

$u_b^{j+1} \geqslant u_b^j.$

Ces deux inégalités respectent d’ailleurs la physique de la diffusion de la chaleur.

Par induction, on en déduit $\| u_.^j\|_\infty \leqslant \| u_.^0\|_\infty$ , et donc la stabilité du schéma.

Schéma de Crank-Nicholson

Ce schéma est défini en choisissant comme second membre la moyenne des seconds membres respectifs des deux schémas précédents, soit

$u_i^{j+1} \, = \, u_i^j + (k/2) \, \{ (u_{i-1}^{j+1} \, - \, 2 \, u_i^{j+1} \, + \, u_{i+1}^{j+1}) + (u_{i-1}^j \, - \, 2 \, u_i^j \, + \, u_{i+1}^j)\}.$

Référence

↑ (en) P. D. Lax, R.D. Richtmyer, « Survey of the stability of linear finite difference equations », dans Comm. Pure Appl. Math. (en), vol. 9, 1956, p. 267-293 [texte intégral, lien DOI]

Voir aussi

Portail de l'analyse

Catégorie :

Analyse numérique

Wikimedia Foundation. 2010.

Contenu soumis à la licence CC-BY-SA. Source : Article Stabilité d'un schéma numérique de Wikipédia en français (auteurs)

Игры ⚽ Нужно сделать НИР?

Regardez d'autres dictionnaires:

Stabilite (d'un schema numerique) — Stabilité d un schéma numérique Un schéma numérique est dit stable lorsque la différence entre la solution numérique obtenue et la solution exacte des équations discrétisées est bornée. La solution numérique doit vérifier l équation d origine… … Wikipédia en Français
Stabilité (d'un schéma numérique) — Stabilité d un schéma numérique Un schéma numérique est dit stable lorsque la différence entre la solution numérique obtenue et la solution exacte des équations discrétisées est bornée. La solution numérique doit vérifier l équation d origine… … Wikipédia en Français
Stabilité numérique — En analyse numérique, une branche des mathématiques, la stabilité numérique est une propriété globale d’un algorithmique numérique, une qualité nécessaire pour espérer obtenir des résultats ayant du sens. Une définition rigoureuse de la stabilité … Wikipédia en Français
Stabilité — Cette page d’homonymie répertorie les différents sujets et articles partageant un même nom. Sur les autres projets Wikimedia : « Stabilité », sur le Wiktionnaire (dictionnaire universel) Stabilité peut se référer à : Stabilité … Wikipédia en Français
Chimie Numérique — La chimie numérique (désignée parfois aussi par l anglicisme chimie computationnelle) est une branche de la chimie et/ou de la physico chimie qui utilise les lois de la chimie théorique exploitées dans des programmes informatiques spécifiques… … Wikipédia en Français
Chimie numerique — Chimie numérique La chimie numérique (désignée parfois aussi par l anglicisme chimie computationnelle) est une branche de la chimie et/ou de la physico chimie qui utilise les lois de la chimie théorique exploitées dans des programmes… … Wikipédia en Français
Chimie numérique — La chimie numérique (désignée parfois aussi par l anglicisme chimie computationnelle) est une branche de la chimie et/ou de la physico chimie qui utilise les lois de la chimie théorique exploitées dans des programmes informatiques spécifiques… … Wikipédia en Français
Méthodes de chimie numérique dans les systèmes solides — Chimie numérique La chimie numérique (désignée parfois aussi par l anglicisme chimie computationnelle) est une branche de la chimie et/ou de la physico chimie qui utilise les lois de la chimie théorique exploitées dans des programmes… … Wikipédia en Français
DÉRIVÉES PARTIELLES (ÉQUATIONS AUX) - Analyse numérique — Plus peut être que tout autre domaine des mathématiques, les équations aux dérivés partielles étaient prédisposées à bénéficier de l’utilisation des ordinateurs, pour de nombreuses raisons. La plus importante est leur intervention dans de… … Encyclopédie Universelle
Méthode des différences finies — En analyse numérique, la méthode des différences finies est une technique courante de recherche de solutions approchées d équations aux dérivées partielles qui consiste à résoudre un système de relations (schéma numérique) liant les valeurs des… … Wikipédia en Français

Dictionnaires et Encyclopédies sur 'Academic'

Stabilité d'un schéma numérique

Sommaire

La stabilité dans le théorème de Lax