Stabilité d'un schéma numérique

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

Par conséquent ${\displaystyle \|C(\Delta t)\|_{\infty }\leqslant 1}$ et ${\displaystyle \|C^{m}(\Delta t)^{n}\|_{\infty }\leqslant (\|C(\Delta t)\|_{\infty })^{m}\leqslant 1.}$

Notons ${\displaystyle {\vec {u}}_{j}}$ le vecteur dont les composantes sont les ${\displaystyle u_{i}^{j}}$ pour ${\displaystyle 0<i<N}$.

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

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

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

${\displaystyle \displaystyle \lambda _{N}=1-2\,k\,(1+cos(\pi /N))=1-4\,k\,cos^{2}(\pi /2N).}$

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

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

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

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

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

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

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

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

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

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

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

Montrons finalement ${\displaystyle \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 \displaystyle \phi (\Delta t)}$ est le produit :

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

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

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

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

${\displaystyle \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 :

${\displaystyle |\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 \displaystyle k>0}$) : ${\displaystyle \lim _{\Delta t\to 0}\phi ^{j}(\Delta t)=1}$ avec une convergence en ${\displaystyle \displaystyle O(\Delta t).}$

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

${\displaystyle u_{a}^{j+1}\leqslant u_{a}^{j},}$

${\displaystyle 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 ${\displaystyle \|u_{.}^{j}\|_{\infty }\leqslant \|u_{.}^{0}\|_{\infty }}$, et donc la stabilité du schéma.