3ème MARINE

M2 MATH
2012/2013

VTISI & TSR

FU
Pro

s
ue
iq
blè

ol
m

rb
es

pe
hy

NO NI
A CCA
G. F i s e -à-jo
ur
m
ière
Dern 2013
m ars
8
di 1
Lun
 Gloria FACCANONI

IMATH Bâtiment U-318 T 0033 (0)4 94 14 23 81

Université du Sud Toulon-Var
Avenue de l’université B [Link]@[Link]
83957 LA GARDE - FRANCE i [Link]

2
 Table des matières

I. Devoir maison 5
1. Étude d’équations hyperboliques 6

2. Étude d’un système hyperbolique 12

II. Rappels Théoriques 16

3. Équations linéaires 17

4. Équations non-linéaires 18
4.1. Exemple : équation de B URGERS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.2. Exercice : trafic dans un tunnel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

5. Systèmes hyperboliques 30
5.1. Exemple : le système de S AINT-V ENANT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
5.2. Exemple : Euler isotherme en coordonnées Lagrangiennes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
5.3. Exemple : le système de la dynamique des gaz isentropique par relaxation . . . . . . . . . . . . . . . . . . . . . . . . 41
5.4. Exemple : mélange isotherme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

III. Rappels Numériques 53

6. Résolution numérique de l’équation de transport 59
6.1. Exemples d’études de stabilité et de consistance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

7. Résolution numérique d’une équation non-linéaire 65

8. Systèmes hyperboliques 67

3
 Première partie .

Devoir maison
Le DM sera évalué. Vous avez jusqu’au 1 mars 2013 pour rendre un rapport.
Vous pouvez écrire le rapport seul ou en binôme avec un autre membre de votre
groupe.

5
 1. Étude d’équations hyperboliques
Cas linéaire On considère des équations linéaires qui s’écrivent sous la forme
(
∂t (u(t , x)) + a(t , x)∂x (u(t , x)) + b(t , x)u(t , x) = f (t , x), t > 0, x ∈ R,
u(x, 0) = g (x) x ∈ R.

Calculer la solution faible u(t , x) pour t ≥ 0 et x ∈ R des problèmes suivants en utilisant la méthode des caractéristiques :
Dans la correction on utilise les notations du chapitre 3.
(
∂t u + x∂x u = 0 x ∈ R, t > 0
1.
u(0, x) = x 4 x ∈ R,

a(t , x) = x, b(t , x) = f (t , x) = 0, g (x) = x 4 donc pour (t̂ , x̂) donné on a x(t ) = x̂e t −t̂ et donc
³ ´4
u(t̂ , x̂) = g (x(0)) = g (x̂e −t̂ ) = x̂e −t̂

(
(t + 1)∂t u + ∂x u = 0 x ∈ R, t > 0
2.
u(0, x) = x x ∈ R,
³ ´
1 t +1
a(t , x) = t +1 , b(t , x) = f (t , x) = 0, g (x) = x donc pour (t̂ , x̂) donné on a x(t ) = x̂ + ln t̂ +1
et donc

u(t̂ , x̂) = g (x(0)) = g (x̂ − ln(t̂ + 1)) = x̂ − ln(t̂ + 1)

(
∂t u + xt ∂x u = t x ∈ R, t > 0
3. 2
u(0, x) = x x ∈ R,
2 −t̂ 2 )/2
a(t , x) = xt , b(t , x) = 0, f (t , x) = t , g (x) = x 2 donc pour (t̂ , x̂) donné on a x(t ) = x̂e (t et donc

t t t2
Z Z
2 /2 2
u(t̂ , x̂) = g (x(0)) + f (s, x(s))d s = g (x̂e −t̂ )+ sd s = x̂ 2 e −t̂ + .
0 0 2

(
∂t u + ∂x u + u = e −t +2x x ∈ R, t > 0
4.
u(0, x) = 0 x ∈ R,
a(t , x) = 1, b(t , x) = 1, f (t , x) = e −t +2x , g (x) = 0 donc pour (t̂ , x̂) donné on a x(t ) = x̂ + t − t̂ et donc
Z t̂ 1 1
u(t̂ , x̂) = g (x(0))e −t̂ + e −s+2x(s) e −(t̂ −s) d s = g (x̂ − t̂ )e −t̂ + e 2x̂−t̂ = e 2x̂−t̂ .
0 2 2

Cas non-linéaire On considère l’équation non-linéaire

u(2−u)
(
∂t (u) + ∂x (q(u)) = 0, avec q(u) = 2 ,
u(x, 0) = g (x).

Donner l’unique solution faible entropique u(t , x) pour t ≥ 0 et x ∈ R des problèmes suivants en utilisant la méthode des
caractéristiques. Justifier clairement la réponse avec tous les dessins nécessaires.
2
³ ´
Dans la correction on utilise les notations du chapitre 4. On a q(u) = u − u2 , q 0 (u) = 1 − u, (q 0 )−1 x−x 0 x−x 0
t −t 0 = 1 − t −t 0 et
q 00 (u) = −1 (le flux est concave). L’équation de la caractéristique de pied (ξ, 0) est donc x(t ) = ξ + q 0 (g (ξ))t = ξ + (1 − g (ξ))t .

6
Dernière mise à jour : Lundi 18 mars 2013 1. Étude d’équations hyperboliques

(
1 si x < 1,
1. g (x) =
0 si x > 1,
(
ξ, ξ < 1,
L’équation de la caractéristique de pied (ξ, 0) est x(t ) =
ξ + t , ξ > 1.
t

1 x

La donnée initiale du problème de Riemann a un saut décroissant et le flux est concave donc l’unique solution
entropique présente une onde de raréfaction centrée en (1, 0).
t

1 x

La solution est donc 

1,
 x ≤ 1,
u(x, t ) = 1 − x−1
t , 1 < x ≤ 1+ t,

0, x > 1+ t.


u
t =0
t = 0.5
t =1
x
(
0 si x < 1,
2. g (x) =
1 si x > 1,
(
ξ + t , ξ < 1,
L’équation de la caractéristique de pied (ξ, 0) est x(t ) =
ξ, ξ > 1.
t

1 x

La donnée initiale du problème de Riemann a un saut croissant et le flux est concave donc l’unique solution entropique
présente une onde de choc qui part en (1, 0). Pour calculer l’équation x = s(t ) de l’onde de choc on utilise les relations
de R ANKINE -H UGONIOT :
q(u R ) − q(u L ) uR + uL 1
s 0 (t ) = = 1− =
uR − uL 2 2
avec la donnée initiale s(0) = 1. On trouve donc x = s(t ) = 1 + 12 t .
t

1 x

La solution est donc (

0, x < 1 + 12 t ,
u(x, t ) =
1, x > 1 + 12 t .

© G. Faccanoni 7
1. Étude d’équations hyperboliques Dernière mise à jour : Lundi 18 mars 2013

u
t =0
t = 0.5
t =1
x

2 si x < 1,

3. g (x) = 1 si x ∈ ]1; 2[,

0 si x > 2,


ξ − t , ξ < 1,


L’équation de la caractéristique de pied (ξ, 0) est x(t ) = ξ, 1 < ξ < 2,

ξ + t , ξ > 2.


1 2 x

La donnée initiale a deux sauts décroissants, un en x = 1 et l’autre en x = 2. Puisque q est concave, on s’attend à ce
que la solution faible entropique présente deux ondes de raréfaction centrées l’une en (1, 0) et l’autre en (2, 0).
t

1 2 x



 2, si x ≤ 1 − t ,
x−1

1 − t , si 1 − t < x ≤ 1,



u(x, t ) = 1, si 1 < x ≤ 2,

x−2
1 − t , si 2 < x ≤ 2 + t ,






0, si x > 2 + t .

u t =0
t = 1.5
t =3

x

0 si x < 1,

4. g (x) = 1 si x ∈ ]1; 2[,

2 si x > 2,


Elle a deux sauts croissants, un en x = 1 et un en x = 2. Puisque q est concave, on s’attend à ce que la solution faible
entropique présente deux ondes de choc qui partent de (1, 0) et de (2, 0).
ξ + t , ξ < 1,


B En effet, l’équation de la caractéristique de pied (ξ, 0) est x(t ) = ξ, 1 < ξ < 2,

ξ − t , ξ > 2.


8 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 1. Étude d’équations hyperboliques

1 2 x
B Pour calculer les équations des deux ondes de choc on utilise les relations de R ANKINE -H UGONIOT :
B choc de pied (1, 0)
q(u R )−q(u L )
= 1 − uR +u
(
1
s 0 (t ) = u R −u L 2
L
= 2
s(0) = 1
1
d’où x = s(t ) = + 1 ;
2t
B choc de pied (2, 0)
q(u R )−q(u L )
= 1 − uR +u
(
s 0 (t ) = u R −u L 2
L
= − 21
s(0) = 2
d’où x = s(t ) = − 21 t + 2.
t

(x ∗ , t ∗ )

1 2 x
∗ ∗
Ces deux chocs vont ensuite interagir à partir d’un (x , t ) qu’on détermine en résolvant le système linéaire
(
x = 12 t + 1,
x = − 12 t + 2.

On obtient ainsi le pied de la nouvelle onde de choc qui est (x ∗ , t ∗ ) = ( 23 , 1). Donc la solution pour 0 < t < 1 est
1

0, si x < 2 t + 1,

u(x, t ) = 1, si 12 t + 1 < x < − 12 t + 2,

2, si x > − 12 t + 2,


B Pour t ≥ 1 seul les données initiales pour x < 1 et x > 2 sont transportées et ce nouveau choc, qui part de ( 32 , 1),
vérifie encore la relation de R ANKINE -H UGONIOT :
q(u )−q(u )
s 0 (t ) = uR −u L = 1 − uR +u
(
L
2 =0
R L
3
s(1) = 2

d’où x = s(t ) = 32 .
t

1 2 x
L’unique solution faible entropique pour tout t > 0 est donc
1 3

0, si x < min{ 2 t + 1, 2 },

u(x, t ) = 1, si 12 t + 1 < x < − 12 t + 2,

2, si x > max{− 12 t + 2, 23 },


© G. Faccanoni 9
1. Étude d’équations hyperboliques Dernière mise à jour : Lundi 18 mars 2013

u
t =0

t = 0.5

t ≥1

x

2 si x < 1,

5. g (x) = 0 si x ∈ ]1; 2[,

1 si x > 2.


ξ − t , ξ < 1,


L’équation de la caractéristique de pied (ξ, 0) est x(t ) = ξ + t , 1 < ξ < 2,

ξ, ξ > 2.


1 2 x
On s’attend donc à ce que la solution faible entropique présente une onde de raréfaction centrée en (1, 0) et une onde
de choc qui part de (2, 0) et que, après un certain temps, les deux ondes interagissent.
B La raréfaction est comprise entre la droite d’équation x(t ) = 1 − t et la droite d’équation x(t ) = 1 + t et a équation
x −1 x −1
µ ¶
u(x, t ) = (q 0 )−1 = 1− .
t −0 t
B Pour déterminer l’équation de l’onde de choc qui a pied en (2, 0) on utilise les relations de R ANKINE -H UGONIOT :
q(u R )−q(u L )
= 1 − uR +u
(
1
s 0 (t ) = u R −u L 2
L
= 2
s(0) = 2

d’où x = s(t ) = 12 t + 2.
Ces considérations sont valables jusqu’à ce que l’onde de choc et l’onde de raréfaction se touchent, c’est-à-dire
jusqu’à t = 2. Donc la solution pour 0 < t < 2 est


 2, si x ≤ 1 − t ,

1 − x−1 , si 1 − t < x ≤ t + 1,

u(x, t ) = t


 0, si t + 1 < x < 12 t + 2,
si x > 12 t + 2.

1,


1 2 3 x
B Pour t > 2 on continue à avoir une onde de choc mais l’état gauche n’est plus constant car il est donné par la
raréfaction i.e. u L = 1 − s(t t)−1 donc le nouveau choc a pied en (3, 2) et vérifie
 ³ ´2
q(u )−q(u ) 0− s(t t)−1
s (t ) = uR −u L = 1 − uR +u = s(t2t)−1

 0 L
=
R L 2 2(0− s(t t)−1 )

s(2) = 3

10 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 1. Étude d’équations hyperboliques

p
d’où l’équation x = 1 + 2t .
t

1 2 3 x
L’unique solution faible entropique pour t > 0 est donc

2,
 si x ≤ 1 − t ,

 x−1
p
1 − t , si 1 − t < x < min{1 + t , 1 + 2t },



u(x, t ) = 0, si 1 + t ≤ x < 12 t + 2,

1, si x > 12 t + 2 et t < 2,




 p
1, si x > 1 + 2t et t ≥ 2.

u
t =0

t =1

t =2

t =3

t =4

© G. Faccanoni 11
 2. Étude d’un système hyperbolique
On considère le système de AW et R ASCLE pour la modélisation macroscopique du trafic routier :
(
∂t ρ + ∂x (ρu) = 0,
avec x ∈ R, t > 0, (2.1)
∂t (ρ(u + p(ρ))) + ∂x (ρu(u + p(ρ))) = 0,

où ρ > 0 est la densité des véhicules, u leur vitesse moyenne et p joue le rôle d’un facteur d’anticipation (on assume p 0 (ρ) > 0 et
p 00 (ρ) ≥ 0).
Sous forme vectorielle ce système s’écrit
∂t U + ∂x F(U) = 0
ayant défini
ρ ρu
µ ¶ µ ¶
U= , F(U) = .
ρ(u + p(ρ)) ρu(u + p(ρ))

1. Soit W = (ρ, u), alors le système (2.1) se réécrit, pour des solutions régulières, sous la forme quasi-linéaire
∂t W + A(W)∂x W = 0. (2.2)
Montrer que la matrice A(W) s’écrit
ρ
µ ¶
u
.
0 u − ρp 0 (ρ)

On développe les dérivées du système (2.1) pour des solutions régulières :

(
∂t ρ + ∂x (ρu) = 0,
∂t (ρ(u + p(ρ))) + ∂x (ρu(u + p(ρ))) = 0,
⇓
(
∂t (ρ) + u∂x (ρ) + ρ∂x (u) = 0,
ρ∂t (u + p(ρ)) + (u + p(ρ))∂t (ρ) + (ρu)∂x (u + p(ρ)) + (u + p(ρ))ρ∂x (u) + (u + p(ρ))u∂x (ρ) = 0,
⇓
(
∂t (ρ) + u∂x (ρ) + ρ∂x (u) = 0,
ρ∂t (u) + ρp 0 (ρ)∂t (ρ) + (u + p(ρ))∂t (ρ) + (ρu)∂x (u) + (ρu)p 0 (ρ)∂x (ρ) + (u + p(ρ))ρ∂x (u) + (u + p(ρ))u∂x (ρ) = 0,
⇓
(
∂t (ρ) + u∂x (ρ) + ρ∂x (u) = 0,
∂t (u) + (u − ρp 0 (ρ))∂x (u) = 0.

2. Calculer les deux valeurs propres λ1 (W) et λ2 (W) de la matrice A(W). Afin de fixer les notations on ordonne les deux valeurs
propres selon λ1 (W) < λ2 (W). Proposer une base associée de vecteurs propres à droite {r1 (W), r2 (W)}. En déduire que le
système (5.1) est strictement hyperbolique. (Que se passe-t-il quand ρ → 0 ?)
On cherche les deux solutions λ1 (W) et λ2 (W) de l’équation det(A(W) − λ(W)Id) = 0, c’est-à-dire de l’équation
(u − λ)(u − ρp 0 (ρ) − λ) = 0.
On obtient
λ1 = u − ρp 0 (ρ), λ2 = u.
On peut alors prendre comme vecteurs propres à droite associés à ces valeurs propres les deux vecteurs
µ ¶ µ ¶
−1 1
r1 = 0 , r2 = .
p (ρ) 0
Puisque p 0 (ρ) > 0, les valeurs propres sont réelles et distinctes donc le système (2.1) est strictement hyperbolique. Si ρ → 0,
A est une matrice diagonale dont les deux valeurs propres sont égaux à u.

12
Dernière mise à jour : Lundi 18 mars 2013 2. Étude d’un système hyperbolique

3. Vérifier que le champs 1 est vraiment non linéaire tandis que le champ 2 est linéairement dégénéré. Quelle est la nature du
2-champ ? Autrement dit, il s’agit d’un choc, d’une détente ou d’une discontinuité de contact ?
Pour déterminer la nature des champs caractéristiques on calcule (∇λk )T · rk pour k = 1, 2 :

(∇λ1 )T · r1 = 2p 0 (ρ) + ρp 00 (ρ) > 0 =⇒ le 1-champ est VNL

(∇λ2 )T · r2 = 0 =⇒ le 2-champ est LD

4. On note I k l’invariant de Riemann du k-ème champ caractéristique. Montrer qu’un choix possible pour I 1 et I 2 est

I 1 = u + p(ρ), I 2 = λ2 .

On vérifie que (∇I k )T · rk = 0 pour k = 1, 2 :

(∇I 1 )T · r1 = −p 0 (ρ) + p 0 (ρ) = 0, (∇I 2 )T · r2 = (∇λ2 )T · r2 = 0.

5. Considérons le 2-champ. Quelle est la vitesse σ̇2 de cette discontinuité ? Expliciter les états W = (ρ, u) qui peuvent être relié
à un état droit WR par une discontinuité de vitesse σ̇2 .
Le 2-champ est un champ linéairement dégénéré, donc on a une discontinuité de contact. La condition sur l’invariant de
Riemann donne I 2 (W) = I 2 (WR ), c’est-à-dire u = u R et la discontinuité se déplace à la vitesse σ̇2 = u = u R .
6. Pour un état gauche WL = (ρ L , u L ) donné on cherche les états W = (ρ, u) qui peuvent être relié à WL par une 1-onde de
détente. Montrer que u > u L et ρ < ρ L . Calculer u en fonction de ρ L , u L et ρ. Plus précisément, montrer que u peut se
mettre sous la forme
u = u L + r (ρ L , ρ)
en explicitant la fonction r . Étudier la fonction ρ 7→ u = u L + r (ρ L , ρ) et tracer son graphe dans le plan (ρ, u).

La condition d’entropie λ1 (WL ) < λ1 (W) s’écrit u L − ρ L p 0 (ρ L ) < u

u −ρp 0 (ρ), donc u > u L −ρ L p 0 (ρ L )+ρp 0 (ρ). Dans une détente les
invariants de Riemann sont conservés, ici I 1 (W) = I 1 (WL ), donc
u + p(ρ) = u L + p(ρ L ) et on obtient la courbe intégrale
WL
uL
u : ρ 7→ u(ρ) ≡ u L + p(ρ L ) − p(ρ)

(fonction strictement décroissante et concave). Par conséquent

ρ L p 0 (ρ L )+p(ρ L ) > ρp 0 (ρ)+p(ρ) ; la fonction f : ρ 7→ ρp 0 (ρ)+p(ρ)
est strictement croissante donc de f (ρ L ) > f (ρ) on déduit ρ L > ρ ρL ρ
et par conséquent u L < u.

7. Pour un état gauche WL = (ρ L , u L ) donné on cherche les états W = (ρ, u) qui peuvent être relié à WL par une 1-onde de choc
entropique de vitesse σ̇1 . Montrer que les relations de R ANKINE -H UGONIOT peuvent s’écrire sous la forme
(
j 1 ≡ ρ(u − σ̇1 ) = ρ L (u L − σ̇1 ),
(u + p(ρ))ρ(u − σ̇1 ) = (u L + p(ρ L ))ρ L (u L − σ̇1 ).

Montrer que u < u L et ρ > ρ L . Calculer u en fonction de ρ L , u L et ρ. Plus précisément, montrer que u peut se mettre sous la
forme
u = u L + d (ρ L , ρ)
en explicitant la fonction d . Étudier la fonction ρ 7→ u = u L + d (ρ L , ρ) et tracer son graphe dans le plan (ρ, u).
Les relations de R ANKINE -H UGONIOT pour un 1-choc de vitesse σ̇1 s’écrivent
ρu−ρ u
σ̇1 = ρ−ρL L
(
L
ρu(u+p(ρ))−ρ L u L (u L +p(ρ L ))
σ̇1 = ρ(u+p(ρ))−ρ L (u L +p(ρ L ))

qu’on peut réécrire comme

( (
j 1 ≡ ρ(u − σ̇1 ) = ρ L (u L − σ̇1 ), j 1 ≡ ρ(u − σ̇1 ) = ρ L (u L − σ̇1 ),
=⇒
(u + p(ρ))ρ(u − σ̇1 ) = (u L + p(ρ L ))ρ L (u L − σ̇1 ). u + p(ρ) = u L + p(ρ L ),

© G. Faccanoni 13
2. Étude d’un système hyperbolique Dernière mise à jour : Lundi 18 mars 2013

et on obtient la courbe intégrale

u : ρ 7→ u(ρ) ≡ u L + p(ρ L ) − p(ρ)
(fonction strictement décroissante et concave).
La condition d’entropie (Lax) demande à ce que la vitesse σ̇1
vérifie ( u
λ1 (W) < σ̇1 < λ2 (W),
σ̇1 < λ1 (WL ),
c’est-à-dire WL
( uL
u − ρp 0 (ρ) < σ̇1 < u,
σ̇1 < u L − ρ L p 0 (ρ L ).

Par conséquent p(ρ) + ρp 0 (ρ) > p(ρ L ) + ρ L p 0 (ρ L ) ; la fonction

f : ρ 7→ ρp 0 (ρ) + p(ρ) est strictement croissante donc de f (ρ) > ρL ρ
f (ρ L ) on déduit ρ L < ρ et par conséquent u L > u.

8. À l’aide du dessin d’ondes dans le plan (ρ, u) résoudre le problème de R IEMANN : pour un état gauche WL et un état droit
WR on construira une solution composée d’une 1-onde et d’une 2-onde séparant un état intermédiaire W∗ . On précisera les
valeurs de cet état intermédiaire ainsi que les vitesses des ondes.
Soit un problème de Riemann avec les deux états donnés suivants :

ρL ρR
µ
¶ µ ¶
WL = , WR = .
uL uR

La solution est constituée de trois états constants

ρL ρ ∗ = p −1 (−u R + u L + p(ρ L )) ρR
µ ¶ µ ¶ µ ¶
WL = , W∗ = , WR = ,
uL u∗ = uR uR

séparés par deux ondes : la première est soit un choc soit une raréfaction, la deuxième est toujours une discontinuité de
contact de vitesse u R . On note que la courbe intégrale du 1-choc et la courbe intégrale de la 1-raréfaction ont la même
expression, on parle alors de système de T EMPLE.
B Si u R < u L on a un 1-choc et la solution est

WL , si x < σ̇1 t ,


ρ ∗ uR − ρ L uL

W(x, t ) = W∗ , si σ̇1 t < x < u R t , avec σ̇1 = .
 ρ∗ − ρL
WR , si x > u R t ,


u u

t 2-ddc
WL WL
uL uL

1-choc W∗
u∗ = uR WR W∗ u∗ = uR W∗ WR
ρR ρL ρ∗ ρ ρL ρ∗ ρR ρ
WL WR

x
B Si u R > u L on a une 1-raréfaction et la solution est


WL , si x < (u L − ρ L p 0 (ρ L ))t ,

W , si (u − ρ p 0 (ρ ))t < x < (u − ρ p 0 (ρ ))t ,

d L L L R ∗ ∗
W(x, t ) =


W∗ , si (u R − ρ ∗ p 0
(ρ ∗ ))t < x < u R t ,

WR , si x > u R t ,


x̂
Pour calculer Wd en un point (x̂, t̂ ) à l’intérieure de l’onde de raréfaction (i.e. λ1 (WL ) < t̂
< λ1 (W∗ )) on considère la
caractéristique qui relie (x̂, t̂ ) à l’origine de l’onde (0, 0). La vitesse caractéristique de l’onde est dx
dt = λ1 (W), c’est-à-
dire x̂t̂ = u − ρp 0 (ρ). De plus, l’invariant de R IEMANN est conservé dans l’onde, c’est-à-dire u + p(ρ) = u L + p(ρ L ). Par
conséquent, ρ(x̂, t̂ ) est la solution de l’équation p(ρ) + ρp 0 (ρ) = u L + p(ρ L ) − x̂t̂ et u(x̂, t̂ ) = −p(ρ(x̂, t̂ )) + u L + p(ρ L ).

14 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 2. Étude d’un système hyperbolique

u u

WR WR 2-ddc
u∗ = uR W∗ u∗ = uR W∗ 1-détente t

WL WL W∗
uL uL

ρR ρ∗ ρL ρ ρ∗ ρL ρR ρ
WL WR

x
B Si u R = u L on a W∗ = WL et la solution est
(
WL , si x < u R t ,
W(x, t ) =
WR , si x > u R t .

t 2-ddc

u u
WL WR

W∗ = WL
u∗ = uR = uL WR W∗ = WL u ∗ = u R = u L WR
ρR ρ∗ = ρL ρ ρ∗ = ρL ρR ρ
x

© G. Faccanoni 15
 Deuxième partie .

Rappels Théoriques

16
 3. Équations linéaires
On cherche

u : R × R+ → R
(x, t ) 7→ u(x, t )

solution faible du problème

(
∂t u(t , x) + a(t , x)∂x u(t , x) + b(t , x)u(t , x) = f (t , x), x ∈ R, t > 0,
u(x, 0) = g (x), x ∈ R.

Soit (t̂ , x̂) ∈]0; +∞[×R donné. On appelle courbe caractéristique passant par (t̂ , x̂) la solution (si elle existe) du problème de
C AUCHY (
x 0 (t ) = a(t , x(t )) t > 0,
x(t̂ ) = x̂.
En particulier, si a(t , x) = A constante, alors x(t ) = x̂ + A(t − t̂ ).

Notons v(t ) ≡ u(t , x(t )) la restriction de u à la courbe caractéristique. Alors u(t̂ , x̂) = v(t̂ ) et l’on a

v 0 (t ) = ∂t u(t , x(t )) + x 0 (t )∂x u(t , x(t )) = ∂t u(t , x(t )) + a(t , x(t ))∂x u(t , x(t )).

B Si b(t , x) = f (t , x) = 0 pour tout (t , x) ∈ [0; +∞[×R alors v 0 (t ) = 0 donc

u(t̂ , x̂) = v(t̂ ) = v(0) = u(0, x(0)) = g (x(0)).

En particulier, si a(t , x) = A constante, alors x(t ) = x̂ + A(t − t̂ ) et donc u(t̂ , x̂) = g (x̂ − A t̂ ).
B Si f (t , x) = 0 pour tout (t , x) ∈ [0; +∞[×R alors v 0 (t ) + b(t , x(t ))v(t ) = 0 donc
Rt Rt Rt
b(s,x(s))d s b(s,x(s))d s b(s,x(s))d s
v(t ) = v(0)e − 0 = u(0, x(0))e − 0 = g (x(0))e − 0 .

On a alors R t̂
b(s,x(s))d s
u(t̂ , x̂) = v(t̂ ) = g (x(0))e − 0 .
En particulier, si b(t , x) = B constante, alors u(t̂ , x̂) = g (x(0))e −B t̂ . Si, de plus, a(t , x) = A constante, alors x(t ) = x̂ + A(t − t̂ ) et donc u(t̂ , x̂) = g (x̂ − A t̂ )e −B t̂ .
B On doit résoudre l’EDO v 0 (t ) + b(t , x(t ))v(t ) = f (t , x(t )). Comme on connait déjà la solution de l’équation homogène, il ne
reste à calculer qu’une solution particulière de l’équation complète et on trouve

R t̂ Z t̂ R t̂ R t̂ Z t̂ R t̂
b(s,x(s))d s b(r,x(r ))d r b(s,x(s))d s b(r,x(r ))d r
u(t̂ , x̂) = v(t̂ ) = v(0)e − 0 + f (s, x(s))e − s d s = g (x(0))e − 0 + f (s, x(s))e − s d s.
0 0

17
 4. Équations non-linéaires
On cherche
u : R × R+ → R
(x, t ) 7→ u(x, t )
solution faible entropique du problème (
∂t u + ∂x q(u) = 0, x ∈ R, t > 0,
u(x, 0) = g (x), x ∈ R,
avec le flux
q: R→R
u 7→ q(u)

Courbes caractéristiques : on appelle courbe caractéristique de l’EDP une courbe x = x(t ) dans le demi-espace t > 0 le long de
laquelle la solution u est constante.
Si pour un point (x, t ) passe une et une seule caractéristique et cette caractéristique a pied en (ξ, 0), alors u(x, t ) = g (ξ).
La caractéristique de pied (ξ, 0) a équation
x(t ) = ξ + q 0 (g (ξ))t .
Pour t suffisamment petit, la solution u est définie implicitement par l’équation
u = g (x − q 0 (g (ξ))t ).
Condition de Rankine-Hugoniot : si deux caractéristiques se croisent elles génèrent une discontinuité dans la solution.
La courbe de discontinuité est appelée onde de choc et si elle est régulière et a équation x = s(t ) alors u doit vérifier les
relations de R ANKINE -H UGONIOT : si on appelle u R et u L respectivement les traces de u à droite et à gauche de la courbe de
choc, on a
q(u R (s(t ), t )) − q(u L (s(t ), t ))
s 0 (t ) = .
u R (s(t ), t ) − u L (s(t ), t )
Si (x 0 ; t 0 ) est le premier point où les caractéristiques se croisent, alors on obtient l’équation de la courbe de choc en résolvant
l’EDO avec la condition initiale s(t 0 ) = x 0 .
Condition d’entropie : une condition qui permet de sélectionner parmi les solutions faibles qui vérifient les relations de saut
de R ANKINE -H UGONIOT la solution «physique» est la condition d’entropie de Lax le long de la courbe de choc x = s(t ) :
B si le flux q est convexe, le choc est entropique si et seulement si
q 0 (u R (s(t ), t )) < s 0 (t ) < q 0 (u L (s(t ), t )),
B si le flux q est concave, le choc est entropique si et seulement si
q 0 (u L (s(t ), t )) < s 0 (t ) < q 0 (u R (s(t ), t )),
En particulier, lorsque l’on considère le générique problème de R IEMANN
(
u L , si x < 0
g (x) =
u R , si x > 0,

B si le flux est convexe (par exemple q(u) = u 2 /2) on a

Bune onde de choc ssi u L > u R
Bune onde de raréfaction ssi u L < u R
B si le flux est concave (par exemple q(u) = u(2 − u)/2) on a
B une onde de choc ssi u L < u R
B une onde de raréfaction ssi u L > u R
Onde de raréfaction ou de détente : dans les régions du demi-espace t > 0 qui ne sont pas rejointes par les caractéristiques
on construit la solution (qui se connecte avec continuité avec les autres régions du demi-espace) par des ondes de
raréfaction.
Une onde de raréfaction centrée en (x 0 , t 0 ) a équation
0 −1 x − x 0
µ ¶
u(x, t ) = (q ) .
t − t0

18
Dernière mise à jour : Lundi 18 mars 2013 4. Équations non-linéaires

4.1. Exemple : équation de Burgers

³ ´
u2
Considérons l’équation de B URGERS ∂t u + ∂x 2 = 0 avec la donnée initiale u(x, 0) = g (x). Il s’agit d’une équation de type
∂t u + ∂x q(u) = 0 avec

u2
q(u) = (convexe) q 0 (u) = u, (q 0 )−1 (r ) = r.
2

on a q(u) = u 2 /2, q 0 (u) = u et (q 0 )−1 (r ) = r . L’équation de la caractéristique de pied (ξ, 0) est alors

ξ + Lt , ξ < −1,


0
x(t ) = ξ + q (g (ξ))t = ξ + (g (ξ))t = ξ + M t , −1 < ξ < 1,

ξ + R t , ξ > 1.


Donnée “ 2 − 1” Considérons la donnée initiale

(
2, si x < 1,
u(x, 0) = g (x) ≡
−1, si x > 1.

L’équation de la caractéristique de pied (ξ, 0) est donc

(
0 ξ + 2t , ξ < 1,
x(t ) = ξ + q (g (ξ))t = ξ + (g (ξ))t =
ξ− t, ξ > 1.

1 x

La donnée initiale du problème de Riemann a un saut décroissant et le flux est convexe donc l’unique solution entropique
présente une onde de choc qui part en (1, 0). Pour calculer l’équation x = s(t ) de l’onde de choc on utilise les relations de
R ANKINE -H UGONIOT :
q(u R ) − q(u L ) u R + u L 1
s 0 (t ) = = =
uR − uL 2 2
avec la donnée initiale s(0) = 1. On trouve donc x = s(t ) = 1 + 21 t .

1 x

La solution est donc (

2, x < 1 + 12 t ,
u(x, t ) =
−1, x > 1 + 21 t .

t =0
t = 0.5
t =1
x

© G. Faccanoni 19
4. Équations non-linéaires Dernière mise à jour : Lundi 18 mars 2013

Donnée “ 012” 
0 si x < 1,

g (x) = 1 si 1 < x < 2,

2 si x > 2.


L’équation de la caractéristique de pied (ξ, 0) est

ξ, ξ < 1,


0
x(t ) = ξ + q (g (ξ))t = ξ + (g (ξ))t = ξ + t , 1 < ξ < 2,

ξ + 2t , ξ > 2.


1 2 x

La donnée initiale a deux sauts croissants, un en x = 1 et l’autre en x = 2. Puisque q est convexe, on s’attend à ce que la solution
faible entropique présente deux ondes de raréfaction centrées l’une en (1, 0) et l’autre en (2, 0).

1 2 x




 0, si x < 1,
 x−1

t , si 1 < x < 1 + t ,



u(x, t ) = 1, si 1 + t < x < 2 + t ,

 x−2
t , si 2 + t < x < 2 + 2t ,





2, si x > 2t + 2,

t =0 u

t = 1.5

t =3

20 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 4. Équations non-linéaires

Donnée “ 210” Considérons la donnée initiale


2 si x < 1,

g (x) = 1 si 1 < x < 2,

0 si x > 2.


Elle a deux sauts décroissants, un en x = 1 et un en x = 2. Puisque q est convexe, on s’attend à ce que la solution faible entropique
présente deux ondes de choc qui partent de (1, 0) et de (2, 0).
B En effet, l’équation de la caractéristique de pied (ξ, 0) est

ξ + 2t , ξ < 1,


x(t ) = ξ + q 0 (g (ξ))t = ξ + (g (ξ))t = ξ + t , 1 < ξ < 2,

ξ, ξ > 2.


1 2 x
B Pour calculer les équations des deux ondes de choc on utilise les relations de R ANKINE -H UGONIOT :
B choc de pied (1, 0)
q(u R )−q(u L ) (u R )2 −(u L )2
(
1−4 3
s 0 (t ) = u R −u L = 2(u R −u L ) = 2(1−2) = 2
s(0) = 1

d’où x = s(t ) = 23 t + 1 ;
B choc de pied (2, 0)
q(u R )−q(u L ) (u R )2 −(u L )2
(
0−1 1
s 0 (t ) = u R −u L = 2(u R −u L ) = 2(0−1) = 2
s(0) = 2
1
d’où x = s(t ) = 2t + 2.
t

(x ∗ , t ∗ )

1 2 x
Ces deux chocs vont ensuite interagir à partir d’un (x ∗ , t ∗ ) qu’on détermine en résolvant le système linéaire
(
x = 12 t + 2,
x = 32 t + 1.

On obtient ainsi le pied de la nouvelle onde de choc qui est (x ∗ , t ∗ ) = ( 25 , 1). Donc la solution pour 0 < t < 1 est

3

2, si x < 2 t + 1,

u(x, t ) = 1, si 2 t + 1 < x < 12 t + 2,
3

0, si x > 21 t + 2,


© G. Faccanoni 21
4. Équations non-linéaires Dernière mise à jour : Lundi 18 mars 2013

B Pour t ≥ 1 seule les données initiales pour x < 1 et x > 2 sont transportées et ce nouveau choc, qui part de ( 52 , 1), vérifie encore
la relation de R ANKINE -H UGONIOT :
q(u R )−q(u L ) (u R )2 −(u L )2
(
0−4
s 0 (t ) = u R −u L = 2(u R −u L ) = 2(0−2) =1
5
s(1) = 2

d’où x = s(t ) = t + 32 .
t

1 2 x
L’unique solution faible entropique pour tout t > 0 est donc

3 3

2, si x < min{ 2 t + 1, t + 2 },

u(x, t ) = 1, si 32 t + 1 < x < 12 t + 2,

0, si x > max{ 12 t + 2, t + 32 },


u
t =0

t = 0.5

t =1

t = 1.5

Donnée “ 010” Considérons la donnée initiale


0 si x < 1,

g (x) = 1 si 1 < x < 2,

0 si x > 2.


L’équation de la caractéristique de pied (ξ, 0) est

ξ, ξ < 1,


x(t ) = ξ + q 0 (g (ξ))t = ξ + (g (ξ))t = ξ + t , 1 < ξ < 2,

ξ, ξ > 2.


1 2 x

On s’attend donc à ce que la solution faible entropique présente une onde de raréfaction centrée en (1, 0) et une onde de choc qui
part de (2, 0) et que, après un certain temps, les deux ondes interagissent.

22 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 4. Équations non-linéaires

B La raréfaction est comprise entre la droite d’équation

x(t ) = 1
et la droite d’équation
x(t ) = 1 + t
et a équation
x −1 x −1
µ ¶
u(x, t ) = (q 0 )−1 = .
t −0 t
B Pour déterminer l’équation de l’onde de choc qui a pied en (2, 0) on utilise les relations de R ANKINE -H UGONIOT :
q(u R )−q(u L ) (u R )2 −(u L )2
(
1−0 1
s 0 (t ) = u R −u L = 2(u R −u L ) = 2(1−0) = 2
s(0) = 2

d’où x = s(t ) = 21 t + 2.
Ces considérations sont valables jusqu’à ce que l’onde de choc et l’onde de raréfaction se touchent, c’est-à-dire jusqu’à t = 2.
Donc la solution pour 0 < t < 2 est 
0,
 si x < 1,

 x−1 , si 1 < x < t + 1,

u(x, t ) = t


 1, si t + 1 < x < 21 t + 2,
si x > 12 t + 2.

0,


1 2 3 x
B Pour t > 2 on continue à avoir une onde de choc mais l’état gauche n’est plus constante car il est donné par la raréfaction i.e.
s(t )−1
uL = t donc le nouveau choc a pied en (3, 2) et vérifie
 ³ ´2
 q(u R )−q(u L ) (u R )2 −(u L )2 0− s(t t)−1 s(t )−1
0
s (t ) = = = =

u R −u L 2(u R −u L ) 2(0− s(t t)−1 ) 2t

s(2) = 3

p
d’où l’équation x = 1 + 2t .
t

1 2 3 x
L’unique solution faible entropique pour t > 0 est donc

 0, si x < 1,
p


x−1

 t , si 1 < x < min{1 + t , 1 + 2t },



u(x, t ) = 1, si 1 + t < x < 12 t + 2,

0, si x > 12 t + 2 et t < 2,





 p
0, si x > 1 + 2t et t ≥ 2.

© G. Faccanoni 23
4. Équations non-linéaires Dernière mise à jour : Lundi 18 mars 2013

t =0
u
t =1

t =2

t =3

t =4

x
q q
2 1
Remarquons que l’amplitude du choc vaut t et le choc se déplace à vitesse 2t . Par conséquent, pour t → ∞, l’amplitude et la
vitesse du choc tendent vers zéro.

Donnée “ 101” 
1 si x < 1,

g (x) = 0 si 1 < x < 2,

1 si x > 2.


L’équation de la caractéristique de pied (ξ, 0) est

ξ + t , ξ < 1,


0
x(t ) = ξ + q (g (ξ))t = ξ + (g (ξ))t = ξ, 1 < ξ < 2,

ξ + t , ξ > 2.


1 2 x

La donnée initiale a un saut décroissant en x = 1 et un saut croissant en x = 2. Puisque q est convexe, on s’attend à ce que
la solution faible entropique présente une onde de choc qui part de (1, 0) et une onde de raréfaction centrée en (2, 0). Pour
déterminer l’équation de l’onde de choc on utilise les relations de R ANKINE -H UGONIOT :

q(u R )−q(u L ) (u R )2 −(u L )2

(
1−0 1
s 0 (t ) = u R −u L = 2(u R −u L ) = 2(1−0) = 2
s(0) = 1

d’où x = 21 t + 1.

1 2 x

24 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 4. Équations non-linéaires

Ces considérations sont valables jusqu’à t = 2, ensuite on continue à avoir une onde de choc mais avec
 ³ ´2
q(u R )−q(u L ) (u R )2 −(u L )2 12 − s(tt )
= 12 + s(t )

 0
s (t ) = u −u = 2(u −u ) =
R L R L 2(1− s(tt ) ) 2t

s(2) = 2

p
d’où l’équation x = t − 2t + 2.

1 2 x

 p

1, si x < min{t /2 + 1, t − 2t + 2},


0, si t + 1 < x < 2 et t < 2,
u(x, t ) = x−2 p
 t ,

 si max{2, t − 2t + 2} < x < t + 2,

1, si x > t + 2.


t =0

t =1
u
t =2

t =3

Donnée “ L < M < R ” 

L
 si x < −1,
g (x) = M si − 1 < x < 1,

R si x > 1.


Comme L < M < R, l’unique solution faible entropique présente deux ondes de raréfaction centrées l’une en (−1, 0) et l’autre en
(1, 0) : 


 L, si x < −1,
x+1

 t , si − 1 < x < −1 + Lt ,



u(x, t ) = M, si − 1 + Lt < x < 1 + M t ,

x−1
, si 1 + M t < x < 1 + R t ,


 t



R, si x > 1 + R t .

Donnée “ L > M > R ” 

L
 si x < −1,
g (x) = M si − 1 < x < 1,

R si x > 2.


Comme L > M > R, l’unique solution faible entropique présente deux ondes de choc qui partent de (−1, 0) et de (1, 0).
Pour calculer les équations des deux ondes de choc on utilise les relations de R ANKINE -H UGONIOT :

© G. Faccanoni 25
4. Équations non-linéaires Dernière mise à jour : Lundi 18 mars 2013

B choc de pied (−1, 0) :

(u R )2 (u L )2

q(u R )−q(u L ) 2 − 2 u R +u L M +L
s 10 (t ) =

u R −u L = u R −u L = 2 = 2 ,
s 1 (0) = −1,


M +L
d’où x = s 1 (t ) = 2 t −1;
B choc de pied (1, 0) :
(u R )2 (u L )2

q(u R )−q(u L ) 2 − 2 u R +u L R+M
s 20 (t ) =

u R −u L = u R −u L = 2 = 2 ,
s 2 (0) = 1,


R+M
d’où x = s 2 (t ) = 2 t + 1.
M +L R+M ∗ ∗
Comme 2 > 2 , ces deux chocs vont interagir à partir d’un (x , t ) qu’on détermine en résolvant le système linéaire
(
M +L
x= 2 t − 1,
R+M
x= 2 t + 1.

On obtient ainsi le pied de la nouvelle onde de choc qui est (x ∗ , t ∗ ) = ( M +R 4

L−R , L−R ). Donc la solution pour 0 < t <
4
L−R est

M +L

L, si x < 2 t − 1,

M +L R+M
u(x, t ) = M , si 2 t − 1 < x < 2 t + 1,

R, si x > R+M
2 t + 1.


4
Pour t ≥ L−R seule les données initiales pour x < −1 et x > 1 sont transportées et ce nouveau choc, qui part de ( M +R 4
L−R , L−R ), vérifie
encore la relation de R ANKINE -H UGONIOT :
(u R )2 (u L )2

q(u )−q(u ) −
s 3 (t ) = uR −u L == 2u −u 2 = uR +u = R+L
 0 L
R L R L 2 2 ,
4 M +R
s 3 ( L−R ) = L−R ,


M +R
d’où x = s 3 (t ) = L−R t − 1.
4
L’unique solution faible entropique pour tout t ≥ L−R est donc
(
M +R
L, si x < L−R t − 1,
u(x, t ) = M +R
R, si x > L−R t − 1.

4.2. Exercice : trafic dans un tunnel

Soit ρ ≡ ρ(x, t ) la densité des voitures et q ≡ q(ρ) le flux des voitures. Un modèle qui décrit la vitesse des voitures dans un tunnel
est le suivant (
∂t ρ + ∂x (q(ρ)) = 0, x ∈ R, t ∈ R+ ,
ρ(x, 0) = g (x), x ∈ R,
avec q(ρ) = ρv(ρ) et v ≡ v(ρ) la vitesse des voitures ainsi définie :

si 0 ≤ ρ ≤ ρ c ,
(
vm ,
v(ρ) = ρm
³ ´
λ ln ρ , si ρ c ≤ ρ ≤ ρ m ,

et λ ≡ (v m )/(ln(ρ m /ρ c )). Les constantes v m et ρ m représentent respectivement la vitesse maximale, la densité maximale et
ρ c la densité critique en dessous de laquelle les voitures peuvent voyager à la vitesse maximale. (Valeurs vraisemblables sont
v m = 90 km h−1 , ρ m = 110 voitures/km et ρ c = 7 voitures/km qui donnent λ ≈ 32,727272.)
Supposons que l’entrée du tunnel est en x = 0 et qu’on ouvre le tunnel à l’instant t = 0 ; supposons de plus que précédemment les
voitures s’étaient accumulées en un bouchon. La donnée initiale est donc
(
ρ m , si x < 0,
g (x) =
0, si x > 0.

1. Calculer la densité et la vitesse du trafic et les dessiner dans les plans (x, ρ) et (x, v) pour t = 0 et pour t̄ > 0.
2. Déterminer et dessiner dans le plan (x, t ) la trajectoire d’une voiture qu’à l’instant initiale se trouve en x = x̄ < 0 et calculer à
quel instant elle rentre dans le tunnel.

26 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 4. Équations non-linéaires

Solution
1. Le graphe de la vitesse en fonction de la densité est le suivant (remarquons que v est continue) :
v

vm

80

70

60

50

40

30

20

10

ρc
10 20 30 40 50 60 70 80 90 100 ρm ρ

Le flux q s’écrit
v m ρ, si 0 ≤ ρ ≤ ρ c ,
(
q(ρ) = ³
ρm
´
λρ ln ρ , si ρ c ≤ ρ ≤ ρ m ,

donc
si 0 ≤ ρ ≤ ρ c ,
(
0
vm ,
q (ρ) = ρm
h ³ ´ i
λ ln ρ − 1 , si ρ c ≤ ρ ≤ ρ m .

Les graphes de q et q 0 en fonction de la densité sont les suivants (remarquons le saut de q 0 en ρ c ) :

q
q0

1300
90 q 0 (ρ Lc ) = v m
1200
80

1100
70

1000
60
q 0 (ρ Rc ) = (v m − λ)
900
50

800
40

700
30

600
20

500
10

400 ρc ρm
00
10 20 30 40 50 60 70 80 90 100 ρ
300
−10

200
−20

100
−30
ρc
10 20 30 40 50 60 70 80 90 100 ρm ρ
−40

Le saut de q 0 en ρ c est
‚q 0 (ρ c )ƒ = q 0 (ρ Rc ) − q 0 (ρ Lc ) = (v m − λ) − v m = −λ.

La caractéristique de pied (ξ, 0) a équation

(
0 ξ − λt , si ξ < 0,
x = ξ + q (g (ξ))t =
ξ + v m t , si ξ > 0.

x = −λt

x = vm t

© G. Faccanoni 27
4. Équations non-linéaires Dernière mise à jour : Lundi 18 mars 2013

On peut donc écrire la solution au moins pour certaines parties du plan :

ρ m , si x < −λt ,


ρ(x, t ) = ?, si − λt < x < v m t ,

0, si x > v m t .


Sachant que q 0 est discontinue en ρ = ρ c , on va donc étudier ρ séparément sur

S 1 ≡ {(x, t ) : −λt < x < v m t , ρ c < ρ(x, t ) < ρ m } = {(x, t ) : −λt < x < (v m − λ)t }

et
S 2 ≡ {(x, t ) : −λt < x < v m t , 0 < ρ(x, t ) < ρ c } = {(x, t ) : (v m − λ)t < x < v m t }.
B Étude dans S 1 : pour ρ c < ρ(x, t ) < ρ m on a q 00 (ρ) = −λ/ρ < 0. Puisque ρ est décroissante on cherche une solution sous
forme d’onde de rarefaction, centrée en (0, 0), qui se colle avec continuité avec la droite x = −λt (de cette manière elle
vérifie aussi la condition d’entropie). Cette onde a équation
³ x ´
ρ(x, t ) = (q 0 )−1 (x/t ) = ρ m exp −1 − .
λt
B Étude dans S 2 : pour ρ(x, t ) < ρ c on a q 0 (ρ) = v m . Par conséquent ρ est constante le long des lignes caractéristiques du
type
x = v m t + k.
Ces caractéristiques sortent de la droite x = (v m − λ)t et sur cette droite, qui appartient à la région S 1 , la densité vaut
(v m − λ)t ³ v ´
µ ¶
m
ρ(x, t ) = ρ((v m − λ)t , t ) = ρ m exp −1 − = ρ m exp − = ρc .
λt λ
t

x = −λt x = (v m − λ)t

x = vm t

En résumé, on a 

 ρm , si x < −λt ,

ρ exp ¡−1 − x ¢ ,

si − λt < x < (v m − λ)t ,
m λt
ρ(x, t ) =


 ρ c , si (v m − λ)t < x < v m t ,

0, si x > v m t .


ρ
ρm
t =0
t = 0.2
t = 0.4
t = 0.6
t = 0.8

ρc
x
Pour un temps fixé, on voit que la densité décroît de la densité maximale (vitesse nulle) graduellement jusqu’à la densité
critique qui permet la vitesse maximale. On remarque la discontinuité de la densité pour x = v m t .
2. Considérons une voiture qu’à l’instant t = 0 se trouve en x = x̄ < 0. La voiture restera immobile jusqu’à l’instant t 0 = |x̄|/λ. À
cet instant elle rentre dans la région S 1 où la vitesse est donnée par
x
v(ρ(x, t )) = λ + .
t

28 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 4. Équations non-linéaires

v
vm

t =0
t = 0.2
t = 0.4
t = 0.6
t = 0.8

x
Si on indique par x = x(t ) la trajectoire de la voiture, on a

x 0 (t ) = λ + x(tt ) ,
(

x(t 0 ) = x̄.

On obtient
λt
µ µ ¶ ¶
x(t ) = λt ln −1 .
|x̄|
e|x̄|
La voiture entre donc dans le tunnel à l’instant T pour lequel x(T ) = 0, c’est-à-dire T = λ .
x

t 0 = |x̄|/λ
T= e|x̄| t
λ

© G. Faccanoni 29
 5. Systèmes hyperboliques
Dans ce dernier chapitre on s’intéresse à la résolution de problèmes de Riemann associés aux systèmes hyperboliques non-
linéaires en une dimension d’espace. Plus précisément, on cherche une fonction

W : R+ × R → R p
(t , x) 7→ W(t , x)

qui vérifie (au sens faible) le système d’EDPs

∂t W + ∂x F(W) x ∈ R, t > 0,


 = 0,
(
WL si x < 0,
W(0, x) =

WR si x > 0,

avec le flux

F : Rp → Rp
W 7→ F(W)

Pour des solutions régulières on peut réécrire ce système sous la forme quasi-linéaire

∂t W + B(W)∂x W = 0

avec B(W) la matrice jacobienne du flux F et on note λ1 (W) < λ1 (W) < · · · < λp (W) ses valeurs propres.

5.1. Exemple : le système de Saint-Venant

On cherche à simuler un écoulement d’eau sous l’hypothèse de faible
profondeur dans le cas 1D et à fond plat. De plus, l’eau sera considé-
rée comme un fluide incompressible non visqueux et les frottements
air/eau et eau/sol seront négligés. Les inconnues, fonction du temps
t ∈ [0; +∞[ et de l’espace x ∈ R, sont h ≡ h(x, t ) > 0 la hauteur de l’eau h(x, t ) : hauteur de l’eau
et u ≡ u(x, t ) ∈ R sa vitesse horizontale (on suppose qu’elle est la même
sur toute la hauteur de l’eau). On note g > 0 la constante de gravité. En u(x, t ) : vitesse de l’eau
dimension un d’espace et pour une topographie plate, on modélise ce
type d’écoulements par le système de Saint Venant (ou des eaux peu x
profondes) (
∂t h + ∂x (hu) = 0,
avec x ∈ R, t > 0. (5.1)
∂t (hu) + ∂x (hu 2 + g h 2 /2) = 0,

Sous forme vectorielle il s’écrit

∂t U + ∂x F(U) = 0
ayant défini
µ ¶ µ ¶
h hu
U= , F(U) = g .
hu hu 2 + 2 h 2

1. Soit W = (h, u), alors le système (5.1) se réécrit, pour des solutions régulières, sous la forme quasi-linéaire

∂t W + A(W)∂x W = 0. (5.2)

Montrer que la matrice A(W) s’écrit µ ¶

u h
.
g u

30
Dernière mise à jour : Lundi 18 mars 2013 5. Systèmes hyperboliques

2. Calculer les deux valeurs propres λ1 (W) et λ2 (W) de la matrice A(W). Afin de fixer les notations on ordonne les deux valeurs
propres selon λ1 (W) < λ2 (W). Proposer une base associée de vecteurs propres à droite {r1 (W), r2 (W)}. En déduire que le
système (5.1) est strictement hyperbolique.
3. Vérifier que les champs 1 et 2 sont vraiment non linéaires.
4. On note I k l’invariant de Riemann du k-ème champ caractéristique. Montrer qu’un choix possible pour I 1 et I 2 est
q
I 1 = u + 2 g h,
q
I 2 = u − 2 g h.

5. Pour un état gauche WL = (h L , u L ) donné on cherche les états W = (h, u) qui peuvent être relié à WL par une onde de détente.
5.1. Considérons le 1-champ. Montrer que u > u L , h < h L . Calculer u en fonction de h L , u L et h. Plus précisément, montrer
que u peut se mettre su la forme
u = u L + r (h L , h)
en explicitant la fonction r . Étudier la fonction h 7→ u = u L + r (h L , h) et tracer son graphe dans le plan (h, u).
5.2. Reprendre la question pour le 2-champ.
6. Pour un état gauche WL = (h L , u L ) donné on cherche les états W = (h, u) qui peuvent être relié à WL par une onde de choc
entropique de vitesse σ̇k .
6.1. Montrer que les relations de R ANKINE -H UGONIOT peuvent s’écrire sous la forme
(
j = h(u − σ̇k ) = h L (u L − σ̇k ),
g g
uh(u − σ̇k ) + 2 h 2 = u L h L (u L − σ̇k ) + 2 h L2 .

En déduire que
2 2
g hL − h
j= .
2 u − uL
6.2. Considérons le 1-champ. Calculer u en fonction de h L , u L et h. Plus précisément, montrer que u peut se mettre sous
la forme
u = u L + d (h L , h)
en explicitant la fonction d . Étudier la fonction h 7→ u = u L + d (h L , h) et tracer son graphe dans le plan (h, u). Montrer
que u < u L , h > h L .
6.3. Reprendre la question pour le 2-champ.
7. À l’aide du dessin d’onde dans le plan (h, u) résoudre le problème de Riemann : pour un état gauche WL et un état droit WR
on construira une solution composée d’une 1-onde et d’une 2-onde séparant un état intermédiaire W∗ . On précisera les
valeurs de cet état intermédiaire ainsi que les vitesses des ondes.
8. Cherchons maintenant une entropie pour le système (5.1). Montrer que la fonction

hu 2 g h 2
η(W) = +
2 2
est une entropie du système avec flux d’entropie

hu 3
Φ(W) = + g h 2 u.
2

Solution
1. On développe les dérivées du système (5.1) pour des solutions régulières :
(
∂t h + u∂x h + h∂x u = 0,
u∂t h + h∂t u + 2hu∂x u + u 2 ∂x h + g h∂x h = 0.

Il se réécrit alors (
∂t h + u∂x h + h∂x u = 0,
u(−u∂x h − h∂x u) + h∂t u + 2hu∂x u + u 2 ∂x h + g h∂x h = 0
et on trouve le système quasi-linéaire suivant :
µ ¶ µ ¶ µ ¶ µ ¶
h u h h 0
∂t + ∂ = .
u g u x u 0

© G. Faccanoni 31
5. Systèmes hyperboliques Dernière mise à jour : Lundi 18 mars 2013

2. On cherche les deux solutions λk (W) de l’équation det(A(W) − λ(W)Id) = 0, i.e. de l’équation

(u − λ)2 − hg = 0.

On obtient
q q
λ1 (W) = u − g h < λ2 (W) = u + g h.

Puisque h 6= 0, les valeurs propres sont réelles et distinctes donc le système (5.1) est strictement hyperbolique. On peut alors
prendre
µ p ¶ µp ¶
− h h
r1 (W) = p , r2 (W) = p .
g g

3. Pour déterminer la nature des deux champs caractéristiques on calcule (∇λk )T · rk pour k = 1, 2 :

∂ u − gh ∂ u − gh p
¡ p ¢ ¡ p ¢
T
p 3p
(∇λ1 (W)) · r1 (W) = (− h) + g= g > 0, 1-champ VNL,
∂h ∂u 2
∂ u + gh p ∂ u + gh p
¡ p ¢ ¡ p ¢
3p
(∇λ2 (W))T · r2 (W) = h+ g= g > 0, 2-champ VNL.
∂h ∂u 2

Étant donné qu’aucun champ n’est LD, il n’y aura aucune discontinuité de contact.
4. On vérifie que (∇I k (W))T · rk (W) = 0 pour k = 1, 2 :

∂ u +2 gh ∂ u +2 gh p
¡ p ¢ ¡ p ¢
T
p
(∇I 1 (W)) · r1 (W) = (− h) + g = 0,
∂h ∂u
∂ u −2 gh p ∂ u −2 gh p
¡ p ¢ ¡ p ¢
(∇I 2 (W))T · r2 (W) = h+ g = 0.
∂h ∂u

5. Étude des détentes : on cherche à déterminer les états droits W = (h, u) qui peuvent être reliés à un état gauche WL = (h L , u L )
par une onde de détente.

5.1. 1-champ :
B Dans une détente les invariants de Riemann sont conservés. Ici k = 1 d’où

I 1 (WL ) = I 1 (W)

donc q q
uL + 2 g hL = u + 2
gh
p p p
d’où h 7→ u = u L + 2 g ( h L − h).
p p
λ1 (WL ) < λ1 (W) se réécrit u L − g h L < u − g h donc u > u L + g ( h − h L ). Comme
p p p
B La condition d’entropie
p p p
u = u L + 2 g ( h L − h), on obtient
h < hL et u > u L .
On obtient ainsi
u = u L + r (h L , h) pour h < h L
p p p
avec r (h L , h) = 2 g ( h L − h). De plus,
r
g
u 0 (h) = − < 0, pour h < h L ,
h
r
1 g
u 00 (h) = > 0, pour h < h L .
2 h3

On a donc les graphes suivants :

32 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 5. Systèmes hyperboliques

t
WL
uL

x=
x=

λ1
WL WR

(W
λ1
(W

R
)t
hL h

L
)t
x
Pour calculer la solution en un point (x̂, t̂ ) à l’intérieure de l’onde on considère la caractéristique qui relie ce point à
x̂
l’origine de l’onde (0, 0). La vitesse caractéristique de l’onde est dx λ
p
dt = 1 , c’est-à-dire t̂
= u − g h. De plus, l’invariant
³ ´
de Riemann est conservé dans l’onde, c’est-à-dire u + 2 g h = u L + 2 g h L . On obtient u(x̂, t̂ ) = 3 u L + 2 g h L + 2 x̂t̂
1
p p p
³ ´2
1
u L + 2 g h L − x̂t̂ .
p
et h(x̂, t̂ ) = 9g
5.2. 2-champ :
B Dans une détente les invariants de Riemann sont conservés. Ici k = 2 d’où
I 2 (WL ) = I 2 (W)

donc q q
uL − 2 g hL = u − 2
gh
p p p
d’où h 7→ u = u L + 2 g ( h − h L ).
p p p
λ2 (WL ) < λ2 (W) se réécrit u L + g h L < u + g h donc u > u L + g ( h L − h). Comme
p p
B La condition d’entropie
p p p
u = u L + 2 g ( h − h L ), on obtient
h > hL et u > u L .
On obtient ainsi
u = u L + r (h L , h) pour h < h L
p p p
avec r (h L , h) = 2 g ( h − h L ). De plus,
r
g
u 0 (h) = > 0, pour h > h L ,
h
r
1 g
u 00 (h) = − < 0, pour h > h L .
2 h3
On a donc les graphes suivants :
u

t
WL
uL
L )t
(W

t
R)

WL WR
W
λ2

2(
x=

hL h
x=

x
Pour calculer la solution en un point (x̂, t̂ ) à l’intérieure de l’onde on considère la caractéristique qui relie ce point à
x̂
l’origine de l’onde (0, 0). La vitesse caractéristique de l’onde est dx λ
p
dt = 2 , c’est-à-dire t̂
= u + g h. De plus, l’invariant
³ ´
de Riemann est conservé dans l’onde, c’est-à-dire u−2 g h = u R −2 g h R . On obtient u(x̂, t̂ ) = 3 −u R + 2 g h R + 2 x̂t̂
1
p p p
³ ´2
1
u R − 2 g h R + x̂t̂ .
p
et h(x̂, t̂ ) = 9g

© G. Faccanoni 33
5. Systèmes hyperboliques Dernière mise à jour : Lundi 18 mars 2013

6. Étude des chocs : on cherche à déterminer les états droits W = (h, u) qui peuvent être reliés à un état gauche WL = (h L , u L )
par une discontinuité de vitesse σ̇k . Les relations de R ANKINE -H UGONIOT pour un choc de vitesse σ̇k s’écrivent
h L u L −hu
σ̇k =
(
h L −h ,
h L u L2 +(g /2)h L2 −hu 2 −(g /2)h 2
σ̇k = h L u L −hu ,

qu’on peut réécrire comme (

hu − σ̇k h = h L u L − σ̇k h L ,
(u − σ̇k )hu − (u L − σ̇k )h L u L = (g /2)(h L2 − h 2 ),
ou encore (
j k ≡ h(u − σ̇k ) = h L (u L − σ̇k ),
j k (u − u L ) = (g /2)(h L2 − h 2 ).
6.1. 1-champ : la condition d’entropie (Lax) pour k = 1 demande à ce que la vitesse σ̇1 du 1-choc vérifie
(
λ1 (W) < σ̇1 < λ2 (W),
σ̇1 < λ1 (WL ),

c’est-à-dire (
g h < σ̇1 < u +
p p
u− g h,
σ̇1 < u L − g h L ,
p

donc
< σ̇1
 p  p
g (h 2 −h 2 )
u − g hp  j1 < h p gh L
(
 
h > hL , j1 = 2(u−u L )
σ̇1 < u + g h, =⇒ j 1 > −h g h, =⇒ =⇒ u < uL
  j 1 > 0,
σ̇1 < u L − g h L ,
p p
j 1 > hL g hL ,
 

et on conclut que
h > hL et u < u L .
En éliminant σ̇1 dans les relations de R ANKINE -H UGONIOT on trouve

h 7→ u = u L + d (h L , h) pour h > h L
q
g h L +h
avec d (h L , h) = (h L − h) 2 hL h et la vitesse du 1-choc est
s
g hL + h
σ̇1 = u L − h .
2 hL h

De plus,
2 2
g h L h + 2h + h L
r
0
u (h) = − p < 0, pour h > h L ,
8 h 3 h L (h L + h)
q
r
g h L3 (5h + 3h L )
u 00 (h) = p > 0. pour h > h L .
32 h 5 (h L + h)3

On a donc les graphes suivants :

u

x = σ̇1 t t
WL
uL

WL WR
hL h

34 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 5. Systèmes hyperboliques

6.2. 2-champ : la condition d’entropie (Lax) pour k = 2 demande à ce que la vitesse σ̇2 du 2-choc vérifie
(
λ2 (W) < σ̇2 ,
λ1 (WL ) < σ̇2 < λ2 (WL ),

c’est-à-dire (
u + g h < σ̇2 ,
p

u L − g h L < σ̇2 < u L + g h L ,

p p

donc
u + g h <pσ̇2
 p  p
g (h 2 −h 2 )
 j 2 < −hp g h L
(
 
h < hL , j2 = 2(u−u L )
σ̇2 > u L − g h L , =⇒ j 2 < h L g h L , =⇒ =⇒ u < uL
  j 2 < 0,
σ̇2 < u L + g h L ,
p p
j 2 > −h L g h L ,
 

et on conclut que
h < hL et u < u L .
En éliminant σ̇2 dans les relations de R ANKINE -H UGONIOT on trouve

h 7→ u = u L + d (h L , h) pour h < h L
q
g h L +h
avec d (h L , h) = (h − h L ) 2 hL h et la vitesse du 2-choc est
s
g hL + h
σ̇2 = u L + h .
2 hL h

De plus,
2 2
g h L h + 2h + h L
r
u 0 (h) =p > 0, pour h < h L ,
8 h 3 h L (h L + h)
q
r
g h L3 (5h + 3h L )
00
u (h) = − p < 0. pour h < h L .
32 h 5 (h L + h)3

On a donc les graphes suivants :

u

t x = σ̇2 t
WL
uL

WL WR
hL h

Récapitulatif : quel qui soit WL = (τL , u L ), le demi-plan R∗+ × R se décompose en quatre zones séparées par les quatre
demi-courbes 1-choc, 2-choc, 1-détente et 2-détente.

© G. Faccanoni 35
5. Systèmes hyperboliques Dernière mise à jour : Lundi 18 mars 2013

1-détente
2-détente
III IV
WL
uL
I

1-choc
II
hL h
2-choc

7. Soit un problème de Riemann avec les deux états constants donnés suivants :
µ ¶ µ ¶
hL hR
WL = , WR = .
uL uR

La solution est constituée de trois états constants séparés par deux ondes.
1-onde t

W∗
2-onde
WL

WR

x
Pour expliciter cette solution on cherche à définir l’inconnue

h∗
µ ¶
W∗ =
u∗

à l’aide de l’étude des ondes précèdent.

On a cinq cas possibles :

Cas 1) 1-choc et 2-choc (il correspond au cas où WR appartient à la zone I I , c’est-à-dire si h ∗ > h L et h ∗ > h R )
x = σ̇1 (WL , W∗ )t t x = σ̇2 (W∗ , WR )t

W∗

WL WR

L’unique solution faible entropique est

WL , si x < σ̇1 (h L , h , u L )t ,

∗


W(x, t ) = W∗ , si σ̇1 (h L , h ∗ , u L )t < x < σ̇2 (h ∗ , h R , u ∗ )t ,

WR , si x > σ̇2 (h ∗ , h R , u ∗ )t ,


36 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 5. Systèmes hyperboliques

avec
s
∗ ∗ g hL + h ∗
σ̇1 (h L , h , u L ) = u L − h
2 hL h ∗
s
∗ ∗ ∗ g h ∗ + hR
σ̇2 (h , h R , u ) = u + h R
2 h ∗ hR

et h ∗ et u ∗ est l’unique solution du système

 q
u ∗ = u + (h − h ∗ ) g hL +h∗∗ ,
L L
q2 h L h
u = u ∗ + (h − h ∗ ) g h ∗ +hR .
R R 2 h ∗ hR

Cas 2) 1-choc et 2-détente (il correspond au cas où WR appartient à la zone I , c’est-à-dire si h ∗ > h L et h ∗ < h R )
x = σ̇1 (WL , W∗ )t t

W∗

)t
∗
2 (W

t
R)
W
λ

2(
x=

λ
=
x
WL WR

L’unique solution faible entropique est



WL , si x < σ̇1 (h L , h ∗ , u L )t ,

W∗ , si σ̇1 (h L , h ∗ , u L )t < x < (u ∗ + g h ∗ )t ,
 p
W(x, t ) =
si (u ∗ + g h ∗ )t < x < (u R + g h R )t ,
p p
W2−det ,


 p
WR , si x > (u R + g h R )t ,


avec s
∗ g hL + h ∗
∗
σ̇1 (h L , h , u L ) = u L − h ,
2 hL h ∗
x 2
à ¡ ¢ !
1
p
−u R + 2 g h R + t¢
W2−det = 9g1 ¡ p x
,
3 uR − 2 g hR + 2 t

et h ∗ et u ∗ l’unique solution du système

q ∗
g
(
u ∗ = u L + (h L − h ∗ ) 2 hhL +h
h ∗ ,
p L
p p
u R = u ∗ + 2 g ( h R − h ∗ ).

Cas 3) 1-détente et 2-choc (il correspond au cas où WR appartient à la zone I I I , c’est-à-dire si h ∗ < h L et h ∗ > h R )

t x = σ̇2 (W∗ , WR )t

W∗
x=
x=

λ1
(W
λ1

∗
(W

)t
L
)t

WL WR

© G. Faccanoni 37
5. Systèmes hyperboliques Dernière mise à jour : Lundi 18 mars 2013

L’unique solution faible entropique est

 p

WL , si x < (u L − g h L )t ,

g h L )t < x < (u ∗ − g h ∗ )t ,
 p p
1−det , si (u L −
W
W(x, t ) =
W∗ , si (u ∗ − g h ∗ )t < x < σ̇2 (h ∗ , h R , u ∗ )t ,
p



si x > σ̇2 (h ∗ , h R , u ∗ )t ,

WR ,


avec s
∗ g h ∗ + hR
∗
σ̇2 (h , h R ) = u + h R ,
2 h ∗ hR
à ¡ ¢2 !
1
u L + 2 g h L − xt
p
9g
W1−det = 1 ¡ ¢ ,
x
p
3 uL + 2 g hL + 2 t

et h ∗ et u ∗ l’unique solution du système

( ∗ p p p
u = u L + 2 g ( h L − h ∗ ),
q ∗
g
u R = u ∗ + (h R − h ∗ ) 2 hh ∗+h R
h . R

Cas 4) 1-détente et 2-détente (il correspond au cas où WR appartient à la zone I V , c’est-à-dire si h ∗ < h L et h ∗ < h R )
t

W∗
x=

)t
∗
2 (W
x=

t
λ1

R)
W
(W
λ1

2(
∗

x=
(W

λ
)t

=
L
)t

WL WR

L’unique solution faible entropique est

 p


WL , si x < (u L − g h L )t ,
si (u L − g h L )t < x < (u ∗ − g h ∗ )t ,
 p p
 1−det ,
W



W(x, t ) = W∗ , si (u ∗ − g h ∗ )t < x < (u ∗ + g h ∗ )t ,
p p

si (u ∗ + g h ∗ )t < x < (u R + g h R )t ,
p p
W2−det ,





 p
W ,
R si x > (u R + g h R )t

avec
x 2
à ! à ¢2 !
1 1
−u R + 2 g h R + xt
¡ p ¢ ¡ p
W1−det = 9g ¡ u L + 2p g h L − t ¢ , W2−det = 9g ¢ ,
x
1 1
u R − 2 g h R + 2 xt
¡ p
3 uL + 2 g hL + 2 t 3

et h ∗ et u ∗ l’unique solution du système

( p p p
u ∗ = u L + 2 g ( h L − h ∗ ),
p p p
u R = u ∗ + 2 g ( h R − h ∗ ).

Remarque : il peut se générer une zone de hauteur h nulle lorsque WR appartient à la zone I V : un état intermédiaire
h∗ > 0 n’est possible que si
p p p
u R − u L < 2 g ( h R − h R ).

La méthode de résolution du problème de Riemann consiste donc à calculer cet état intermédiaire (h ∗ , u ∗ ). On a vu que cet
état peut être relié à l’état de gauche (h L , u L ) par une 1-onde :
p
u L − (h ∗ − h L ) p 2 gp

∗
si h ∗ < h L (1-détente)
∗
u = qhL +∗ h
∗ g h +h L
u L − (h − h L ) 2 h ∗ h si h ∗ > h L (1-choc)

L

38 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 5. Systèmes hyperboliques

et à l’état de droite (h R , u R ) par une 2-onde :

p
u R + (h ∗ − h R ) p 2 gp

∗
si h ∗ < h R (2-détente)
u∗ = qhR +∗ h
∗ g h +h R
u R + (h − h R ) 2 h ∗ h si h ∗ > h R (2-choc)

R

Pour calculer h ∗ on doit donc résoudre

u R + (h ∗ − h R )z(h ∗ , h R ) = u L − (h ∗ − h L )z(h ∗ , h L )

avec
 2p g
p p si h ≤ ψ,
z(h, ψ) = qh+ ψ
 g h+ψ sinon.
2 hψ

ce qui revient à chercher le zéro de la fonction

f (h) = u R − u L + (h − h L )z(h, h L ) + (h − h R )z(h, h R )

et u ∗ est déduit tout simplement par

u ∗ = u R + (h ∗ − h R )z(h ∗ , h R ).

Remarque : Notons que pour h L = 0 (resp. h R = 0), la 1-onde (resp. la 2-onde) de choc n’est pas définie, dans ce cas la seule
solution possible est une onde de détente.
8. Pour vérifier que η est une entropie du système (5.1) avec le flux d’entropie Φ, on montre que

∇W Φ(W) = ∇W η(W) · A(W).

Or,
hu 3
µ ¶
h
W= , Φ(W) = + g h 2 u,
u 2
donc
u3
à !
∇W Φ(W) = 2 + 2g hu,
3hu 2 2
2 +gh ,
et
hu 2 g h 2
η(W) = +
2 2
d’où Ã !
u2
+gh
∇W η(W) = 2 .
hu
Donc !T µ
Ã
u2
¶ Ã u3 !
+gh u h + 2g hu
∇W η(W) · A(W) = 2 · 2
= 3hu 2 = ∇W Φ(W).
hu g u + g h2
2

Il ne reste à prouver que la convexité de l’entropie en calculant la matrice hessienne :

µ ¶
g u
d2 η(W) = .
u h

Si g h − u 2 > 0 l’entropie est convexe.

5.2. Exemple : Euler isotherme en coordonnées Lagrangiennes

On considère le système d’équations aux dérivées partielles
(
∂t τ − ∂x u = 0,
(5.3)
∂t u + ∂x τ1 = 0,

(ce système est celui des équations d’Euler en régime isotherme en coordonnées Lagrangiennes).

© G. Faccanoni 39
5. Systèmes hyperboliques Dernière mise à jour : Lundi 18 mars 2013

1. Trouver les vecteurs V : R∗ × R → R2 et F(V) : R2 → R2 tels que le système (5.3) s’écrit

∂t V + ∂x F(V) = 0.

2. En se plaçant dans les variables W = (τ, u), le système (5.3) se réécrit, pour des solutions régulières, sous la forma quasili-
néaire
∂t W + A(W)∂x W = 0. (5.4)
Montrer que la matrice A(W) s’écrit
µ ¶
u −1
.
− τ12 0
3. Calculer les deux valeurs propres λ1 et λ2 . Afin de fixer les notations on ordonne les deux valeurs propres selon λ1 < λ2 .
Proposer une base associée de vecteurs propres à droite {r1 , r2 } de la matrice A. En déduire que le système (5.4) est
hyperbolique pour (τ, u) ∈ Ω où Ω est un ouvert de R2 à précises.
4. Sachant que τ est l’inverse d’une densité spécifique et u une vitesse, déterminer un domaine physique d’hyperbolicité.
5. Vérifier que les champs 1 et 2 sont vraiment non linéaires.
6. On note I k les invariants de Riemann du k-ème champ caractéristique. Montrer qu’un choix possible pour I k est
1
I1 = + u,
τ
1
I 2 = − + u.
τ
7. Écrire les relations de R ANKINE -H UGONIOT pour ce système.
8. Soit σ la vitesse d’une discontinuité qui vérifie les relations de R ANKINE -H UGONIOT qui relie un état gauche (τL , u L ) à un
état droit (τR , u R ).
8.1. Exprimer σ en fonction de τL et u L .
8.2. Exprimer u R en fonction de τL , τR et u L .
8.3. Parmi ces discontinuités, chercher celles vérifiant

σ > 0,


u R = 0,

τR = 1.


Tracer, dans un diagramme (τ, u), l’allure de l’ensemble des états à droite (τR , u R ) admissibles.

Solution
1. On a
τ
µ ¶ µ ¶
−u
V= , F(V) = 1 .
u τ

2. On développe les dérivées du système (5.9) pour des solutions régulières :

(
∂t τ − ∂x u = 0,
∂t u − τ12 ∂t τ = 0.

On a le système quasilinéaire suivant :

τ τ
µ ¶ µ ¶ µ ¶ µ ¶
0 −1 0
∂t + ∂ = .
u − τ12 0 x u 0
3. On cherche les 2 solutions λk (W) de l’équation det(A(W) − λ(W)Id) = 0, i.e. de l’équation
1
(−λ)2 − = 0.
τ2
On obtient
1 1
λ1 = − , λ2 = .
τ2 τ2
On peut alors prendre
µ 2¶ µ 2¶
τ −τ
r1 = , r2 = .
1 1
Pour τ 6= 0, les valeurs propres sont réelles et distinctes : le système (5.3) est strictement hyperbolique. Donc Ω = R∗ × R.

40 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 5. Systèmes hyperboliques

4. Physiquement τ > 0 donc le domaine physique d’hyperbolicité est R∗+ × R.

5. Pour déterminer la nature des trois champs on calcul ∇λk · rk pour k = 1, 2 :
∂(−1/τ2 ) 2 ∂(−1/τ2 ) 2
∇λ1 · r1 = τ + 1= 1-champ VNL
∂τ ∂u τ
∂(1/τ2 ) ∂(1/τ2 ) 2
∇λ2 · r2 = (−τ2 ) + 1=− 2-champ VNL
∂τ ∂u τ
6. On vérifie que ∇I k · rk = 0 pour k = 1, 2 :
∂( τ1 + u) ∂( τ1 + u)
∇I 1 · r1 = τ2 + 1 = 0,
∂τ ∂u
∂(− τ1 + u) ∂(− τ1 + u)
∇I 2 · r2 = (−τ2 ) + 1 = 0.
∂τ ∂u
7. Les relations de R ANKINE -H UGONIOT sont (
σ(τL − τR ) = −(u L − u R ),
1
σ(u L − u R ) = τL − τ1 .
R

8. Discontinuité.
8.1. Les relations de R ANKINE -H UGONIOT donnent
1 1
σ2 (τL − τR ) = −
τL τR
soit s
1
σ=± .
τL τR
8.2. La première relation de R ANKINE -H UGONIOT permet d’écrire
τL − τR
u R = u L + σ(τL − τR ) = u L ± p .
τL τR
q
8.3. On prend ici u L = 0, τR = 1 et on cherche une discontinuité se propageant à une vitesse positive σ = + τ1 . L’expression
R
de u R est
1 − τR
uR = p
τR
et son graphe est remprésenté sur la figure suivante
u

0
1 2 3 4 5 6 7 8 9 τ

−1

−2

−3

5.3. Exemple : le système de la dynamique des gaz isentropique par

relaxation
On considère le système de la dynamique des gaz isentropique en une dimension d’espace
(
∂t ρ + ∂x (ρu) = 0,
(5.5)
∂t (ρu) + ∂x (ρu 2 + p) = 0,
avec une loi d’état p ≡ p(ρ). On note ρ ≡ ρ(x, t ) > 0 la densité et u ≡ u(x, t ) ∈ R la vitesse du fluide.

© G. Faccanoni 41
5. Systèmes hyperboliques Dernière mise à jour : Lundi 18 mars 2013

1. Dans un premier temps on suppose que la loi d’état est celle d’un gaz parfait :

p(ρ) = ρ γ , γ > 1.

1.1. Trouver les vecteurs V : R+ × R → R2 et F(V) : R2 → R2 tels que le système (5.5) s’écrit

∂t V + ∂x F(V) = 0.

1.2. En se plaçant dans les variables W ≡ (ρ, u), le système (5.5) se réécrit, pour des solutions régulières, sous la forma
quasi-linéaire
∂t W + A(W)∂x W = 0. (5.6)
Montrer que la matrice A(W) s’écrit
ρ
à !
u
A(W) = p0 .
ρ u

1.3. Calculer les valeurs propres λ1 (W) et λ2 (W) de la matrice A(W). Afin de fixer les notations on ordonne les valeurs
propres selon λ1 (W) < λ2 (W). Proposer une base associée de vecteurs propres à droite {r1 (W), r2 (W)}. En déduire que
le système (5.6) est strictement hyperbolique.
1.4. Vérifier que les deux champs caractéristiques sont vraiment non linéaires.
2. Lorsque le fluide n’est pas un gaz parfait, le caractère vraiment non linéaire des champs caractéristiques de (5.5) rend la
résolution du problème de Riemann assez compliquée. C’est pourquoi il a été proposé de remplacer (5.5) par un système
approché, dit système de relaxation, qui s’écrit

∂t ρ + ∂x (ρu) = 0,


∂t (ρu) + ∂x (ρu 2 + π) = 0, (5.7)
³ ´
∂ (ρV ) + ∂ (ρuV ) = µρ 1 − V ,


t x ρ

avec µ > 0 un paramètre constant, V ≡ V (x, t ) une nouvelle variable dite volume de relaxation et π ≡ π(ρ,V ) une nouvelle
loi d’état µ ¶
2 1
π(ρ,V ) = p(1/V ) − a −V
ρ
où a > 0 est une constante et p(1/V ) est la «vraie» loi d’état du fluide (sur laquelle on ne fait aucune hypothèse).
2.1. Montrer formellement que, lorsque le paramètre µ tend vers +∞, alors on retrouve le système (5.5) à partir du
système (5.7) (et on aura V = 1/ρ).
2.2. Désormais et pour tout ce qui suit on ne s’intéresse qu’à la structure différentielle de (5.7), c’est-à-dire qu’on ne tient
pas compte du terme source. Autrement dit, on pose µ = 0. On suppose aussi que ρ et V sont strictement positifs.
2.2.1. Trouver les vecteurs U : R+ × R → R3 et H(U) : R3 → R3 tels que le système (5.7) s’écrive

∂t U + ∂x H(U) = 0.

2.2.2. En se plaçant dans les variables Y ≡ (ρ, u,V ), le système (5.7) se réécrit, pour des solutions régulières, sous la
forme quasi-linéaire
∂t Y + B(Y)∂x Y = 0. (5.8)
Montrer que la matrice B(Y) s’écrit
ρ
 
u 0
 2 1 ∂π 
B(Y) =  ρa 3 u ρ ∂V  .
0 0 u
2.2.3. Calculer les valeurs propres λ1 (Y), λ2 (Y) et λ3 (Y) de la matrice B(Y). Afin de fixer les notations on ordonne les va-
leurs propres selon λ1 (Y) < λ2 (Y) < λ3 (Y). Proposer une base associée de vecteurs propres à droite {r1 (Y), r2 (Y), r3 (Y)}.
En déduire que le système (5.8) est strictement hyperbolique.
2.2.4. Vérifier que les trois champs caractéristiques sont linéairement dégénérés. (Par conséquent, les ondes associées
à chaque champ sont des discontinuités de contact). Les ondes associées sont-elles entropiques ?
2.2.5. Pour chaque k = 1, 2, 3, on note I k et J k deux invariants de Riemann du k-ème champ caractéristique. Montrer
qu’un choix possible pour I k et J k est
a
I 1 = V, J1 = u − ,
ρ
I 2 = u, J 2 = π,
a
I 3 = V, J3 = u + .
ρ

42 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 5. Systèmes hyperboliques

2.2.6. En utilisant les invariants de Riemann et les relations de Rankine-Hugoniot, montrer que u et π sont constants
dans une 2-onde, tandis que les courbes de 1-onde et de 3-onde sont des droites dans le plan (u, π). Vérifier
en particulier que les courbes de détente et de choc coïncident (i.e. les détentes vérifient les conditions de
Rankine-Hugoniot).
2.2.7. À l’aide du dessin d’ondes dans le plan (u, π) résoudre le problème de Riemann : pour un état gauche YL et un
état droit YR on construira une solution composée d’une 1-onde, d’une 2-onde et d’une 3-onde séparant deux
états intermédiaires Y1 et Y2 . On précisera les valeurs de ces états intermédiaires ainsi que les vitesses des ondes.
2.2.8. Montrer que la fonction
1
Z ρ π(r,V )
η(Y) = ρu 2 + ρ dr
2 ρ0 r2
est une entropie du système (5.7). Pour cela on vérifiera que
· ¸
1 2
∂t η + ∂x (Φ) = 0 avec Φ(Y) = u ρu +G(ρ,V )
2
c’est-à-dire que
∇Y (Φ(Y))T = ∇Y (η(Y))T · B(Y)
(on ne vérifiera pas que η est convexe).
2.2.9. Expliciter la fonction G dans le flux d’entropie.
2.2.10. Montrer qu’une discontinuité de contact de la deuxième famille reliant un état gauche YL = (ρ L , u L ,VL ) à un
état droit Y = (ρ, u,V ) avec vitesse σ = λ2 (YL ) = λ2 (Y) vérifie

Φ(Y) − Φ(YL ) = σ(η(Y) − η(YL )),

autrement dit on a conservation de l’entropie.

Solution
1. Système (5.5).
1.1. On a
ρ ρu
µ ¶ µ ¶
V= , F(V) = .
ρu ρu 2 + p(ρ)
1.2. On développe les dérivées du système (5.5) pour des solutions régulières :
(
∂t ρ + ρ∂x u + u∂x ρ = 0,
dp
ρ∂t u + u∂t ρ + u 2 ∂x ρ + 2ρu∂x u + dρ ∂x ρ = 0.

Il se réécrit alors (
∂t ρ + ρ∂x u + u∂x ρ = 0,
dp
ρ∂t u − u ρ∂x u + u∂x ρ + u 2 ∂x ρ + 2ρu∂x u + dρ ∂x ρ = 0,
£ ¤

d’où (
∂t ρ + u∂x ρ + ρ∂x u = 0,
dp
ρ∂t u + uρ∂x u + dρ ∂x ρ = 0,
et donc finalement on a le système quasi-linéaire suivant :

u ρ
µ ¶ Ã ! µ ¶ µ ¶
ρ ρ 0
∂t + p0 ∂x = .
u ρ u u 0

1.3. On cherche les deux solutions λk (W) de l’équation det(A(W) − λk (W)Id) = 0, i.e. de l’équation

λ2 − 2uλ + u 2 − p 0 = 0.

On obtient q q
λ1 = u − p 0, λ2 = u + p0

et p 0 = γρ γ−1 > 0. On peut alors prendre

 r  r 
ρ γ−1 ρ γ−1
ρ/ p
p ¶ µ p 0¶
0
µ
−ρ/ p − γ , γ .
r1 = = r2 = =
1 1
1 1

Puisque p 0 > 0, les valeurs propres sont réelles et distinctes donc le système (5.5) est strictement hyperbolique.

© G. Faccanoni 43
5. Systèmes hyperboliques Dernière mise à jour : Lundi 18 mars 2013

1.4. Pour déterminer la nature des deux champs on calcule ∇λk · rk pour k = 1, 2 :
s
∂(u − ρ γ−1 ∂(u − p 0 ) γ−1
p p
p 0) q
∇λ1 · r1 = − + (1) = γρ γ−3 + 1 > 0, 1-champ VNL,
∂ρ γ ∂u 2
p s
∂(u + p 0 ) ρ γ−1 ∂(u + p 0 ) γ−1
p q
∇λ2 · r2 = + (1) = γρ γ−3 + 1 > 0, 2-champ VNL.
∂ρ γ ∂u 2

2. Considérons maintenant le système de relaxation (5.7).

2.1. La troisième équation du système (5.7) se réécrit

µ ¶
1£ 1
∂t (ρV ) + ∂x (ρuV ) = ρ
¤
−V .
µ ρ
³ ´
1
Si µ → +∞ alors ρ ρ − V → 0, c’est-à-dire V → ρ1 , ce qui implique π(ρ, 1/ρ) = p(ρ).
2.2. Désormais et pour tout ce qui suit on pose µ = 0. On suppose aussi que ρ et V sont strictement positifs.
2.2.1. On a

ρ ρu
   
2
U =  ρu  , H(V) = ρu + π .
ρV ρuV

2.2.2. On développe les dérivées du système (5.7) pour des solutions régulières :

∂t ρ + u∂x ρ + ρ∂x u = 0,


ρ∂t u + u∂t ρ + 2uρ∂x u + u 2 ∂x ρ + ∂π ∂π
∂ρ ∂x ρ + ∂V ∂x V = 0,

ρ∂ V + V ∂ ρ + uV ρ∂ ρ + ρV ∂ u + ρu∂ V = 0.

t t x x x

Il se réécrit alors 
∂t ρ + u∂x£ρ + ρ∂x u = 0,¤


ρ∂t u − u u∂x ρ + ρ∂x u + 2uρ∂x u + u 2 ∂x ρ + ∂π ∂π
∂ρ ∂x ρ + ∂V ∂x V = 0,

ρ∂ V − V u∂ ρ + ρ∂ u + uV ρ∂ ρ + ρV ∂ u + ρu∂ V = 0,
 £ ¤
t x x x x x

d’où 
∂t ρ + u∂x ρ + ρ∂x u2 = 0,


∂π
ρ∂t u + uρ∂x u + ρa 2 ∂x ρ + ∂V ∂x V = 0,

∂ V + u∂ V = 0,

t x

et donc finalement on a le système quasi-linéaire suivant :

  u ρ

ρ 0 ρ
   
0
 2 1 ∂π 
∂t  u  +  a 3 ρ
u ρ ∂V  ∂x u = 0 .
   
V 0 0 u V 0

2.2.3. On cherche les 3 solutions λk (Y) de l’équation det(B(Y) − λ(Y)Id) = 0, i.e. de l’équation
³ a2 ´
(u − λ) λ2 − 2uλ + u 2 − 2 = 0.
ρ

On obtient
a a
λ1 = u − , λ2 = u, λ3 = u + .
ρ ρ

On peut alors prendre

−ρ 2 ∂π
     
1 1
 a 2 ∂V 
r1 = −a/ρ 2  , r2 =  0 , r3 = a/ρ 2  .
0 1 0

Puisque par hypothèse a > 0, les valeurs propres sont réelles et distinctes donc le système (5.7) est strictement
hyperbolique.

44 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 5. Systèmes hyperboliques

2.2.4. Pour déterminer la nature des trois champs caractéristiques on calcule ∇λk · rk pour k = 1, 2, 3 :
∂(u − ρa ) ∂(u − ρa ) −a ∂(u − ρa )
∇λ1 · r1 = 1+ + 0=0 1-champ LD
∂ρ ∂u ρ2 ∂V
∂u ∂π −ρ 2 ∂u ∂u
∇λ2 · r2 = + 0+ 1=0 2-champ LD
∂ρ ∂V a 2 ∂u ∂V
∂(u + ρa ) ∂(u + ρa ) a ∂(u + ρa )
∇λ3 · r3 = 1+ + 0=0 3-champ LD
∂ρ ∂u ρ2 ∂u
2.2.5. On vérifie que ∇I k · rk = 0 et ∇J k · rk = 0 pour k = 1, 2, 3 :
∂V ∂V −a ∂V
∇I 1 · r1 = 1+ + 0 = 0,
∂ρ ∂u ρ 2 ∂V
∂(u − a/ρ) ∂(u − a/ρ) −a ∂(u − a/ρ)
∇J 1 · r1 = 1+ + 0 = 0,
∂ρ ∂u ρ2 ∂V
∂u ∂π −ρ 2 ∂u ∂u
∇I 2 · r2 = + 0+ 1 = 0,
∂ρ ∂V a 2 ∂u ∂V
∂π ∂π −ρ 2 ∂π ∂π ∂π ∂π
∇J 2 · r2 = + 0+ 1=− + = 0,
∂ρ ∂V a 2 ∂u ∂V ∂V ∂V
∂V ∂V a ∂V
∇I 3 · r3 = 1+ + 0 = 0,
∂ρ ∂u ρ 2 ∂V
∂(u + a/ρ) ∂(u + a/ρ) a ∂(u + a/ρ)
∇J 3 · r3 = 1+ + 0 = 0.
∂ρ ∂u ρ2 ∂V

2.2.6. Étude des trois ondes :

B k = 1 : la condition sur les invariants de Riemann donne
(
I 1 (Y) = I 1 (YL ),
J 1 (Y) = J 1 (YL ),

c’est-à-dire (
V = VL ,
u − ρa = u L − ρa ,
L

ce qui donne (
V = VL ,
π(ρ,V ) = −a(u − u L ) + π(ρ L ,VL );
B k = 2 : la condition sur les invariants de Riemann donne
(
I 2 (Y) = I 2 (YL ),
J 2 (Y) = J 2 (YL ),

c’est-à-dire (
u = uL ,
π(ρ,V ) = π(ρ L ,VL ).
Les conditions de R ANKINE -H UGONIOT s’écrivent
2 2
(ρu) − (ρ L u L ) (ρu + π(ρ,V )) − (ρ L u L + π(ρ L ,VL )) (ρuV ) − (ρ L u L VL )
σ= = =
ρ − ρL ρu − ρ L u L ρV − ρ L VL
et, étant donné que u = u L et π(ρ,V ) = π(ρ L ,VL ), on a σ = u L ;
B k = 3 : la condition sur les invariants de Riemann donne
(
I 3 (Y) = I 3 (YL ),
J 3 (Y) = J 3 (YL ),

c’est-à-dire (
V = VL ,
u + ρa = u L + ρa ,
L

ce qui donne (
V = VL ,
π(ρ,V ) = a(u − u L ) + π(ρ L ,VL ).

© G. Faccanoni 45
5. Systèmes hyperboliques Dernière mise à jour : Lundi 18 mars 2013

2.2.7. Soit un problème de Riemann avec les deux états donnés suivants :

ρL ρR
   

Y L = u L  , YR = u R  .
VL VR

La solution est constituée de quatre états constants



 YL , si x < (u − a/ρ)t ,


Y ,
1 si (u − a/ρ)t < x < u L t ,
Y(x, t ) =


 Y 2, si u L t < x < (u + a/ρ)t ,

YR , si x > (u + a/ρ)t .


t 2-onde

1-onde Y1
Y2 3-onde

YL YR

Pour expliciter cette solution on cherche à définir les six inconnues

ρ1 ρ2

  

Y1 = u 1  , Y 2 = u 2 
V1 V2

à l’aide des invariants de Riemann qu’ici donnent




V1 = VL ,
u 1 − a = u L − a ,




 ρ1 ρL

u = u ,
1 2


π(ρ 1 ,V1 ) = π(ρ 2 ,V2 ),

V 2 = VR ,




u 2 + ρa2 = u R + ρa .



R

Grâce à la première, à la troisième et à la dernière relations on peut tout de suite écrire

ρ1 ρ2

  

Y 1 = u ∗  , Y2 =  u ∗ 
VL VR

et il ne reste à calculer que les trois inconnues ρ 1 , ρ 2 et u ∗ . Le point u ∗ est donné par l’intersection, dans le plan
(π, u), des deux courbes (
π = −au + au L + π(ρ L ,VL ),
π = au − au R + π(ρ R ,VR ),
donc
π(ρ L ,VL ) − π(ρ R ,VR )
µ ¶
1
∗
u = uL + uR +
2 a
et
a a
ρ1 = a , ρ2 = .
u∗ − u L+ ρL u R − u ∗ + ρa
R

2.2.8. Pour vérifier que η est une entropie du système (5.7) avec le flux d’entropie Φ, on montre que

∇Y Φ(Y) = ∇Y η(Y) · B(Y).

46 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 5. Systèmes hyperboliques

Or,
ρ
 
ρu 2
· ¸
Y = u  , Φ(Y) = u +G(ρ, u)
2
V
donc  u3 ∂G

2 + u ∂ρ
2
∇Y Φ(Y) =  3ρu
 
2 +G

∂G
u ∂V
et ρ
1 π(r,V )
Z
η(Y) = ρu 2 + ρ dr
2 ρ0 r2
d’où  u2 Rρ π(r,V ) π(ρ,V ) 
2 + ρ0 r2
dr + ρ
∇Y η(Y) =  uρ .
 
Rρ ∂V π(r,V )
ρ ρ0 r2
dr
Donc

∇Y η(Y) · B(Y) =
 u2 Rρ π(r,V ) π(ρ,V ) T
ρ
 
+ ρ0 dr + ρ
u 0
2 r2
  a2 1 ∂π 
= uρ  ·  ρ 3 u ρ ∂V  =

Rρ
ρ ρ 0 ∂V π(r,V
r2
)
dr 0 0 u
ρ π(r,V ) π(ρ,V
h 2 R i 2

u )
+ dr + ρ u + ua
h 22 R ρ 0 r 2 ρ2 
ρ π(r,V ) π(ρ,V )
i
 u
=  2 + ρ 0 r 2 dr + ρ ρ + u2ρ =

ρ
+ ρ ρ 0 ∂V π(r,V
 h R i 
∂π )
u ∂V r 2 dr u
R ρ π(r,V ) π(ρ,V )
 3 2 
u
2 + u ρ 0 r 2 dr + u ρ + ua
ρ 2
 2
ρu R ρ π(r,V ) 
2 
=  2 + ρ ρ 0 r 2 dr + π(ρ,V ) + u ρ 
Rρ
∂π
u ∂V + ρu ρ 0 ∂V π(r,V
r2
)
dr

qui donne
u3
R ρ π(r,V ) π(ρ,V ) 2
+ ua
 u3 ∂G
 
2 + u ρ 0 r 2 dr + u ρ

2 + u ∂ρ ρ 2

3ρu 2   ρu 2 R ρ π(r,V ) 
2 
 2 + ρ ρ 0 r 2 dr + π(ρ,V ) + u ρ  .
=

2 +G

∂G Rρ
u ∂V ∂π
u ∂V + ρu ρ 0 ∂V π(r,V
r2
)
dr

Identifiant terme à terme on obtient les conditions suivantes

∂G
R ρ π(r,V ) π(ρ,V ) a2

 ∂ρ = Rρ 0 r 2 dr + ρ + ρ 2 ,


ρ
G = ρ ρ 0 π(r,V
r2
)
dr + π(ρ,V ),
 ∂G = ∂π + ρ R ρ ∂V π(r,V ) dr.


∂V ∂V ρ0 r2

2.2.9. On obtient Z ρ π(r,V )

G =ρ dr + π(ρ,V ).
ρ0 r2

2.2.10. Pour k = 2 on a (
u = uL ,
π(ρ,V ) = π(ρ L ,VL ).
ce qui donne

ρ 3
Z ρ π(r,V )
Φ(Y) = u + uρ dr + π(ρ,V ),
2 ρ0 r2
Z ρL
ρL u3 π(r,V )
Φ(YL ) = + uρ L dr + π(ρ L ,VL ),
2 ρ0 r2
Z ρ
1 π(r,V )
η(Y) = ρu 2 + ρ dr,
2 ρ0 r2

© G. Faccanoni 47
5. Systèmes hyperboliques Dernière mise à jour : Lundi 18 mars 2013

1
Z ρL π(r,V )
η(YL ) = ρ L u 2 + ρ L dr,
2 ρ0 r2
σ = u.

Donc on a bien
Z ρ Z ρL
ρ − ρL 3 π(r,V ) π(r,V )
Φ(Y) − Φ(YL ) = u + uρ 2
dr − uρ L dr,
2 ρ0 r ρ0 r2
Z ρ Z ρL =
ρ − ρL 3 π(r,V ) π(r,V )
σ(η(Y) − η(YL )) = u + uρ 2
dr − uρ L dr.
2 ρ0 r ρ0 r2

5.4. Exemple : mélange isotherme

On considère un mélange homogène et isentrope de deux fluides compressibles. On note ρ ≡ ρ(x, t ) > 0 la densité du mé-
lange, u ≡ u(x, t ) la vitesse du mélange, c ≡ c(x, t ) la concentration massique du premier fluide (on a toujours 0 ≤ c ≤ 1 et
1 − c est la concentration massique du deuxième fluide), p ≡ p(ρ, c) > 0 la pression du mélange qu’on supposera vérifier les
conditions
∂p ¯¯
 ¯
 = a 2 > 0,
∂ρ ¯c¯


∂p ¯¯ ∂2 p ¯¯
¯
2 +ρ > 0,


∂ρ ¯c ∂ρ 2 ¯c
où a > 0 désigne la vitesse du son du mélange.
En dimension un d’espace, on modélise ce type d’écoulements bifluide par le système suivant

∂t ρ + ∂x (ρu) = 0,


∂t (ρu) + ∂x (ρu 2 + p) = 0, (5.9)

∂t (ρc) + ∂x (ρuc) = 0.


1. Trouver les vecteurs V : R+ × R → R3 et F(V) : R3 → R3 tels que le système (5.9) s’écrit

∂t V + ∂x F(V) = 0.

2. En se plaçant dans les variables W = (ρ, u, c), le système (5.9) se réécrit, pour des solutions régulières, sous la forma
quasilinéaire
∂t W + A(W)∂x W = 0. (5.10)
Montrer que la matrice A(W) s’écrit
u ρ 0¯
 
a2 1 ∂p ¯ 

 ρ u ρ ∂c ¯ρ  .
0 0 u
3. Calculer les valeurs propres λ1 , λ2 et λ3 . Afin de fixer les notations on ordonne les trois valeurs propres selon λ1 < λ2 < λ3 .
Proposer une base associée de vecteurs propres à droite {r1 , r2 , r3 } de la matrice A. En déduire que le système (5.10) est
strictement hyperbolique.
4. Vérifier que les champs 1 et 3 sont vraiment non linéaires tandis que le champ 2 est linéairement dégénéré.
5. On note I k et J k les deux invariants de Riemann du k-ème champ caractéristique. Montrer qu’un choix possible pour I k et
J k est

I 1 = c, J 1 = u + a ln(ρ),
I 2 = u, J 2 = c,
I 3 = c, J 3 = u − a ln(ρ).

6. Soit une détente de la première famille reliant un état gauche WL = (ρ L , u L , c L ) à un état droit W = (ρ, u, c).
6.1. Montrer que c = c L , u > u L , ρ < ρ L .
6.2. Calculer u en fonction de ρ L , u L et ρ. Plus précisement, montrer que u peut se mettre su la forme

u = u L + g (ρ L , ρ)

en explicitant la fonction g .

48 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 5. Systèmes hyperboliques

7. Soit une onde de la deuxième famille reliant un état gauche WL = (ρ L , u L , c L ) à un état droit W = (ρ, u, c) (il s’agit d’une
discontinuité de contact).
7.1. Montrer que u = u L et p(ρ, c) = p(ρ L , c L ).
7.2. Quelle est la vitesse de déplacement σ de cette discontinuité ?
8. Cherchons maintenant une entropie pour le système (5.9).
8.1. Montrer que la fonction
1
Z ρ p(r, c)
η(W) = ρu 2 + ρ dr
2 ρ0 r2
est une entropie du système avec flux d’entropie
· Z ρ p(r, c)
¸
Φ(W) = u p(ρ, c) + dr +G(ρ, u)
ρ0 r2

(on ne vérifiera pas que η est convexe).

8.2. Expliciter la fonction G.
8.3. Montrer qu’une discontinuité de contact reliant un état gauche WL = (ρ L , u L , c L ) à un état droit W = (ρ, u, c) avec
vitesse σ vérifie
Φ(W) − Φ(WL ) = σ(η(W) − η(WL )),
autrement dit on a conservation de l’entropie.

Correction
1. On a

ρ ρu
   

V = ρu  , F(V) = ρu 2 + p  .

ρc ρuc

2. On développe les dérivées du système (5.9) pour des solutions régulières :


∂ ρ + u∂x ρ + ρ∂x u = 0,
 t


∂p ¯ ∂p ¯
¯ ¯
ρ∂t u + u∂t ρ + 2uρ∂x u + u 2 ∂x ρ + ∂ρ ¯ ∂x ρ + ∂c ¯ρ ∂x c = 0,

 c
ρ∂t c + c∂t ρ + ucρ∂x ρ + ρc∂x u + ρu∂x c = 0.


Il se réécrit alors 
∂t ρ + u∂x ρ + ρ∂x u = 0,

∂p ¯
 ¯
ρ∂t u − u u∂x ρ + ρ∂x u + 2uρ∂x u + u 2 ∂x ρ + a 2 ∂x ρ + ∂c ¯ ∂x c = 0,
£ ¤

 ρ
ρ∂t c − c u∂x ρ + ρ∂x u + uc∂x ρ + ρc∂x u + ρu∂x c = 0,
 £ ¤

d’où 
∂ ρ + u∂x ρ + ρ∂x u = 0,
 t


∂p ¯
2
¯
∂t u + aρ ∂x ρ + u∂x u + ρ1 ∂c ¯ρ ∂x c = 0,


ρ∂t c + ρu∂x c = 0,


et donc finalement on a le système quasilinéaire suivant :

  u ρ 0¯
    
ρ ρ 0
 a2 1 ∂p ¯ 
∂ t u  +  ρ u ρ ∂
∂c ρ  x
¯ u  =  0 .
c 0 0 u c 0

3. On cherche les 3 solutions λk (W) de l’équation det(A(W) − λ(W)Id) = 0, i.e. de l’équation

u 3 − 3u 2 λ + 3uλ2 − λ3 − a 2 u + a 2 λ = 0.

On obtient

λ1 = u − a, λ2 = u, λ3 = u + a.

© G. Faccanoni 49
5. Systèmes hyperboliques Dernière mise à jour : Lundi 18 mars 2013

On peut alors prendre

1 ∂p ¯
 ¯

ρ/a
   
−ρ/a a 2 ∂c ρ
¯
r1 =  1  , r2 =  , r3 =  1  .
 
0
0 1 0

Puisque par hypothèse a > 0, les valeurs propres sont réelles et distinctes donc le système (5.9) est strictement hyperbolique.
4. Pour déterminer la nature des trois champs on calcul ∇λk · rk pour k = 1, 2, 3 :

∂(u − a) −ρ ∂(u − a) ∂(u − a)

∇λ1 · r1 = + 1+ 0=1 1-champ VNL
∂ρ a ∂u ∂c
∂p ¯
¯
∂u ∂c ρ ∂u ∂u
¯
∇λ2 · r2 = + 0+ 1=0 2-champ LD
∂ρ a 2 ∂u ∂c
∂(u + a) ρ ∂(u + a) ∂(u + a)
∇λ3 · r3 = + 1+ 0=1 3-champ VNL
∂ρ a ∂u ∂c

5. On vérifie que ∇I k · rk = 0 et ∇J k · rk = 0 pour k = 1, 2, 3 :

∂c −ρ ∂c ∂c
∇I 1 · r1 = + 1 + 0 = 0,
∂ρ a ∂u ∂c
∂(u + a ln(ρ)) −ρ ∂(u + a ln(ρ)) ∂(u + a ln(ρ))
∇J 1 · r1 = + 1+ 0 = 0,
∂ρ a ∂u ∂c
∂p ¯
¯
∂u ∂c ρ ∂u ∂u
¯
∇I 2 · r2 = + 0+ 1 = 0,
∂ρ a 2 ∂u ∂c
∂p ¯
¯
∂p ∂c ρ ∂p ∂p
¯
∇J 2 · r2 = + 0+ 1 = 0,
∂ρ a 2 ∂u ∂c
∂c ρ ∂c ∂c
∇I 3 · r3 = + 1 + 0 = 0,
∂ρ a ∂u ∂c
∂(u − a ln(ρ)) ρ ∂(u − a ln(ρ)) ∂(u − a ln(ρ))
∇J 3 · r3 = + 1+ 0 = 0.
∂ρ a ∂u ∂c

6. Soit le 1-champ une détente.

6.1. La condition sur les vitesse entraîne λ1 (WL ) < λ1 (W) d’où

u > uL .

La condition sur les invariants de Riemann pour k = 1 donne

(
I 1 (W) = I 1 (WL ),
J 1 (W) = J 1 (WL ),

c’est-à-dire (
c = cL ,
u − a ln(ρ) = u L − a ln(ρ L ),
d’où (
c = cL ,
u − u L = a ln(ρ L /ρ) > 0.
On conclut donc que

c = c L ,

u > uL ,

ρ < ρL .


6.2. La condition sur l’invariant de Riemann J 1 (W) = J 1 (WL ) donne g (ρ L , ρ) = a ln(ρ L /ρ).
7. Considérons le 2-champ qui est une discontinuité de contact.

50 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 5. Systèmes hyperboliques

7.1. La condition sur les invariants de Riemann pour k = 2 donne

(
I 2 (W) = I 2 (WL ),
J 2 (W) = J 2 (WL ),

c’est-à-dire (
u = uL ,
p(ρ, c) = p(ρ L , c L ).

7.2. La condition sur la vitesse σ de la discontinuité de contact découle des conditions de R ANKINE -H UGONIOT et notam-
ment de la relation
σ(ρ L u L − ρu) = (ρ L u L2 + p(ρ L , c L ) − ρu 2 − p(ρ, c))
qui donne
σ(ρ L − ρ)u L = (ρ L − ρ)u L2
c’est-à-dire σ = u L .
8. Entropie.
8.1. Pour vérifier que η est une entropie du système (5.9) avec le flux d’entropie Φ, on montre que

∇W Φ(W) = ∇W η(W) · A(W).

Or,
ρ
 
· Z ρ p(r, c)
¸
W = u  , Φ(W) = u p(ρ, c) + dr +G(ρ, u)
ρ0 r2
c
donc R ρ p(r,c)
ua + u ρ 0 r 2 dr + u ρ + ∂G
 2 p(ρ,c) 
∂ρ
R ρ p(r,c)
p(ρ, c) + ρ 0 r 2 dr + ∂G
 
∇W Φ(W) =  ∂u
 

∂p(ρ,c)
· ¸
∂p(ρ,c) ρ
 
∂c
R
u ∂c + ρ 0 r 2 dr

et ρ
1 p(r, c)
Z
η(W) = ρu 2 + ρ dr
2 ρ0 r2
d’où Rρ
u2 p(r,c) p(r,c)
 
2 + ρ0 r2
dr + ρ
∇W η(W) =  uρ
 
.

Rρ ∂p(ρ,c) 
ρ ρ0
∂c
r2
dr
Donc

∇W η(W) · A(W) =
u2
Rρ p(r,c) p(r,c) T
 
+ dr + u ρ 0¯
 
2 ρ0 r2 ρ
  a2 ∂p ¯
uρ u ρ1 ∂c ¯ 

=  ·
ρ =

R ρ ∂p(ρ,c)
 ρ
ρ ρ 0 r 2 dr
∂c 0 0 u
ρ p(r,c)
h 2 R i 
u p(r,c)
+ dr + ρ u + ua 2
h 22 Rρ 0 r 2
ρ p(r,c)
i 
 u p(r,c) 2 
= 2
 + ρ0 r 2 dr + ρ ρ + u ρ =
∂p(ρ,c)

∂p ¯ ρ
 ¯ 
u ∂c ¯ + uρ ρ 0 r∂c2 dr
R
ρ
 3
u
R ρ p(r,c) p(r,c) 2

2 + u ρ 0 r 2 dr + u ρ + ua
 3 2 R ρ p(r,c)
=  2 ρu + ρ ρ 0 r 2 dr + p(r, c) 

 
∂p(ρ,c)
∂p ¯ ρ
 ¯ 
u ∂c ¯ + uρ ρ 0 r∂c2 dr
R
ρ

qui donne
R ρ p(r,c) R ρ p(r,c)
ua + u ρ 0 r 2 dr + u ρ + ∂G
 2 p(ρ,c)   3
u p(r,c) 2

∂ρ 2 + u ρ 0 r 2 dr + u ρ + ua
R ρ p(r,c)
p(ρ, c) + ρ 0 r 2 dr + ∂G   ρu + ρ ρ p(r,c)
   3 2 R 
∂u = ρ 2 dr + p(ρ, c)  .
 
   2 0 r
∂p(ρ,c) ∂p(ρ,c)
· ¸
∂p(ρ,c) ρ ∂p ¯ ρ
   ¯ 
+ ρ 0 r∂c2 dr u ∂c ¯ + uρ ρ 0 ∂c2 dr
R R
u ∂c ρ r

© G. Faccanoni 51
5. Systèmes hyperboliques Dernière mise à jour : Lundi 18 mars 2013

Identifiant terme par terme on obtient les conditions suivantes

u3
( ∂G
∂ρ = 2 ,
∂G 3u 2 ρ
∂u = 2 .

8.2. On peut prendre G(ρ, u) = 12 u 3 ρ.

8.3. Sur une discontinuité de contact on a (
u = u L = σ,
p(ρ, c) = p(ρ L , c L ).
donc
· Z ρ
p(r, c)
¸
1 3
· Z ρL
p(r, c)
¸
1
Φ(W) − Φ(WL ) = u p(ρ, c) + dr + u ρ − u L p(ρ, c)L + dr − u L3 ρ L =
ρ0 r2 2 ρ0 r2 2
·Z ρ
p(r, c)
¸
1
= uL dr + (u 3 ρ − u L3 ρ L ) =
ρL r2 2
Z ρ Z ρL
p(r, c) p(r, c)
· ¸
1 2 1 2
= uL ρu + ρ dr − ρu − ρ L dr =
2 ρ0 r2 2 L ρ0 r2
= u L (η(W) − η(WL )).

52 © G. Faccanoni
 Troisième partie .

Rappels Numériques

53
 Rappels sur les différences finies
À part dans quelques cas très particulier, il est impossible de calculer explicitement des solutions de modèles issus de la physique.
Il est donc nécessaire d’avoir recours au calcul numérique sur ordinateur pour estimer qualitativement et quantitativement
ces solutions. Le principe de toutes les méthodes de résolution numérique des équations aux dérivées partielles est d’obtenir
des valeurs numériques discrètes (c’est-à-dire en nombre fini) qui approchent (en un sens convenable à préciser) la solution
exacte. Dans ce procédé il faut bien être conscient de deux points fondamentaux : premièrement, on ne calcul pas des solutions
exactes mais approchées ; deuxièmement, on discrétise le problème en représentant des fonctions par un nombre finis de valeurs,
c’est-à-dire que l’on passe du continu au discret.
Il existe de nombreuses méthodes d’approximation numérique des solutions d’équations aux dérivées partielles. Pour simplifier
la présentation, nous nous limiterons dans ce chapitre à la dimension un d’espace.

Principe de la méthode des différences finies

Soit f : R → R une fonction de classe C 1 (R). y
Comme
f (x 0 + h) − f (x 0 )
f 0 (x 0 ) = lim
h→0 h
il est naturel d’introduire les approximations f (x 0 + ∆x)
f (x 0 )
f (x 0 + ∆x) − f (x 0 )
f 0 (x 0 ) ≈ , (5.11)
∆x
f (x 0 ) − f (x 0 − ∆x)
f 0 (x 0 ) ≈ , (5.12)
∆x f (x)
f (x 0 + ∆x) − f (x 0 − ∆x)
f 0 (x 0 ) ≈ . (5.13)
2∆x
f (x 0 − ∆x)
De manière analogue, la dérivée seconde peut être approchée
∆x ∆x
par
f (x 0 + ∆x) − 2 f (x 0 ) + f (x 0 − ∆x) x 0 − ∆x x 0 x 0 + ∆x x
f 00 (x 0 ) ≈ .
(∆x)2

Si on fait un développement de Taylor en x autour du point x 0

f (x 0 ± ∆x) = f (x 0 ) ± ∆x f 0 (x 0 ) + (∆x)2 f 00 (x 0 ) + O((∆x)3 ),

on a
f (x 0 + ∆x) − f (x 0 ) f (x 0 ) + ∆x f 0 (x 0 ) + (∆x)2 f 00 (x 0 ) + O((∆x)3 ) − f (x 0 )
= = f 0 (x 0 ) + O(∆x),
∆x ∆x
f (x 0 ) − f (x 0 − ∆x) f (x 0 ) − f (x 0 ) + ∆x f 0 (x 0 ) − (∆x)2 f 00 (x 0 ) + O((∆x)3 )
= = f 0 (x 0 ) + O(∆x),
∆x ∆x
f (x 0 + ∆x) − f (x 0 − ∆x) f (x 0 ) + ∆x f 0 (x 0 ) + (∆x)2 f 00 (x 0 ) + O((∆x)3 ) − f (x 0 ) + ∆x f 0 (x 0 ) − (∆x)2 f 00 (x 0 )
=
2∆x 2∆x
= f 0 (x 0 ) + O((∆x)2 ),

et pour l’approximation de la dérivée seconde on a

f (x 0 + ∆x) − 2 f (x 0 ) + f (x 0 − ∆x)
(∆x)2
f (x 0 ) + ∆x f 0 (x 0 ) + (∆x)2 f 00 (x 0 ) + O((∆x)3 ) − 2 f (x 0 ) + f (x 0 ) − ∆x f 0 (x 0 ) + (∆x)2 f 00 (x 0 )
=
(∆x)2
= f 00 (x 0 ) + O((∆x)2 ).

Si ∆x est «petit», ces formules sont des «bonnes» approximations.

54
Dernière mise à jour : Lundi 18 mars 2013 5. Systèmes hyperboliques

Généralisation aux équations aux dérivées partielles

Nous nous limitons pour le moment à la dimension un d’espace et considérons une équation aux dérivées partielles F (u) = 0
définie pour (x, t ) ∈ R × R+ avec une condition initiale u(x, 0) = g (x) pour x ∈ R (remarquons que F (u) est une notation pour une
fonction de u et de ses dérivées partielles en tout point). Pour discrétiser le domaine R × R+ , on introduit un pas d’espace ∆x > 0
et un pas de temps ∆t > 0 et on définit les nœuds d’un maillage régulier

(x j , t n ) = ( j ∆x, n∆t ) pour j ∈ Z, n ∈ N.

On note u nj la valeur d’une solution discrète approchée au point (x j , t n ) et u(x, t ) la solution exacte (inconnue). Le principe de la
méthode des différences finies est de remplacer les dérivées par des différences finies en utilisant des formules de Taylor dans
lesquelles on néglige les restes.
Dans tous schéma il y a bien sûr une donnée initiale pour démarrer les itération en n : les valeurs initiales (u 0j ) j ∈Z sont définies
par exemple par u 0j = g (x j ) où g est la donnée initiale de l’équation.

S’il y a un second membre f (x, t ) dans l’équation aux dérivées partielles, alors les schémas se modifient en remplaçant zéro au
second membre par une approximation consistante de f (x, t ) au point (x j , t n ).
Si l’équation est définie sur un domaine borné, par exemple x ∈ [α; β], le maillage spatiale sera restreint à cet intervalle c’est-à-dire
j ∈ {0, 1, . . . , N } avec x 0 = α et x N = β et ∆x = (β − α)/(N + 1). Il faut de plus ajouter des conditions aux limites qui peuvent être de
plusieurs types.
Par exemple, si on a des conditions aux limites de Dirichlet

u(α, t ) = L, u(β, t ) = R, pour t ∈ R+

∗,

elles se traduisent au niveau discret en

u 0n = L, n
uN +1 = R, pour n ∈ N.

Si on a des conditions de Neumann

∂x u(α, t ) = L, ∂x u(β, t ) = R, pour t ∈ R+
∗,

elles se traduisent au niveau discret en

n n
u 1n − u 0n uN − uN −1
= L, = R, pour n ∈ N.
∆x ∆x

Si on a des conditions périodiques

u(x + β, t ) = u(x + α, t ), pour x ∈ [α; β], t ∈ R+
∗

elles se traduisent au niveau discret en

u 0n = u N
n
, pour n ∈ N,
et plus généralement u nj = u N
n
+j
.

Niveaux et stencil Un schéma est dit à m niveaux s’il ne fait intervenir que m indices de temps. Les schémas les plus populaires
sont des schémas à deux ou trois niveaux.
La collection des couples ( j 0 , n 0 ) qui interviennent dans l’équation discrète au point ( j , n) est appelé stencil du schéma (qu’on
peut essayer de traduire par support). En général, plus le stencil est large, plus le schéma est coûteux et difficile à program-
mer.

Consistance et ordre de précision

De manière générale, un schéma aux différences finies est défini, pour tous les indices possibles j ∈ Z et n ∈ N, par la for-
mule à !
n o
n+m
F ∆x,∆t u j +k k − ≤k≤k + = 0
m − ≤m≤m +
− + − +
où les entiers k , k , m et m définissent la largeur du stencil du schéma.
Un des buts de l’analyse numérique est de comparer et de sélectionner les meilleurs schémas suivant des critères de précision, de
coût ou de robustesse.

© G. Faccanoni 55
5. Systèmes hyperboliques Dernière mise à jour : Lundi 18 mars 2013

n + m+
n + m+ − 1

n − m− + 1
n − m−

j − k− + 1

j + k+ − 1
j − k−

j + k+
j
F IGURE 5.1.: Exemple de stencil pour un schéma à m − + m + + 1 niveaux et k − + k + + 1 points.

Erreur de troncature Considérons le schéma aux différence finies F ∆x,∆t ({u n+m j +k
}) pour l’approximation de l’équation aux dérivées
partielles F (u) = 0. Soit u(x, t ) une solution suffisamment régulière de cette équation. On appelle erreur de troncature du schéma
la quantité Ã !
n o
n
τ j ≡ F ∆x,∆t u(x + k∆x, t + m∆t ) k − ≤k≤k + .
m − ≤m≤m +

Concrètement on calcule l’erreur de troncature d’un schéma en remplaçant u n+m

j +k
dans la formule par u(x+k∆x, t +m∆t ).

Consistance Le schéma aux différence finies F ∆x,∆t ({u n+m

j +k
}) est dit consistant avec l’équation aux dérivées partielles F (u) = 0
si l’erreur de troncature du schéma tend vers zéro, uniformément par rapport à (x, t ), lorsque ∆x et ∆t tendent vers zéro
indépendamment.

Ordre de consistance Le schéma aux différence finies F ∆x,∆t ({u n+m

j +k
}) est précis à l’ordre p en espace et à l’ordre q en temps avec
l’équation aux dérivées partielles F (u) = 0 si l’erreur de troncature du schéma tend vers zéro comme O((∆x)p + (∆t )q ) lorsque ∆x
et ∆t tendent vers zéro.

Convergence
La convergence d’un schéma aux différences finies est une propriété naturelle qui assure que, pour des valeurs suffisamment
petites des pas d’espace et de temps, la solution numérique calculée sera proche de la solution exacte du problème de dé-
part.

Convergence Le schéma aux différence finies F ∆x,∆t ({u n+m

j +k
}) utilisé pour la résolution numérique de l’équation aux dérivées
partielles F (u) = 0 est convergent si, pour toute solution u de l’équation F (u) = 0, la suite u nj converge vers u(x j , t n ) avec
(∆x, ∆t ) → (0, 0).

Malheureusement la notion de consistance ne suffit pas à garantir que le schéma soit convergente comme on verra sur des
exemples. Pour introduire un critère (très pratique) qui permet de voir si un schéma donné est convergent nous allons introduire
la notion de stabilité.

Stabilité
Autre les outils qui permettent de comparer les performances des différents schémas, on doit également choisir les pas ∆x
et ∆t de sorte que le schéma correspondant donnera une solution approchée correcte, au sens où une petite perturbation
de la donnée initiale g n’induira par une perturbation trop grande sur la solution calculée au temps final. Cette idée, déjà
rencontrée pour la définition de problème bien posé, est à la base du concept de stabilité pour les schémas aux différences
finies.
Soit u n ≡ (u nj )1≤ j ≤N −1 la solution numérique d’un schéma.

56 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 5. Systèmes hyperboliques

Stabilité Un schéma aux différences finies est dit stable pour la norme k·k s’il existe une constante K > 0 indépendante de ∆x et
∆t (lorsque ces valeurs tendent vers zéro) telle que

ku n k ≤ K ku 0 k pour tout n ≥ 0,

quelle que soit la donnée initiale u 0 . Si cette inégalité n’a lieu que pour des pas ∆x et ∆t astreints à certaines inégalités, on dit que
le schéma est conditionnellement stable.

On définit les normes classiques

à !1
NX
−1 p
ku n kp = ∆x |u nj |p pour 1 ≤ p < +∞,
j =1

ku n k∞ = max |u nj |.
1≤ j ≤N −1

Schéma linéaire Un schéma aux différences finies est dit linéaire si la formule F ∆x,∆t ({u n+m
j +k
}) = 0 qui le définit est linéaire par
rapport à ses arguments u n+m
j +k
.

La stabilité d’un schéma linéaire à deux niveaux est facile à interpréter. En effet, par linéarité tout schéma linéaire à deux niveaux
peut s’écrire sous la forme condensée
Au n = u n+1 ,
où A est une matrice (dite d’itération) et on obtient An u 0 = u n+1 (attention, la notation An désigne ici la puissance n-ème de A)
et par conséquent la stabilité du schéma est équivalente à

kAn u 0 k ≤ K ku 0 k, ∀n ≥ 0, ∀u 0 ∈ RN −1 .

Introduisant la norme matricielle subordonnée kMk = supu∈RN −1 ,u6=0 kMuk

kuk , la stabilité du schéma est équivalente à

kAn k ≤ K ∀n ≥ 0

qui veut dire que la suite des puissances de A est bornée.

Principe du maximum discret - stabilité L ∞ Un schéma aux différences finies vérifie le principe du maximum discret si pour
tout n ≥ 0 et tout 1 ≤ j ≤ N − 1 on a ³ ´ ³ ´
min u 0j ≤ u nj ≤ max u 0j
0≤ j ≤N 0≤ j ≤N
0
quelle que soit la donné initiale u .

Stabilité L 2
La norme L 2 se prête bien à l’étude de la stabilité grâce à l’outil très puissant de l’analyse de Fourier. Supposons désormais que les
conditions aux limites pour l’équation aux dérivées partielles sont des conditions aux limites de périodicité. À chaque vecteur
u n ≡ (u nj )1≤ j ≤N −1 on associe une fonction u n (x), constante par morceaux, périodique, définie sur [α; β] par

u n (x) = u nj si x j −1/2 < x < x j +1/2

avec x j +1/2 = α + ( j + 1/2)∆x pour 0 ≤ j ≤ N , x −1/2 = α et x N +1/2 = β. Ainsi définie, la fonction u n (x) appartient à L 2 ([α; β]), elle
peut donc se décomposer en la somme de Fourier

u n (x) = û (k)e 2i πkx

X n
k∈Z
Rβ
avec û n (k) = α u n (x)e −2i πkx dx et la formule de Plancherel
Z β
|u n (x)|2 dx = |û n (k)|2 .
X
α k∈Z

Remarquons que même si u n (x) est une fonction réelle, les coefficients û n (k) de la série de Fourier sont complexes. Une
propriété importante pour l’étude de stabilité de la transformée de Fourier des fonctions périodiques est la suivantes : si on note
v n (x) = u n (x + ∆x) alors v̂ n (k) = û n (k)e 2i πk∆x .

«Recette» pour un schéma à deux niveaux On injecte dans le schéma un mode de Fourier, on obtient ainsi

u nj = A(k)n e 2i πkx j

© G. Faccanoni 57
5. Systèmes hyperboliques Dernière mise à jour : Lundi 18 mars 2013

et on en déduit la valeur du facteur d’amplification A(k). Rappelons que pour l’instant nous nous sommes limité au cas scalaire,
c’est-à-dire que A(k) est un nombre complexe. On appelle condition de stabilité de Von Neumann l’inégalité

|A(k)| ≤ 1 pour tout mode k ∈ Z.

Si la condition de stabilité de Von Neumann est satisfaite (avec éventuellement des restrictions sur ∆x et ∆t ), alors le schéma est
stable pour la norme L 2 , sinon il est instable.

Dans la plupart des cas, on va trouver des restrictions sur ∆x et ∆t pour obtenir la stabilité au sens L 2 du schéma. Comme ∆x
est initialement fixé, ceci nous oblige à nous donner un pas de temps ∆t petit. Plus cette condition de stabilité est restrictive,
plus le schéma sera couteux à utiliser d’un point de vue du temps de calcul. Au contrario, les schémas inconditionnellement
stables ne nécessitent aucune restriction particulière et donc peuvent être à priori utilisés pour une valeur quelconque de ∆t . Ceci
ne signifie pas pour autant qu’ils seront des «bons» schémas, et notamment que la solution calculée sera proche de la solution
exacte. En effet, un choix trop grand de ∆t donne une mauvaise approximation de la dérivée partielle par rapport au temps. En
pratique, un schéma instable est inutilisable car, même si on part d’une donnée initiale spécialement préparée de manière à ce
qu’aucun des modes de Fourier instables ne soit excité par elle, les inévitables erreurs d’arrondi vont créer des composantes non
nulles (bien que très petites) de la solution sur ces modes instables. La croissance exponentielle de ces modes instables entraîne
qu’après seulement quelque pas en temps ces modes deviennent énormes et polluent complètement le reste de la solution
numérique.

Stabilité + Consistance = Convergence

Théorème de Lax Soit u(x, t ) la solution suffisamment régulière de l’équation aux dérivée partielles F (u) = 0 avec des conditions
aux limites appropriées. Soit u nj la solution numérique discrète obtenue par un schéma aux différences finies avec la donnée
initiale u 0j = g (x j ). Si le schéma est linéaire, à deux niveaux, consistant et stable pour une norme k·k, alors le schéma est convergent
au sens où µ ¶
∀T > 0, lim sup ku nj − u(x j , t n )k = 0.
∆x,∆t →0 t n ≤T

De plus, si le schéma est précis à l’ordre p en espace et à l’ordre q en temps, alors pour tout T > 0 il existe une constante C T > 0
telle que
sup ku nj − u(x j , t n )k ≤ C T (∆x)p + (∆t )q .
¡ ¢
t n ≤T

D’un point de vue pratique ce théorème est très rassurant : si l’on utilise un schéma consistant (ils sont construit pour cela en
général) et que l’on n’observe pas d’oscillations numériques (c’est-à-dire qu’il est stable), alors la solution numérique est proche
de la solution exacte (le schéma converge).

Équation équivalente : diffusion et dispersion

Pour comparer divers schémas consistants et stables (donc convergents) d’un point de vue pratique, un concept pertinent
(quoique formel) est celui d’équation équivalente.

Équation équivalente On appelle équation équivalente d’un schéma l’équation obtenue en ajoutant au modèle étudié la partie
principale (c’est-à-dire le terme d’ordre le plus bas) de l’erreur de troncature.

Tous les schéma qu’on va voir sont consistants. Cependant, si on ajoute à l’équation la partie principale de l’erreur de troncature
d’un schéma, alors ces schémas non seulement sont encore consistants avec cette nouvelle équation «équivalente», mais sont
même strictement plus précis pour cette équation équivalente. En d’autre termes, les schémas sont «plus consistants» avec
l’équation équivalente qu’avec l’équation qu’on veut approcher.
Cette équation va nous donner des renseignements précieux sur le comportement numérique du schéma. Le coefficient de
diffusion (c’est-à-dire le coefficient de la dérivée seconde) de l’équation équivalente est appelé diffusion numérique. S’il est
grand on dit que le schéma est diffusif (ou dissipatif). Le comportement typique d’un schéma diffusif est sa tendance à étaler
artificiellement les données initiales au cours du temps. Si le schéma est précis d’ordre 2 alors l’équation équivalente ne contient
pas de terme de diffusion mais un terme du troisième ordre, dit dispersif. Le comportement typique d’un schéma dispersif est
qu’il produit des oscillation lorsque la solution est discontinue. En effet, le terme dispersif modifie la vitesse de propagation des
modes de Fourier de la solution (particulièrement des modes de fréquence élevée), alors qu’un terme diffusif ne fait qu’atténuer
son amplitude.

58 © G. Faccanoni
 6. Résolution numérique de l’équation de transport
On considère le problème de transport en une dimension d’espace dans le domaine borné [0; L] avec une vitesse c constante non
nulle et des conditions aux limites de périodicité

∂t u(x, t ) + c∂x u(x, t ) = 0 pour (x, t ) ∈ [0; L] × [0; T ],



u(x + L, t ) = u(x, t ) pour (x, t ) ∈ [0; L] × [0; T ],

u(x, 0) = g (x) pour x ∈ [0; L].


On souhaite calculer la valeur de la solution u en un ensemble discret de points en espace et en temps. Plus précisément, en
fixant un pas d’espace ∆x = L/N > 0 (N entier positif) et un pas de temps ∆t > 0, on cherche à calculer u nj ≈ u( j ∆x, n∆t ) la valeur
d’une solution discrète approchée au point (x j , t n ). On sait que u(x, t ) = g (x − c t ) est la solution exacte. Les conditions aux limites
de périodicité conduisent aux égalités u 1n = u N
n
+1 pour tout n ≥ 0, par conséquent l’inconnue discrète à chaque pas de temps est
n n
un vecteur u = (u j )1≤ j ≤N .

Notre stratégie consiste à remplacer des opérateurs différentiels par des quotients aux différences finies. En utilisant différentes
façon d’évaluer les dérivées partielles, beaucoup de choix de schémas sont possibles. Nous voulons en étudier ici quelques
uns.
∆t
Soit α := c ∆x . On considère les schémas aux différences finies suivants :
¶ le schéma décentré à gauche

u n+1
j
− u nj u nj − u nj−1
+c =0 i.e. u n+1
j = u nj − α(u nj − u nj−1 )
∆t ∆x

· le schéma décentré à droite

u n+1
j
− u nj u nj+1 − u nj
+c =0 i.e. u n+1
j = u nj − α(u nj+1 − u nj )
∆t ∆x

¸ le schéma centré
u n+1
j
− u nj u nj+1 − u nj−1 u nj+1 − u nj−1
+c =0 i.e. u n+1
j = u nj − α
∆t 2∆x 2
¹ le schéma upwind (décentré amont)

u n+1
j
− u nj c + |c| n c − |c| n
+ (u j − u nj−1 ) + (u j +1 − u nj ) = 0
∆t 2∆x 2∆x
i.e.
¶ ( n
α + |α| n α − |α| n u j − α(u nj − u nj−1 ) si α > 0,
µ
u n+1
j = u nj − n
(u j − u j −1 ) + n
(u j +1 − u j ) =
2 2 u nj − α(u nj+1 − u nj ) si α < 0.

º le schéma de Lax-Friedrichs
u nj+1 +u nj−1
u n+1
j
− 2 u nj+1 − u nj−1 1−α n 1+α n
+c =0 i.e. u n+1
j = u j +1 + u j −1
∆t 2∆x 2 2

» le schéma de Lax-Wendroff
u n+1
j
− u nj u nj+1 − u nj−1 n n n
c 2 ∆t u j +1 − 2u j + u j −1
+c − =0
∆t 2∆x 2 (∆x)2
i.e.
u nj+1 − u nj−1 u nj+1 − 2u nj + u nj−1
u n+1
j = u nj − α + α2
2 2

59
6. Résolution numérique de l’équation de transport Dernière mise à jour : Lundi 18 mars 2013

¼ le schéma de Beam-Warmimg (pour c > 0)

³ ´ α−1 3−α
u n+1
j = u nj − α g (u nj−1 , u nj ) − g (u nj−2 , u nj−1 ) avec g (A, B ) = A+ B
2 2
i.e.
(1 − α) n
µ ¶
u n+1
j = u n
j − α (u n
j − u n
j −1 ) + (u j − 2u n
j −1 + u n
j −2 )
2
ou encore
α(α − 1) n (α − 1)(α − 2) n
u n+1
j = u j −2 + α(2 − α)u nj−1 + uj
2 2
½ le schéma de Fromm (pour c > 0)

u n+1
j
− u nj g (u nj−1 , u nj , u nj+1 ) − g (u nj−2 , u nj−1 , u nj ) α−1 1−α
+c = 0 avec g (A, B,C ) = A +B + C
∆t ∆x 4 4
i.e.
α(α − 1) n α(5 − α) n (1 − α)(α + 4) n α(α − 1) n
u n+1
j = u j −2 + u j −1 + uj + u j +1
4 4 4 4
¾ le schéma anti-diffusif de Després-Lagoutière (pour c > 0)

u n+1
j
− u nj g (u nj−1 , u nj , u nj+1 ) − g (u nj−2 , u nj−1 , u nj )
+c =0
∆t ∆x
i.e.
u n+1
j = u nj − α(g (u nj−1 , u nj , u nj+1 ) − g (u nj−2 , u nj−1 , u nj ))
avec 
 A,
 si R ≤ A,
g (L,C , R) = B, si R ≥ B,

R, sinon,


où
C − max(L,C ) C − min(L,C )
A = max(L,C ) + , B = min(L,C ) + .
α α

6.1. Exemples d’études de stabilité et de consistance

Étude du schéma ¶ (décentré à gauche)

u n+1
j
− u nj u nj − u nj−1
+c =0 i.e. u n+1
j = u nj − α(u nj − u nj−1 )
∆t ∆x
Stencil

n +1

j −1 j

Stabilité L 2 . On utilise l’analyse de Fourier : pour k ∈ Z, le coefficient de Fourier û n (k) de la solution du schéma vérifie
h i
û n+1 (k) = 1 − α + αe −i 2πk∆x û n (k).

En notant ξ ≡ 2πk∆x, on a
h i
û n+1 (k) = 1 − α + αe −i ξ û n (k) =
= 1 − α + α cos(−ξ) + i sin(−ξ) û n (k) =
£ ¡ ¢¤

60 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 6. Résolution numérique de l’équation de transport

= [1 − α + α cos(ξ) − i α sin(ξ)] û n (k).

Après simplification on obtient

|û n+1 (k)|2 = |A(k)|2 |û n (k)|2
avec

|A(k)|2 ≡ [1 − α + α cos(ξ)]2 + α2 sin2 ξ =

= 1 + 2α(α − 1)(1 − cos(ξ)).

On a
|A(k)| ≤ 1 ∀k ∈ Z ⇐⇒ 2α(α − 1)(1 − cos(ξ)) ≥ 0 ∀ξ ∈ R ⇐⇒ α(α − 1) ≥ 0.
Pour 0 ≤ α ≤ 1 on a |A(k)| ≤ 1 pour toute fréquence k ∈ Z, ce qui prouve que le schéma est stable en norme L 2 sous la
condition CFL 0 ≤ α ≤ 1.
Ordre de consistance. On remplace u im par u(x i , t m ) où u est une fonction régulière, i = j − 1, j et m = n, n + 1. On définit
l’erreur de troncature par

u(x j , t n+1 ) − u(x j , t n ) u(x j , t n ) − u(x j −1 , t n )

τnj ≡ +c .
∆t ∆x
On fait un développement de Taylor en x autour du point x j et en t autour du point t n et, comme u est solution de
l’équation ∂t u = −c∂x u, on a

∂u (∆x)2 ∂2 u
u(x j −1 , t n ) = u(x j , t n ) − ∆x (x j , t n ) + (x j , t n ) + O((∆x)3 ),
∂x 2 ∂x 2
∂u (∆t )2 ∂2 u
u(x j , t n+1 ) = u(x j , t n ) + ∆t (x j , t n ) + (x j , t n ) + O((∆t )3 )
∂t 2 ∂t 2
∂u (∆t )2 ∂ ∂u
µ ¶ µ ¶
n n n
= u(x j , t ) + ∆t −c (x j , t ) + −c (x j , t ) + O((∆t )3 )
∂x 2 ∂t ∂x
∂u (∆t )2 ∂ ∂u
µ ¶
n n n
= u(x j , t ) − c∆t (x j , t ) − c (x j , t ) + O((∆t )3 )
∂x 2 ∂x ∂t
∂u (∆t )2 ∂2 u
= u(x j , t n ) − c∆t (x j , t n ) + c 2 (x j , t n ) + O((∆t )3 ).
∂x 2 ∂x 2
Par conséquent l’erreur de troncature se réécrit

u(x j , t n+1 ) − u(x j , t n ) u(x j , t n ) − u(x j −1 , t n )

τnj ≡ +c
∆t ∆x
¶ 2
∆t ∆x ∂ u
µ
= c2 −c (x j , t n ) + O((∆x)2 + (∆t )2 )
2 2 ∂x 2
= O((∆x) + (∆t )).

Le schéma est donc d’ordre 1 en temps et en espace.

On aurait pu s’arrêter un ordre avant dans les développements de Taylor, mais vous voyez bien que même si on va plus loin on obtient le bon résultat !
Cependant, dans la consigne de l’exercice je vous ai suggéré l’ordre pour éviter de faire des calculs inutiles. Dans les exemples qui suivent on va s’arrêter au
minimum nécessaire.

Exemple : étude du schéma ¸ (centré)

u n+1
j
− u nj u nj+1 − u nj−1 α n
+c =0 i.e. u n+1
j = u nj − (u − u nj−1 )
∆t 2∆x 2 j +1
Stencil

n +1

j −1 j j +1

© G. Faccanoni 61
6. Résolution numérique de l’équation de transport Dernière mise à jour : Lundi 18 mars 2013

Stabilité L 2 . On utilise l’analyse de Fourier : pour k ∈ Z, le coefficient de Fourier û n (k) de la solution du schéma vérifie
h α ³ i 2πk∆x ´i
û n+1 (k) = 1 − e − e −i 2πk∆x û n (k).
2
En notant ξ ≡ 2πk∆x, on a
h α ³ iξ ´i
û n+1 (k) = 1 − e − e −i ξ û n (k) =
2
h α¡ ¢i
= 1− cos(ξ) + i sin(ξ) − cos(−ξ) − i sin(−xi ) û n (k) =
2
= [1 − i α sin(ξ)] û n (k).

Après simplification on obtient

|û n+1 (k)|2 = |A(k)|2 |û n (k)|2
avec

|A(k)|2 ≡ 1 + α2 sin2 (ξ).

On a
|A(k)| ≤ 1 ∀k ∈ Z ⇐⇒ α2 sin2 (ξ) ≤ 0 ∀ξ ∈ R ⇐⇒ 6 ∃α ∈ R.
Ce qui prouve que le schéma est inconditionnellement instable en norme L 2 .
Ordre de consistance. On remplace u im par u(x i , t m ) où u est une fonction régulière, i = j − 1, j , j + 1 et m = n, n + 1. On définit
l’erreur de troncature par

u(x j , t n+1 ) − u(x j , t n ) u(x j +1 , t n ) − u(x j −1 , t n )

τnj ≡ +c .
∆t 2∆x
On fait un développement de Taylor en x autour du point x j et en t autour du point t n et, comme u est solution de
l’équation ∂t u = −c∂x u, on a

∂u (∆x)2 ∂2 u
u(x j −1 , t n ) = u(x j , t n ) − ∆x (x j , t n ) + (x j , t n ) + O((∆x)3 ),
∂x 2 ∂x 2
∂u (∆x)2 ∂2 u
u(x j +1 , t n ) = u(x j , t n ) + ∆x (x j , t n ) + (x j , t n ) + O((∆x)3 ),
∂x 2 ∂x 2
∂u
u(x j , t n+1 ) = u(x j , t n ) + ∆t (x j , t n ) + O((∆t )2 )
∂t
∂u
µ ¶
= u(x j , t n ) + ∆t −c (x j , t n ) + O((∆t )2 )
∂x
∂u
= u(x j , t n ) − c∆t (x j , t n ) + O((∆t )2 ).
∂x
Par conséquent l’erreur de troncature se réécrit

u(x j , t n+1 ) − u(x j , t n ) u(x j +1 , t n ) − u(x j −1 , t n )

τnj ≡ +c = O((∆x)2 + (∆t )).
∆t 2∆x
Le schéma est donc d’ordre 1 en temps et 2 en espace.

Exemple : étude du schéma º (Lax-Friedrichs)

u nj+1 +u nj−1
u n+1
j
− 2 u nj+1 − u nj−1 1−α n 1+α n
+c =0 i.e. u n+1
j = u j +1 + u j −1
∆t 2∆x 2 2
Stencil

n +1

j −1 j j +1

62 © G. Faccanoni
Dernière mise à jour : Lundi 18 mars 2013 6. Résolution numérique de l’équation de transport

Stabilité L 2 . On utilise l’analyse de Fourier : pour k ∈ Z, le coefficient de Fourier û n (k) de la solution du schéma vérifie
1 − α i 2πk∆x 1 + α −i 2πk∆x n
· ¸
n+1
û (k) = e − e û (k).
2 2

En notant ξ ≡ 2πk∆x, on a

1 − α i ξ 1 + α −i ξ n
· ¸
û n+1 (k) = e + e û (k) = [cos(ξ) − i α sin(ξ)] û n (k).
2 2

Après simplification on obtient

|û n+1 (k)|2 = |A(k)|2 |û n (k)|2
avec

|A(k)|2 ≡ cos2 (ξ) + α2 sin2 (ξ).

On a
|A(k)| ≤ 1 ∀k ∈ Z ⇐⇒ cos2 (ξ) + α2 sin2 (ξ) ≥ 1 ∀ξ ∈ R ⇐⇒ α2 ≤ 1.
Ce qui prouve que le schéma est stable en norme L 2 sous la condition CFL |α| ≤ 1.
Ordre de consistance. On remplace u im par u(x i , t m ) où u est une fonction régulière, i = j − 1, j , j + 1 et m = n, n + 1. On définit
l’erreur de troncature par
u(x j +1 ,t n )+u(x j −1 ,t n )
u(x j , t n+1 ) − 2 u(x j +1 , t n ) − u(x j −1 , t n )
τnj ≡ +c .
∆t 2∆x
On fait un développement de Taylor en x autour du point x j et en t autour du point t n et, comme u est solution de
l’équation ∂t u = −c∂x u, on a

∂u
u(x j −1 , t n ) = u(x j , t n ) − ∆x (x j , t n ) + O((∆x)2 ),
∂x
∂u
u(x j +1 , t n ) = u(x j , t n ) + ∆x (x j , t n ) + O((∆x)2 ),
∂x
∂u
u(x j , t n+1 ) = u(x j , t n ) + ∆t (x j , t n ) + O((∆t )2 )
∂t
∂u
µ ¶
= u(x j , t n ) + ∆t −c (x j , t n ) + O((∆t )2 )
∂x
∂u
= u(x j , t n ) − c∆t (x j , t n ) + O((∆t )2 ).
∂x
Par conséquent l’erreur de troncature se réécrit

u(x j , t n+1 ) − u(x j , t n ) u(x j +1 , t n ) − u(x j −1 , t n ) (∆x)2

µ ¶
τnj ≡ +c = O (∆x) + (∆t ) + .
∆t 2∆x ∆t
(∆x)2
³ ´
Étant donné que, sous la condition CFL calculée précédemment, O ∆t = O(∆x), le schéma est donc d’ordre 1 en temps
et en espace.

Exemple : étude du schéma » (Lax-Wendroff)

Pour j = 0, . . . , N − 1 et n ∈ N, le schéma de Lax-Wendroff s’écrit

u n+1
j
− u nj u nj+1 − u nj−1 u nj+1 − 2u nj + u nj−1 u nj+1 − u nj−1 u nj+1 − 2u nj + u nj−1
+c − c 2 ∆t =0 i.e. u n+1
j = u nj − α + α2
∆t 2∆x 2(∆x)2 2 2

Stencil

n +1

j −1 j j +1

© G. Faccanoni 63
6. Résolution numérique de l’équation de transport Dernière mise à jour : Lundi 18 mars 2013

Stabilité L 2 . On utilise l’analyse de Fourier : pour k ∈ Z, le coefficient de Fourier û n (k) de la solution du schéma vérifie

−α + α2 i 2πk∆x α + α2 −i 2πk∆x n
· ¸
û n+1 (k) = e + 1 − 2α2 + e û (k).
2 2

En notant ξ ≡ 2πk∆x, on a

−α + α2 i ξ α + α2 −i ξ n
· ¸
n+1 2
e +1−α + û (k) = 1 − α2 (1 − cos(ξ)) − i α sin(ξ) û n (k).
£ ¤
û (k) = e
2 2

Après simplification on obtient

|û n+1 (k)|2 = |A(k)|2 |û n (k)|2
avec

|A(k)|2 ≡ α2 (α2 − 1)(cos(ξ) − 1)2 + 1.

On a
|A(k)| ≤ 1 ∀k ∈ Z ⇐⇒ α2 ≤ 1.
Ce qui prouve que le schéma est stable en norme L 2 sous la condition CFL |α| ≤ 1.
Ordre de consistance. On remplace u im par u(x i , t m ) où u est une fonction régulière, i = j − 1, j , j + 1 et m = n, n + 1. On définit
l’erreur de troncature par

u(x j , t n+1 ) − u(x j , t n ) u(x j +1 , t n ) − u(x j −1 , t n ) u(x j +1 , t n ) − 2u(x j , t n ) + u(x j −1 , t n )

τnj ≡ +c − c 2 ∆t = 0.
∆t 2∆x 2(∆x)2

On fait un développement de Taylor en x autour du point x j et en t autour du point t n et, comme u est solution de
l’équation ∂t u = −c∂x u, on a

∂u (∆x)2 ∂2 u
u(x j −1 , t n ) = u(x j , t n ) − ∆x (x j , t n ) + (x j , t n ) + O((∆x)3 ),
∂x 2 ∂x 2
∂u (∆x)2 ∂2 u
u(x j +1 , t n ) = u(x j , t n ) + ∆x (x j , t n ) + (x j , t n ) + O((∆x)3 ),
∂x 2 ∂x 2
∂u (∆t )2 ∂2 u
u(x j , t n+1 ) = u(x j , t n ) + ∆t (x j , t n ) + (x j , t n ) + O((∆t )3 )
∂t 2 ∂t 2
∂u (∆t )2 2 ∂2 u
µ ¶ µ ¶
= u(x j , t n ) + ∆t −c (x j , t n ) + c (x j , t n
) + O((∆t )3 )
∂x 2 ∂x 2
∂u (∆t )2 ∂2 u
= u(x j , t n ) − c∆t (x j , t n ) + c 2 (x j , t n ) + O((∆t )3 ).
∂x 2 ∂x 2
Par conséquent l’erreur de troncature se réécrit

u(x j , t n+1 ) − u(x j , t n ) u(x j +1 , t n ) − u(x j −1 , t n ) u(x j +1 , t n ) − 2u(x j , t n ) + u(x j −1 , t n )

τnj ≡ +c − c 2 ∆t = O((∆t )2 ) + O((∆x)2 ).
∆t 2∆x 2(∆x)2
¡ ∆t ¢
car, sous la condition CFL calculée précédemment, O c ∆x = O(1), le schéma est donc d’ordre 2 en temps et en espace.

64 © G. Faccanoni
 7. Résolution numérique d’une équation non-linéaire
On cherche

u : R × R+ → R
(x, t ) 7→ u(x, t )

solution faible entropique du problème (

∂t u + ∂x q(u) = 0, x ∈ R, t > 0,
u(x, 0) = g (x), x ∈ R,
avec le flux

q: R→R
u 7→ q(u)

∆t
Soit α := ∆x et cfl une constante réelle positive et posons

1
Z x j +1/2
u nj ≈ u(x, n∆t )dx
∆x x j −1/2

la solution approchée dans la mailles de centre ( j ∆x, n∆t ).

¶ On considère tout d’abord le schéma upwind, qui se base sur la forme non-conservative

∂t u + q 0 (u)∂x u = 0,

tel qu’on l’a écrit pour l’équation de transport :

u n+1
j
− u nj q 0 (u nj ) + |q 0 (u nj )| q 0 (u nj ) − |q 0 (u nj )|
+ (u nj − u nj−1 ) + (u nj+1 − u nj ) = 0
∆t 2∆x 2∆x
i.e.
q 0 (u nj ) + |q 0 (u nj )| q 0 (u nj ) − |q 0 (u nj )|
( n
u j − αq 0 (u nj )(u nj − u nj−1 ) si q 0 (u nj ) > 0,
à !
u n+1
j = u nj −α (u nj − u nj−1 ) + (u nj+1 − u nj ) =
2 2 u nj − αq 0 (u nj )(u nj+1 − u nj ) si q 0 (u nj ) < 0.

On considère ensuite des schémas de type Volumes Finis qui se basent sur la forme conservative de l’équation non-linéaire

∂t u + ∂x q(u) = 0.

Principe : on intègre l’EDP entre x j −1/2 ≡ ( j − 1/2)∆x et x j +1/2 ≡ ( j + 1/2)∆x pour tout t > 0 :
ÃZ ! ÃZ !
x j +1/2 x j +1/2
∂t u(x, t )d x + ∂x q(u(x, t ))d x = 0.
x j −1/2 x j −1/2

On pose
1
Z x j +1/2
u j (t ) = u(x, t )d x
∆x x j −1/2

ce qui donne
q(u j +1/2 (t )) − q(u j −1/2 (t ))
µ ¶
∂t (u j (t )) + ∂x = 0.
∆x
Idée : approcher q(u j +1/2 (t )) ≈ g (u j (t ), u j +1 (t )) où g est le flux numérique du schéma qu’il faut définir :
∆x
u n+1
j = u nj − α(g (u nj−1 , u nj ) − g (u nj , u nj+1 )) avec ∆t ≤ cfl ,
sup|q 0 (u nj )|
j

65
7. Résolution numérique d’une équation non-linéaire Dernière mise à jour : Lundi 18 mars 2013

Les différents schémas se caractérisent par la définition du flux numérique g qui est une approximation du flux q :
· le schéma de Lax-Friedrichs :
q(L) + q(R) + L−R
α
g (L, R) =
2
¸ le schéma de Lax-Wendroff : ¡ L+R ¢
q(L) + q(R) − α(q(R) − q(L))q 0 2
g (L, R) =
2
¹ le schéma de Godunov : si on dénote w( xt , L, R) la solution exacte du problème de Riemann à deux états donnés par L
(gauche) et R (droite), le flux numérique s’écrit

g (L, R) = q(w(0, L, R))

º le schéma de Murman-Roe :
( (
q(L) si a(L, R) ∈ {0, 1}, signe(q(L) − q(R)) · signe(L − R) si L 6= R,
g (L, R) = où a(L, R) =
q(R) si a(L, R) = −1 signe(q 0 (L)) sinon.

66 © G. Faccanoni
 8. Systèmes hyperboliques
Dans ce dernier chapitre on s’intéresse à la résolution de problèmes de Riemann associés aux systèmes hyperboliques non-
linéaires en une dimension d’espace. Plus précisément, on cherche une fonction

W : R+ × R → R p
(t , x) 7→ W(t , x)

qui vérifie (au sens faible) le système d’EDPs

∂t W + ∂x F(W) x ∈ R, t > 0,


 = 0,
(
WL si x < 0,
W(0, x) =

WR si x > 0,

avec le flux

F : Rp → Rp
W 7→ F(W)

Pour des solutions régulières on peut réécrire ce système sous la forme quasi-linéaire

∂t W + B(W)∂x W = 0

avec B(W) la matrice jacobienne du flux F et on note λ1 (W) < λ1 (W) < · · · < λp (W) ses valeurs propres.
Pour j ∈ Z et n ∈ N on pose

1
Z x j +1/2
Wnj ≈ W(x, t n )d x et G(Wnj , Wnj+1 ) ≈ F(W(x j +1/2 , t n )).
∆x x j −1/2

On considère des schémas Volumes Finis :

∆t n h i
Wn+1
j = Wnj − G(Wnj , Wnj+1 ) − G(Wnj−1 , Wnj ) ,
∆x
à savoir
∆x
F(a) + F(b) − ∆t n (b − a)
le schéma de Lax-Friedrichs G(a, b) =
2
F(a) + F(b) b−a
le schéma de Rusanov G(a, b) = − max {|λk (a)|, |λk (b)|}
2 k=1,...,p 2

On associe à ces schémas la condition de stabilité de Courant-Friedrichs-Lewy

∆x
∆t n ≤ n o.
2 max |λ1 (W j )|, · · · , |λp (Wnj )|
n
j ∈Z

67