1.5 Étude de la convergence du schéma explicite

Pour étudier la convergence de la solution numérique $u_{i}^{n}$ du schéma aux différences finies (1.34) vers la solution exacte $u(x,t)$ de l'équation (1.1), on peut soit essayer d'étudier directement la convergence, ce qui peut être compliqué et pas toujours possible, soit utiliser un résultat d'analyse du à Lax (Richtmyer et Norton 1967):

Théorème de Lax

pour un problème linéaire aux valeurs initiales, la solution numérique d'un schéma itératif en temps aux différences finies converge vers la solution exacte si le schéma est consistant et stable.

La consistance et la stabilité d'un schéma sont en général beaucoup plus facile à étudier que sa convergence.

1.5.1 Étude de la consistance du schéma explicite

La consistance caractérise la façon dont l'équation aux différences finies (EDF) approche l'équation aux dérivées partielles (EDP).

consistance

un schéma aux différences finies est consistant à l'équation exacte si l'EDF tends vers l'EDP lorsque les pas de discrétisation en temps et en espace tendent vers zéro indépendamment.

Pour formaliser cette approche, notons $u(x,t)$ la solution de l'EDP, i.e $u(x,t)$ vérifie exactement l'équation: $EDP[u]=0$ , et $u_{i}^{n}$ la solution de l'EDF, i.e. $u_{i}^{n}$ vérifie exactement l'équation: $EDF[u]=0$ . Nous allons calculer la différence $E_{t}$ entre l'EDP et l'EDF aux noeuds du maillage et pour la solution exacte $u$ :

$$E_{t}=EDF[u(i\Delta x,n\Delta t)]-EDP[u]_{(i\Delta x,n\Delta t)}$$

ce qui donne, compte tenu du fait que $EDP[u]=0$

$$E_{t}=EDF[u(i\Delta x,n\Delta t)]$$

Cette différence est l'erreur de troncature du schéma aux différences finies, et corresponds donc à l'erreur commise lorsque l'on remplace la solution approchée $u_{i}^{n}$ par la solution exacte aux noeuds du maillage $u(i\Delta x,n\Delta t)$ dans l'équation aux différences EDF.

Le schéma EDF est dit consistant à l'équation EDP si cette erreur de troncature tends vers zéro lorsque le pas de discrétisation en temps $\Delta t$ et le pas de discrétisation en espace $\Delta x$ tendent vers zéro indépendamment.

$$\lim_{\Delta x\rightarrow0, \Delta t\rightarrow0}E_{t}=EDF[u(i\Delta x,n\Delta t)]=0$$

Appliquons cette définition au schéma explicite (1.34), en notant $u(i\Delta x,n\Delta t)=U_{i}^{n}$ la solution exacte aux noeuds du maillage:

$$E_{t}=\frac{U_{i}^{n+1}-U_{i}^{n}}{\Delta t}-\kappa \frac{U_{i+1}^{n}-2U_{i}^{n}+U_{i-1}^{n}}{\Delta x^{2}}$$

Pour calculer cette erreur de troncature, on exprime toutes les quantités en fonction de $U_{i}^{n}$ en effectuant des développement en série de Taylor au voisinage du point $(i\Delta x,n\Delta t)$ :

$$U_{i}^{n+1}=U_{i}^{n}+\left.\frac{\Delta t}{1!}\frac{\partial U}{... ...c{\partial^{2}U}{\partial t^{2}}\right\vert _{\left(i,n\right)}+O(\Delta t^{3})$$

$$U_{i+1}^{n}=U_{i}^{n}+\left.\frac{\Delta x}{1!}\frac{\partial U}{... ...c{\partial^{4}U}{\partial x^{4}}\right\vert _{\left(i,n\right)}+O(\Delta x^{5})$$

$$U_{i-1}^{n}=U_{i}^{n}-\left.\frac{\Delta x}{1!}\frac{\partial U}{... ...l^{4}U}{\partial x^{4}}\right\vert _{\left(i,n\right)}+O(\Delta x^{5})\textrm{}$$

En reportant ces développements dans l'erreur de troncature , il vient

$$E_{t}=\left.\frac{\partial U}{\partial t}\rig... ...^{4}}\right\vert _{\left(i,n\right)}\right)+O(\Delta t^{2},\Delta x^{4})$$

On utilise ensuite le fait que $U_{i}^{n}$ vérifie l'équation exacte au point $(i,n):$

$$\left.\frac{\partial U}{\partial t}\right\ver... ...}-\left.\kappa\frac{\partial^{2}U}{\partial x^{2}}\right\vert _{(i,n)}=0$$

ce qui donne l'expression de l'erreur de troncature du schéma explicite:

$$E_{t}=\left.\frac{\Delta t}{2}\frac{\partial^{2}U}{\partial t^{2}... ...\partial^{4}U}{\partial x^{4}}\right\vert _{(i,n)}+O(\Delta t^{2},\Delta x^{4})$$ (1.38)

Cette erreur de troncature tends bien vers zéro, lorsque $\Delta t$ et $\Delta x$ tendent vers zéro indépendamment.

Le schéma explicite 1.34 est donc consistant à l'équation de la chaleur (1.1) et l'erreur de troncature est en $O(\Delta t,\Delta x^{2})$ , i.e. d'ordre 1 en temps et d'ordre 2 en espace.

L'ordre de troncature est aussi la précision de la solution numérique, i.e l'erreur entre la solution numérique $u_{i}^{n}$ et la solution exacte $U_{i}^{n}$ est en $O(\Delta t,\Delta x^{2})$ .

Remarque: on peut transformer l'erreur de troncature (1.37) en dérivant l'équation exacte par rapport à t pour calculer $\frac{\partial u^{2}}{\partial t^{2}}$ :

$$\frac{\partial^{2}U}{\partial t^{2}}=\kappa\f... ...}{\partial x^{2}}\left(\kappa\frac{\partial^{2}U}{\partial x^{2}}\right)$$

En remplaçant dans (1.37) il vient:

$$E_{t}=\kappa\left(\kappa\frac{\Delta t}{2}-\f... ...c{\partial^{4}U}{\partial x^{4}}\right\vert+O(\Delta t^{2},\Delta x^{4})$$

on remarque que pour la valeur particulière de $\Delta t$ telle que:

$$\kappa\frac{\Delta t}{\Delta x^{2}}=\frac{1}{6}$$ (1.39)

le premier terme de l'erreur de troncature s'annulle, et la précision du schéma augmente pour être en $O(\Delta t^{2},\Delta x^{4})$ .

1.5.2 Etude de la stabilité du schèma explicite

La notion de stabilité s'applique à des schèmas, pour lesquels on calcule des solutions de façon itérative. Les calculs s'éffectuent sur des ordinateurs avec une précision finie, et donc sont sujet à des erreurs d'arrondis. Lors d'un calcul itératif, ces erreurs peuvent être amplifiées par le schèma numérique. Le but de l'étude de stabilité est donc de déterminer quelle est l'amplification des erreurs (ou perturbations) par le schèma.

stabilité

un schèma itératif est dit stable, si les perturbations de la solution numérique ne sont pas amplifiées au cours des itérations.

Soit $u_{i}^{n}$ la solution à l'étape $n$ de l'équation aux différences (1.34). la solution $u_{i}^{n+1}$ à l'étape $n+1$ vérifie:

$$u_{i}^{n+1}=u_{i}^{n}+r\left(u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}\right)$$ (1.40)

Soit $\epsilon_{i}^{n}$ une perturbation de la solution à l'étape $n$ . La solution perturbée $u_{i}^{n+1}+\epsilon_{i}^{n+1}$ à l'étape $n+1$ est solution de l'équation aux différences:

$$u_{i}^{n+1}+\epsilon_{i}^{n}=u_{i}^{n}+\epsilon_{i}^{n}+r\left(u_... ...n}-2\left(u_{i}^{n}+\epsilon_{i}^{n}\right)+u_{i-1}^{n}+\epsilon_{i}^{n}\right)$$ (1.41)

L'équation sur la perturbation est obtenue en effectuant la différence (1.40)-(1.39):

$$\epsilon_{i}^{n+1}=\epsilon_{i}^{n}+r\left(\epsilon_{i+1}^{n}-2\epsilon_{i}^{n}+\epsilon_{i-1}^{n}\right)$$ (1.42)

On décompose cette perturbation en tout point $x_{i}=i\Delta x$ du maillage et à tout instant $t_{n}$ sous la forme d'une série de modes de Fourier (en notant $I=\sqrt{-1}$ ):

$$\varepsilon_{i}^{n}=\sum\limits _{m=1}^{\infty}\psi_{m}^{n}e^{I\omega_{m}i\Delta x}$$ (1.43)

Le problème étant linéaire, chacun des modes vérifie l'équation (1.41), qui s'écrit pour un mode $m$ :

$$\psi_{m}^{n+1}e^{I\omega_{m}i\Delta x}=\psi_{m}^{n}e^{I\omega_{m}... ...{m}(i+1)\Delta x}-2e^{I\omega_{m}i\Delta x}+e^{I\omega_{m}(i-1)\Delta x}\right)$$ (1.44)

En simplifiant par $e^{I\omega_{m}i\Delta x}$ l'équation (1.43) s'écrit :

$$\psi_{m}^{n+1}=\psi_{m}^{n}\left[1+r\left(e^{I\omega_{m}\Delta x}-2+e^{-I\omega_{m}\Delta x}\right)\right]$$ (1.45)

En posant $\beta=\omega_{m}\Delta x$ et en utilisant les relations trigonométriques $\cos\beta=\frac{e^{I\beta}+e^{-I\beta}}{2}$ et $sin^{2}\frac{\beta}{2}=\frac{1}{2}\left(1-\cos\beta\right),$ l'équation (1.44) devient :

$$\psi_{m}^{n+1}=\psi_{m}^{n}\left(1-4r\sin^{2}\frac{\beta}{2}\right)$$ (1.46)

Il est clair que si l'amplitude $\psi_{m}^{n}$ de chaque mode de la perturbation diminue d'une itération à l'autre, l'erreur décroît et la perturbation finit par disparaître. Ceci se traduit par la condition:

$$\left\vert\frac{\psi_{m}^{n+1}}{\psi_{m}^{n}}\right\vert=\left\vert 1-4r\sin^{2}\frac{\beta}{2}\right\vert\leq1$$ (1.47)

qui doit être vérifiée pour chaque mode $m$ , i.e. $\forall\omega_{m}$ c'est à dire $\forall\beta$ .

Dans l'équation (1.46) on note par $G=\left\vert\frac{\psi_{m}^{n+1}}{\psi_{m}^{n}}\right\vert$ le facteur d'amplification du schèma.

Remarque: cette analyse de stabilité n'a pas prise en compte les conditions aux limites. Le type des conditions aux limites et leurs discrétisations sélectionnent certains modes de Fourier. Dans notre étude, la condition aux limites sur la perturbation est une condition homogéne de Dirichlet: $\varepsilon_{0}^{n}=\varepsilon_{N}^{n}=0$ , et elle impose des modes de Fourier en $\sin(m\pi\frac{i\Delta x}{L})$ pour vérifier ces conditions aux limites ( i.e. $\omega_{m}=\frac{m\pi}{L} m\in\mathbb{Z}$ ). On verra à la fin du chapitre une autre façon de faire l'analyse de stabilité: l'approche matricielle qui prend en compte directement les conditions aux limites.

L'inégalité (1.46), s'écrit:

$$-1\le1-4r\sin^{2}\beta\le1 \forall\beta$$

soit

$$\frac{1}{2\sin^{2}\beta}\ge r\ge0 \forall\beta$$

Cette inégalité impose une valeur maximale à $r$ :

$$r\le\frac{1}{2}$$ (1.48)

Le schèma explicite est donc conditionnellement stable, avec un critère de stabilité donné par la condition:

$$\kappa\frac{\Delta t}{\Delta x^{2}} \leq \frac{1}{2}$$ (1.49)

Cette condition s'utilise de la façon suivante: on fixe le maillage en espace, i.e. le pas $\Delta x$ et le pas en temps d'intégration en temps $\Delta t$ doit alors être choisi tel que:

$$\Delta t\le\frac{\Delta x^{2}}{2\kappa}$$

On remarque que plus on augmente le nombre de points en espace (i.e. plus $\Delta x$ diminue), plus il faut prendre un pas en temps $\Delta t$ petit. C'est l'un des principales inconvénients de ce schéma explicite.

1.5.3 Étude de la convergence du schéma explicite

Pour ce schéma explicite, il est possible d'étudier directement la convergence de la solution numérique vers la solution exacte. Pour cela nous allons déterminer la solution exacte de l'équation aux différences finies (1.35):

$$u_{i}^{n+1}=u_{i}^{n}+r\left(u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}\right)$$

Par analogie avec la solution exacte (1.17), nous allons chercher une solution $u_{i}^{n}$ sous la forme d'une série:

$$u_{i}^{n}=\sum_{m=1}^{M}C_{m}\Lambda_{m}^{n}\sin\left(\frac{m\pi}{L}i\Delta x\right)$$

où les coefficients $\Lambda_{m}^{n}$ donnent la dépendance temporelle et les $C_{m}$ prennent en compte la condition initiale.

En notant $N$ le nombre d'intervalle en espace, on a $\Delta x=L/N$ et la série s'écrit:

$$u_{i}^{n}=\sum_{m=1}^{M}C_{m}\Lambda_{m}^{n}\sin\left(m\pi\frac{i}{N}\right)$$

Cette série discréte ne comprends que $M=N-1$ modes indépendants. En effet pour $m>N$ , le mode $m$ est confondu aux noeuds du maillage avec un mode $q<N$ tel que $m=p N+q$ , puisque:

$$\sin\left((pN+q)\pi\frac{i}{N}\right)=sin(q\pi\frac{i}{N})$$

et pour $m=N$ on a un mode identiquement nul. Ce nombre de modes $M=N-1$ est en fait le nombre de degrés de liberté de la solution numérique: $\left\{ u_{i}^{n}\right\}$ est définie par $N+1$ valeurs aux points de maillage, mais on impose 2 conditions aux limites $u_{0}^{n}=0$ et $u_{N}^{n}=0$ , ce qui donne en définitive $N-1$ degrés de liberté.

La solution $u_{i}^{n}$ s'écrit donc:

$$u_{i}^{n}=\sum_{m=1}^{N-1}C_{m}\Lambda_{m}^{n}\sin\frac{m\pi i}{N}$$

Le problème étant linéaire, chacun des modes est solution de l'équation (1.35), ce qui donne en simplifiant par $C_{m}$ qui ne dépend ni de $i$ ni de $n$ :

$$\Lambda_{m}^{n+1}\sin\frac{m\pi i}{N} = \Lambda_{m}^{n}\sin\frac{m\pi i}{N}+r\left(\Lambda_{m}^{n}\sin\fr... ...ambda_{m}^{n}\sin\frac{m\pi i}{N}+\Lambda_{m}^{n}\sin\frac{m\pi(i-1)}{N}\right)$$

$$= \Lambda_{m}^{n}\sin(\frac{m\pi i}{N})\left(1+2r(\cos\frac{m\pi}{N}-1)\right)$$

$$= \Lambda_{m}^{n}\sin(\frac{m\pi i}{N})\left(1-4r\sin^{2}\frac{m\pi}{2N}\right)$$

d'où le coefficient $\Lambda_{m}^{n}$

$$\Lambda_{m}^{n}=\left(1-4r\sin^{2}\frac{m\pi}{2N}\right)^{n}$$

Les coefficients $C_{m}$ sont déterminés de façon à vérifier la condition initiale aux noeuds du maillage $u_{i}^{0}=C(i\Delta x)$ :

$$C(i\Delta x)=\sum_{m=1}^{N-1}C_{m}\sin\frac{m\pi i}{N} \forall i=1,N-1$$

Ce sont donc les coefficients de la transformée de Fourier de la condition initiale $C(x)$ calculée aux noeuds du maillage:

$$C_{m}=\frac{2}{N}\sum_{i=1}^{N-1}C(i\Delta x)\sin\frac{m\pi i}{N}$$

La solution numérique s'écrit donc:

$$u_{i}^{n} = \sum_{m=1}^{N-1}\left(C_{m}\left(1-4r\sin^{2}\frac{m\pi}{2N}\right)^{n}\sin\frac{m\pi i}{N}\right)$$ (1.50)

$$C_{m} = \frac{2}{N}\sum_{i=1}^{N-1}C(i\Delta x)\sin\frac{m\pi i}{N}$$

Comparons cette solution à la solution exacte (1.17) aux noeuds du maillage. La solution exacte au point $x=i\Delta x$ et à $t_{0}=n\Delta t$ s'écrit, en remplaçant $\kappa$ et $\Delta t$ en fonction de $r$ et $N$ :

$$u(i\Delta x,n\Delta t)=\sum\limits _{m=1}^{\infty}C_{m}e^{-rn\left(\frac{m\pi}{N}\right)^{2}}\sin\left(\frac{m\pi i}{N}\right)$$ (1.51)

En comparant les deux expressions (1.49) et (1.50), on voit que pour montrer la convergence de $u_{i}^{n}$ vers $u(x,t_{0})$ , il faut montrer que $\Lambda_{m}^{n}$ tends vers $e^{-rn\left(\frac{m\pi}{N}\right)^{2}}$ . Attention, on se place en un point $x=i\Delta x$ fixe et à un instant $t_{0}=n\Delta t$ fixé. Donc lorsque l'on fait tendre $\Delta x$ et $\Delta t$ vers zéro, les indices $i$ et $n$ augmentent en fonction de $\Delta x$ et $\Delta t$ . On introduit alors un petit paramètre $\delta=\frac{m\pi}{N}=\frac{m\pi h}{L}$ , et on effectue les développements limités par rapport à $\delta$ de $\Lambda_{m}^{n}$ :

$$\Lambda_{m}^{n} = \left(1-4r\sin^{2}\frac{\delta}{2}\right)^{n}$$ (1.52)

$$= 1-nr\delta^{2}+\left(\frac{n}{12}r+\frac{n^{2}-n}{2}r^{2}\right)\delta^{4}+O(\delta^{6})$$

et de $e^{-rn\left(\frac{m\pi}{N}\right)^{2}}$ :

$$e^{-rn\left(\frac{m\pi}{N}\right)^{2}}=1-nr\delta^{2}+\frac{n^{2}}{2}r^{2}\delta^{4}+O(\delta^{6})$$

L'écart entre ces 2 développements s'écrit en remplaçant $r$ et $\delta$ en fonction de $\Delta x$ et $\Delta t$ et en notant que $n=t_{0}/\Delta t$ :

$$\Lambda_{m}^{n}-e^{-rn\left(\frac{m\pi}{N}\right)^{2}} \approx \frac{nr}{12}\delta^{4}-\frac{nr^{2}}{2}\delta^{4}$$

$$\approx \frac{\kappa t_{0}}{12}\left(\frac{m\pi}{L}\right)^{4}\Delta x^{2}-\frac{\kappa^{2}t_{0}}{2}\left(\frac{m\pi}{L}\right)^{4}\Delta t$$

Cet écart tend bien vers zéro lorsque $\Delta t$ et $\Delta x$ tendent vers zéro, et est d'ordre 2 en espace et d'ordre 1 en temps.

A partir de l'expression (1.49), on peut aussi en déduire le comportement de la solution numérique. Pour chaque mode $m$ de cette solution, on approxime la dépendance temporelle en exponentielle décroissante de la solution exacte par une fonction puissance $\Lambda_{m}^{n}$ . Pour qu'il y ait convergence, il faut donc que $\left\vert\Lambda_{m}^{n}\right\vert\le1$ (c'est aussi la condition pour que le développement limité 1.51 converge). Cette condition implique:

$$\left\vert 1-4r\sin^{2}\frac{\delta}{2}\right\vert\le1$$

qui est justement la condition de stabilité (1.47) du schéma explicite:

$$r\le\frac{1}{2}$$

Nous avons tracé sur la figure (1.5) le facteur d'amplification $G_{m}$ du mode $m$

$$G_{m}=1-4r\sin^{2}\frac{m\pi}{2N}$$

en fonction de $\frac{m}{N}$ et pour différentes valeurs de $r$ .

Figure 1.5: Facteur d'amplification $G_{m}$

On note que ce facteur d'amplification est positif et inférieure à 1 pour les petits nombres d'ondes $m$ , et donc les grandes longueurs d'ondes $\frac{L}{m\pi}$ sont toujours amorties comme pour la solution exacte. Par contre le comportement des petites longueurs d'ondes (i.e. des grandes nombres d'ondes $m$ ) dépend de $r$ .

Pour $r>\frac{1}{2}$ , le facteur d'amplification des grands nombres d'ondes $m$ est négatif et inférieure à $-1$ , donc les petites longueurs d'ondes se mettent à osciller avec une amplitude croissante et le schéma diverge. La divergence du schéma explicite se traduit donc par l'apparition d'oscillations à hautes fréquences d'amplitudes croissantes.

Pour $\frac{1}{4}<r<\frac{1}{2}$ , le facteur d'amplification des grands nombres d'ondes $m$ est négatif, mais supérieure à $-1$ , donc les petites longueurs d'ondes se mettent aussi à osciller mais avec une amplitude décroissante et le schéma ne diverge pas. La solution numérique peut présenter des oscillations, mais qui décroissent au cours du temps.

Pour $r<\frac{1}{4}$ , le facteur d'amplification des grands nombres d'ondes $m$ reste toujours positif, et les petites longueurs d'ondes décroissent de façon monotone, comme dans la solution exacte.

Cette étude nous a donc fournit un renseignement supplémentaire par rapport à l'étude de stabilité: si on veut une convergence monotone de la solution numérique vers la solution exacte, il faut choisir $r\le\frac{1}{4}$ .