1.6 Validation avec Matlab

Pour valider les analyses précédentes, nous allons utiliser Matlab pour simuler les problèmes $(P_{1})$ et $(P_{2})$ .

1.6.1 Simulation du problème $(P_{1})$

La fonction Matlab expliciteP1, donnée dans le programme 1.6.1, calcule l'erreur $max(\vert u_{e}-u_{i}^{n}\vert)$ au temps $t=t_{0}$ entre la solution exacte $u_{e}(x,t)$ et la solution numérique $u_{i}^{n}$ calculée avec $N$ intervalles, un paramètre $r$ fixé, et pour une condition initiale modale:

$$C(x)=\sin\frac{m\pi}{L}x$$

La programmation Matlab est très simple. Avec Matlab, les vecteurs ont un indice qui commence à 1. La solution est donc stockée dans un vecteur $U=[U(1),U(2),...U(N+1)]$ de dimension $N+1$ . On effectue une boucle en temps (ligne 16 à 20), pour calculer la solution $u_{i}^{n}$ . En utilisant les notations matricielles de Matlab, la solution aux noeuds internes $i=2,N-1$ s'écrit (ligne 18):

$$U(2:N-1)=(1-2*r)*U(2:N-1)+r*(U(1:N-2)+U(3:N+1))$$

Dans cette expression, à droite du signe =, $U$ représente le vecteur solution connu à l'instant $t-dt$ . La nouvelle valeur de $U$ ( à gauche du signe =) représente la nouvelle solution calculée à l'instant $t$ .

programme Matlab 1.6.1: Calcul de la solution numérique explicite du problème ( $P_1$ )

function [ErrMax,U]=expliciteP1(N,r,t0,m)
% schema explicite pble P1
% calcul la solution U et l'erreur ErrMax
% a t=t0 sur N points
% pour une valeur du parametre r fixe
% pour une condition initiale sin(m*Pi*X)
% parametres
L=1;
h=1.0/N; dt=r*h^2; Ndt=round(t0/dt)
% liste des nds internes et pts du maillage
Ni=[2:N];
X=[0:h:L];
% condition initiale C
C=sin(m*pi*X);
t=0; U=C;
for it=1:Ndt
        t=t+dt;
        U(Ni)=(1-2*r)*U(Ni)+r*(U(Ni-1)+U(Ni+1));
        U(1)=0; U(N+1)=0;
end;
% calcul de l'erreur
Ue=C*exp(-m^2*pi^2*t);
ErrMax=norm(Ue-U,inf);

englishANIMATION:Calcul de la solution numérique explicite du problème ( $P_1$ )

On a tout d'abord tracé la solution numérique calculée à un instant $t=0.014$ pour $m=4$ , $\kappa=1$ , $L=1$ et deux valeurs de $r$ (figure 1.6).

Pour $r=\frac{1}{2}$ (figure 1.6a), la solution numérique décroit comme la solution exacte, i.e. l'amplitude de la solution numérique décroit de $0,327$ , alors que la solution exacte décroit de $e^{-\kappa\frac{m^{2}\pi^{2}}{L^{2}}t}\approx0,331$ .

Par contre pour $r=1$ (figure 1.6b), la solution numérique diverge avec une apparition d'oscillations à hautes fréquences d'amplitudes croissantes, comme prévue dans la thèorie.

Figure 1.6: Solution numérique du problème $P_{1}$ (mode $m=4$ )

[ $r=\frac {1}{2}$ et $N=50$ ]

[$r=1$ et $N=50$ ]

Pour les mêmes paramêtres , on a ensuite tracé sur la figure (1.7) l'évolution temporelle de la valeur maximale de la solution numérique calculée pour différentes valeurs du pas en temps $dt$ et pour deux maillages $N=50$ et $N=100$ . On les a comparé à la valeur de exacte déduite de (1.18):

$$U_{max}=e^{-\kappa\frac{m^{2}\pi^{2}}{L^{2}}t}$$

Figure 1.7: Evolution temporelle de $U_{max}$ entre $t=0$ et $t=\frac {4}{m^{2}\pi ^{2}}$

[N=50]

[N=100]

Pour $N=50$ , on constate que la solution numérique et la solution exacte sont très proches tant que $dt\leq2 10^{-4}$ , et que la solution numérique diverge pour $dt>2 10^{-4}$ . De même pour $N=100$ , la solution numérique converge tant que $dt\leq0,5 10^{-4}$ et diverge sinon. Ces valeurs du pas en temps correspondent justement à la valeur limite de stabilité du paramêtre $r$ (1.36):

$$\kappa\frac{\Delta t}{\Delta x^{2}}=\frac{1}{2}$$

De façon à analyser la précision de la solution numérique, on a tracé sur la figure (1.8) pour $m=4$ et $t_{0}=\frac{1}{m^{2}\pi^{2}}$ , l'évolution de l'erreur $max(\vert u_{e}-u_{i}^{n}\vert)$ pour des valeurs de $N$ de 8 à 256, et différentes valeurs de $r$ .

Figure 1.8: Erreur du schéma explicite en fonction de $N$ (pble $P_{1}$ )

On vérifie sur cette figure (1.8) la décroissance quadratique de l'erreur, à $r$ fixé, en fonction de $N$ , puisque l'erreur théorique est en $O(\Delta t,\Delta x^{2})$ . A $r$ fixé, lorsque l'on divise $\Delta x$ par 2, on divise $\Delta t$ par 4 et l'erreur est divisée d'un facteur 4. On note aussi que pour $r=\frac{1}{6}$ , la convergence est plus rapide car la valeur du pas en temps $\Delta t$ correspond au cas particulier c1eq28b pour lequel le premier terme de l'erreur de troncature est nul. L'erreur de troncature est dans ce cas en $O(\Delta t^{2},\Delta x^{4})$ et l'erreur décroit en $O(N^{-4})$ .

On constate enfin que pour $r=1$ , le schéma diverge, mais uniquement pour les grandes valeurs de $N$ . En effet pour $N$ grand, on effectue beaucoup plus d'itérations que pour $N$ petit et les erreurs d'arrondis sont suffisament amplifiées pour faire diverger la solution.

On note enfin que l'erreur diminue à $N$ fixé avec $r$ , c'est à dire avec $\Delta t$ (avec un cas particulier pour $r=\frac{1}{6}$ ).

1.6.2 Simulation du problème $(P_{2})$

La discrétisation explicite du problème $(P_{2})$ conduit à la même équation aux différences (1.35), que pour le problème $(P_{1})$ . La seule différence provient de la condition aux limites en $x=0$ . Pour le problème $(P_{2})$ , la condition aux limites est une condition de symétrie: $\frac{\partial u}{\partial x}=0$ , que l'on doit discrétiser pour obtenir la valeur numérique $u_{0}^{n+1}$ au premier point du maillage. Nous étudierons deux façons d'imposer cette condition aux limites de type Neumann:

1. discrétisation décentrée de la condition aux limites:
   on approxime la condition aux limites par un schéma aux différences décentré faisant intervenir la valeur $u_{0}$ :
   

   $$\left.\frac{\partial u}{\partial x}\right\vert _{0}^{n+1}\approx\frac{u_{1}^{n+1}-u_{0}^{n+1}}{\Delta x}$$ (1.53)

   
   
   d'où l'on déduit l'équation au noeud $i=0$
   

   $$u_{0}^{n+1}=u_{1}^{n+1}$$ (1.54)

   
   

2. condition miroir:
   on utilise l'équation aux différences (1.35) au noeud $i=0$ . Dans cette équation apparaît la valeur inconnue $u_{-1}^{n}$ , que l'on calcule avec une condition miroir, traduction de la condition de symétrie:
   

   $$\left.\frac{\partial u}{\partial x}\right\vert _{0}^{n}\approx\frac{u_{1}^{n}-u_{-1}^{n}}{2\Delta x}$$ (1.55)

   
   
   ce qui conduit à $u_{-1}^{n}=u_{1}^{n}$ et à l'équation suivante au noeud $i=0$
   

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

   
   

La fonction Matlab expliciteP2, donnée dans le programme 1.6.2, calcule l'erreur $max(\vert u_{e}-u_{i}^{n}\vert)$ au temps $t=t_{0}$ entre la solution exacte $u_{e}(x,t)$ et la solution numérique $u_{i}^{n}$ calculée avec $N$ intervalles, un paramètre $r$ fixé, et pour une condition initiale modale si $m\ge0$ :

$$C(x)=\cos\frac{(2m+1)\pi}{2L}x$$

ou une condition initiale d'Heaviside si $m<0$ (dont on a calculé une solution exacte au paragraphe c1par1):

$$C(x)=Heaviside(\delta-x)$$

programme Matlab 1.6.2: Calcul de la solution numérique explicite du problème ( $P_2$ )

function [ErrMax,U]=expliciteP2(N,r,t0,m)
% schema explicite pble P2
% calcul la solution U et l'erreur ErrMax
% a t=t0 sur N points
% pour une valeur du parametre r fixe
% pour une condition initiale cos((m+1)*Pi/2*X) si m>=0
% ou fonction porte si m<0
% parametres
L=1;
h=1.0/N; dt=r*h^2; Ndt=round(t0/dt)
% liste des nds internes et pts du maillage
Ni=[2:N];
X=[0:h:L];
% condition initiale C
if (m>=0)
   C=cos((2*m+1)*pi/2*X);
else
   % coefficient de Fourier
   delta=1/10;
   M=[0:20];
   CM=4*(sin(pi*M*delta)*cos(pi*delta/2)+...
         cos(pi*M*delta)*sin(pi*delta/2))./((2*M+1)*pi);
   C=CM*cos(pi/2*(2*M+1)'*X);
end;
t=0; U=C;
for it=1:Ndt
        t=t+dt;
	Un=U;
        U(Ni)=(1-2*r)*Un(Ni)+r*(Un(Ni-1)+Un(Ni+1));
 	% conditions aux limites
        U(1)=(1-2*r)*Un(1)+r*2*Un(2);
	U(N+1)=0;
end;
% calcul de l'erreur
if (m>=0)
   Ue=C*exp(-((2*m+1)*pi/2)^2*t);
else
   Ue=(CM.*exp(-(2*M+1).^2*pi^2/4*t))*cos(pi/2*(2*M+1)'*X);
end;
ErrMax=norm(Ue-U,inf);

La programmation Matlab est identique à celle du programme (1.6.1). On notera simplement l'utilisation de la notation matricielle pour le calcul de la solution analytique Ue (équation 1.20) aux noeuds du maillage (ligne 21 et 38) sur $M$ modes à l'instant $t=n\Delta t$ :

$$U_{e}(i)=\sum_{k=0}^{M}C_{k}e^{-\frac{1}{4}\l... ...^{2}\pi^{2}n\Delta t}\cos(\frac{\left(2 k+1\right)}{2}\pi i\Delta x)$$

où les $C_{k}$ sont données par l'équation 1.21:

$$C_{k}=4 \frac{\sin(1/10 \pi k)\cos(1/20 \pi)+\cos(1/10 \pi k)\sin(1/20 \pi)}{\left(2 k+1\right)\pi}$$

En notant $M=[0:20]$ les modes utilisés pour la calcul de la solution, le tableau $CM$ des coefficients $C_{k}$ est calculé directement à la ligne 21. Le vecteur solution exacte $U_{e}$ aux noeuds du maillage est calculé à la ligne 38 en multipliant terme à terme le vecteur $CM$ par le vecteur $\left[e^{-1/4 \left(2 M+1\right)^{2}\pi^{2}t}\right]$ (opérateur $.*$ ), puis en effectuant le produit scalaire (opérateur $*$ ) avec le vecteur colonne $\cos(\frac{\left(2 M+1\right)}{2}\pi X)$ pour effectuer la somme.

La boucle en temps correspond aux lignes 26 à 33 , et la condition aux limites en $x=0$ est imposée à la ligne 31 (condition miroir 1.55 dans le programme 1.6.2). Pour l'autre condition (1.53), il suffit de remplacer la ligne 31 par:

$$U(1)=U(2);$$

On a tracé sur la figure (1.9) la solution calculée pour $N=160$ et $r=0.1$ à différents instants ( $t_{0}$ est le temps caractéristique de diffusion $t_{0}=\frac{1}{\kappa(\frac{\pi}{2})^{2}}$ ). En comparant avec la solution exacte de la figure (1.2b), on constate que les solutions sont très proches.

Figure 1.9: Solution numérique du problème $(P_{2})$ avec $N=160$ et $r=0.1$

Pour le vérifier, on a tracé sur la figure (1.10) l'évolution de l'erreur entre la solution exacte (calculée avec 20 modes) et la solution numérique à $t=\frac{t_{0}}{10}$ pour des valeurs de $N$ de 20 à 640, différentes valeurs de $r$ et les deux types de prise en compte de la condition aux limites en $x=0$ . Avec la condition miroir (1.55), on constate que l'erreur, à $r$ fixé, décroit de façon quadratique en $O(N^{-2})$ , ce que prévoit bien la théorie du schéma explicite avec une précision en $O(\Delta t,\Delta x^{2})$ . Par contre si on impose la condition (1.53), l'erreur décroit linéairement en $O(N^{-1})$ .

Figure 1.10: Erreur du schéma explicite (pble $P_{2}$ )

[condition miroir 1.55]

[condition décentrée 1.53]

On constate ici un point important sur la précision des schémas numériques.

La précision d'un schéma numérique dépend à la fois de la précision du schéma interne et de la précision de la discrétisation des conditions aux limites.

Avec la condition (1.53), on calcule la valeur $u_{0}^{n+1}$ avec une approximation d'ordre $O(\Delta x)$ (discrétisation décentrée), alors que avec la condition (1.55), la valeur de $u_{-1}^{n}$ est calculée à l'ordre $O(\Delta x^{2})$ (discrétisation centrée), et on conserve une approximation en $O(\Delta x^{2})$ pour $u_{0}^{n+1}$ . Donc avec une condition au limite d'ordre 1 en $O(\Delta x)$ , même avec un schéma interne d'ordre 2 en $O(\Delta x^{2})$ , la solution a une précision d'ordre 1 en $O(\Delta x)$ .

Figure 1.11: Solution numérique du problème $(P_{2})$ avec $N=160$ et $r=0.5$

Pour terminer cette analyse, nous avons calculé la solution numérique avec $r=0.5$ et $N=160$ pour une condition initiale $C(x)=Heaviside(\delta-x)$ . Le résultat est tracé sur la figure (1.11). On constate sur cette figure l'apparition d'oscillations numériques dès les premières itérations en temps, qui sont ensuite amorties au cours du temps. Ces oscillations numériques n'apparaisent que pour des valeurs de $r$ comprises entre $\frac{1}{4}$ et $\frac{1}{2}$ , et uniquement si la condition initiale possède des modes hautes fréquences. C'est ce que prévoit l'analyse théorique du paragraphe c1convge.

englishANIMATION:Calcul de la solution numérique explicite du problème ( $P_1$ )