3.5 Équation des ondes

3.5.1 Problème physique: propagation d'une onde sonore

On considère la propagation d'une onde sonore plane dans tube de longueur $2L$ contenant un milieu au repos de densité $\rho_{0}$ et de pression $p_{0}$ . Le milieu est perturbé à l'instant initial par une fluctuation de pression $w(x)$ ne dépendant que de la direction spatiale $x$ et sans vitesse initiale (i.e. $\frac{\partial p}{\partial t}(x,0)=0$ ). La densité $\rho(x,t)+\rho_{0}$ , la pression $p(x,t)+p_{0}$ et la vitesse $u(x,t)$ du milieu perturbée sont solutions des équations de conservation d'Euler, qui décrivent la dynamique d'un gaz non visqueux. En supposant que les perturbations sont faibles (hypothèse de l'acoustique), l'équation de conservation de la masse s'écrit au premier ordre:

$$\frac{\partial\rho}{\partial t}+\rho_{0}\frac{\partial u}{\partial x}=0$$

De même l'équation de conservation de la quantité de mouvement s'écrit au premier ordre:

$$\rho_{0}\frac{\partial u}{\partial t}+\frac{\partial p}{\partial x}=0$$

Les fluctuations étant faibles, on peut supposer l'écoulement isentropique, i.e.:

$$\frac{dp}{p_{0}}-\gamma\frac{d\rho}{\rho_{0}}=0 \Rightarrow \frac{p}{p_{0}}=\gamma\frac{\rho}{\rho_{0}}$$

De ces relations on en déduit un système d'équations hyperbolique sur $p$ et $u$ :

$$\frac{\rho_{0}}{\gamma p_{0}}\frac{\partial p}{\partial t}+\rho_{0}\frac{\partial u}{\partial x} =$$ 0

$$\rho_{0}\frac{\partial u}{\partial t}+\frac{\partial p}{\partial x} =$$ 0

En dérivant la première équation par rapport à t et la seconde par rapport à x, on obtiens par différence l'équation des ondes pour la fluctutation de pression $p$ :

$$\frac{\partial^{2}p}{\partial t^{2}}+\frac{\gamma p_{0}}{\rho_{0}}\frac{\partial^{2}p}{\partial x^{2}}=0$$

Si les deux extrémités du tube sont ouvertes, on peut considérer que la pression y est constante et égale à $p_{0}$ . La fluctuation de pression vérifie alors une condition de Dirichlet homogène $p=0$ en $x=-L$ et $x=L$ . Dans le cas d'extrémités fermés (parois solides), la fluctuation de vitesse y est nulle, et la fluctuation de pression vérifie une condition de Neumann $\frac{\partial p}{\partial x}=0$ en $x=-L$ et $x=L$ (car $\frac{\partial p}{\partial x}=-\rho_{0}\frac{\partial u}{\partial t}=0$ si $u=0$ ).

Pour un tube ouvert, la fluctuation de pression $p(x,t)$ est solution de l'équation des ondes, avec les conditions aux limites et initiales suivantes:

$$\frac{\partial^{2}p}{\partial t^{2}}-\frac{\gamma p_{0}}{\rho_{0}}\frac{\partial^{2}p}{\partial x^{2}}=0\textrm{}$$

$$p(-L,t)=0, p(L,t)=0$$

$$p(x,0)=w(x), \frac{\partial p}{\partial t}(x,0)=0$$

Cette équation décrit la propagation d'ondes de pression avec une célérité $c_{0}=\sqrt{\frac{\gamma p_{0}}{\rho_{0}}}$ . Le problème modèle associé s'écrit:

$$(P_{1})\left\{ \begin{array}{ll} \frac{\partial^{2}u}{\partial t^... ...{\partial u}{\partial t}(x,0)=0 & \mbox{conditions initiales}\end{array}\right.$$ (3.45)

3.5.2 Étude de la solution analytique

En effectuant le changement de variables $X=x-c_{0}t$ et $Y=x+c_{0}t$ , l'équation (3.45) s'écrit:

$$\frac{\partial^{2}u(X,Y)}{\partial X\partial Y}=0$$

dont la solution est $u(X,Y)=F(X)+G(Y)$ . La solution générale $u(x,t)$ de l'équation des ondes (3.45) s'écrit:

$$u(x,t)=F(x-c_{0}t)+G(x+c_{0}t)$$

qui est la somme d'une onde progressive $F$ qui se propage avec la célérité $c_{0}$ et d'une onde régressive $G$ qui se propage avec la célérité $-c_{0}$ . Les fonctions $F$ et $G$ sont déterminées par les conditions initiales et aux limites du problème. Nous étudierons plus particulièrement les deux cas suivants:

3.5.2.1 cas 1: propagation d'un créneau

On suppose que la perturbation à l'instant initiale est un créneau de largeur $2\delta$ centré en $x=0$ , qui peut être définie mathématiquement sous la forme:

$$w(x)=Heaviside(x+\delta)*Heaviside(\delta-x)$$

Les fonctions $F$ et $G$ doivent donc vérifier à l'instant initial les deux équations:

$$F(x)+G(x)=w(x)$$  et $$-c_{0}F'(x)+c_{0}G'(x)=0$$

dont les solutions, à une constante additive près, s'écrivent:

$$F(x)=\frac{1}{2}w(x)+K, G(x)=\frac{1}{2}w(x)-K$$

La solution $u_{1}(x,t)$ :

$$u_{1}(x,t)=\frac{1}{2}(w(x-c_{0}t)+w(x+c_{0}t))$$ (3.46)

décrit la propagation dans des directions opposées de deux créneaux de hauteur égale à la moitié de la hauteur initiale. Cette solution vérifie les conditions aux limites tant que les créneaux n'ont pas atteint les frontières, i.e. pour $t<\frac{L}{c_{0}}$ . Lorsque le créneau sort du domaine à $t_{1}=\frac{L}{c_{0}}$ , la solution (3.46) ne vérifie plus les conditions aux limites puisque $u_{1}(-L,t_{1})=\frac{1}{2}$ et $u_{1}(+L,t_{1})=\frac{1}{2}$ . Pour que ces conditions soient vérifiées, il faut qu'une autre onde $u_{2}(x,t)$ de signe opposé rentre dans le domaine en $x=-L$ et $x=L$ , de telle sorte que $u_{1}(-L,t)+u_{2}(-L,t)=0$ et $u_{1}(L,t)+u_{2}(L,t)=0$ . Cela correspond à la réflexion de l'onde initiale sur la frontière. Cette onde entrante s'écrit:

$$u_{2}(x,t)=-\frac{1}{2}(w(x+2L-c_{0}t)+w(x-2L+c_{0}t))$$

Elle correspond à la propagation de deux ondes d'amplitude $-\frac{1}{2}$ qui traversent le domaine et ressortent à $t_{2}=3t_{1}=\frac{3L}{c_{0}}$ . A nouveau pour que la solution vérifie les conditions aux limites, il faut qu'une onde $u_{3}(x,t)$ de signe opposé rentre dans le domaine. Cette onde s'écrit:

$$u_{3}(x,t)=+\frac{1}{2}(w(x+4L-c_{0}t)+w(x-4L+c_{0}t))$$

Elle correspond à la propagation de 2 ondes d'amplitude $\frac{1}{2}$ , qui traversent le domaine et qui à $t_{3}=4t_{1}=\frac{4L}{c_{0}}$ se rencontrent pour reformer l'onde initiale $w(x)$ . Le processus est ensuite périodique.

La solution $u(x,t)$ est donc une solution périodique, de période $4t_{1}=\frac{4L}{c_{0}}$ , qui correspond à la somme de ces 3 ondes pour $0<t<\frac{4L}{c_{0}}$ :

$$u(x,t) = \frac{1}{2}(w(x-c_{0}t)+w(x+c_{0}t)-w(x+2L-c_{0}t)-w(x-2L+c_{0}t)$$

$$+ w(x+4L-c_{0}t)+w(x-4L+c_{0}t)) 0<t<\frac{4L}{c_{0}}$$ (3.47)

On a tracé sur la figure (3.23) l'évolution temporelle de cette solution sur une période. On observe bien la réflexion des ondes sur les parois du domaine.

Figure 3.23: Solution exacte de l'équation des ondes (cas 1)

[][][]

Le programme Maple (3.5.2) permet de tracer l'animation temporelle de cette solution en utilisant la fonction Maple animate.

programme Maple 3.5.2: Solutions exactes de l'équation des ondes 3.45 (cas 1 et 2)

> restart: with(plots):
# Solution exacte equation des ondes: cas 1
> delta:=1/10;
> w:=x->Heaviside(x+delta)*Heaviside(delta-x);
> u(t,x)=1/2*(w(x-c0*t)+w(x+c0*t))-1/2*(w(x+2*L-c0*t)+
  w(x-2*L+c0*t))+1/2*(w(x+4*L-c0*t)+w(x-4*L+c0*t));sol:=rhs(%):
> animate(subs(c0=1,L=1,sol),x=-1..1,t=0..4,numpoints=200,frames=100);
# Solution exacte equation des ondes: cas 2
> m:=3;Tf:=4*L/c0/(2*m+1);
> u(t,x)=cos((2*m+1)*Pi/2*x/L)*cos((2*m+1)*Pi/2*c0/L*t);sol:=rhs(%):
> animate(subs(c0=1,L=1,sol),x=-1..1,t=0..4/7,numpoints=200,frames=100);

3.5.2.2 cas 2: onde stationnaire

Dans le cas d'une solution initiale harmonique, i.e. de période multiple de $2L$ :

$$w(x)=\cos\left(\frac{(2k+1)\pi}{2}\frac{x}{L}\right)$$

la solution s'écrit d'après (3.46):

$$u(x,t)=\frac{1}{2}\left(\cos\left(\frac{(2k+1... ...}\right)+\cos\left(\frac{(2k+1)\pi}{2} \frac{x+c_{0}t}{L}\right)\right)$$

qui après transformation trigonométrique donne:

$$u(x,t)=\cos\left(\frac{(2k+1)\pi}{2}\frac{x}{L}\right)\cos\left(\frac{(2k+1)\pi}{2}\frac{c_{0}t}{L}\right)$$ (3.48)

C'est une onde stationnaire de période $T=\frac{4L}{(2k+1)c_{0}}$ , dont on a tracé la représentation sur la figure (3.24) pour $m=3$ . Le programme Maple (3.5.2) permet de tracer l'animation temporelle de cette solution.

Figure 3.24: Solution exacte de l'équation des ondes (cas 2)

3.5.2.3 cas général

On peut montrer, en utilisant la méthode de séparation de variables, que la solution générale de l'équation des ondes (3.45) est une combinaison linéaire d'ondes stationnaires:

$$u(x,t)=\sum_{k=0}^{\infty}U_{k}\cos\left(\frac{(2k+1)\pi}{2}\frac{x}{L}\right)\cos\left(\frac{(2k+1)\pi}{2}\frac{c_{0}t}{L}\right)$$ (3.49)

Les coefficients $U_{k}$ sont alors les coefficients de Fourier de la solution initiale $w(x)$ . Ainsi la solution du créneau (3.47) peut aussi s'écrire sous la forme d'une combinaison d'ondes stationnaires (3.49), mais dans ce cas le nombre d'ondes stationnaires à considérer est très grands (infinie en théorie), car la série de Fourier associée converge très lentement.

3.5.3 Discrétisation avec un schéma explicite centré

La discrétisation de l'équation des ondes (3.45) par un schéma de différences finies explicite centré s'écrit:

$$\frac{u_{i}^{n+1}-2u_{i}^{n}+u_{i}^{n-1}}{dt^{2}}=c_{0}^{2}\left(\frac{u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}}{dx^{2}}\right)$$ (3.50)

C'est un schéma explicite qui donne la valeur inconnue $u_{i}^{n+1}$ à l'étape $t^{n+1}$ en fonction des valeurs connues $\{u_{i}^{n},u_{i+1}^{n},u_{i-1}^{n},u_{i}^{n-1}\}$ à l'étape $t^{n}$ et $t^{n-1}$ comme indiqué sur le diagramme (3.25).

Figure 3.25: schéma explicite centrée pour l'équation des ondes

3.5.3.1 Étude de la stabilité

L'étude de la stabilité se fait classiquement avec la méthode de perturbation de Neumann. L'évolution temporelle d'un mode de Fourier d'une perturbation $\varepsilon_{i}^{n}=\psi^{n}e^{I\omega idx}$ vérifie l'équation aux différences, soit après simplification par $e^{I\omega idx}$ :

$$\psi^{n+1}-2\psi^{n}+\psi^{n-1}=2\psi^{n}\frac{dt^{2}c_{0}^{2}(\cos(\omega dx)-1)}{dx^{2}}$$

En notant $G=\frac{\psi^{n+1}}{\Psi^{n}}$ le facteur d'amplification du schéma, la relation de récurrence à 2 niveaux précédente conduit à l'équation du second degré suivante:

$$G^{2}-2(1-2C_{FL}^{2}\sin(\frac{\omega dx}{2})^{2})G+1=0$$ (3.51)

dans laquelle on a introduit le nombre de Courant $C_{FL}$ :

$$C_{FL}=\frac{c_{0}dt}{dx}$$ (3.52)

Le mode de Fourier $\varepsilon_{i}^{n}$ est donc une combinaison linéaire des racines $G_{1}$ et $G_{2}$ de cette équation (3.51):

$$\varepsilon_{i}^{n}=(\alpha_{1}(G_{1})^{n}+\alpha_{2}(G_{2})^{n})e^{I\omega idx}$$ (3.53)

Si cette perturbation ne croît pas au cours du temps, le schéma est stable. Il faut donc:

$$\left\vert G_{1}\right\vert\le1$$   et  $$\left\vert G_{2}\right\vert\le1 \forall\omega$$

D'après l'équation (3.51), le produit des racines $G_{1}G_{2}$ est égal à $1$ . Donc si les racines sont réelles, l'une des racines est en module supérieure à $1$ et le schéma est instable. Par contre si les racines sont complexes conjuguées, leur module est égale à 1 et le schéma est stable. Pour que l'équation (3.51) ait des racines complexes, il faut que son discriminant soit négatif, i.e. que:

$$(1-2C_{FL}^{2}\sin(\frac{\omega dx}{2})^{2})^{2}<1$$

soit:

$$2C_{FL}^{2}\sin(\frac{\omega dx}{2})^{2}<2 \forall\omega$$

ce qui conduit à la condition classique de Courant:

$$C_{FL}<1$$ (3.54)

Le schèma explicite est donc conditionnellement stable avec une condition de stabilité donnée par la condition de Courant (3.54).

3.5.3.2 Étude de la consistance

Pour étudier la consistance de ce schéma, nous utiliserons le programme Maple (3.5.3).

programme Maple 3.5.3: Etude de la consistance du schèma explicite 3.50

> restart;with(plots):
# Etude de l'equation des ondes
> diff(U(t,x),t$2)-c0^2*diff(U(t,x),x$2)=0;eq:=%:
# Schema discret D.F.
> (U[n+1,i]-2*U[n,i]+U[n-1,i])/dt^2=
  c0^2*(U[n,i+1]-2*U[n,i]+U[n,i-1])/dx^2; eqh:=%:
# Consistance
> Uex:=(p,q)->U(t+(p-n)*dt,x+(q-i)*dx);
> subs(U[n,i]=Uex(n,i),U[n-1,i]=Uex(n-1,i),
      U[n,i+1]=Uex(n,i+1),U[n,i-1]=Uex(n,i-1),
      U[n+1,i]=Uex(n+1,i),lhs(eqh)-rhs(eqh)); rel3:=%:
> expand(simplify(rel3-lhs(eq)));rel4:=%:
# Develt en serie de Taylor
> U(t+dt,x)=convert(mtaylor(U(t+dt,x),[dt],6),diff);S1:=%:
> U(t-dt,x)=convert(mtaylor(U(t-dt,x),[dt],6),diff):S2:=%:
> U(t,x+dx)=convert(mtaylor(U(t,x+dx),[dx],6),diff):S3:=%:
> U(t,x-dx)=convert(mtaylor(U(t,x-dx),[dx],6),diff):S4:=%:
# Erreur de tronacture
> simplify(subs(S1,S2,S3,S4,rel4));
> ErrT:=collect(%,dt);
# transformation a l'aide de l'equation
> eq;diff(U(t,x),t$4)=c0^4*diff(U(t,x),x$4);
> simplify(subs(%,ErrT));

Le programme fournit l'erreur de troncature suivante:

$$ErrT=\left(-\frac{1}{12}c_{0}^{2} \frac{\par... ...1}{12}\frac{\partial^{4}u}{\partial t^{4}}dt^{2}\right)+O(dx^{4},dt^{4})$$

que l'on peut transformer en dérivant 2 fois l'équation (3.45), et en introduisant le nombre de Courant (3.52):

$$ErrT=\frac{c_{0}^{2}dx^{2}}{12}\left(C_{FL}^{2}-1\right)\frac{\partial^{4}u}{\partial x^{4}}+O(dx^{4},dt^{4})$$

Le schéma explicite est donc consistant à l'équation des ondes. Il est d'ordre 2 en espace et en temps.

On peut montrer que pour la valeur particulière $C_{FL}=1$ , l'erreur de troncature est exactement nulle, et la solution exacte est alors solution de l'équation discrète. Pour $C_{FL}=1$ (limite de stabilité), le schéma explicite fournit donc la solution exacte de l'équation des ondes.

D'après ce calcul d'erreur de troncature, l'équation aux différences finies (3.50) est donc équivalente à l'équation suivante:

$$\frac{\partial^{2}u}{\partial t^{2}}-c_{0}^{2... ...}u}{\partial x^{2}}+c_{0}^{2}\beta\frac{\partial^{4}u}{\partial x^{4}}=0$$

avec un coefficient $\beta$ égal à :

$$\beta=\frac{dx^{2}}{12}\left(C_{FL}^{2}-1\right)$$

C'est une équation d'ondes avec dispersion dont la solution élémentaire s'écrit:

$$u(x,t)=e^{I\omega x}e^{I\omega c_{0}\sqrt{1+\beta\omega^{2}}t}$$

Cette solution correspond à la propagation d'une onde avec une célérité $c_{0}\sqrt{1+\beta\omega^{2}}$ fonction de la pulsation $\omega$ . Dans notre cas $\beta$ est négatif (car $C_{FL}<1$ ) et les ondes hautes fréquences sont donc ralenties par le schéma numérique.

3.5.3.3 Propriétés de dispersion

Pour étudier les propriétés de la solution approchée, nous allons utiliser la même démarche que dans le paragraphe c3disp. Pour cela nous allons considérer le problème de la convection d'une onde $u_{0}(x)=e^{I\omega x}$ , vérifiant les conditions aux limites (i.e. avec $\omega=\frac{(2k+1)\pi}{2L}$ ). La solution exacte de l'équation des ondes (3.45) s'écrit:

$$u(x,t)=\frac{1}{2}e^{I \omega x}\left(e^{-I \omega c_{0} t}+e^{+I \omega c_{0} t}\right)$$

Cette solution initiale correspond justement à la perturbation $\epsilon_{i}^{n}$ dans la méthode de stabilité de Neumann. La solution numérique $u_{i}^{n}$ peut dans ce cas se calculer à partir de la relation (3.53). Elle s'écrit en fonction de la solution initiale $u_{0}$ au noeud $i$ du maillage :

$$u_{i}^{n}=\frac{1}{2}e^{I\omega idx}\left((G_{1})^{n}+(G_{2})^{n}\right)$$

$G_{1}$ et $G_{2}$ sont les 2 racines complexes de l'équation 3.51:

$$G_{1}=1-2y^{2}+2Iy\sqrt{1-y^{2}}, G_{2}=1-2y^{2}-2Iy\sqrt{1-y^{2}}$$

en notant $y=C_{FL}\sin(\frac{\omega dx}{2})$ . Ce sont deux nombres complexes conjugués de module unité et de phase $\psi$ :

$$G_{1}=e^{I\psi}, G_{2}=e^{-I\psi}$$     avec   $$\psi=\arctan\frac{2y\sqrt{1-y^{2}}}{1-2y^{2}}$$

La solution discrète $u_{i}^{n}$ s'écrit sous la forme:

$$u_{i}^{n}=\frac{1}{2}e^{I\omega idx}\left(e^{-In\psi}+e^{+In\psi}\right)$$

qu'il faut comparer à la solution exacte:

$$u(idx,ndt)=\frac{1}{2}e^{I\omega idx}\left(e^{-I\omega nc_{0}dt}+e^{+I\omega nc_{0}dt}\right)$$

Ces deux expressions décrivent la propagation de l'onde initiale: la première avec la célérité $\frac{\psi}{\omega dt}$ et la seconde avec la célérité $c_{0}$ .

En effectuant un développement limité de $\psi$ par rapport à $dx$ , on montre que:

$$\psi=C_{FL}\omega dx+O(dx^{3})\approx c_{0}\omega dt$$

La solution approchée converge vers la solution exacte à l'ordre 2. L'erreur entre la solution exacte et la solution approchée provient essentiellement d'une erreur de dispersion: i.e. la solution numérique se propage sans dissipation, mais avec une célérité légèrement différente de la célérité exacte $c_{0}$ . Nous avons tracé sur la figure (3.26) la célérité (multipliée par $\omega dt$ ) de l'onde pour la solution exacte et la solution approchée en fonction de la pulsation $\omega dx=\frac{(2k+1)\pi}{2N}$ . On note que pour les petits nombres d'ondes (i.e. les basses fréquences) l'écart est minime, et que cet écart croit avec le nombre d'onde (i.e avec la fréquence).

La dispersion numérique du schéma augmente donc avec la fréquence. Les ondes hautes fréquences sont ralenties par le schéma numérique puisque leur célérité est plus petite que $c_{0}$ , comme prédit par l'analyse de consistance.

Figure 3.26: Célérité des ondes $c$ multipliée par $\omega dt$ en fonction de la pulsation $\omega dx$ pour $C_{FL}=\frac {1}{2}$

3.5.4 Expérimentation numérique avec Matlab

La résolution numérique du schéma explicite (3.50) nécessite l'initialisation de la solution $u^{0}$ à $t=0$ et $u^{1}$ à $t=dt$ . Disposant de deux conditions initiales (3.45), la valeur $u^{0}$ est donnée par la solution initiale $w(x)$ :

$$u_{i}^{0}=w(idx)$$ (3.55)

La valeur de $u^{1}$ est obtenue à partir d'un développement limité en temps à l'ordre 2 de la solution au premier pas en temps $u(idx,dt)$ :

$$u_{i}^{1}=u_{i}^{0}+\left(\frac{\partial u}{\... ...+\left(\frac{\partial^{2}u}{\partial t^{2}}\right)_{t=0}\frac{dt^{2}}{2}$$

La valeur de $\frac{\partial u}{\partial t}$ est fournie par la seconde condition initiale, et on utilise l'équation exacte pour calculer $\frac{\partial^{2}u}{\partial t^{2}}$ en fonction de $\frac{\partial^{2}w}{\partial x^{2}}$ , que l'on discrétise ensuite par différences finies centrées:

$$\left(\frac{\partial^{2}u}{\partial t^{2}}\ri... ...^{2}}\approx c_{0}^{2}\left(\frac{w_{i+1}-2w_{i}+w_{i-1}}{dx^{2}}\right)$$

On obtiens ainsi la valeur de $u^{1}$ avec une précision $O(dt^{2},dx^{2})$ , identique à celle du schéma:

$$u_{i}^{1}=w_{i}+\frac{c_{0}^{2}dt^{2}}{2dx^{2}}(w_{i+1}-2w_{i}+w_{i-1})$$ (3.56)

Le programme Matlab (3.5.4) implémente ce calcul pour les 2 conditions initiales étudiées: le créneau (lignes 14 à 15) et l'onde stationnaire (lignes 18 à 19). L'initialisation des champs utilise les relations précédentes (3.55) et (3.56) (lignes 22 à 24). Dans les itérations en temps (lignes 27 à 30), on utilise l'équation aux différences (3.50) pour les noeuds internes $1<i<N$ en utilisant la programmation matricielle Matlab (ligne 27) . Les conditions aux limites fournissent la valeur aux noeuds frontières $i=1$ et $i=N$ (ligne 28).

programme Matlab 3.5.4: Résolution par différences finies de l'équation des ondes 3.45

% resolution equation des ondes
clear;
L=1; c0=2; N=201; dx=2*L/(N-1);
CFL=0.8
dt=CFL*dx/c0; Tf=4*L/c0;
nit=round(Tf/dt);
X=[-L:dx:L];
% noeuds internes
I=[2:N-1]; Ip1=[3:N]; Im1=[1:N-2];
% C.I.
cas=1;
if (cas==1) 
% C.I. cas 1
 delta=0.05;
 W=(X<=delta).*(X>=-delta);
else
% C.I. cas 2
 m=20;
 W=cos((2*m+1)*pi/2*X/L); Ue=W;
end
% initialisation du calcul
Un0=W;
Un(I)=W(I)+CFL^2/2*(W(Im1)-2*W(I)+W(Ip1));
Un(1)=0; Un(N)=0;
% iteration en temps
for it=2:nit
  Un1(I)=2*Un(I)-Un0(I)+CFL^2*(Un(Im1)-2*Un(I)+Un(Ip1));
  Un1(1)=0; Un1(N)=0;
  Un0=Un; Un=Un1;
end;

On a tracé sur la figure (3.27) la solution numérique de la propagation d'un créneau (cas 1) pour les deux valeurs $1.0$ et $0,9$ du nombre de Courant $C_{FL}$ . En comparant ces solutions avec la solution exacte sur la figure (3.23), on constate que pour $C_{FL}=1$ on obtiens la solution exacte. Pour $C_{FL}<1$ , des oscillations hautes fréquences apparaissent à l'arrière du créneau pour $t<t_{1}$ (i.e. avant réflexion sur la frontière). Cela confirme la dispersion numérique des hautes fréquences par le schéma pour $C_{FL}<1$ . En effet la solution initiale exacte possède un nombre infinie de modes de Fourier. La solution initiale discrète contient donc des grands nombres d'onde, qui se propagent numériquement avec une vitesse plus faible que $c_{0}$ . Cela explique l'apparition d'oscillations hautes fréquences à l'arrière du créneau avant la réflexion (figure 3.27a). Après la réflexion sur les frontières (figure 3.27b), ces oscillations ont été réfléchies et ont polluées tout le domaine.

Figure 3.27: propagation d'un créneau avec le schéma explicite (N=401)

[$t=3/4 t_1$ et $N=401$ ]

[$t=5/4 t_1$ et $N=401$ ]

Sur la figure (3.28), on a tracé la solution numérique avec $C_{FL}=0,9$ pour une onde stationnaire de petit nombre d'onde ( $m=3$ ) comparée à la solution exacte à deux instants sur une période. On constate que la solution numérique est très proche de la solution exacte (les courbes coïncident), ce qui montre que la dispersion numérique est très faible pour les ondes à base fréquence.

Figure 3.28: onde stationnaire ($m=3$ ) avec le schéma explicite ( $C_{FL}=0,9 N=101)$

Figure 3.29: onde stationnaire ($m=10$ ) avec le schéma explicite ( $C_{FL}=0,9 N=101)$

Par contre, pour une onde stationnaire avec un plus grand nombre d'onde ( $m=10$ ), on a comparé sur la figure (3.29) la solution numérique et la solution exacte à un instant $t$ , calculée avec les mêmes paramêtres que précédemment. On note sur cette figure le déphasage introduit par la dispersion de la solution numérique, caractérisé pour une onde stationnaire par une amplitude déphasée.