Introduction aux équations aux

dérivées partielles et à leur

approximation numérique

Anne-Sophie Bonnet-Ben Dhia, Sonia Fliss, Patrick Joly, Philippe Moireau

April 24, 2018

2
Table des matières

1 Introduction aux problèmes hyperboliques linéaires 7

1.1 Exemples de modèles hyperboliques . . . . . . . . . . . . . . . . . . . . 7
1.1.1 Equation de transport linéaire . . . . . . . . . . . . . . . . . . . . 7
1.1.2 L’équation des ondes et ses variantes . . . . . . . . . . . . . . . . 9
1.2 La méthode des caractéristiques pour l’équation de transport . . . . . 12
1.2.1 Le problème de Cauchy pour l’équation de transport linéaire . . 12
1.2.2 Résolution de l’équation de transport à coefficients constants . 13
1.2.3 Résolution de l’équation de transport à coefficients variables . . 17
1.3 Systèmes hyperboliques linéaires . . . . . . . . . . . . . . . . . . . . . . 23
1.3.1 Systèmes hyperboliques linéaires . . . . . . . . . . . . . . . . . . 23
1.3.2 L’équation des ondes en dimension 1 . . . . . . . . . . . . . . . 25

2 Problèmes hyperboliques non linéaires : 29

2.1 Lois de conservation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.1.2 Construction d’une équation de conservation . . . . . . . . . . 30
2.1.3 Exemples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.1.4 Ajout d’un terme de viscosité . . . . . . . . . . . . . . . . . . . . 33
2.2 Solutions classiques : méthode des caractéristiques . . . . . . . . . . . 34
2.2.1 Solutions classiques au problème de Cauchy . . . . . . . . . . . 34
2.2.2 Solutions classiques locales et naissance d’un choc . . . . . . . 39
2.3 Solutions faibles - Solutions entropiques . . . . . . . . . . . . . . . . . . 42
2.3.1 Solutions faibles . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
2.3.2 Solutions C 1 par morceaux. Condition de choc . . . . . . . . . . 44
2.3.3 Solutions entropiques . . . . . . . . . . . . . . . . . . . . . . . . 49
2.3.4 Existence et unicité de la solution entropique . . . . . . . . . . . 57
2.4 Le problème de Riemann . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
2.4.1 Présentation du problème . . . . . . . . . . . . . . . . . . . . . . 58
2.4.2 Une solution auto-semblable . . . . . . . . . . . . . . . . . . . . 59
2.4.3 Le cas d’une onde de détente . . . . . . . . . . . . . . . . . . . . 60
2.4.4 Le cas d’une onde de choc . . . . . . . . . . . . . . . . . . . . . . 61
2.4.5 Solution du problème de Riemann à 2 états . . . . . . . . . . . . 61
2.4.6 Le problème de Riemann à trois états . . . . . . . . . . . . . . . 62

3
4 TABLE DES MATIÈRES

3 Approximation par différences finies 67

3.1 Présentation de la méthode de différences finies . . . . . . . . . . . . . 67
3.1.1 Approximation de quelques opérateurs différentiels . . . . . . . 67
3.1.2 La démarche : Etude de la consistance et la stabilité . . . . . . . 71
3.1.3 Application aux équations hyperboliques linéaires . . . . . . . . 72
3.2 Le schéma explicite centré pour l’équation d’advection . . . . . . . . . 74
3.2.1 Consistance et ordre du schéma . . . . . . . . . . . . . . . . . . . 76
3.2.2 Propagation numérique et condition nécessaire de convergence 78
3.2.3 Analyse de la stabilité d’un schéma : principes généraux . . . . . 79
3.2.4 Analyse de stabilité L2 du schéma explicite centré : la méthode
de Fourier-Von Neumann . . . . . . . . . . . . . . . . . . . . . . 82
3.3 Le schéma de Lax-Friedrichs pour l’équation d’advection . . . . . . . . 87
3.3.1 Consistance et ordre du schéma . . . . . . . . . . . . . . . . . . . 88
3.3.2 Vitesse de propagation numérique . . . . . . . . . . . . . . . . . 89
3.3.3 Analyse de stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . 90
3.3.4 Convergence du schéma . . . . . . . . . . . . . . . . . . . . . . . 92
3.4 Le schéma saute-mouton pour l’équation des ondes . . . . . . . . . . . 97
3.4.1 Présentation du schéma . . . . . . . . . . . . . . . . . . . . . . . 97
3.4.2 Le schéma de démarrage : choix des conditions initiales appro-
chées . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
3.4.3 Ordre du schéma : consistance et erreur de troncature . . . . . . 99
3.4.4 Analyse de stabilité L2 par techniques énergétiques . . . . . . . 99

4 Discrétisation des équations hyperboliques non linéaires 1D 103

4.1 Présentation du problème . . . . . . . . . . . . . . . . . . . . . . . . . . 103
4.2 Propriétés des schémas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
4.2.1 Schémas sous forme conservative . . . . . . . . . . . . . . . . . 104
4.2.2 Consistance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
4.2.3 Linéarisation du schéma - Stabilité L2 . . . . . . . . . . . . . . . 106
4.2.4 Vers la convergence : le théorème de Lax-Wendroff . . . . . . . . 108
4.2.5 Schémas entropiques . . . . . . . . . . . . . . . . . . . . . . . . . 110
4.2.6 Schémas monotones . . . . . . . . . . . . . . . . . . . . . . . . . 112
4.3 Exemples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
4.3.1 Quelques schémas d’ordre 1 . . . . . . . . . . . . . . . . . . . . . 117
4.3.2 Un exemple de schéma d’ordre 2 : le schéma de Lax-Wendroff . 121
4.3.3 Un exemple fondamental - Le schéma de Godounov . . . . . . . 122
Préambule

Les équations aux dérivées partielles sont aux fonctions de plusieurs variables ce
que les équations différentielles sont aux fonctions d’une variable réelle. Elles ex-
priment, sous forme d’égalités, des relations que doivent satisfaire les dérivées par-
tielles d’une certaine fonction inconnue u de plusieurs variables afin de décrire un
phénomène physique, satisfaire une propriété prescrite, ...etc ... Quand on s’intéresse
à une application donnée, on est très vite amené à faire la dictinction entre les EDPs
stationnaires (où toutes les variables jouent un rôle équivalent) et les EDPs d’évolu-
tion (où une des variables joue un rôle privilégié : le temps). On rencontre de telles
équations dès qu’on s’intéresse à des questions de modélisation : en physique, en
électromagnétisme, en mécanique du solide et des fluides bien sûr, mais aussi en bio-
logie, en chimie, en économie, en finance... On y est naturellement confronté quand
on s’intéresse à des problèmes de calcul des variations ou plus généralement d’op-
timisation. Savoir manipuler des équations aux dérivées partielles fait de nos jours
partie du quotidien de l’ingénieur.

L’objectif de ce cours est d’introduire à l’étude et à la résolution approchée sur ordi-

nateur des équations aux dérivées partielles les plus simples qu’on rencontre dans les
applications. Ce cours est un cours de mathématiques appliquées qui prétend à un
certain niveau de rigueur et de généralité. Toutefois, compte tenu de la complexité
technique du sujet, nous serons parfois amenés à sacrifier quelque peu cette rigueur
et cette généralité : la théorie des équations aux dérivées partielles et de leur approxi-
mation numérique est une science difficile, au coeur de la recherche moderne en
mathématiques, dont seuls les rudiments sont accessibles en première année d’Ecole
d’ingénieurs. C’est pourquoi nous avons fait le choix pédagogique d’aborder le sujet
par le « petit bout de la lorgnette » en choisissant de présenter, plutôt qu’une grande
théorie, des exemples et des problèmes modèles simples pour lesquels on peut aller
assez loin dans la description des propriétés des solutions, souvent par le biais de cal-
culs explicites. Qu’on ne s’y trompe pas pourtant : les propriétés que nous mettrons
en évidence sont le plus souvent représentatives de ce qui se passe dans le cas géné-
ral (cela sera d’ailleurs mis en lumière de temps à autre). De plus, la compréhension
de ces exemples simples est très utile pour aborder des théories plus élaborées.

Nous avons l’habitude de classer les équations aux dérivées partielles en trois grandes
classes fondamentales d’équations : les équations elliptiques (qui servent typique-

5
6 TABLE DES MATIÈRES

ment à décrire des phénomènes d’équilibre en physique) pour les problèmes sta-
tionnaires, les équations paraboliques (qui permettent de décrire des phénomènes
de diffusion) et les équations hyperboliques (qui permettent de décrire les phéno-
mènes de propagation) pour les problèmes d’évolution. C’est cette dernière classe
d’équations qui fera l’objet de ce cours. Les deux autres catégories seront abordés
dans deux cours de deuxième année sur la méthode des éléments finis : les équations
elliptiques dans le cours MA201 et les équations paraboliques dans le cours MA206.

Dans ce cours, nous abordons les problèmes linéaires et non linéaires que nous trai-
terons par le biais d’une équation modèle. Pour les équations hyperboliques linéaires,
ce sera l’équation de transport (Chapitre 1). Pour les équations hyperboliques non li-
néaires, enfin, ce seront les lois de conservation scalaires dont le prototype est l’équa-
tion de Burgers (Chapitre 3). Par souci de simplicité, nous nous limiterons la plupart
du temps à la dimension 1 d’espace.

Rares sont les situations réalistes pour lesquelles on sait calculer une solution ana-
lytique à la main. Pour aller plus loin dans la connaissance fine et « quantitative »
des solutions d’une EDPs (dont l’analyse mathématique ne fournira en général que
l’existence d’une solution, voire l’unicité et quelques propriétes qualitatives), l’ingé-
nieur sera amené à se tourner vers l’utilisation de méthodes numériques pour cal-
culer une solution approchée. Le second objectif du cours est d’introduire le lecteur
à la méthode de discrétisation la plus classique (et la plus simple) : la méthode des
différences finies qui est la « mère » des méthodes de discrétisation plus élaborées
telles que les méthodes d’éléments finis, de volumes finis ou de Galerkin discontinu
(qui seront abordées plus tard dans le cursus). La mise en oeuvre de ces méthodes
est en général simple. Leur théorie ne l’est pas forcément et un des buts du cours est
d’aborder - dans le cas des différences finies - les aspects théoriques des méthodes
numériques, ce que l’on appelle l’analyse numérique. Nous introduirons en parti-
culier les concepts de consistance, de stabilité et de convergence et les liens étroits
entre ces trois notions. Le lecteur exigeant pourra s’offusquer de voir que nous avons
choisi la plupart du temps de présenter et analyser les méthodes numériques sur les
cas simples pour lesquels on saurait par exemple calculer la solution à la main. A
nouveau, le but est pédagogique : il est bien clair qu’une méthode numérique doit
bien marcher dans ces cas simples si on veut avoir une chance qu’elle marche dans
les cas compliqués. Par ailleurs, sur ces cas simples, on arrive assez facilement à ob-
tenir des résultats d’analyse numérique précis.

Signalons qu’il y a des connections fortes entre ce cours et d’autres cours de Mathé-
matiques de 1ère année, en particulier le cours MA102 sur les outils fondamentaux
d’analyse et le cours AO102 sur la théorie des équations différentielles ordinaires,
dans lesquels le lecteur trouvera les prérequis nécessaires pour appréhender aisé-
ment le contenu de ce cours.
Chapitre 1

Introduction à la théorie et
l’approximation des problèmes
hyperboliques linéaires

1.1 Exemples de modèles hyperboliques

Commençons par donner quelques exemples d’équations intervenant dans des ap-
plications classiques issues de la physique, de la biologie ou des sciences de l’ingé-
nieur. Elles sont toutes de nature hyperbolique, notion qui sera définie plus loin, ce
qui signifie concrètement que la vitesse de propagation est finie, contrairement à
l’équation de la chaleur.

1.1.1 Equation de transport linéaire

L’équation hyperbolique la plus simple est l’équation de transport linéaire. Elle consiste
à trouver la solution u(x, t) ∈ R de l’équation aux dérivées partielles suivante :

 ∂u ∂u
 + c(x, t) = 0, t ≥ 0, x ∈ R,
∂t ∂x (1.1)

u(x, t = 0) = u0 (x) x ∈ R.

Ici le scalaire c(x, t) est donné : c’est par définition la vitesse de propagation associée
à l’équation, au point x à l’instant t. On comprendra plus loin cette dénomination.

Cette équation qui peut sembler très simple pose des difficultés numériques consi-
dérables, elle est toujours l’objet de recherches actuellement (il s’agit notamment de
savoir calculer la solution sur des temps très grands). Par ailleurs, couplée à d’autres
équations, elle pose des difficultés théoriques également.

7
8 C HAPITRE 1 : Introduction aux problèmes hyperboliques linéaires

On rencontre cette équation dans un grand nombre d’applications. Citons en quelques

unes.

E XEMPLE 1.1.1 (L A CIRCULATION AUTOMOBILE )

On étudie la circulation automobile sur une route. La variable u(x, t) représente la
quantité de voitures présentes entre les bornes x et x + ∆x à l’instant t et on appelle
F (x, t) le flux de voitures par minute qui passent à l’instant t devant la borne x. On
suppose que chaque conducteur ajuste la vitesse de sa voiture en fonction unique-
ment de la vitesse de la voiture qui le précède. Alors la conservation de la quantité
de voitures (il n’y a ni station-service, ni itinéraires de délestage) se traduit par :
∂u ∂F
(x, t) + (x, t) = 0.
∂t ∂x
Si toutes les voitures roulaient à vitesse constante donnée c, on aurait :

F (x, t) = c u(x, t)

et donc l’équation de transport linéaire (1.1).

E XEMPLE 1.1.2 (L ES ÉQUATIONS CINÉTIQUES )
La physique cinétique décrit les plasma ou gaz dilués et fournit un ensemble impor-
tant d’équations de transport dont les variables sont : le temps t, la position x ∈ Ω des
particules avec Ω le domaine d’étude et leur vitesse v ∈ V avec V l’ensemble des vi-
tesses admissibles. L’exemple le plus classique est l’équation de scattering décrivant
l’évolution d’une densité f (x, t, v) de particules (neutrons, amibes ou bactéries en ce
qui concerne les applications à la biologie)

 ∂ ∂
 f (x, t, v) + v f (x, t, v) = K[f ],
∂t ∂x

f (x, t = 0, v) = f (x, v) donnée,
0

où K est donnée par :

Z Z
′ ′ ′
K[f ] = k(v, v )f (x, t, v ) dv − k(v ′ , v)f (x, t, v ′) dv ′ ,
V V

où k(v, v ′ ) représente une probabilité de « tourner » d’une vitesse v ′ à une vitesse v .

Cette équation découle encore d’une relation de conservation de la quantité de par-

ticules.
E XEMPLE 1.1.3 (M ODÈLE DE DÉMOGRAPHIE / RENOUVELLEMENT CELLULAIRE )
Ici u(x, t) représente la densité d’individus d’âge x ≥ 0 à l’instant t. Le taux de morta-
lité est noté d(x) et les individus peuvent donner naissance à des nouveaux nés d’âge
1.1 Exemples de modèles hyperboliques 9

x = 0 avec un taux de fécondité b(x) :

∂ ∂

 ∂t u(x, t) + ∂x u(x, t) + d(x)u(x, t) = 0,

Z +∞


u(0, t) = b(x′ ) u(x′ , t)dx′
0

Dans le cas de la mitose cellulaire, il est naturel de choisir d(x) =χ {x>x0 } 1 (disparition
des cellules qui se divisent) et b(x) = 2χ {x>x0 } (la mitose donne naissance à deux cel-
lules identiques).

1.1.2 L’équation des ondes et ses variantes

L’équation des ondes fournit un des exemples les plus simples de système hyperbo-
lique. Voyons cet aspect en détail. Les ondes scalaires se propagent dans un milieu
monodimensionnel selon l’équation :
 2 2
 ∂ u 2∂ u

 − c = 0, t ∈ R, x ∈ R

 ∂t2 ∂x2


u(x, t = 0) = u0 (x), (1.2)





 ∂u (x, t = 0) = u1 (x).

∂t
Ici la constante c désigne la vitesse de l’onde. Cette équation du deuxième ordre peut
aussi être écrite comme un système d’équations hyperboliques du premier ordre. En
effet, posons :
∂u ∂u
v= , w=c .
∂t ∂x
On a alors : 
 ∂w ∂v

 −c = 0,
∂t ∂x
 ∂v
 ∂w
 −c = 0.
∂t ∂x
Si on introduit le vecteur U et la matrice C définis par :
   
w 0 −c
U= et C =
v −c 0
l’équation (1.2) se met sous la forme :
∂U ∂U
+C =0
∂t ∂x
1
On désigne par χA la fonction caractéristique de A
10 C HAPITRE 1 : Introduction aux problèmes hyperboliques linéaires

qu’on peut voir comme une équation de transport linéaire vectorielle où la solution
est à valeurs dans R2 . C’est en fait un premier exemple de système hyperbolique li-
néaire tel que nous les définirons plus loin.

Nous donnons dans la suite des exemples de modèles de propagation d’ondes en

physique. La liste est évidemment non exhaustive : nous aurions pu citer également
la propagation des ondes élastiques issue des équations de l’élastodynamique, ...
E XEMPLE 1.1.4 (C ORDE VIBRANTE )
Le problème est de calculer le mouvement d’une corde, de longueur L, fixée en ses
extrémités et qui est, soit écartée de sa position d’équilibre et lâchée (corde de gui-
tare), soit frappée (corde de piano) de façon à lui imprimer, en ses différents points,
des vitesses de déplacement vertical.

Ainsi au repos, la corde occupe un intervalle [a, b]. Sous l’action d’une force normale
d’intensité f (dans le cas de la corde de piano) ou d’un écart de sa position d’équi-
libre u0 (écart de position) et u1 (écart de vitesse), la corde se déforme. Si on appelle
u(x, t) le déplacement latéral de la corde au point x à l’instant t, les variations de u
sont décrites par les équations :


 ∂2u ∂2u

 − 2 = f, pour t ∈ R+ , x ∈ [a, b],

 ∂t2 ∂x



u = 0 en x = a et x = b

 u(x, 0) = u0 (x), pour x ∈ [a, b],





 ∂

 u(x, 0) = u1 (x), pour x ∈ [a, b].
∂t
On montre que dans le cas d’une peau de tambour par exemple, le déplacement

f (x, t)
a b x
u(x, t)

F IGURE 1.1 : Déplacement d’une corde vibrante

normal de la peau vérifie également une équation des ondes, la seule différence est
que la position x est dans un domaine de R2 .
E XEMPLE 1.1.5 (ACOUSTIQUE LINÉAIRE )
Le son est en fait une conséquence d’un mouvement matériel d’oscillation : une
corde qui vibre ou la membrane d’un haut-parleur par exemple. Cette vibration pro-
voque un mouvement des atomes l’avoisinant qui va se déplacer de proche en proche
1.1 Exemples de modèles hyperboliques 11

sous forme d’onde de pression. La vélocité du son varie suivant le milieu dans lequel
il se propage. Le principal facteur est la densité de ce milieu : dans un gaz, sa vitesse
est plus faible que dans un liquide. Par exemple, le son se propage approximative-
ment à 340 m/s dans l’air à 15◦ C, et à 1500 m/s (5400 km/h) dans l’eau.

On appelle p(x, t) la pression du fluide au point x à l’instant t et c la vitesse du son

dans le milieu qu’on suppose constante ici (le milieu est homogène). La fonction p
est alors solution de l’équation des ondes dans R3 :

∂2p  2 ∂2p ∂2p 

2 ∂ p
− c + + = f (x, y, z, t)
∂t2 ∂x2 ∂y 2 ∂z 2

où f est la source à l’origine du son.

E XEMPLE 1.1.6 (L ES ÉQUATIONS DE M AXWELL )

Dans cet exemple, la solution recherchée n’est plus un scalaire mais un vecteur. Les
inconnues de ce problème sont le champ électrique E(x, t) ∈ R3 , le champ magné-
tique H(x, t) ∈ R3 , l’induction électrique D(x, t) ∈ R3 et l’induction magnétique
B(x, t) ∈ R3 .
Ces champs obéissent, en l’absence de charges et de courants, aux équations de Max-
well : 
 ∂B

 = −rot E,
∂t

 ∂D = rot H,

∂t
et sont par ailleurs reliés par les lois de comportement :

B(x, t) = µ H(x, t),

D(x, t) = ε E(x, t),

où µ est la perméabilité magnétique et ε la permittivité diélectrique du milieu (constantes

strictement positives dans le cas où le milieu est homogène). Ils caractérisent le com-
portement électromagnétique du matériau dans lequel l’onde se propage.

On montre facilement qu’on peut éliminer B , D et H des équations pour aboutir à

une équation où seul le champ électrique E apparaît :

∂2E 1
ε + rot (rot E) = 0
∂t2 µ

ce qui correspond à une équation des ondes vectorielle dans R3 .

Si on suppose que E est de direction constante, on montre facilement qu’on retombe

sur une équation d’onde scalaire.
12 C HAPITRE 1 : Introduction aux problèmes hyperboliques linéaires

1.2 La méthode des caractéristiques pour l’équation de

transport
1.2.1 Le problème de Cauchy pour l’équation de transport linéaire
On considère ici le problème suivant : trouver u(x, t) : R × R+ → R tel que

 ∂u + c(x, t) ∂u = 0, x ∈ R, t > 0, (i)
∂t ∂x (1.3)

u(x, 0) = u0 (x), x ∈ R, (ii)

où, au moins pour définir une solution classique (voir plus loin), on suppose que
la donnée initiale u0 est une fonction de classe C 1 (on peut admettre une régularité
moindre, par exemple u0 ∈ L∞ ou u0 ∈ L2 , avec le concept de solution faible - voir
remarque 1.2.3).

L’hypothèse sur la régularité du coefficient c(x, t) est plus fondamentale. Nous sup-
poserons que c(·, ·) est une fonction continue de x, t ∈ R × R i. e. c ∈ C 0 (R × R) (nous
autorisons le temps à être négatif ) lipchitzienne en x, uniformément par rapport à t :

∃L>0 / ∀ t ∈ R, ∀ (x1 , x2 ) ∈ R × R, | c(x1 , t) − c(x2 , t) | ≤ L |x1 − x2 | . (1.4)

Pour simplifier un peu les aspects techniques des démonstrations, nous travaillerons
avec l’hypothèse légèrement plus forte suivante :

c(·, ·) est continûment dérivable en x

∂c (1.5)
et ∃L>0 / ∀ t ∈ R, ∀ x ∈ R, | (x, t) | ≤ L .
∂x
R EMARQUE 1.2.1
Contrairement à ce qui se passait pour le coefficient de diffusion dans l’équation de
la chaleur, le signe du coefficient c n’a aucune importance fondamentale, en tous cas
aucune incidence sur l’existence ou non d’une solution. Il en est du même du signe
ou du sens du temps. L’équation de transport peut s’intégrer indifféremment pour
t < 0 ou t > 0. C’est une équation réversible.

R EMARQUE 1.2.2
Rappelons que, le fait que x → c(x, t) soit lipschitzienne entraîne (c’est un résultat
non trivial de théorie de la mesure) que la dérivée au sens des distributions

∂c
,
∂x
est une fonction L∞ . En ce sens les hypothèses (1.4) et (1.5) sont très proches.
1.2 La méthode des caractéristiques pour l’équation de transport 13

Solutions classiques du problème de Cauchy

D ÉFINITION 1.2.1
On dit que u est une solution classique de l’équation (1.3) si c’est une fonction C 1
de x et de t et si elle satisfait (1.3) point par point.

Dans ce qui suit, nous allons montrer que le problème (1.3) admet une unique so-
lution classique et que celle-ci peut être construite par la technique des caractéris-
tiques.

R EMARQUE 1.2.3 (S UR LA RÉGULARITÉ DES DONNÉES )

Le lecteur notera que, pour espérer l’existence d’une solution classique, il est néces-
saire que la donnée initiale u0 soit C 1 et que la fonction c(·, ·) soit continue en ses deux
variables. L’utilité de la condition de Lipschitz (1.4) (et même de (1.5)) (qui est tou-
tefois essentielle) n’apparaîtra qu’au moment où nous exposerons la méthode des
caractéristiques. Lorsque u0 est moins régulière (u0 ∈ Lp par exemple) on peut tou-
jours donner un sens à grâce à la notion de solution faible. Ceci sera abordé plus en
détail dans l’étude des équations hyperboliques non linéaires.

1.2.2 Résolution de l’équation de transport à coefficients constants

Dans le cas où c(·, ·) est une fonction constante, c(x, t) = c, l’opérateur de transport
peut s’interpréter comme un opérateur de dérivation dans une direction oblique du
plan (x, t) :  
∂ ∂  ∂ ∂ t
c
+c = ∇x,t · , ∇x,t = , .
∂t ∂x 1 ∂x ∂t
Par conséquent, les solutions classiques de (1.3) sont constantes le long de la famille
de droites parallèles au vecteur (c, 1)t du plan (x, t)
n o
Dλ = (y, s) ∈ R × R / y − cs = λ . (1.6)

Par définition, les droites Dλ sont la famille des droites caractéristiques de (1.3). Par
tout point (x, t), passe une droite caractéristique et une seule
n o
D(x, t) = (y, s) ∈ R × R / (y − x) − c(s − t) = 0 .

Bien entendu, on observe que

D(x, t) = D(x′ , t′ ) ⇐⇒ x − x′ = c(t − t′ ) , (1.7)

14 C HAPITRE 1 : Introduction aux problèmes hyperboliques linéaires

alors que si on fixe t par exemple, l’ensemble des droites caractéristiques D(x, t) forme,
quand x décrit R, un ensemble de droites parallèles disjointes qui remplissent le plan.
[
D(x, t) ∩ D(x′ , t) = ∅ pour x 6= x′ , D(x, t) = R2 . (1.8)
x∈R

Voyons maintenant comment construire u la solution de l’équation de transport (1.3)

à l’aide des droites caractéristiques. Il nous faut bien évidemment utiliser la condi-
tion initiale . Fixons t > 0, pour calculer la valeur de la solution en (x, t) nous remar-
quons que D(x, t) coupe l’axe t = 0 en un point unique, appelé pied de la caractéris-
tique issue de (x, t) à savoir :

(x0 , 0), x0 = x − ct.

La solution u de (1.3) étant, si elle existe, constante le long de D(x, t), cette solution
est nécessairement donnée par :

u(x, t) = u(x − ct, 0) = u0 (x − ct). (1.9)

A posteriori, on vérifie que la solution donnée par (1.9) est bien solution de (1.3).

u0 (x) u(x, t)
t=0 t>0

x x

Droites caractéristiques
(a) Condition initiale (b) Construction de la solution à l’instant t
F IGURE 1.2 : Méthode des caractéristiques pour l’équation du transport

Nous avons établi le

T HÉORÈME 1.2.1
Le problème admet une unique solution classique u ∈ C 1 (R × R+ ) donnée par
(1.9).

Le graphe de la solution à l’instant t se déduit du graphe de la donnée initiale u0

par une translation de longueur ct selon l’axe des abscisses : la formule (1.9) repré-
sente une fonction qui se propage (ou se transporte) sans déformation à la vitesse
1.2 La méthode des caractéristiques pour l’équation de transport 15

constante c (voir figure (1.2)). Ceci justifie le fait que l’on appelle le coefficient c vi-
tesse de l’équation de transport.

Il est immédiat à partir de (1.9) de vérifier les propriétés suivantes.

• Propagation à vitesse finie. De façon évidente

supp u0 ⊂ [α, β] =⇒ supp u(·, t) ⊂ [α + ct, β + ct].

• Préservation de la régularité. De façon tout à fait évidente, la régularité en la

variable x de la solution u(·, t) est pour tout t la même que celle de u0 .

• Conservation de la norme Lp . Du fait de l’invariance de l’intégrale par transla-

tion, on a évidemment

∀ p ∈ [1, +∞], ku(·, t)kLp = ku0 kLp .

ce qui signifie que l’application linéaire S(t) : u0 7→ u(·, t) est une isométrie dans
tous les espaces Lp . On a donc aucun effet de dissipation.

Quand on évoque la conservation de la norme L2 , on parle souvent de conservation

de l’énergie. Ce résultat peut s’obtenir directement. En effet
Z +∞ Z +∞
∂u ∂u
u dx + c u dx = 0.
−∞ ∂t −∞ ∂x
Z +∞ Z +∞  
∂u ∂ u2
Or u dx = = 0 et par conséquent
−∞ ∂x −∞ ∂x 2
Z +∞
d
u2 dx = 0.
dt −∞

On va utiliser la transformation de Fourier en espace, dont nous rappelons les prin-

cipales propriétés.

Rappels sur la transformation de Fourier

Nous considérons la transformation de Fourier définie par :
F
→ fˆ ∈ C 0 (R) ∩ L∞ (R),
f ∈ L1 (R) −

où fˆ est donnée par Z ∞

fˆ(k) := f (x) e−ikx dx.
−∞
16 C HAPITRE 1 : Introduction aux problèmes hyperboliques linéaires

La transformation F est injective et si fˆ ∈ L1 (R), on a la formule d’ inversion :

Z ∞
−1 ˆ 1
(F f )(x) = f (x) = fˆ(k) eikx dk.
2π −∞
On rappelle enfin que F s’étend aux distributions en un isomorphisme de S ′ (R) et
que  
∂
F (δ) = 1, F = i k F.
∂x
Par ailleurs, il est bien connu que F envoie L2 (R) dans lui même et qu’on a le théo-
rème de Plancherel :
Z +∞ Z +∞
2
∀ f ∈ L (R), ˆ 2
|f (k)| dk = 2π |f (x)|2 dx,
−∞ −∞

On rappelle enfin que , si on définit le produit de convolution de deux fonctions dans

L2 (R) par
Z ∞ Z ∞
(f ∗ g)(x) := f (x − y) g(y) dy = f (y) g(x − y) dy, ∈ L∞ (R),
∞ ∞

alors
F (f ∗ g) = fˆ ĝ.
Enfin, nous rappelons les propriétés de la transformation de Fourier vis à vis des
gaussiennes et d’un changement d’échelle
2
− k2 F −1 2
− x2 F −1 1 x
e −−→ e et b(ak)
u −−→ u (a ∈ R),
a a
et finalement des translations (a ∈ R)


∀ f ∈ L2 (R), τa f (x) = f (x + a) =⇒ F τa f (k) = eika F f (k). (1.10)

Résolution de l’équation de transport à coefficients constants par Fourier

Il est instructif de chercher à résoudre l’équation de transport (1.3) en utilisant la
transformation de Fourier en espace :
u(x, t) → û(k, t).
En appliquant la transformation de Fourier spatiale à (1.3) on obtient l’équation dif-
férentielle ordinaire
dû
+ ick û = 0,
dt
ce qui donne, compte tenu de la condition initiale
û(k, t) = e−ickt û0 (k).
Le lecteur notera alors que
1.2 La méthode des caractéristiques pour l’équation de transport 17

• L’expression (1.9) de la solution exacte se retrouve en utilisant les propriétés de

la transformation de Fourier vis à vis des translations.

• On a en particulier
u(k, t)| = |û0(k)|.
|b
Grâce au théorème de Plancherel, on retrouve la conservation de la norme L2 .

• On a également

b ∆t) u
b(k, t + ∆t) = S(k,
u b(k, t), b ∆t) = e−ick∆t .
S(k, (1.11)

b ∆t)
Le coefficient d’amplification complexe S(k,

– est de module 1 : c’est la manifestation du caractère non dissipatif de l’équa-

tion de transport à coefficients constants,
– a une phase proportionnelle à k : c’est la manifestation du caractère non
dispersif de l’équation de transport à coefficients constants.

1.2.3 Résolution de l’équation de transport à coefficients variables

Cette fois, l’opérateur de transport à coefficients variables
∂ ∂
+ c(x, t)
∂t ∂x
n’est plus vraiment un opérateur de dérivation oblique. Par contre, il peut s’interpré-
ter, comme nous allons le détailler, comme un opérateur de dérivation tangente le
long d’une famille de courbes, appelées courbes caractéristiques, qui vont jouer le
rôle des droites Dλ dans le cas où c était constant. Géométriquement, ces courbes
sont les courbes enveloppes du champ de vecteur 2D
 
c(x, t)
(x, t) =⇒ .
1

Pour les introduire de façon plus analytique,

 cherchons à priori les courbes du plan
(x, t) paramétrées par le temps, soit (X(s), s) s ∈ R où X(s) est une fonction -
supposée dérivable - à déterminer, telles que toute solution classique u de
∂u ∂u
+ c(x, t) =0 (1.12)
∂t ∂x
soit constante le long de ces courbes. Nous étudions les variations de U(s) = u(X(s), s)
le long d’une telle courbe. Il vient

d
∂u
∂u

U ′ (s) = u X(s), s = X(s), s + X(s), s X ′ (s).
ds ∂t ∂x
18 C HAPITRE 1 : Introduction aux problèmes hyperboliques linéaires

On voit, grâce à (1.12), que si on impose à X d’être solution de l´équation différen-

tielle ordinaire :

X ′ (s) = c X(s), s , (1.13)
alors on aura

U ′ (s) = 0 =⇒ u X(s), s = constante.
Les courbes cherchées sont donc les courbes intégrales (ou trajectoires) de l’équation
différentielle (1.13), ce qui établit un lien intime entre la résolution des équations de
transport linéaires et la théorie des équations différentielles. Comme on a affaire à
une équation différentielle ordinaire du premier ordre, on s’attend à ce qu’il existe
une trajectoire et une seule passant par un point (x, t) : cela revient à ajouter à (1.13)
une condition « initiale » en s = t. Cela nous amène naturellement à introduire une
fonction X non plus de la seule variable s mais des trois variables (s; x, t) de la façon
suivante : (x, t) étant fixés et jouant le rôle de paramètres, s → X(s; x, t) est la solution
du problème de Cauchy :

 ∂X (s; x, t) = c(X(s; x, t), s), s ∈ R, (i)
∂s (1.14)

X(t ; x, t) = x. (ii)
Notons que nous cherchons à intégrer (1.13) aussi bien pour s < t que s > t.

L’hypothèse (1.4) est précisément l’hypothèse qui nous permet d’appliquer le célèbre
théorème de Cauchy-Lipschitz relatif à la résolution des équations différentielles or-
dinaires :

L EMME 1.2.2
Pour tout (x0 , t0 ) ∈ R × R, le problème (1.14) admet une unique solution

s 7→ X(s ; x, t) ∈ C 1 (R).

L’application X(·; ·, ·) ainsi définie est appelée flot de l’équation différentielle (1.13).
Par ailleurs, sous l’hypothèse (1.5), on peut démontrer (mais ce n’est pas si évident
que cela) que X est également une fonction de classe C 1 par rapport aux variables x
et t et que les fonctions

 ∂X

 s → Ẋx (s ; x, t) := (s ; x, t),
∂x

 s → Ẋt (s ; x, t) := ∂X (s ; x, t),

∂t
sont respectivement les solutions des équations différentielles linéaires


 ∂ Ẋx (s ; x, t) − ∂c (X(s ; x, t), s) Ẋ (s ; x, t) = 0, s ∈ R,
x
∂s ∂x (1.15)


Ẋx (t ; x, t) = 1,
1.2 La méthode des caractéristiques pour l’équation de transport 19

et 

 ∂ Ẋt (s ; x, t) − ∂c (X(s ; x, t), s) Ẋ (s ; x, t) = 0,
t s ∈ R,
∂s ∂x (1.16)


Ẋt (t ; x, t) = − c(x, t).
A titre d’exercice instructif, nous invitons le lecteur à établir, au moins formellement,
les équations (1.15) et (1.16). Ces équations s’intègrent « à la main » et on obtient :
 Z s 
 ∂X ∂c


 (s ; x, t) = exp X(τ ; x, t), τ dτ , (i)
 ∂x t ∂x
Z s  (1.17)

 ∂X ∂c


 (s ; x, t) = − c(x, t) exp X(τ ; x, t), τ dτ . (ii)
∂t t ∂x

En particulier on voit que que s étant fixé la fonction (x, t) → X(s; x, t) vérifie une
équation aux dérivées partielles qui n’est autre que l’équation de transport,

∂X ∂X
∀x ∈ R, ∀ (s, t) ∈ R × R, (s; x, t) + c(x, t) (s; x, t) = 0. (1.18)
∂t ∂x
R EMARQUE 1.2.4
L’expression (1.17) montre que

∂X
(s ; x, t) > 0,
∂x
et par conséquent que, pour tout (s, t) ∈ R × R, la fonction x −→ X(s; x, t) est stric-
tement croissante.
R EMARQUE 1.2.5
Si c ne vérifie que la condition (1.4), on peut montrer que X admet presque partout
des dérivées partielles par rapport aux variables x et t qui sont toujours caractérisées
par les équations (1.15) et (1.16), la dérivée en x de c étant prise au sens des distribu-
tions.
Nous appellerons courbes caractéristiques de (1.12) les courbes intégrales de (1.14),
c’est à dire les courbes du plan (y, s) (paramétrées par s) définies par :
n o
C(x, t) = (y, s) ∈ R × R / y = X(s ; x, t), s ∈ R .

La courbe C(x, t) est la caractéristique de (1.12) qui passe par le point (x, t). C’est
l’équivalent de la droite D(x, t) du paragraphe précédent (cf. (1.6)).

Les propriétés (1.7) et (1.8) sont donc les équivalents pour les courbes C(x, t) de celles
que nous avons rencontrés pour les des droites D(x, t). Elles se traduisent au travers
20 C HAPITRE 1 : Introduction aux problèmes hyperboliques linéaires

de celles de la fonction X et leur démonstration est simplement un exercice d’ap-

plication (répétée) du théorème de Cauchy-Lipschitz et principalement du résultat
d’unicité associé.

La première propriété précise ce qui est presque une tautologie, à savoir que deux
courbes caractéristiques C(x, t) et C(x′ , t′ ) coïncident si et seulement si (x, t) et (x′ , t′ )
appartiennent à une même trajectoire de l´équation différentielle (1.13) :
C(x, t) = C(x′ , t′ ) ⇐⇒ x′ = X(t′ ; x, t) ⇐⇒ x = X(t; x′ , t′ ). (1.19)
En termes du flot X, cela traduit la relation fondamentale :
∀ (t, t′ ) ∈ R2 , ∀ s ∈ R, X(s; x, t) = X(s; X(t′ ; x, t), t′ ). (1.20)
Cette propriété exprime une propriété (assez intuitive d’ailleurs) dite de semi-groupe
de l’équation différentielle (1.13) qui exprime le fait que intégrer (1.13) de s à t en
partant de la position x équivaut à intégrer (1.13) d’abord entre t et t′ pour calculer
X(t′ ; x, t), puis repartir de t′ pour aller jusqu’au temps s à partir de la nouvelle posi-
tion « initiale » X(t′ ; x, t) au temps t′ .

Par ailleurs, à un instant t donné, l’ensemble des courbes caractéristiques C(x, t) forme,
quand x décrit R un ensemble de courbes disjointes qui remplissent le plan (on dit
qu’elles forment un fibrage du plan).
[
C(x, t) ∩ C(x′ , t) = ∅ pour x 6= x′ , C(x, t) = R2 . (1.21)
x∈R

Cela se traduit en termes de l’application X par la propriété que

Pour tout (s, t), l’application X(s; ·, t) est une bijection de R dans R. (1.22)
On en connait en outre l’inverse :
∀ (s, t) ∈ R × R, X(s; ·, t)−1 = X(t; ·, s), (1.23)
ce qui exprime que pour revenir de x à x, on peut tout d’abord aller au point X(s; x, t)
en intégrant de t à s en partant du point x à l’instant t, puis revenir en x en intégrant
en marche arrière de s à t, en partant du point X(s; x, t) à l’instant s.

P REUVE DES PROPRIÉTÉS (1.19) À (1.23):

Démontrer (1.20) est essentiellement un jeu d’écriture. Fixons t, t′ et x et introduisons les

deux fonctions de la variable réelle s

Y1 (s) = X(s; x, t), Y2 (s) = X(s; X(t′ ; x, t), t′ ).

Par définition de X(·; ·, ·) (voir 1.14), Y1 et Y2 satisfont la même équation différentielle du

premier ordre :
dY1 dY2
(s) = c(Y1 (s), s), (s) = c(Y2 (s), s).
ds ds
1.2 La méthode des caractéristiques pour l’équation de transport 21

De plus si on regarde en s = t′ , il vient grâce à (1.14-(ii))

Y2 (t′ ) = X(t′ ; X(t′ ; x, t), t′ ) = X(t′ ; x, t) = Y1 (t′ ).

Grâce au théorème de Cauchy-Lipschitz (et plus précisément le résultat d’unicité), nous sa-
vons que Y1 et Y2 coïncident pour tout s, ce que nous voulions démontrer.

Démontrons la première équivalence de (1.19). La seconde se déduit de (1.23).

L’égalité C(x, t) = C(x′ , t′ ) entraine en particulier que (x′ , t′ ) appartient à C(x, t) ce qui veut
dire que x′ = X(t′ ; x, t).

La réciproque repose sur (1.20). En effet, si x′ = X(t′ ; x, t)

n
o
C(x′ , t′ ) = s, X(s; x′ , t′ ) , s ∈ R , (par définition)
n
o
= s, X s; X(t′ ; x, t), t′ ) , s ∈ R , (par hypothèse)
n
o
= s, X(s; x, t) , s ∈ R , (par (1.20))

= C(x, t). (par définition)

L’injectivité de X(s; ·, t) résulte du fait que cette fonction est strictement croissante (cf. re-
marque 1.2.4). On peut aussi la voir comme une conséquence de Cauchy-Lipschitz : si X(s; x, t) =
X(s; x′ , t), les fonctions τ → X(τ ; x, t) et τ → X(τ ; x′ , t) seraient deux solutions de (1.13) qui
coïncideraient en τ = s. D’après le théorème d’unicité on en déduit qu’elles sont égales pour
tout τ et donc pour τ = t ce qui donne x = x′ compte tenu de (1.14-(ii)).

Si nous choisissons t = s et t′ = t dans (1.20) nous obtenons, compte tenu de (1.14-(ii)),

x = X(s; X t; x, s), t),

ce qui démontre à la fois la surjectivité de X(s; ·, t) et la formule (1.23).

Finalement, la propriété (1.21) ne fait que traduire géométriquement la bijectivité des appli-
cations X(s; ·, t). 

On peut retrouver (1.18) à partir de (1.20). En effet, dérivons (1.20) par rapport à t′ ,
(x, s, t) étant fixés. En utilisant le théorème de dérivée composée, il vient :
∂X
∂X
∂X ′
0= s; X(t′ ; x, t), t′ + s; X(t′ ; x, t), t′ (t ; x, t),
∂t ∂x ∂t
c’est dire , compte tenu de (1.14-(i))
∂X
∂X

s; X(t′ ; x, t), t′ + c X(t′ ; x, t), t′ ) s; X(t′ ; x, t), t′ = 0.
∂t ∂x
En choisissant t′ = t et en utilisant (1.14-(ii)), nous aboutissons à (1.18).
22 C HAPITRE 1 : Introduction aux problèmes hyperboliques linéaires

Nous pouvons maintenant mettre en place la méthode des caractéristiques. Suppo-

sons que (1.3) ait une solution classique. Pour calculer sa valeur en (x, t), nous utili-
sons le fait que u est constante le long de la caractéristique C(x, t). Autrement dit :

u(x, t) = u(X(s; x, t), s) ∀ s ∈ R.

En particulier, en prenant s = 0, on
introduit naturellement le pied de la caractéris-

tique C(x, t), à savoir X(0; x, t), 0 et on obtient en utilisant la condition initiale :

u(x, t) = u0 X(0; x, t) . (1.24)

Il reste pour conclure à vérifier que la fonction u donnée par (1.24) est bien solu-
tion. C’est bien une fonction de classe C 1 en tant que composée de deux fonctions de
classe C 1 . Par ailleurs,
   
∂u ∂u 0′

∂X
∂X
+c (x, t) = u X(0; x, t) (0; x, t) + c X(0; x, t) (0; x, t) = 0,
∂t ∂x ∂t ∂x

grâce à (1.18). Nous avons donc démontré le

T HÉORÈME 1.2.3
Sous les hypothèses (1.4), le problème admet une unique solution classique don-
née par (1.24).

On peut démontrer aisément à partir de la formule que certaines propriétés de l’équa-

tion de transport à coefficients constants sont préservées.

C OROLLAIRE 1.2.4
La solution u se propage à vitesse finie

supp u0 ⊂ [α, β] =⇒ supp u(·, t) ⊂ [ X(t; α, 0), X(t; β, 0) ] (1.25)

Elle vérifie la propriété de conservation

∀ t ≥ 0, sup u(x, t) = sup u0 (x), inf u(x, t) = inf u0 (x), ku(·, t)kL∞ = ku0 kL∞
x∈R x∈R x∈R x∈R
(1.26)

En revanche certaines propriétés sont perdues. Ainsi, comme les caractéristiques

peuvent s’écarter ou se rapprocher au cours du temps, la solution ne se propage pas
à vitesse constante (elle peut « ralentir » ou « accélérer ») et le profil (ou graphe) de la
solution u(·, t) va se déformer au cours du temps. De ce fait, on n’a plus en général
conservation de la norme Lp de la solution pour p < +∞.
1.3 Systèmes hyperboliques linéaires 23

1.3 Systèmes hyperboliques linéaires

1.3.1 Systèmes hyperboliques linéaires
Pour simplifier l’exposé, nous nous limiterons dans ce qui suit aux systèmes hyper-
boliques linéaires à coefficients constants.

Considérons un système n × n du premier ordre sous la forme suivante :

d n
X X (k) ∂
∂
ui + aij uj = 0, i ∈ {1, 2, . . . , n}. (1.27)
∂t j=1
∂xk
k=1

On cherche ici la solution (u1(x, t), u2 (x, t), . . . , un (x, t)), avec x qui, dans les cas les
plus simples, est un scalaire (x ∈ R, d = 1), un vecteur du plan (x ∈ R2 , d = 2) ou de
l’espace (x ∈ R3 , d = 3). La donnée est une famille de matrices :
(k)
A(k) = (aij ) ∈ Mn×n , ∀k ∈ {1, . . . , d}.

Par exemple, pour le cas n = 2 et d = 1, on retombe sur l’équation des ondes 1D (1.2)
avec :  
0 −c
A=
−c 0
La complexité de cette structure à trois indices peut s’éclaircir en introduisant la no-
tion d’onde plane.

D ÉFINITION 1.3.1
Une onde plane v(y, t) de direction ν = (ν1 , . . . , νd ) ∈ Sd−1 (ensemble des vecteurs
de Rd de norme 1) est une solution ne dépendant que d’une variable d’espace, de
direction ν :
u(x, t) = v(x · ν, t) ∈ Rd

L EMME 1.3.1
Une onde plane solution du système (1.27) satisfait, pour 1 ≤ i ≤ n, l’équation
monodimensionnelle
X n
∂ ∂
vi (y, t) + aij (ν) vj (y, t) = 0,
∂t j=1
∂y

avec la notation
d
X (k)
aij (ν) = νk aij .
k=1
24 C HAPITRE 1 : Introduction aux problèmes hyperboliques linéaires

Par la suite, on notera également

d
X
A(ν) = νk A(k) ∈ Mn×n .
k=1

D ÉFINITION 1.3.2
Considérons le système linéaire du premier ordre (1.27),

1. ce système est dit hyperbolique si ∀ν ∈ Sd−1 , la matrice A(ν) est diagonali-

sable dans R. On note λl (ν), 1 ≤ l ≤ n ses valeurs propres ordonnées,

2. ce système est dit strictement hyperbolique si les n valeurs propres sont dis-
tinctes, on a alors λ1 (ν) < . . . < λn (ν).

R EMARQUE 1.3.1
Pour un système hyperbolique n × n en dimension un (d = 1), il existe une matrice Λ
diagonale et une matrice de passage P telle que A = P ΛP −1, avec :
 
λ1
 λ2 
 
Λ= ... .
 
λn

Posons v(x, t) = P −1 u(x, t) et v 0 (x) = P −1u0 (x), on peut alors écrire n équations de
transport découplées dont les inconnues vj (x, t) sont les composantes de v(x, t) :

 ∂ ∂
 vj + λj vj = 0, j = 1, 2, · · · , n,
∂t ∂x

 vj (x, 0) = vj0 (x).

Pour tout u0 ∈ C 1 (R)n , l’unique solution classique de ce problème de Cauchy sera

donc donnée par :


 u(x, t) = P v(x, t),

vj (x, t) = vj0 (x − λj t),


 v 0 (x) = P −1 u0 (x).

Ainsi, l’étude des systèmes hyperboliques linéaires en dimension un se ramène fa-

cilement à l’étude de plusieurs équations de transport découplées.
1.3 Systèmes hyperboliques linéaires 25

R EMARQUE 1.3.2
En dimension d d’espace, cette base dépend de la direction ν de l’onde plane et il
n’est que très rarement vrai que l’on puisse trouver une base commune. Considérer
par exemple l’équation des ondes en dimension d > 1.

1.3.2 L’équation des ondes en dimension 1

Nous allons donner ici des résultats spécifiques à l’équation des ondes 1D :
 2 2
 ∂ u 2∂ u

 − c = 0, t ∈ R, x ∈ R,

 ∂t2 ∂x2


u(x, 0) = u0 (x), (1.28)





 ∂u (x, 0) = u1 (x).

∂t

La formule de D’Alembert
T HÉORÈME 1.3.2
Toute solution de l’équation des ondes homogène (1.28) se décompose sous la forme

u(x, t) = u+ (x, t) + u− (x, t).

• u+ (x, t) = f (x − ct) est une onde progressive se propageant à la vitesse c dans la

direction x > 0 :
L
u+ (x + L, t) = u+ (x, t − )
c
(le graphe de x → u (x, t + T ) se déduit de celui de x → u+ (x, t) par une simple
+

translation de L = cT vers la droite). Une telle solution est constante sur les
droites x − ct = Cte. Cette famille de droites constitue ce que l’on appelle la pre-
mière famille de courbes caractéristiques associée à l’équation des ondes (1.28).

• u− (x, t) = g(x + ct) est une onde progressive se propageant à la vitesse c dans la
direction x < 0 :
L
u− (x − L, t) = u− (x, t − )
c
(le graphe de x → u (x, t + T ) se déduit de celui de x → u− (x, t) par une simple
−

translation de L = cT vers la gauche). Une telle solution est constante sur les
droites x + ct = Cte qui constitue la seconde famille de courbes caractéristiques
associée à l’équation des ondes.
P REUVE : Nous avons vu au début de se chapitre que l’équation des ondes (1.28) se met sous
forme hyperbolique :
∂U ∂U
+A =0
∂t ∂x
26 C HAPITRE 1 : Introduction aux problèmes hyperboliques linéaires

où : 
∂u  
 ∂t  0 −c
U =
 ∂u 
 et A=
−c 0
c
∂x
On diagonalise aisément sous la forme A = P ΛP −1 , où :
     
c 0 1 1 1/2 −1/2
Λ= P = P −1 =
0 −c −1 1 1/2 1/2

Avec V (x, t) = P −1 U (x, t) et V0 (x) = P −1 U0 (x), on a de manière évidente :

 
∂u0
−c u1 (x)
(x) 
1 ∂x
V0 (x) =  
2 1 ∂u0 
u (x) + c (x)
∂x
et par application de la méthode des caractéristiques :

 
∂u0
− ct) − c u1 (x
(x − ct) 
1 ∂x
V (x, t) =  
2 1 ∂u0 
u (x + ct) + c (x + ct)
∂x
d’où l’on déduit immédiatemment que :
 0 0 
u1 (x − ct) + u1 (x + ct) − c ∂u (x − ct) + c ∂u (x + ct)
1 ∂x ∂x 
U (x, t) =  0 0


2 ∂u ∂u
−u1 (x − ct) + u1 (x + ct) + c (x − ct) + c (x + ct)
∂x ∂x

Alors, en intégrant par rapport au temps la première composante du vecteur U ci-dessus, on

obtient :
Z t Z t Z t Z t
0 1 1 1 1 c ∂u0 c ∂u0
u(x, t)−u (x) = u (x+cs)ds+ u (x−cs)ds+ (x+cs)ds− (x−cs)ds,
2 0 2 0 2 0 ∂x 2 0 ∂x
soit, après deux changements de variable évidents :
Z x+ct Z x−ct
1 1 1 1 1
u(x, t) = u (s)ds − u1 (s)ds + u0 (x + ct) + u0 (x − ct),
2c x 2c x 2 2
et donc la décomposition demandée avec :
Z x Z x+ct
1 1 1 1
+
u (x, t) = u0 (x − ct) + 1
u (s)ds et u (x, t) = u0 (x + ct) +
−
u1 (s)ds,
2 2c x−ct 2 2c x

ce qui donne la célèbre formule de d’Alembert :

1.3 Systèmes hyperboliques linéaires 27

Z x+ct
1 0
1
u(x, t) = u (x + ct) + u0 (x − ct) + u1 (s)ds. (1.29)
2 2c x−ct

Propagation à vitesse finie

La formule (1.29) montre que la valeur de la solution u au point x à l’instant t ne

dépend que des valeurs des données initiales dans l’intervalle [x − ct; x + ct] qui est
aussi la base du cône caractéristique (ou cône de dépendance) D(x, t) issu du point
(x, t) c’est à dire le cône délimité par les deux droites caractéristiques de l’équation
des ondes qui passent par le point (x, t) (voir figure 1.3) :
n o
D(x, t) = (y, s) , 0 ≤ s ≤ t, |y − x| ≤ c(t − s) .

(x, t)

D(x, t)

x − ct x + ct

F IGURE 1.3 : Le cône de dépendance

On en déduit que la solution se propage à vitesse finie au sens où, si on part des
données à support compact, la solution reste à support compact à tout instant. Plus
précisément (voir la figure 1.4) :
[
Supp u0 Supp u1 ⊂ [a, b] ⇒ Supp u(·, t) ⊂ [a − ct; b + ct]

Conservation de l’énergie

Il est naturel d’associer à toute solution u de (1.28) la densité d’énergie :

1  ∂u 2 ∂u 
e(x, t) = | | (x, t) + c2 | |2 (x, t) .
2 ∂t ∂x
28 C HAPITRE 1 : Introduction aux problèmes hyperboliques linéaires

x = a − ct x = b + ct

u=0 u=0

(x2 , t2 )
(x1 , t1 )

D(x1 , t1 ) D(x2 , t2 )
a b
S
Supp u 0
Supp u1

F IGURE 1.4 : Propagation à vitesse finie

On dit que la solution u est d’énergie finie (les solutions qui ont un sens physique
sont généralement d’énergie finie) dès que :
Z
E(t) = e(x, t)dx < +∞,
R

où E(t) est par définition l’énergie totale de la solution.

T HÉORÈME 1.3.3
Il y a conservation d’énergie au cours du temps :
∂
E(t) = 0.
∂t
∂u
P REUVE : En multipliant l’équation des ondes par ∂t , on obtient :

∂ 2 u ∂u 2
2 ∂ u ∂u
− c = 0.
∂t2 ∂t ∂x2 ∂t
Nous observons ensuite que :
Z Z
∂ 2 u ∂u 1 d ∂u 2 
dx = | | dx ,
R ∂t2 ∂t 2 dt R ∂t
alors que, moyennant une intégration par parties :
Z 2 Z
2 ∂ u ∂u c2 d  ∂u 2 
−c 2
dx = | | dx ,
R ∂x ∂t 2 dt R ∂x
Ainsi en ajoutant ces deux dernières égalités, on obtient :
∂
E(t) = 0.
∂t

Chapitre 2

Problèmes hyperboliques non

linéaires : étude des lois de
conservation scalaires

2.1 Lois de conservation

2.1.1 Introduction
Dans ce chapitre, nous nous intéressons aux équations scalaires de la forme :

∂u ∂

+ f (u) = 0, (2.1)
∂t ∂x
où f est une fonction régulière de u. Comme nous l’expliquons au paragraphe 2.1.2,
cette équation traduit la conservation de la quantité u. On dit qu’il s’agit d’une loi de
conservation.

Plus particulièrement, nous nous intéresserons à l’équation aux dérivées partielles

de Bürgers :  
∂u ∂ u2
+ = 0.
∂t ∂x 2
Notre objectif est de résoudre le problème de Cauchy, qui consiste à trouver une fonc-
tion u(x, t) telle que l’on ait :


 ∂u + ∂ f (u) = 0, ∀x ∈ R, ∀t > 0,
∂t ∂x

u(x, 0) = u0 (x), ∀x ∈ R,

où la condition initiale u0 est une fonction donnée.

29
30 C HAPITRE 2 : Problèmes hyperboliques non linéaires :

Notre but dans ce chapitre est d’ébaucher la théorie de ce problème. Nous montre-
rons tout d’abord comment déterminer la solution classique, quand elle existe, par
la méthode des caractéristiques (section 2.2.1). Cependant, de par le caractère non-
linéaire de l’équation, de telles solutions classiques ne peuvent exister globalement
en temps et, typiquement, des discontinuités peuvent apparaître et disparaître au
cours du temps : ce sont les chocs. La fonction u est alors une solution faible du
problème, au sens des distributions (section 2.3.1). Nous montrons que la vitesse
de propagation du choc doit vérifier une condition, appelée condition de Rankine-
Hugoniot (section 2.3.2).

Mais on vérifie aisément qu’il n’y a pas unicité de la solution faible du problème de
Cauchy. L’unique solution physique est celle qui vérifie une condition supplémen-
taire dite condition d’entropie (section 2.3.3). Pour conclure (section 2.4), nous cal-
culons explicitement la solution dans le cas du problème de Riemann où la condition
initiale u0 est de la forme :
(
ug si x < 0,
u0 (x) = .
ud si x > 0.

Cette étude théorique du problème continu est un préalable indispensable à l’étude

des schémas d’approximation par différences finies que nous présenterons dans le
cours suivant. En effet, une méthode numérique de discrétisation devra tenir compte
au mieux des propriétés de la solution du problème continu.

Nous considérerons la plupart du temps une fonction f régulière mais quelconque.

Toutefois, on verra que, lorsque la fonction f est strictement convexe (dont le pro-
totype f (u) = u2 /2 donne l’équation de Bürgers), un certain nombre de résultats
prennent une forme plus simple. En particulier, nous nous limiterons à ce cas pour
l’étude du problème de Riemann.

2.1.2 Construction d’une équation de conservation

L’équation (2.1) résulte en fait de l’écriture d’une simple loi de conservation.

En effet, considérons par exemple une répartition linéique de matière. Notons u(x, t)
la densité de matière en x (x ∈ R) à l’instant t (t ≥ 0) et soit F (x, t) le flux de ma-
tière passant par le point x à l’instant t. Alors, la variation de la quantité de matière
dans l’intervalle [x1 , x2 ] entre les temps t1 et t2 est égale à la différence entre le flux de
matière entrant en x1 et le flux sortant en x2 pendant cette période. Autrement dit :
Z x2 Z t2
{u(ξ, t2) − u(ξ, t1)} dξ = {F (x1 , s) − F (x2 , s)} ds. (2.2)
x1 t1
2.1 Lois de conservation 31

Ceci s’écrit également :

Z x2 Z t2
1 u(ξ, t2) − u(ξ, t1) −1 F (x1 , s) − F (x2 , s)
dξ = ds.
x2 − x1 x1 t2 − t1 t2 − t1 t1 x1 − x2

Lorsque (x2 −x1 ) et (t2 −t1 ) tendent vers 0, si les fonctions u et F sont assez régulières,
cette identité intégrale devient l’équation aux dérivées partielles suivante :

∂u ∂F
(x, t) + (x, t) = 0. (2.3)
∂t ∂x
Supposons de plus que le flux F soit lui-même une fonction de la densité u :

F (x, t) = f (u(x, t)). (2.4)

Alors l’équation (2.3) prend la forme suivante :

∂u ∂
(x, t) + f (u(x, t)) = 0. (2.5)
∂t ∂x
On dit que l’équation (2.5) est sous forme conservative, par opposition à sa forme
non conservative :
∂u ∂u
(x, t) + c(u(x, t)) (x, t) = 0 (2.6)
∂t ∂x
où c(u) = f ′ (u).
u2
Dans le cas particulier où f (u) = , les équations (2.5) et (2.6) s’écrivent :
2
 
∂u ∂ u2
+ = 0, (2.7)
∂t ∂x 2

∂u ∂u
+u = 0.
∂t ∂x
Il s’agit alors de l’équation de Bürgers.

2.1.3 Exemples
La circulation automobile
L’équation (2.5) fournit un modèle simple pour la circulation automobile sur une
route à une seule voie : u représente la quantité de voitures par centaine de mètres et
F le flux de voitures par minute qui passent à l’heure t devant la borne x. On suppose
que chaque conducteur ajuste la vitesse de sa voiture en fonction uniquement de la
vitesse de la voiture qui le précède. Alors la conservation de la quantité de voitures (il
n’y a ni station-service, ni itinéraire de délestage) se traduit par l’équation (2.5).
32 C HAPITRE 2 : Problèmes hyperboliques non linéaires :

Si toutes les voitures roulaient à une vitesse constante donnée c, on aurait :


 f (u) = c u
 ∂u + c ∂u = 0
∂t ∂x
et on retrouve l’équation de transport qui est linéaire.

En réalité, le flux de voitures F est une fonction f non-linéaire de la densité de voi-

tures u. Essayons d’imaginer à quoi peut ressembler cette fonction...

Bien sûr, le flux de voitures est nul si la densité est elle-même nulle, mais il est égale-
ment nul si la densité est trop importante (u voisin de la densité de saturation umax ) :
on a alors un embouteillage. Enfin, on observe un flux maximum au voisinage d’une
densité intermédiaire idéale ū.

u
0 u umax

Supposons que la fonction f est strictement concave sur l’intervalle des densités ad-
missibles [0, umax ]. Notons v l’écart entre la densité de voitures idéale ū et la densité
effective u :
v(x, t) = ū − u(x, t).
On a alors, sur l’intervalle [ū − umax , ū] :
∂v ∂
+ (g(v)) = 0 (2.8)
∂t ∂x
où g est la fonction strictement convexe donnée par g(v) = −f (ū − v).

La dynamique des gaz

En réalité, la loi de conservation scalaire (2.5) a surtout été étudiée parce qu’elle
constitue un modèle simple pour l’étude des systèmes de lois de conservation qui
2.1 Lois de conservation 33

modélisent de nombreux phénomènes physiques. L’exemple fondamental est celui

de la dynamique des gaz.

En effet, considérons un fluide parfait compressible. Alors un écoulement unidimen-

sionnel est modélisé, en coordonnées Eulériennes, par :
• L’équation de conservation de la masse
∂ρ ∂
+ (ρu) = 0 ;
∂t ∂x

• L’équation de conservation de la quantité de mouvement

∂ ∂
(ρu) + (ρu2 + p) = 0 ;
∂t ∂x

• L’équation de conservation de l’énergie

∂ ∂

(ρe) + (ρe + p)u = 0 ;
∂t ∂x

où ρ représente la densité du fluide, u la vitesse et e l’énergie totale spécifique. La

pression p est donnée par une équation d’état en fonction des quantités ρ, u et e.

L’étude de ce système de trois équations et trois inconnues présente de nombreuses

difficultés théoriques dont certaines n’ont pas encore été résolues à ce jour. C’est
pourquoi l’on s’intéresse à l’équation scalaire (2.5) qui, bien que beaucoup plus simple
à étudier, présente les mêmes caractères fondamentaux que le système précédent.

2.1.4 Ajout d’un terme de viscosité

Jusque là, nous avons fait l’hypothèse que le flux F (x, t) avait la forme (2.4). Pour
beaucoup de problèmes issus de la physique, il est plus réaliste de supposer que l’on
a:
∂u
F (x, t) = f (u(x, t)) − ε (x, t) (2.9)
∂x
où ε est un petit paramètre strictement positif. L’équation (2.5) devient alors :
∂u ∂ ∂2u
+ (f (u)) − ε 2 = 0. (2.10)
∂t ∂x ∂x
Le terme
∂2u
−ε
∂x2
permet de tenir compte des effets de dissipation de l’énergie et de dispersion. Par
analogie avec la dynamique des gaz, on dit qu’il s’agit d’un terme de viscosité.
34 C HAPITRE 2 : Problèmes hyperboliques non linéaires :

Pour l’exemple de la circulation automobile, l’expression du flux (2.9) traduit le fait

que le conducteur ne tient pas seulement compte de la densité de voitures u là où il
se trouve mais qu’il anticipe en fonction de ce qu’il aperçoit devant lui.

L’équation (2.5) est en fait la limite de l’équation visqueuse (2.10) lorsque le para-
mètre de viscosité ε tend vers 0. Or, contrairement à l’équation (2.5), l’équation (2.10)
admet une solution unique et très régulière. Cette remarque, très importante, nous
permettra de repérer la solution physique parmi toutes les solutions faibles de l’équa-
tion (2.5).

2.2 Solutions classiques : méthode des caractéristiques

On considère le problème suivant :

∂u ∂

+ f (u) = 0 t > 0, x∈R (2.11)
∂t ∂x

u(x, 0) = u0 (x) x∈R (2.12)

où la condition initiale u0 est une fonction C 1 de x. Rappelons que la fonction f est

supposée être très régulière (disons au moins de classe C 2 ) et posons :

c(u) = f ′ (u), (2.13)

de sorte que a est une fonction croissante dès que f est convexe.

R EMARQUE 2.2.1
Dans le cas linéaire, f (u) = c u, c(u) = f ′ (u) = c n’est autre que la vitesse de propa-
gation (constante) des solutions de l’équation de transport. En quelque sorte, dans le
cas non linéaire, on peut dire que la vitesse de propagation dépend de la solution !

2.2.1 Solutions classiques au problème de Cauchy

D ÉFINITION 2.2.1 (S OLUTION CLASSIQUE )

On dit que u est une solution classique de l’équation (2.11) dans un domaine
ouvert Q de R × R+ si c’est une fonction C 1 de x et de t, dans Q et si elle satisfait
(2.11) point par point dans Q.

Nous allons montrer que le problème (2.11-2.12) admet, au moins pour t < T ≤ +∞,
une solution classique (on parle de solution classique locale si T < +∞) et que celle-
ci peut être construite par la technique des droites caractéristiques.
2.2 Solutions classiques : méthode des caractéristiques 35

R EMARQUE 2.2.2 (C AS D ’ UNE DONNÉE INITIALE CONTINUE ET C 1 PAR MORCEAUX )

Dans tout ce paragraphe, nous supposerons que la donnée initiale u0 est C 1 et nous
nous intéresserons aux solutions classiques du problème (2.11-2.12). Cependant, tous
les résultats que nous allons établir peuvent se généraliser à une donnée initiale u0
continue qui est seulement C 1 par morceaux. La méthode des caractéristiques permet
alors de construire une fonction u continue et C 1 par morceaux qui vérifie l’équation
(2.11) point par point dans tout ouvert où elle est C 1 , et telle que (2.12) soit satisfait.
Nous verrons au paragraphe suivant qu’une telle fonction est une solution faible du
problème.

Droites caractéristiques

Soit u une solution classique de l’équation (2.11). Alors u satisfait l’équation sous
forme non conservative :
∂u ∂u
+ c(u) =0 (2.14)
∂t ∂x
où a(u) est donné par (2.13).

Considérons maintenant la courbe (x(t), t) du plan (x, t) définie par :


 dx (t) = c(u(x(t), t))
dt (2.15)

x (t̄) = x̄

où (x̄, t̄) est un point quelconque du plan.

Alors on peut dériver u par rapport au temps le long de cette courbe et l’on obtient :

d ∂u dx ∂u
u(x(t), t) = (x(t), t) (t) + (x(t), t).
dt ∂x dt ∂t
D’après (2.14) et (2.15), cette quantité est nulle. Autrement dit, u est constant le long
de la courbe (2.15) :
u(x(t), t) = u(x̄, t̄) ∀t .
On en déduit que la courbe en question est une droite du plan (x, t) d’équation :

x = x̄ + c(u(x̄, t̄))(t − t̄). (2.16)

D ÉFINITION 2.2.2 (D ROITES CARACTÉRISTIQUES )

Soit une solution classique u de l’équation (2.11). Alors les droites de la forme
(2.15) ou (2.16) sont appelées droites caractéristiques pour u de l’équation (2.11).
De plus, u est constante le long de chacune de ces droites.
36 C HAPITRE 2 : Problèmes hyperboliques non linéaires :

Construction d’une solution classique du problème de Cauchy

Supposons que u est, au moins pour t < T, une solution classique du problème (2.11-
2.12). Considérons alors un point (x̄, t̄) de la forme (ξ, 0) . D’après (2.6), la droite ca-
ractéristique passant par (ξ, 0) a pour équation :


x = ξ + c u0 (ξ) t (2.17)
et le long de cette droite, on a :
u(x, t) = u0 (ξ). (2.18)
Ceci permet donc de calculer la valeur de u en tout point (x, t) par lequel passe une
et une seule droite caractéristique.

Considérons tout d’abord le cas de l’équation de transport (la vitesse c est constante)
∂u ∂u
+c = 0.
∂t ∂x
Les droites caractéristiques ont pour équation
x = ξ + ct
et sont donc parallèles. D’après (2.18), la solution u est alors donnée, ∀t ≥ 0, par :
u(x, t) = u0 (x − ct).
Il s’agit d’une onde qui se propage « vers l’avant » à la vitesse c sans se déformer.

u(x,t) t
t O

u0(x) x x
Equation de transport
t=0

x u(x,t) t
t O
Condition initiale

x x
Equation de Burgers
2.2 Solutions classiques : méthode des caractéristiques 37

En revanche, dans le cas non-linéaire de l’équation (2.11), les droites caractéristiques

ne sont pas parallèles. Toutes les valeurs de u0 ne se propagent pas à la même vitesse.
La solution se déforme au cours du temps.

R EMARQUE 2.2.3 (ATTENTION !)

On représente toujours les droites caractéristiques dans le plan (x, t)
et non dans le
0
plan (t, x). La pente de la droite d’équation (2.17) est donc 1/c u (ξ) . En particulier,
si c(u0 (ξ)) = 0, la droite caractéristique issue de ξ est verticale.

Pour que la méthode des caractéristiques définisse sans ambiguité la valeur de la

solution au point (x, t), il faut qu’il n’y ait qu’une droite caractéristique qui passe par
ce point. Ceci nous amène a introduire la fonction


c0 (x) = c u0 (x) . (2.19)
Notons que c0 (x) peut s’interpréter comme la vitesse de propagation initiale de la
solution au point x : nous l’appellerons vitesse initiale. La caractéristique associée
passant par le point (ξ, 0) a pour équation :
x = ξ + c0 (ξ) t

Le cas d’une vitesse initiale croissante : existence d’une unique solution globale

T HÉORÈME 2.2.1

(E XISTENCE D ’ UNE SOLUTION CLASSIQUE GLOBALE )
Si c0 = c u est croissante sur R, le problème de Cauchy (2.11-2.12) admet une
0

unique solution classique globale, construite par la méthode des caractéristiques.

D ÉMONSTRATION : On fixe t ∈ R+ . Si c0 est croissante, la fonction continue :

ξ −→ F (ξ) = ξ + c0 (ξ) t

définit une bijection de R dans R. En effet :

dF
(ξ) = 1 + c′0 (ξ) t. (2.20)
dξ
De plus, comme c0 (ξ) ≥ c0 (0) pour ξ > 0 et c0 (ξ) ≤ c0 (0) pour ξ < 0 :

lim F (ξ) = ±∞.

ξ−→±∞

Pour tout couple (x, t), il existe donc un unique ξ(x, t) tel que :


x = ξ(x, t) + c0 ξ(x, t) t. (2.21)

Si la solution classique existe, elle est donc donnée par :

u(x, t) = u0 ξ(x, t) .
38 C HAPITRE 2 : Problèmes hyperboliques non linéaires :

Vérifions que u est effectivement solution de l’équation (2.14). On a :


∂ξ
 ∂u
 (x, t) = u′0 ξ(x, t) (x, t)
∂x ∂x
 ∂u (x, t) = u′ ξ(x, t)
∂ξ (x, t)

0
∂t ∂t
et donc
 
∂u
∂u ′

∂ξ
∂ξ
(x, t) + c u(x, t) (x, t) = u0 ξ(x, t) (x, t) + c0 ξ(x, t) (x, t) (2.22)
∂t ∂x ∂t ∂x
Or, en dérivant respectivement par rapport à x et t, il vient :
  
 ∂ξ
 (x, t) 1 + c′0 (ξ) t = 1
∂x (2.23)

 ∂ξ  

(x, t) 1 + c′0 (ξ) t = −c0 ξ(x, t)
∂t
Mais par hypothèse :
c′0 (ξ(x, t)) ≥ 0
donc la quantité entre crochets dans (2.23) reste strictement positive pour tout t ≥ 0 et on
déduit de (2.22) que l’équation (2.14) est donc vérifiée pour tout temps positif. 

R EMARQUE 2.2.4 (C AS CONVEXE )

Lorsque f est convexe, c est croissante et la condition « c0 est croissante » devient
simplement « u0 est croissante »

Exemple
Considérons par exemple la condition initiale suivante :


 0 si x ≤ 0


x
u0 (x) = si 0 ≤ x ≤ α

 α

 1 si x ≥ α

où α est un réel strictement positif. Nous nous permettons de considérer une donnée
initiale qui n’est que C 1 par morceaux en vertu de la remarque 2.2.2. On veut résoudre

le problème de Cauchy (2.11-2.12) pour l’équation de Bürgers f (u) = u2 /2 et c(u) = u.
D’après (2.17), la droite caractéristique issue de (ξ, 0) a pour équation :


 x=ξ si ξ ≤ 0


ξ
 x = ξ + t si 0 ≤ ξ ≤ α .

 α

x=ξ+t si ξ ≥ α
2.2 Solutions classiques : méthode des caractéristiques 39

La solution u est alors donnée par :



 0 si x ≤ 0


x
u(x, t) = u0 (ξ) = si 0 ≤ x ≤ α + t .

 α+t

 1 si x ≥ α + t
On remarque que pour tout t fixé, la fonction u reste croissante en x.

On a représenté ci-dessous la solution u en fonction de x à différents temps :

u
t=
t=
t=

2.2.2 Solutions classiques locales et naissance d’un choc

Nous abandonnons maintenant l’hypothèse que vitesse initiale c0 définie par (2.19)est
une fonction croissante. Dans ce cas, les solutions classiques peuvent ne pas exister
pour tout temps mais seulement localement.

Dans cette section, pour des raisons purement techniques, nous ferons l’hypothèse
40 C HAPITRE 2 : Problèmes hyperboliques non linéaires :

supplémentaire, peu restrictive en pratique :

La fonction x 7→ u0 (x) est bornée sur R. (2.24)

On en déduit bien sur que la vitesse initiale c0 est elle même bornée.

Existence d’une solution classique locale

T HÉORÈME 2.2.2 (T EMPS D ’ APPARITION D ’ UN CHOC )

Si la dérivée de c0 prend des valeurs strictement négatives, alors il n’existe pas de
solution classique pour tout temps t > 0. Le temps maximal t∗ tel qu’il existe une
solution classique pour 0 ≤ t < t∗ est donné par :
1  ′ −1
t∗ = = inf −c0 (ξ) .
sup {−c′0 (ξ)} ξ∈R
ξ∈R

R EMARQUE 2.2.5
L’hypothèse assure que 0 < t∗ < +∞.
D ÉMONSTRATION : Nous commençons par montrer l’existence d’une unique solution clas-
sique locale pour 0 < t < t∗ . On reprend en fait la démonstration du théorème 2.2.1. Pour
tout t < t∗ , la fonction ξ −→ F (ξ) définie dans la preuve du lemme 2.2.1 est strictement crois-
sante (sa dérivée - cf. (2.20) - est supérieure à 1 − t/t∗ ) et tend vers ±∞ quand ξ → ±∞ (le
lecteur notera que c’est l’hypothèse technique (2.24) qui permet de le garantir). Elle définit
donc à nouveau une bijection de R dans R. L’équation (2.21) admet alors une solution ξ(x, t)
unique, ce qui permet d’établir l’existence d’un solution classique u pour t ∈ [0, t∗ [.

Notons que, d’après les formules (2.23), les dérivées de la solution classique tendent vers l’in-
fini lorsque t tend vers t∗ . Démontrons directement que la solution classique ne peut exister
au delà du temps t∗ .

Soit tmax , le temps maximum tmax d’existence d’une solution classique, défini par

tmax = sup{ t > 0 / il existe une solution classique sur R+ × [0, t[ }

Nous savons déjà que tmax ≥ t∗ . Supposons qu’il existe ξ2 > ξ1 tels que c0 (ξ1 ) > c0 (ξ2 ). Alors,
les droites caractéristiques issues des points ξ1 et ξ2 se croisent à l’instant t̄ et au point x̄
donnés par : 

 t̄ = ξ2 − ξ1
,
c0 (ξ1 ) − c0 (ξ2 )

 x̄ = ξ + c (ξ ) t̄ = ξ + c (ξ ) t̄
1 0 1 2 0 2

S’il existait une solution classique u jusqu’à l’instant t̄, elle devrait donc vérifier simultané-
ment u(x̄, t̄) = u0 (ξ1 ) et u(x̄, t̄) = u0 (ξ2 ), ce qui est impossible (comme c0 (ξ1 ) > c0 (ξ2 )). Donc
2.2 Solutions classiques : méthode des caractéristiques 41

t̄ ≥ tmax .

Soit ξ tel que c′0 (ξ) < 0. Pour h > 0 assez petit on a nécessairement c0 (ξ + h) < c0 (ξ) et donc
d’après ce qui précède avec ξ1 = ξ et ξ2 = ξ + h
h 1
tmax ≤ =⇒ tmax ≤ − (en faisant tendre h vers 0).
c0 (ξ) − c0 (ξ + h) c′0 (ξ)

Ceci étant vrai pour tout ξ, on en déduit que tmax ≤ t∗ 

x
x

R EMARQUE 2.2.6
La naissance d’une discontinuité est un phénomène purement non-linéaire, qui peut
se produire même si u0 est très régulière.

Exemple
Considérons cette fois la condition initiale suivante :


 1 si x ≤ 0,


x
u0 (x) = 1− si 0 ≤ x ≤ α,

 α

 0 si x ≥ α.

où α est un réel strictement positif.

Pour l’équation de Bürgers, la droite caractéristique issue de ξ a pour équation :



 x=ξ+t si ξ ≤ 0,

  
ξ
 x=ξ+ 1− t si 0 ≤ ξ ≤ α,

 α

x=ξ si ξ ≥ α.
42 C HAPITRE 2 : Problèmes hyperboliques non linéaires :

Le temps critique t∗ est égal à α. La solution u, pour 0 ≤ t < α, vaut :



 1 si x ≤ t,


u(x, t) = x−α
 si t ≤ x ≤ α,

 t−α

0 si x ≥ α.

t=

On peut représenter la solution pour différents temps :

u
t=
t=
t=
1

On remarque que la solution reste décroissante en x pour tout t, 0 ≤ t < α. Lorsque

t augmente de 0 à α, la pente de la solution devient de plus en plus raide. A l’instant
t = α, la dérivée en x de u est infinie en x = α. La solution u est discontinue en ce
point.

2.3 Solutions faibles - Solutions entropiques

On a vu au paragraphe précédent que, même si la donnée initiale u0 est très régulière,
il peut ne pas exister de solution classique au problème de Cauchy (2.11-2.12) pour
2.3 Solutions faibles - Solutions entropiques 43

tout temps t > 0. A partir d’un temps critique t∗ , la solution u devient discontinue. En
quel sens, peut-on prolonger cette solution au delà du temps t∗ ? C’est la question à
laquelle nous allons chercher à répondre.

2.3.1 Solutions faibles

Rappelons-nous que l’équation (2.11) a été obtenue comme conséquence de la loi de
conservation intégrale (2.2). Or cette dernière a un sens pour les fonctions disconti-
nues. Nous allons donc repartir de cette idée en reformulant cette notion à l’aide de
fonctions tests.

Soit u une solution classique de (2.11-2.12) et ϕ une fonction C 1 (R×[0, +∞[) à support
compact (c’est notre fonction test). Alors, on a, par intégration par parties :
Z +∞ Z +∞  
∂u ∂
0 = + (f (u)) ϕ dt dx
x=−∞ t=0 ∂t ∂x
Z +∞ Z +∞   Z +∞
∂ϕ ∂ϕ
= −u − f (u) dt dx − u0 (x) ϕ(x, 0) dx.
x=−∞ t=0 ∂t ∂x x=−∞

Cette égalité ne fait plus appel explicitement aux dérivées de u et a donc un sens
même si u n’est pas a priori régulière : il suffit par exemple que u soit localement in-
tégrable ! Cette observation nous conduit à introduire la notion de solution faible qui
va seulement demander à la solution d’être localement bornée (donc en particulier
localement intégrable) sur R × [0, +∞[.

D ÉFINITION 2.3.1 (S OLUTION FAIBLE )

Soit u0 ∈ L∞
ℓoc (R). Une fonction u ∈ Lℓoc (R×[0, +∞[), est appelée une solution faible
∞

du problème de Cauchy (2.11-2.12) si elle satisfait

Z +∞ Z +∞   Z +∞
∂ϕ ∂ϕ
u + f (u) dx dt + u0 (x)ϕ(x, 0)dx = 0 (2.25)
x=−∞ t=0 ∂t ∂x x=−∞

pour tout ϕ de classe C 1 à support compact dans R × [0, +∞[.

Le lemme suivant assure la cohérence des définitions 2.2.1 et 2.3.1.

L EMME 2.3.1 (L IEN SOLUTION FAIBLE - SOLUTION CLASSIQUE )

Si u est une solution classique du problème de Cauchy (2.11-2.12), alors c’est une
solution faible. Si u ∈ C 1 (R × [0, +∞[) est une solution faible du problème de Cau-
chy (2.11-2.12), alors c’est une solution classique.

D ÉMONSTRATION :
La première partie a été démontrée. Réciproquement, soit u ∈ C 1 (R × [0, +∞[) une solution
44 C HAPITRE 2 : Problèmes hyperboliques non linéaires :

faible de (2.11-2.12). Alors, d’après (2.25), on a :

Z +∞ Z +∞  
∂ϕ ∂ϕ
u + f (u) dt dx = 0
x=−∞ t=0 ∂t ∂x

pour tout ϕ ∈ D(R×]0, +∞[). Autrement dit, on a :

∂u ∂
+ (f (u)) = 0
∂t ∂x
au sens des distributions. Mais comme u est régulière, ceci est en fait vérifié au sens classique,
point par point.

Soit alors une fonction ϕ ∈ C 1 (R × [0, +∞[) à support compact dans R × [0, +∞[. En intégrant
(2.25) par parties, on obtient :
Z +∞ Z +∞   Z +∞
∂u ∂

0 = − + (f (u)) ϕ dt dx + u0 (x) − u(x, 0) ϕ(x, 0) dx
x=−∞ t=0 ∂t ∂x x=−∞
Z +∞


= u0 (x) − u(x, 0) ϕ(x, 0) dx.
x=−∞

Ceci étant vrai pour tout ϕ, on en déduit (2.12). 

2.3.2 Solutions C 1 par morceaux. Condition de choc

Dans cette section, nous allons nous intéresser à une classe particulière de solutions
faibles, les solutions C 1 par morceaux.

D ÉFINITION 2.3.2 (F ONCTION C 1 PAR MORCEAUX )

Soit O un ouvert borné de R×[0, +∞[. Une fonction v(x, t) est dite C 1 par morceaux
sur O s’il existe un nombre fini de courbes Σ1 , . . . , Σp de la forme :

 x = ξi (t), t ∈ [t(1) (2)
i , ti ],
Σi :
 où ξ est une fonction C 1 de t
i

telles que v soit égale à la restriction d’une fonction de C 1 (R2 ) dans chaque
composante connexe de O\Σ1 ∪ Σ2 · · ∪ Σp .

Une fonction v(x, t) est dite C 1 par morceaux sur R × [0, +∞[ si elle est C 1 par
morceaux sur tout ouvert borné de R × [0, +∞[.

Notons que si une fonction v, C 1 par morceaux, peut admettre, à travers un nombre
2.3 Solutions faibles - Solutions entropiques 45

fini de courbes Σj , des discontinuités en valeur ou en dérivée cette fonction, ainsi

que ses dérivées en temps et en espace admettent des limites de part et d’autre de
ces courbes.

Nous pouvons maintenant considérer des solutions du problème qui soient seule-
ment C 1 par morceaux. Nous allons voir cependant que les discontinuités d’une telle
solution ne sont pas arbitraires.

Pour fixer les notations, si u est une fonction C 1 par morceaux dans le plan (x, t) et
Σ une courbe de discontinuité de u de la forme (ξ(t), t), on note :

• n la normale unitaire à Σ en (ξ(t), t) orientée vers les x positifs, de sorte que :

! !
nx 1 1
n= =√ où s = ξ ′ (t). (2.26)
1 + s 2
nt −s

• u− et u+ les valeurs de u « à gauche » et « à droite » de Σ :

 −

lim u (x, t) − εn ,
 u = ε→0


 u+ = lim u (x, t) + εn .
ε→0

u-
+
u n

x
46 C HAPITRE 2 : Problèmes hyperboliques non linéaires :

Nous pouvons maintenant démontrer le théorème suivant

T HÉORÈME 2.3.2 (E QUATIONS DES DISCONTINUITÉS )

Soit u une fonction C 1 par morceaux dans R × [0, +∞[. Alors u est une solution
faible du problème de Cauchy (2.11-2.12) si et seulement si :

(i) u est une solution classique de (2.11-2.12) dans tout domaine où elle est C 1

(ii) u satisfait la condition :

ξ ′ (t) (u+ − u− ) = f (u+ ) − f (u− ) (2.27)

le long de toute courbe de discontinuité Σ.

D ÉMONSTRATION : Soit u une fonction C 1 par morceaux dans le plan (x, t) et Σ une courbe de
discontinuité de u. Soit alors K un ouvert borné que la courbe Σ partage en deux compo-
santes connexes K − et K + . Pour simplifier, on supposera que K est contenu dans le demi-
plan t > 0 mais ceci n’est pas restrictif. Enfin, on suppose que u est C 1 dans l’adhérence de
chacun des ouverts K − et K + .

-
K
K
K+

Soit ϕ ∈ D(K). Comme u est régulière sur K − et K + :

Z   Z  
∂ϕ ∂ϕ ∂u ∂
u + f (u) dx dt = − + (f (u)) ϕ dx dt
K− ∂t ∂x K− ∂t ∂x
Z


+ u− nt + f (u− )nx ϕ dγ,
Σ

et de même :
Z   Z  
∂ϕ ∂ϕ ∂u ∂
u + f (u) dx dt = − + (f (u)) ϕ dx dt
K+ ∂t ∂x K+ ∂t ∂x
Z


− u+ nt + f (u+ )nx ϕ dγ.
Σ
2.3 Solutions faibles - Solutions entropiques 47

D’où finalement :
Z   Z  
∂ϕ ∂ϕ ∂u ∂
u + f (u) dx dt = − + (f (u) ϕdx dt
K ∂x ∂x K − ∪K + ∂t ∂x
Z (2.28)
 +
− (u − u− )nt + (f (u+ ) − f (u− ))nx ϕ dγ.
Σ

Le théorème s’en déduit. En effet, si u est une solution faible du problème (2.11-2.12), elle
satisfait l’équation (2.11) au sens des distributions dans K et au sens classique dans K − et
K + . Autrement dit :
Z   Z  
∂ϕ ∂ϕ ∂u ∂
0= u + f (u) dx dt = + (f (u)) dx dt.
K ∂t ∂x K − ∪K + ∂t ∂x

D’après (2.28), on a donc :

Z


(u+ − u− ) nt + f (u+ ) − f (u− ) nx ϕ dγ = 0.
Σ

Comme ceci est vrai pour tout ϕ ∈ D(K), on en déduit en utilisant l’expression exacte de la
normale (2.26) :
(u+ − u− )(−ξ ′ (t)) + f (u+ ) − f (u− ) = 0,
le long de Σ ∩ K.

Réciproquement, si u satisfait (i) et (ii), on a d’après (2.28) :

Z  
∂ϕ ∂ϕ
u + f (u) dx dt = 0
K ∂t ∂x

pour tout ϕ ∈ D(K), donc u est solution faible du problème. 

On écrira souvent la relation (2.27) sous la forme synthétique :

s[u] = [f (u)] où s = ξ ′ (t).

Par analogie à la dynamique des gaz, cette condition de choc est généralement appe-
lée condition de Rankine Hugoniot.

R EMARQUE 2.3.1 (E QUATION DE TRANSPORT - ÉQUATION DE BÜRGERS )

Considérons deux cas particuliers.

• Dans le cas particulier de l’équation de transport, on a : f (u) = c u. La condition

de choc s’écrit donc : s = c. Autrement dit, dans le cas linéaire, les discontinui-
tés se propagent nécessairement le long des caractéristiques.
48 C HAPITRE 2 : Problèmes hyperboliques non linéaires :

• Dans le cas de l’équation de Bürgers, la condition de Rankine-Hugoniot s’écrit :

u+ + u−
s= .
2

Autrement dit, la vitesse du choc est la moyenne de la vitesse des caractéris-

tiques de part et d’autre.

Nous sommes maintenant en mesure de justifier rigoureusement la remarque 2.2.2.

C OROLLAIRE 2.3.3 (C AS DES FONCTIONS C 1 PAR MORCEAUX ET CONTINUES )

Soit u une fonction C 1 par morceaux et continue dans R × [0, +∞[. Si u est une
solution classique du problème (2.11-2.12) dans tout domaine où elle est C 1 , alors
c’est une solution faible du problème (2.11-2.12).

D ÉMONSTRATION : En effet, la condition de choc (2.27) est automatiquement satisfaite car on

a toujours u+ ≡ u− . 

Exemple

Nous allons maintenant construire une solution faible globale pour la condition ini-
tiale de l’exemple section 2.2.2. On rappelle qu’il existait une solution continue jus-
qu’au temps t = α. La fonction u vaut alors :
(
1 si x < α,
u(x, α) =
0 si x > α.

Elle est discontinue en x = α. Il nous faut maintenant calculer la vitesse de propaga-

tion de la discontinuité. D’après (2.27), on doit avoir : ξ ′(t) = 12 . La fonction u donnée
par :

 t−α
 1 si x < α + ,
u(x, t) = 2
 0 si x > α + t − α .

2
est bien solution du problème pour t > α.

On a représenté ci-dessous les caractéristiques, et la solution à différents temps.

2.3 Solutions faibles - Solutions entropiques 49

t ligne de choc u
t=
t=
t=

t=

x x

2.3.3 Solutions entropiques

Non unicité de la solution faible - un exemple

Grâce à la notion de solution faible, nous avons calculé une solution globale au pro-
blème de Cauchy (2.11-2.12) dans des situations où il n’existait pas de solution clas-
sique. La question qui se pose alors est de savoir si cette solution faible globale est
unique.

L’exemple suivant montre que ce n’est malheureusement pas le cas.

Considérons l’équation de Bürgers pour la donnée initiale suivante :

(
0 si x < 0,
u0 (x) =
1 si x > 0.

Nous pouvons alors calculer une solution continue au problème de Cauchy : pour
cela, il suffit de considérer l’exemple section 2.2.1 et de faire tendre α vers 0. On ob-
tient : 

 0 si x ≤ 0,


x
u(x, t) = si 0 ≤ x ≤ t,

 t

 1 si x ≥ t,

On vérifie aisément a postériori que u est solution du problème, car c’est une solu-
tion classique de l’équation de Bürgers là où elle est C 1 .
50 C HAPITRE 2 : Problèmes hyperboliques non linéaires :

Mais en suivant la démarche de l’exemple section 2.3.2, nous pouvons aussi calculer
une solution discontinue. La vitesse de propagation du choc est donnée par la rela-
tion ξ ′ (t) = 12 . La solution qui en résulte est donc :

 t
 0 si x < ,
u(x, t) = 2
 1 si x > t .

2

t t
choc

x x

Solution continue Solution discontinue

R EMARQUE 2.3.2
Il est intéressant de noter qu’une condition initiale discontinue peut donner nais-
sance à une solution continue. C’est encore un phénomène purement non-linéaire.

Nous avons pu construire une deuxième solution au problème en introduisant une

discontinuité. En fait, on peut même construire une infinité de solutions si l’on ad-
met trois courbes de discontinuité.

En effet, pour tout γ < 1, la fonction suivante est solution du problème :

 γ

 0 si x < t,

 2



 γ t
 γ si t<x< ,
u(x, t) = 2 2

 t γ

 1 − γ si < x < (1 − ) t,

 2 2

 γ
 1 si x > (1 − ) t.
2
En fait, la solution physique est la solution continue.

Passage à la limite dans l’équation dissipative

Nous allons maintenant présenter un critère qui permettra de trier les solutions faibles
pour extraire l’unique solution physique du problème.
2.3 Solutions faibles - Solutions entropiques 51

Comme nous l’avions signalé, l’équation

∂u ∂

+ f (u) = 0 (2.29)
∂t ∂x

est généralement obtenue comme simplification de l’équation visqueuse ou dissi-

pative :
∂uε ∂
∂ 2 uε
+ f (uε ) − ε 2 = 0, (2.30)
∂t ∂x ∂x
où ε est un petit paramètre strictement positif.

Nous admettrons que le problème de Cauchy constitué de l’équation (2.30) et de la

condition initiale
u(x, 0) = u0 (x) p.p. x ∈ R (2.31)

où u0 ∈ L∞ (R), possède une solution unique uε qui appartient à L∞ (R × [0, +∞[). De

plus, uε est très régulière sur R×]0, +∞[.

Nous souhaitons donc caractériser la solution physique de (2.29-2.31), comme limite

de uε lorsque ε tend vers 0.

L EMME 2.3.4
Soit (uε )ε>0 une suite de solutions de (2.30) telles que :

kuε kL∞ (R×[0,+∞[) ≤ C ∀ε > 0 (2.32)

et
uε −→ u quand ε −→ 0 presque partout dans R × [0, +∞[, (2.33)
alors u est solution au sens des distributions de l’équation (2.29).

D ÉMONSTRATION : En utilisant (2.32), (2.33) et le théorème de convergence dominée de Le-

besgue, on vérifie aisément que uε → u et que f (uε ) → f (u), au sens des distributions, c’est-
à-dire dans D ′ (R × [0, +∞[).

Comme la dérivation est une opération continue dans l’espace des distributions, on en dé-
duit que :
∂uε ∂u ∂ 2 uε ∂2u ∂
∂

→ , 2
→ 2
et f (uε ) → f (u)
∂t ∂t ∂x ∂x ∂x ∂x
dans D ′ (R × [0, +∞[). Le lemme s’en déduit. 

Pour caractériser la limite u des solutions de l’équation visqueuse, on utilise un cri-

tère qui repose sur la notion mathématique d’entropie. Précisons cela.
52 C HAPITRE 2 : Problèmes hyperboliques non linéaires :

D ÉFINITION 2.3.3 (E NTROPIE )

On appelle entropie pour l’équation (2.29) tout couple (U, F ) de fonctions C 1 de R
dans R telles que :

(i) U est une fonction strictement convexe,

(ii) F ′ (u) = U ′ (u) f ′ (u) ∀u ∈ R.

Dès que U(·) et F (·) sont deux fonctions régulières relièes par la relation (ii) ci-dessus,
si u est une solution classique de l’équation (2.29) on a :
∂u ∂u
U ′ (u) + U ′ (u)f ′ (u) =0
∂t ∂x
d’où :
∂ ∂
U(u) + F (u) = 0. (2.34)
∂t ∂x
En revanche, si u est une solution faible de (2.29), par exemple C 1 par morceaux, elle
ne satisfait pas en général l’équation (2.34) au sens des distributions. En effet, il fau-
drait pour cela que la relation s[U(u)] = [F (u)] soit satisfaite le long de toute courbe de
discontinuité de u. Cela est généralement incompatible avec la relation de Rankine-
Hugoniot.

Par contre, lorsque U est convexe, on peut montrer que, si la solution faible u est
construite par passage à la limite sur l’équation avec viscosité, l’égalité (2.34) devient
une inégalité au sens des distributions. Plus précisément, on peut démontrer le théo-
rème suivant.

T HÉORÈME 2.3.5
Soit (uε )ε>0 une suite de solutions régulières de (2.30) vérifiant (2.32) et (2.33), et
soit (U, F ) une entropie pour l’équation (2.29). Alors u vérifie
∂ ∂
U(u) + F (u) ≤ 0 (2.35)
∂t ∂x
au sens des distributions.

D ÉMONSTRATION :
Comme uε est une solution régulière de (2.30), on a :
∂uε ∂uε ∂ 2 uε
U ′ (uε ) + U ′ (uε )f ′ (uε ) − εU ′ (uε ) 2 = 0,
∂t ∂x ∂x
qui s’écrit aussi :
 2
∂ ∂ ∂2 ∂uε
U (uε ) + F (uε ) − ε 2 U (uε ) = −εU ′′ (uε ) .
∂t ∂x ∂x ∂x
2.3 Solutions faibles - Solutions entropiques 53

Par des arguments similaires à ceux de la démonstration du lemme 2.3.4, on a :

∂ ∂ ∂ ∂ ∂2 ∂2
U (uε ) → U (u), F (uε ) → F (u) et U (u ε ) → U (u) dans D ′ (R×[0, +∞[).
∂t ∂t ∂x ∂x ∂x2 ∂x2

 2
′′ ∂uε
En revanche, le terme εU (uε ) ne peut être traité de façon similaire.
∂x

Malgré la présence du facteur ǫ, ce terme ne tend pas nécessairement vers 0 au sens des dis-
tributions. Cela est dû au fait que les hypothèses n’empêchent pas ∂u
∂x de tendre vers +∞.
ε

Cependant, comme U est convexe, on a :

D  2 E
′′ ∂uε
− εU (uε ) ,ϕ ≤0
∂x

pour toute fonction ϕ positive de D(R × [0, +∞[). A la limite, on obtient donc :

D ∂ ∂ E
U (u) + F (u), ϕ ≤ 0
∂t ∂x

pour toute fonction ϕ positive de

D(R × [0, +∞[).

La condition (2.35) est appelée condition d’entropie.

On peut maintenant définir la notion de solution entropique.

D ÉFINITION 2.3.4 (S OLUTION ENTROPIQUE )

Une solution faible u du problème de Cauchy (2.29-2.31) est appelée solution
entropique si elle satisfait la condition d’entropie (2.35) pour toute entropie (U, F )
de l’équation (2.29).

Solutions entropiques C 1 par morceaux. Chocs entropiques

Pour les solutions C 1 par morceaux etudiées au paragraphe 2.3.2, la recherche de so-
lutions entropiques va limiter les chocs admissibles : ce sont les chocs entropiques.
54 C HAPITRE 2 : Problèmes hyperboliques non linéaires :

T HÉORÈME 2.3.6 (E NTROPIE POUR LES FONCTIONS C 1 PAR MORCEAUX )

Soit u une fonction C 1 par morceaux dans R × [0, +∞[ solution faible de (2.29-2.31),
c’est à dire qui vérifie :

• u est une solution classique du problème de Cauchy (2.29-2.31) là où elle est

C 1.

• u vérifie le long de toute courbe de discontinuité Σ :

ξ ′ (t)(u+ − u− ) = f (u+ ) − f (u− ).

Alors u est une solution entropique si et seulement si, pour toute entropie (U, F ),
le long de toute courbe de discontinuité , Σ on a





ξ ′ (t) U u+ (t) − U u− (t) ≥ F u+ (t) − F u− (t) . (2.36)

D ÉMONSTRATION : Soit u une fonction C 1 par morceaux et (U, F ) une entropie pour l’équation
(2.29). En reprenant les notations et la démarche du théorème 2.3.2, on vérifie aisément que,
pour tout ϕ ∈ D(K) :
Z   Z  
∂ϕ ∂ϕ ∂ ∂
U (u) + F (u) dx dt = − U (u) + F (u) ϕ dx dt
K ∂t ∂x K − ∪K + ∂t ∂x
Z (2.37)
 + −

+ −


− U (u ) − U (u ) nt + F (u ) − F (u ) nx ϕ dγ .
Σ
Soit u une solution faible du problème satisfaisant la condition d’entropie (2.35). On a :
Z   Z  
∂ϕ ∂ϕ ∂ ∂
U (u) + F (u) dx dt ≥ 0 et U (u) + F (u) ϕ dx dt = 0
K ∂t ∂x K − ∪K + ∂t ∂x
pour toute fonction ϕ positive de D(K). D’après (2.37), on a, en utilisant l’expression exacte
de la normale,


U (u+ ) − U (u− ) (−ξ ′ (t)) + F (u+ ) − F (u− ) ≤ 0
au sens des distributions le long de Σ. Autrement dit, on a :



ξ ′ (t) U (u+ ) − U (u− ) ≥ F (u+ ) − F (u− ) dans D ′ .


Le critère (2.36) est relativement peu pratique pour vérifier qu’un choc est entro-
pique ou non. Heureusement, on peut en donner une forme équivalente plus simple,
à commencer par le cas où f est convexe.

T HÉORÈME 2.3.7 (C HOC ENTROPIQUE . C AS CONVEXE )

On suppose que f est strictement convexe. Alors la condition de choc entropique
est équivalente à

u− > u+ le long de toute courbe de discontinuité Σ. (2.38)

2.3 Solutions faibles - Solutions entropiques 55

Avant de démontrer ce théorème, nous allons montrer que la condition (2.38) signifie
que lorsqu’il y a choc, les caractéristiques doivent converger vers la ligne de choc et
non s’en écarter. On rappelle que c désigne la dérivée de f.

L EMME 2.3.8
Soit u une solution entropique du problème de Cauchy (2.29-2.31), qui soit C 1 par
morceaux dans R ×[0, +∞[. Alors, u vérifie le long de toute courbe de discontinuité
Σ:
c(u+ ) < ξ ′ (t) < c(u− ).

D ÉMONSTRATION : Comme f est strictement convexe, la fonction :

f (u) − f (u+ )
u 7→
u − u+
est continue et strictement croissante sur [u+ , u− ]. Il en résulte que :

f (u− ) − f (u+ )
c(u+ ) < ·
u− − u+
Mais, par la relation de Rankine-Hugoniot, ceci s’écrit : c(u+ ) < ξ ′ (t). On démontre de même
la seconde inégalité. 

R EMARQUE 2.3.3 (C AS PARTICULIER DE L’ ÉQUATION DE BÜRGERS )

Dans le cas de l’équation de Bürgers, on a tout simplement :

u− + u+
u+ < s = < u− .
2
La condition d’entropie permet donc d’éliminer les solutions discontinues non phy-
siques de l’exemple section 2.3.3.

D ÉMONSTRATION DU THÉORÈME 2.3.7 : Nous partons de la fin de la démonstration du théo-

rème 2.3.6. Considérons l’application G suivante :

f (u) − f −


u 7→ G(u) = −
U (u) − U (u− ) − F (u) − F (u− ) ·
u−u
Nous allons montrer que G est une fonction décroissante de u. Admettons-le pour l’instant.
Comme on a G(u− ) = 0 et que (2.36) s’écrit G(u+ ) ≥ 0, il résulte de la décroissance de G que
les conditions (2.36) et (2.38) sont équivalentes, ce qui démontre le théorème.

Prouvons maintenant que G est décroissante. On a :

′ f ′ (u) (u − u− ) − f (u) − f − − f (u) − f − ′
G (u) = (U (u) − U (u )) + U (u) − F ′ (u).
(u − u− ) u − u−
56 C HAPITRE 2 : Problèmes hyperboliques non linéaires :

Mais comme F ′ (u) = U ′ (u)f ′ (u), ceci s’écrit aussi

′
f ′ (u) (u − u− ) − f (u) − f − − f (u) − f − − f ′ (u) (u − u− ) ′
G (u) = (U (u) − U (u )) + U (u)
(u − u− )2 u − u−



f ′ (u)(u − u− ) − f (u) + f − U (u) − U (u− ) − U ′ (u)(u − u− )
= .
(u − u− )2

Or f et U sont convexes, d’où :

f ′ (u)(u − u− ) − (f (u) − f − ) < 0 et U ′ (u)(u − u− ) − (U (u) − u− ) < 0, ∀u ∈ R,

ce qui achève la démonstration. 

R EMARQUE 2.3.4 (E NTROPIE ET DYNAMIQUE DES GAZ )

En dynamique des gaz, la condition d’entropie se déduit du second principe de la
thermodynamique. En effet, ce principe implique que l’entropie augmente à travers
un choc, de l’amont vers l’aval de l’écoulement. Il en résulte que la vitesse normale
relative du fluide diminue de l’amont vers l’aval. C’est ce qui justifie l’appellation
« condition d’entropie ».

Nous donnons maintenant la version générale du théorème 2.3.7. Sa démonstration

est toutefois plus délicate que celle du théorème 2.3.7 et sera admise ici.

T HÉORÈME 2.3.9 (C HOC ENTROPIQUE . C AS GÉNÉRAL )

La condition de choc entropique (2.36) est équivalente à ∀α ∈ [0, 1]


 f αu− + (1 − α)u+ ≥ αf (u− ) + (1 − α)f (u+ ) si u+ > u− ,
(2.39)
 f αu− + (1 − α)u+
≤ αf (u− ) + (1 − α)f (u+ ) si u+ < u− .

R EMARQUE 2.3.5 (I NTERPRÉTATION GÉOMÉTRIQUE )

La condition s’interprète facilement d’un point de vue géométrique. Un choc est en-
tropique si et seulement si une des deux conditions suivantes est réalisée :

(i) u+ > u− et le graphe de u 7→ f (u) est au dessus de sa corde sur le segment

[u− , u+ ],

(ii) u+ < u− et le graphe de u 7→ f (u) est en dessous de sa corde sur le segment

[u− , u+ ].

Le graphe d’une fonction convexe étant toujours en dessous de sa corde et on re-

tombe bien sur la condition u+ < u− .
2.3 Solutions faibles - Solutions entropiques 57

2.3.4 Existence et unicité de la solution entropique

D’un point de vue mathématique, le gros avantage de la notion de solution entro-
pique est qu’elle permet d’obtenir un résultat d’unicité. Pour conclure , nous allons
énoncer, sans le démontrer, le résultat fondamental de la théorie des lois de conser-
vation scalaire qui assure l’existence et l’unicité d’une solution entropique au pro-
blème de Cauchy.

T HÉORÈME 2.3.10 (E XISTENCE ET UNICITÉ DE LA SOLUTION ENTROPIQUE )

Soit u0 ∈ L∞ (R). Alors le problème de Cauchy (2.11-2.12) admet une solution en-
tropique unique u ∈ L∞ (R × [0, +∞[) qui vérifie l’estimation

ku(·, t)kL∞(R) ≤ ku0 kL∞ (R) p.p. t ≥ 0.

Si u et v sont les solutions entropiques associées aux conditions initiales u0 et v 0 ,

on a Z Z
x2 x2 +At
|u(x, t) − v(x, t)|dx ≤ |u0 (x) − v 0 (x)|dx (2.40)
x1 x1 −At

pour presque tout t ≥ 0, ∀x1 < x2 , où



A = max |f ′ (ξ)| / |ξ| ≤ max ku0 kL∞ (R) , kv 0 kL∞ (R) .

R EMARQUE 2.3.6 (V ITESSE DE PROPAGATION )

L’inégalité (2.40) signifie que la valeur de u au point (x, t) ne dépend que des valeurs
de u0 dans l’intervalle [x − At, x + At]. Autrement dit, la solution entropique a une
vitesse de propagation finie.

R EMARQUE 2.3.7 (A PROPOS DE LA DÉMONSTRATION DU THÉORÈME )

La démonstration de ce théorème est longue et difficile et fait appel à des outils et
des arguments mathématiques qui dépassent largement le cadre de ce cours. Nous
renvoyons à GODLEWSKI-RAVIART1 pour la démonstration détaillée de ce théorème.
Signalons toutefois que

• La démonstration du résultat d’existence est constructive et fait appel à la tech-

nique de viscosité. Les principales étapes sont les suivantes (pour une donnée
u0 assez régulière) :

– On démontre tout dabord, pour tout ǫ > 0, l’exisence et l’unicité de la so-

lution (régulière) de (2.11-2.12).
– A l’aide d’estimation a priori et de techniques de compacité, on montre
que cette solution uǫ satisfait les hypothèses (2.32) et (2.33) du lemme 2.3.4.
1
E. GODLEWSKI, P.A. RAVIART : Hyperbolic systems of conservation laws, Mathématiques et Ap-
plications n◦ 3 − 4, Publication de la SMAI, Ed. Ellipses, 1991.
58 C HAPITRE 2 : Problèmes hyperboliques non linéaires :

– On utilise le lemme 2.3.4 et le théorème pour passer à la limite et montrer

que la limite u est bien solution entropique de (2.30-2.12).
• Le résultat d’unicité est bien entendu une conséquence immédiate du résultat
de « stabilité L1 » (2.40). Il n’est pas très difficile de voir que (2.40) découle de
l’inégalité au sens des distributions (pour tout couple (u, v) de solutions entro-
piques, sgn(·) désigne la fonction « signe »)
∂ ∂ 

|u − v| + sgn(u − v) f (u) − f (v) ≤ 0, dans D ′ (R×]0, +∞[)
∂t ∂x
dont la démonstration, qui est elle beaucoup plus délicate, fait appel aux entro-
pies de Kruzkov Uk (x) = |u − k|, k ∈ R.

T HÉORÈME 2.3.11 (M ONOTONIE )

Si u et v sont les solutions entropiques de (2.11-2.12) associées aux conditions ini-
tiales u0 et v 0 , on a
u0 ≥ v 0 ⇒ u(·, t) ≥ v(·, t) p.p. t ≥ 0 . (2.41)
On en déduit en particulier que

 x → u0 (x) croissante ⇔ x → u(x, t) croissante, p. p. t > 0,
(2.42)
 a ≤ u0(x) ≤ b p. p. x ∈ R ⇔ a ≤ u(x, t) ≤ b p. p. (x, t) ∈ R × R+ .

Nous invitons le lecteur à démontrer (2.42) à partir de (2.41).

2.4 Le problème de Riemann

Dans cette section, nous nous limitons pour simplifier au cas où f est strictement
convexe. Dans les cinq premières parties de cette section, nous considérons le pro-
blème de Riemann à 2 états. Dans la dernière partie, nous nous intéresserons au pro-
blème de Riemann à 3 états et on utilisera les résultats établis dans les sections pré-
cédentes.

2.4.1 Présentation du problème

Pour illustrer les résultats obtenus dans ce chapitre, nous allons calculer la solution
entropique du problème de Cauchy suivant :

 ∂u ∂

 + f (u) = 0 x ∈ R, t > 0
 ∂t ∂x (
ug si x < 0, (2.43)


 u(x, 0) =

ud si x > 0,
2.4 Le problème de Riemann 59

où ug et ud sont deux constantes données.

R EMARQUE 2.4.1 (L IEN AVEC LA DYNAMIQUE DES GAZ )

On dit qu’il s’agit d’un problème de Riemann par analogie à un problème étudié par
ce dernier en dynamique des gaz : on considère un gaz dans un tube cylindrique très
fin et très long. On suppose que le tube est séparé en deux par une membrane et que
le gaz qui est au repos n’a pas les mêmes pression et densité des deux côtés de la
membrane. A l’instant t = 0, la membrane est déchirée et l’on s’intéresse à l’écoule-
ment qui s’en suit.

Nous avons déjà calculé la solution du problème (2.43) pour l’équation de Bürgers
dans les cas suivants :

• ug = 0 et ud = 1 : c’est l’exemple de la section 2.2.1,

• ug = 1 et ud = 0 : c’est l’exemple de la section 2.3.2,

et nous avons pu remarquer que les solutions obtenues étaient de natures très diffé-
rentes.
Pour l’étude du problème général également, il nous faudra distinguer deux cas sui-
vant que ug < ud ou ug > ud .

2.4.2 Une solution auto-semblable

Avant tout, nous allons établir le lemme suivant.

L EMME 2.4.1 (S OLUTION AUTO - SEMBLABLE )

La solution u du problème de Riemann (2.43) est auto-semblable, c’est-à-dire
qu’elle est de la forme
u(x, t) = v (x/t) .
De plus, la fonction v doit vérifier, là où elle est C 1 ,


v ′ (ξ) = 0 ou c v(ξ) = ξ. (2.44)

D ÉMONSTRATION : Soit α > 0. Posons : x̄ = αx, t̄ = αt et ū(x̄, t̄) = u(x, t). Il est clair que ū est
solution tout comme u du problème (2.43). Autrement dit, pour tout α > 0 :

u(x, t) = u(αx, αt).

Ceci signifie exactement que u est auto-semblable. De plus, si u est C 1 , on doit avoir :

∂u ∂u x 1 
+ c(u) = − 2 v ′ (x/t) + c v (x/t) v ′ (x/t) = 0
∂t ∂x t t
60 C HAPITRE 2 : Problèmes hyperboliques non linéaires :

D’où :  



v ′ (x/t) c v x/t − x/t = 0.

Nous pouvons maintenant résoudre explicitement le problème.

2.4.3 Le cas d’une onde de détente

On se place dans le cas où ug < ud . Dans ce cas, tout comme pour l’exemple de la
section 2.2.1, la solution entropique est continue. La méthode des caractéristiques
permet de construire la solution dans les domaines x ≤ c(ug )t et x ≥ c(ud )t. On a :
 x
 ug si ≤ c(ug )
u(x, t) = t .
 u si x ≥ c(u )
d d
t

Pour calculer u dans le domaine c(ug )t ≤ x ≤ c(ud )t, on résout l’équation (2.44). C’est
toujours possible car f est strictement croissante, donc a est inversible. On obtient :
x x
−1
u(x, t) = c si c(ug ) ≤ ≤ c(ud ).
t t

La fonction u ainsi construite est effectivement continue.

t
x x
- = a(ug) - = a(ud) u(.,t)
t t

ud

ug
x x
a(ug)t a(ud)t

On dit qu’il s’agit d’une onde de détente.

2.4 Le problème de Riemann 61

2.4.4 Le cas d’une onde de choc

On se place dans le cas où ug > ud . Dans ce cas, tout comme pour l’exemple de la sec-
tion 2.3.2, la solution est constante de part et d’autre d’un choc. La vitesse de propa-
gation du choc est donnée par la relation de Rankine-Hugoniot (elle est constante) :
f (ug ) − f (ud)
s= .
ug − ud
Le choc est entropique puisque ug > ud . La solution est
(
ug si x < s t,
u(x, t) =
ud si x > s t.

t x=st
u(.,t)

ud

ug
x x
st

2.4.5 Solution du problème de Riemann à 2 états

En résumé, la solution entropique du problème de Riemann (2.43) peut être notée :
x 
wR ; ug , ud
t
- Si u < v : 

 u si ξ ≤ c(u)

wR (ξ; u, v) = c−1 (ξ) si c(u) ≤ ξ ≤ c(v)



v si ξ ≥ c(v)
- Si u > v : (
u si ξ < s, f (u) − f (v)
wR (ξ; u, v) = où s= .
v si ξ > s, u−v
62 C HAPITRE 2 : Problèmes hyperboliques non linéaires :

2.4.6 Le problème de Riemann à trois états

Nous nous proposons de résoudre explicitement le problème de Cauchy suivant (contrai-
rement aux paragraphes précédents, nous nous restreignons ici au cas de l’équation
de Bürgers)   2

 ∂u ∂ u

 + = 0 x ∈ R, t > 0,

 ∂t ∂x 2

 

 u si x < 0,
  1

 u(x, 0) = u2 si 0 < x < 1,



 

  u si x > 1,
3

où u1 , u2 et u3 sont trois constantes données.

Cas d’une condition initiale croissante

On suppose u1 < u2 < u3 . Dans ce cas, la donnée initiale étant croissante, il existe
une solution continue et croissante pour tout temps t > 0. Les deux détentes issues
des points x = 0 et x = 1 n’interfèreront jamais. La solution est donnée par :

 u1 si x < u1 t,





 x/t si u1 t < x < u2 t,


u(x, t) = u2 si u2 t < x < u2 t + 1,





 (x − 1)/t si u2 t + 1 < x < u3 t + 1,



u3 si u3 t + 1 < x.
Application numérique : u1 = −1, u2 = 0, u3 = 1.

t
u(.,t)

1 t=0

t=1

x
1
1

-1
x
1
2.4 Le problème de Riemann 63

Cas d’une condition initiale décroissante

On suppose u1 > u2 > u3 . A l’instant t = 0 naissent deux chocs entropiques issus

des points x = 0 et x = 1. Ils se rencontreront à un temps ultérieur t∗ qu’il faudra
déterminer. La trajectoire du premier choc x = ξ1 (t) est donnée, pour t proche de 0,
par la relation de Rankine-Hugoniot :

1
ξ1 (t) = (u1 + u2 ) t
2

et celle du second choc x = ξ2 (t) par :

1
ξ2 (t) = (u2 + u3 ) t + 1.
2

Le temps t∗ satisfait donc l’équation suivante :

1 1
(u1 + u2 ) t∗ = (u2 + u3 ) t∗ + 1
2 2

d’où :
2
t∗ = .
u1 − u3

A cet instant, on a :
u1 + u2
ξ1 (t) = ξ2 (t) = .
u1 − u3

Au delà du temps t∗ , il ne reste qu’un seul choc dont la trajectoire est donnée par :

1 u2 − u3
x = ξ(t) = (u1 + u3 ) t + .
2 u1 − u3

La solution est finalement donnée par



 u si x < ξ1 (t),
 1
si t < t∗ , u(x, t) = u2 si ξ1 (t) < x < ξ2 (t),


 u si ξ (t) < x,
3 2
(
u1 si x < ξ(t),
si t > t∗ , u(x, t) =
u3 si ξ(t) < x,
64 C HAPITRE 2 : Problèmes hyperboliques non linéaires :

Application numérique : u1 = 1, u2 = 0, u3 = −2.

t
u(.,t)
t=0
t=1/3
t=2/3
1

1 1 x

-2
x
1

Cas d’une condition initiale non monotone

On suppose u1 < u2 et u2 > u3 . Le dernier cas (u1 > u2 et u2 < u3 ) est laissé en
exercice au lecteur ! ! ! On observe donc une détente issue du point x = 0 et un choc
entropique issu du point x = 1. Ce choc et cette détente vont interférer à partir d’un
temps t∗ que l’on déterminera. Pour t < t∗ , la trajectoire du choc est donnée par :
1
x = ξ(t) où ξ(t) = (u2 + u3 ) t + 1.
2
Par conséquent, on a, pour t < t∗ :


 u1 si x < u1 t,



 x/t si u1 t < x < u2 t,
u(x, t) =

 u2 si u2 t < x < ξ(t),




u3 si x > ξ(t).
Le temps t∗ satisfait donc l’équation :
u2 t∗ = ξ (t∗ )
d’où :
2
t∗ = .
u2 − u3
Pour t > t∗ , la trajectoire du choc x = ξ(t) est influencée par la détente. En effet, la
valeur de la solution à gauche du choc est :
ξ(t)
u(ξ(t), t) = ,
t
2.4 Le problème de Riemann 65

et par conséquent, la relation de Rankine-Hugoniot s’écrit :

 
′ 1 ξ(t)
ξ (t) = + u3 .
2 t

La solution générale de cette équation différentielle est donnée par :

√
ξ(t) = α t + u3 t.

Comme de plus, on sait que :

2u2
ξ (t∗ ) =
u2 − u3
on trouve : p
α= 2 (u2 − u3 ).
Pour t > t∗ , la solution se prolonge donc par :


 u1 si x < u1 t,

u(x, t) = x/t si u1 t < x < ξ(t),



u3 si x > ξ(t).

Par la suite, de deux choses l’une. Soit u1 < u3 et les formules ci-dessus restent va-
lables pour tout temps t > t∗ . Soit u1 > u3 et ces formules restent valables tant que

u1 t < ξ(t),

c’est-à-dire pour t < t∗∗ où t∗∗ est solution de

p
u1 t∗∗ = 2 (u2 − u3 ) t∗∗ + u3 t∗∗ .

On a donc :
2 (u2 − u3 )
t∗∗ = .
(u1 − u3 )2
Pour t > t∗∗ , la trajectoire du choc est donnée par :

1 u2 − u3
ξ(t) = (u1 + u3 ) t +
2 u1 − u3

et la solution par :
(
u1 si x < ξ(t),
u(x, t) =
u3 si x > ξ(t).
66 C HAPITRE 2 : Problèmes hyperboliques non linéaires :

Application numérique : u1 = 0, u2 = 1, u3 = −1.

t
u(.,t)

x
1 2
-1
1 t=0
t=1
t=3/2
x
1

Le cas u1 = −1, u2 = 2, u3 = 0

t
u(.,t)
t=0
t=1 1
t=2
t=5
x
1

1 -1

x
1

Le cas u1 = 0, u2 = 1, u3 = −1
Chapitre 3

Approximation par différences finies -

Application aux problèmes
hyperboliques linéaires

3.1 Présentation de la méthode de différences finies

3.1.1 Approximation de quelques opérateurs différentiels
La méthode de différences finies permet de calculer une approximation de la solu-
tion d’une EDP en des points qui sont distribués sur une grille. L’objectif est alors de
construire des approximations des dérivées des fonctions intervenant dans l’EDP à
l’aide de valeurs discrètes de celle-ci, par le biais de formules de Taylor.

Nous nous plaçons, sauf mention contraire, dans le cadre d’une discrétisation par
différences finies en une dimension d’espace. Ceci implique que nous discrétisons le
continuum spatio-temporel par une grille régulière de pas ∆t en temps et h en espace
tel que les coordonnées discrètes soient :
(xj , tn ) = (jh, n∆t).
La solution discrète sera calculée en ces points (voir la figure 3.1). Le principe des
différences finies est donc de remplacer les dérivées en ces points par des différences
finies en utilisant des formules de Taylor dont on néglige les restes.

Pour commencer, reprenons la définition usuelle de la dérivée

u(x + h) − u(x)
u′ (x) := lim ,
h→0 h
afin d’introduire l’approximation décentrée à droite
u(x + h) − u(x)
u′ (x) ≃ Dh+ u(x) := .
h
67
68 C HAPITRE 3 : Approximation par différences finies

bc tn+1
bc × bc tn
bc tn−1
∆t

h xj−1 xj xj+1

F IGURE 3.1 : Grille espace temps

Evidemment, il est, a priori, tout aussi légitime de définir une approximation décen-
trée à gauche
u(x) − u(x − h)
u′ (x) ≃ Dh− u(x) := ,
h
ou encore une approximation centrée
u(x + h) − u(x − h)
u′ (x) ≃ Dh0 u(x) := .
2h
Nous allons voir que ce choix n’est en fait pas anodin et que certains choix a priori
raisonnables ne permettent pas une bonne approximation (dans un sens qu’il nous
faudra définir) de la solution de l’EDP initiale.

La première obligation pour notre approximation est la consistance. Elle exprime

le fait que l’opérateur aux différences finies converge vers l’opérateur différentiel
lorsque le pas de maillage tend vers 0. Ainsi, pour un maillage de pas h, on note Dh un
opérateur d’approximation aux différences finies de dx d
et on appelle erreur de tron-
cature pour cet opérateur la quantité
ε(x, h) := Dh u(x) − u′ (x).
L’approximation est consistante pour la norme k·k, si pour u suffisamment régulière
lim kε(x, h)k = 0.
h→0

On peut même aller plus loin en disant que l’approximation est d’ordre k ∈ N si
ε(x, h) = O(hk ).
En pratique, l’analyse de consistance s’effectue en considérant les développements
de Taylor de la fonction u. Par exemple pour l’approximation décentrée à droite on
sait pour u ∈ C 2 (I), qu’il existe θ ∈]0, 1[ tel que
h2 ′′
u(x + h) = u(x) + h u′(x) + u (x + θh).
2
3.1 Présentation de la méthode de différences finies 69

On obtient ainsi
u(x + h) − u(x) h
− u′ (x) = u′′ (x + θh),
h 2
d’où l’estimation
h
ε(x, h) ≤ sup |u′′|. (3.1)
2
On aurait obtenu le même type d’inégalité pour Dh− donc les approximation décen-
trées à gauche et à droite sont d’ordre 1. En pratique, on voit que l’inégalité (3.1) nous
indique même que le choix de h doit être déterminé par la valeur de la dérivée se-
conde.

Pour l’approximation centrée, les calculs sont un peu plus longs mais pour un bé-
néfice sur l’ordre de l’approximation. On considère u ∈ C 3 (I). Il existe θ+ ∈]0, 1[ tel
que
h2 ′′ h3 (3)
u(x + h) = u(x) + h u′(x) + u (x + θh) + u (x + θ+ h)
2 6
et θ− ∈]0, 1[ tel que

h2 ′′ h3 (3)
u(x − h) = u(x) − h u′ (x) + u (x + θh) − u (x + θ− h)
2 6
ce qui donne par soustraction

u(x + h) − u(x − h) ′ h2  (3) + (3) −


− u (x) = u (x + θ h) + u (x + θ h) .
2h 3
D’après le théorème des valeurs intermédiaires, il existe θ ∈] − 1, 1[ tel que

u(x + h) − u(x − h) h2 (3)

− u′ (x) = u (x + θh),
2h 6
et l’approximation centrée est donc d’ordre 2.

Ces trois approximations sont loin d’être les seules et il existe en fait une infinité de
possibilités d’approximation dont on jugera la pertinence notamment en fonction de
l’erreur de consistance. Il est en particulier possible de faire intervenir plus de points
dans Dh afin d’obtenir une approximation plus précise. Par exemple si nous combi-
nons les approximations suivantes (obtenus comme précédemment mais en choi-
sissant désormais les formules de Taylor-Young au lieu de celles de Taylor-Lagrange),
nous obtenons
u(x + h) − u(x − h) h2 (3)
= u′(x) + u (x) + O(h4 ),
2h 6
et
u(x + 2h) − u(x − 2h) h2
= u′ (x) + 4 u(3) (x) + O(h4 ),
4h 6
70 C HAPITRE 3 : Approximation par différences finies

ce qui donne pour tout coefficient de pondération α ∈ [0, 1]

(4) u(x + h) − u(x − h) u(x + 2h) − u(x − 2h)

Dh u = α + (1 − α)
2h 4h
′

h2 (3)
= u (x) + α + 4(1 − α) u (x) + O(h4 ).
6
Ainsi pour un choix adéquate α = 43 de pondération entre les deux approximations
initiales, on obtient une approximation d’ordre 4.

Afin de présenter les différents choix d’opérateurs d’approximation Dh , on adoptera

une représentation symbolique beaucoup plus visuelle comme utilisée Figure 3.2
pour les exemples d’approximation que nous venons de presenter.

Dh+ : 1
-1 1
h

Dh− : 1
-1 1
h

Dh0 : 1
-1 1
2h

Dh
(4)
: 1
1 -8 8 -1
12h

F IGURE 3.2 : Représentation symbolique des opérateurs aux différences finies

De tout évidence la construction d’opérateurs d’approximation ne se limite pas aux

dérivées d’ordre 1. On utilise le même principe pour les dérivées d’ordre supérieur.
Notamment, on a
u′ (x + h/2) − u′ (x − h/2)
u′′ (x) ≃ ,
h
avec par exemple
u′(x + h) − u′ (x − h)
u′ (x + h/2) ≃ ,
h
et
u′ (x + h) − u′(x − h)
u′ (x − h/2) ≃ .
h
On peut ainsi définir un opérateur Dh2 d’approximation aux différences finies de la
dérivée seconde sous la forme
u(x + h) − 2u(x) + u(x − h)
u′′ (x) ≃ Dh2 u(x) := .
h2
3.1 Présentation de la méthode de différences finies 71

La consistance de cette opérateur se détermine toujours à partir de formule de Taylor

(ici Taylor-Lagrange). Il existe θ+ ∈]0, 1[ tel que

h2 ′′ h2 (3) h4 (4)
u(x + h) = u(x) + h u′(x) + u (x) + u (x) + u (x + θ+ h),
2 6 24
et θ+ ∈]0, 1[ tel que

h2 ′′ h2 (3) h4 (4)
u(x − h) = u(x) − h u′ (x) + u (x) − u (x) + u (x + θ+ h).
2 6 24
On a donc l’existence de θ ∈] − 1, 1[ tel que

u(x + h) − 2u(x) + u(x − h) h2 (4)

Dh2 u(x) = ′′
= u (x) + u (x + θh), (3.2)
h2 12
et l’approximation est donc d’ordre 2.

3.1.2 La démarche : Etude de la consistance et la stabilité

Pour une équation aux dérivées partielles donnée, l’idée est donc de construire une
équation discrète où les opérateurs différentiels sont remplacés par des opérateurs
d’approximation. La solution continue n’a a priori pas de raison d’être solution de
cette équation discrète.

La question qui se pose donc naturellement est de savoir si, pour une équation dis-
crète particulière (on parlera également de schéma) obtenue en utilisant la méthode
de différences finies décrite ci-dessus, la solution de l’équation discrète obtenue ap-
proche de la solution du problème continu et d’autre part quand il y a convergence,
de quantifier la vitesse de convergence.

Cette analyse se fait toujours en deux étapes fondamentales qui seront donc le dé-
nominateur commun de chaque étude de schéma de cet ouvrage :

• La consistance qui correspond à l’erreur d’approximation liées lorsqu’on rem-

place l’EDP par des différences finies. Il est fondamental de comprendre que,
bien que très intuitive en tant que condition nécessaire, la consistance d’une
méthode ne suffit pas à obtenir la convergence.

• La stabilité qui assure que l’opérateur discret est bien inversible et la norme de
son inverse est bornée indépendamment du pas de discrétisation.

On montrera alors que dès qu’un schéma est consistant et stable alors il est convergent,
c’est à dire que la solution numérique calculée à partir de ce schéma tend, quand les
paramètres de discrétisation tendent vers 0, vers la solution de l’équation aux déri-
vées partielles initiale.
72 C HAPITRE 3 : Approximation par différences finies

L’objectif de ce chapitre étant de présenter et d’analyser quelques schémas classiques

pour les équations hyperboliques linéaires, nous étudierons donc par différentes ap-
proches la consistance et la stabilité des schémas avant d’obtenir le résultat final de
convergence.

3.1.3 Application aux équations hyperboliques linéaires

Avant de passer à la construction de schémas de discrétisation particuliers, indi-
quons dès à présent une particularité de la résolution des problèmes d’évolution :
avec la plupart des méthodes numériques couramment utilisées, la résolution ne
s’effectue pas de la même manière en temps et en espace. En effet il existe classi-
quement deux situations bien distinctes suivant que l’EDP fait intervenir ou non le
temps :

• Si on résout par exemple l’équation de Laplace en dimension 2

∂2u ∂2u
− − 2 =f dans R2 ,
∂x2 ∂y

on peut montrer que la méthode des différences finies conduit à inverser un

système linéaire dans lequel les deux variables x et y jouent le même rôle : d’une
certaine façon, toutes les valeurs de la solution approchée sont calculées simul-
tanément.

• Avec un problème d’évolution en dimension 1 d’espace, bien que la solution

soit à nouveau une fonction de deux variables, le temps et l’espace sont traités
différemment. Cela se comprend intuitivement puisque le temps est naturel-
lement orienté et que nous disposons donc de conditions initiales et pas de
conditions finales alors que les directions x et −x jouent le même rôle. De ce
fait, la solution approchée se calculera pas à pas en temps : la solution à l’ins-
tant tn est déterminée à partir de la solution aux instants précédents. Les condi-
tions initiales sont alors utilisées pour déterminer la solution aux premiers pas
de temps. En pratique, pour les problèmes du premier ordre en temps, comme
c’est le cas de l’équation d’advection par exemple, on pourra se contenter d’un
schéma à un pas de temps ; la solution à l’instant tn est calculée à partir de la
solution à l’instant précédent tn−1 . A un instant donné tn :

– soit les valeurs unj , j ∈ Z sont calculées indépendamment les unes des
autres : on parle alors de schéma explicite,
– soit les valeurs unj , j ∈ Z sont calculées en résolvant un système d’équa-
tions couplées. On parle alors de schéma implicite.
3.1 Présentation de la méthode de différences finies 73

Afin de présenter les différents schémas numériques possibles pour les équations hy-
perboliques, nous allons considérer les deux exemples principaux que sont :

• l’équation d’advection sur R avec une vitesse constante c


 ∂u ∂u
+c = 0, (x, t) ∈ R×]0, T [,
∂t ∂x (3.3)
u(0, x) = u0 (x), x ∈ R.

• l’équation des ondes sur R avec une vitesse constante c

 2 2
 ∂ u 2∂ u

 − c = 0, pour t ≥ 0, x ∈ R

 ∂t2 ∂x 2


u(t = 0, x) = u0 (x), pour x ∈ R





 ∂u (t = 0, x) = u1 (x),

pour x ∈ R
∂t
Dans toute la suite et sauf mention contraire, nous supposons la vitesse de propaga-
tion c positive.

Comme annoncé précédemment, nous approchons la solution u de ces EDPs par une
solution discrète unj . On cherche alors à calculer pour chaque n ( n ∈ N, n∆t ≤ T ) un
vecteur un = (unj )j∈Z tel que unj soit une approximation de la valeur de la solution u du
problème précédent en xj à l’instant tn .

R EMARQUE 3.1.1
Il peut sembler surprenant de recourir à une résolution numérique approchée pour
ces problèmes où, la vitesse étant constant, on connaît analytiquement les solutions
exactes. L’intérêt est de calculer dans des cas simples les propriétés des différents
schémas qui seront ensuite généralisés au cas non linéaire. Un schéma qui ne permet
pas de résoudre le cas linéaire sera automatiquement à proscrire dans le cas non
linéaire.
Nous commencerons par l’équation d’advection et présenterons le schéma expli-
cite centré qui semble le plus naturel. À partir de celui ci nous introduirons les no-
tions fondamentales pour l’étude des schémas aux différences finies (schéma expli-
cite/implicite, consistance et erreur de troncature et stabilité). Nous verrons alors
que ce schéma est inconditionnellement instable et donc qu’a priori la solution dis-
crète ne peut converger vers la solution du problème continu. Cet exemple aura pour
but de sensibiliser le lecteur sur le fait qu’un schéma consistant n’est pas nécessaire-
ment convergent. Nous étudierons ensuite le schéma de Lax-Friedrichs pour lequel
74 C HAPITRE 3 : Approximation par différences finies

nous introduirons la notion de convergence. Le schéma saute-mouton sera enfin étu-

dié comme approximation de l’équation des ondes.

3.2 Le schéma explicite centré pour l’équation d’advec-

tion
Nous allons écrire pour tout instant tn , n > 0 et tout point xj , une première équation
discrète qui va remplacer (et approcher) l’équation (3.3) en (xj , tn ). On utilise

∂u
• une approximation décentrée à droite (vers le futur), d’ordre 1, pour
∂t
∂u un+1
j − unj
(xj , tn ) ≃ ,
∂t ∆t

∂u
• une approximation centrée, d’ordre 2, pour , écrite à l’instant tn (voir Section
∂x
3.1) :
∂u unj+1 − unj−1
(xj , tn ) ≃ Dh0 u(xj , tn ) = .
∂x 2h
Le schéma numérique au point (xj , tn ) s’écrit donc

un+1
j − unj unj+1 − unj−1
+c = 0, j ∈ Z, n ≥ 0. (3.4)
∆t 2h
Il s’agit d’un schéma à un pas (ou deux points) en temps et à trois points en espace. Ce
schéma est bien explicite puisqu’on calcule explicitement la solution à l’instant n + 1
à partir de la solution à l’instant n sans recourir à l’inversion d’un système linéaire.
On a en effet
α
un+1
j = unj + (unj+1 − unj−1), j ∈ Z,
2
où
c∆t
α= .
h
Evidemment, pour que (3.4) définisse entièrement la solution approchée, il suffit de
connaître les u0j , j ∈ Z, ce pour quoi on va utiliser la condition initiale et on choisira
typiquement (noter que le choix (i) nécessite que u0 soit une fonction continue)
Z xj − h2
1
(i) u0j 0
= u (xj ) ou (ii) u0j = u0 (x) dx. (3.5)
h xj − h2

Pour l’analyse, il est commode d’associer à une suite (uj )j∈Z à valeurs réelles, une
fonction de la variable x. Nous choisissons dans la suite l’interpolation linéaire par
3.2 Le schéma explicite centré pour l’équation d’advection 75

morceaux (P1 désigne ci-après l’espace des polynômes d’une variable de degré infé-
rieur ou égal à 1). Ainsi pour une suite (uj )j∈Z on associe

uh ∈ Vh = vh ∈ C 0 (R) / ∀ j ∈ Z, vh |[xj ,xj +1] ∈ P1 ,

telle que uh est définie par

∀ j ∈ Z, uh (xj ) = uj .
On peut démontrer que, pour tout réel p > 1,
 X
(uj )j∈Z ∈ ℓp = (uj )j∈Z / |uj |p < +∞ ⇐⇒ uh ∈ Lph = Vh ∩ Lp (R).
j∈Z

et qu’il existe des constantes C∗ (p) et C ∗ (p), indépendantes de h, telles que

X  p1 X  p1
∀ (uj ) ∈ ℓp , C∗ (p) |uj |p h ≤ kuh kLp ≤ C ∗ (p) |uj |p h .
j∈Z j∈Z

De même

(uj )j∈Z ∈ ℓ∞ = { (uj ) / sup |uj |p < +∞} ⇐⇒ uh ∈ L∞ ∞

h = Vh ∩ L (R).
j∈Z

et on a l’identité
kuh kL∞ = sup |uj |.
j∈Z

Nous allons pouvoir réécrire (3.4) de façon abstraite en introduisant l’opérateur li-
néaire Ah tel que

v(x + h) − v(x − h)
v(x) : R → R, 7→ Ah v(x) = c .
2h
Puisque Ah commute avec les opérateurs de translation par un multiple de h, il envoie
linéairement Vh dans lui-même donc

Ah ∈ L(Vh ).

Mieux, on démontre facilement que

kA1 kL(Lp )
∀ p ∈ [1, +∞], Ah ∈ L(Lp ) et kAh kL(Lp ) = .
h2
On voit alors que, avec des notations évidentes, le schéma (3.4, 3.5-(i)) équivaut à

Trouver unh ∈ Vh , n = 0, 1, 2, · · · tels que

(3.6)
un+1 − unh
h
+ Ah unh = 0, n = 0, 1, 2, · · · ,
∆t
76 C HAPITRE 3 : Approximation par différences finies

et satisfaisant la condition initiale

u0h = Πh u0 , (3.7)
où Πh désigne l’opérateur d’interpolation linéaire de C 0 (R) dans Vh .

Il est bien évident que, compte tenu des propriétés de Ah ,

u0h ∈ Lp =⇒ unh ∈ Lp , n = 1, 2, · · ·
La solution de (3.6) s’ écrit donc

n
unh = Sh (∆t) u0h avec Sh (∆t) = I − ∆t Ah . (3.8)

3.2.1 Consistance et ordre du schéma

Rappelons que l’étude de la consistance d’un schéma vise à justifier (et quantifier)
la façon dont les équations discrètes approchent les équations continues. Elle ne dit
rien sur la façon dont la solution discrète approche la solution continue. Pour cela
il faut aussi avoir la stabilité du schéma, concept qui sera abordé dans la prochaine
section. Il est néanmoins bien évidemment nécessaire de vérifier la consistance d’un
schéma pour espérer sa convergence.

Pour ce faire, il faut définir l’erreur de troncature du schéma numérique. Le principe

consiste à partir de la solution exacte u(x, t) de l’équation d’advection (3.3), supposée
régulière, d’introduire la suite des valeurs discrètes
Ujn = u(xj , tn ). (3.9)
Si on remplace unj par Ujn dans le membre de gauche de (3.4), on n’a a priori aucune
raison de trouver 0 comme résultat (sinon le schéma serait « exact »), mais on est en
droit d’exiger que cela soit « petit ». Ce raisonnement justifie la
D ÉFINITION 3.2.1
On appelle erreur de troncature du schéma (3.4) la suite

Ujn+1 − Ujn n
Uj+1 n
− Uj−1
εnj = +c ,
∆t 2h
où Ujn est définie par (3.9), u étant une solution régulière de (3.3).

On dira que le schéma est consistant si, pour tout segment borné I et tout temps
T >0
lim εnj = 0.
(h,∆t)→0

Étant donnés deux entiers strictement positifs k et p, on dit que le schéma est
d’ordre p en temps et k en espace (ce qui entraine a fortiori la consistance) si


εnj = O ∆tp + hk .
3.2 Le schéma explicite centré pour l’équation d’advection 77

T HÉORÈME 3.2.1
Le schéma (3.4) est d’ordre 1 en temps et d’ordre 2 en espace.

D ÉMONSTRATION : C’est une application des développements de Taylor. Il existe θt ∈]0, 1[ tel
que :
 ∂u n ∆t2 ∂ 2 u
Ujn+1 = Ujn + ∆t + (xj , t + θt ∆t),
∂t j 2 ∂t2
où nous avons adopté la notation
 ∂u n ∂u  ∂u n ∂u
= (xj , tn ), = (xj , tn ), ··· (3.10)
∂t j ∂t ∂x j ∂x

Par conséquent
Ujn+1 − Ujn  ∂u n ∆t ∂ 2 u
= + (xj , t + θt ∆t). (3.11)
∆t ∂t j 2 ∂t2

De même, en développant par rapport à la variable x et avec la notation (3.10), il existe θ + ∈

]0, 1[ et θ − ∈]0, 1[ tels que
  ∂u n h2  ∂ 2 u n h3 ∂ 3 u

 n n
 U
 j+1 = U j + h + + (xj + θ + h, tn ),
∂x j 2 ∂x2 j 6 ∂x3

  ∂u n h2  ∂ 2 u n h3 ∂ 3 u

 Uj−1n
= Ujn − h + − (xj − θ − h, tn ).
∂x j 2 ∂x2 j 6 ∂x3

Après soustraction nous obtenons

n − Un
Uj+1  ∂u n h2 h ∂ 3 u ∂3u i
j−1 + n − n
= + (x j + θ h, t ) + (x j − θ h, t ) .
2h ∂x j 12 ∂x3 ∂x3

D’après le théorème des valeurs intermédiaires, il existe donc θ ∈] − 1, 1[ tel que

n − Un
Uj+1  ∂u n h2 ∂ 3 u
j−1
= + (xj + θh, tn ) (3.12)
2h ∂x j 6 ∂x3

De (3.11) et (3.12), nous déduisons, par définition de la solution u, que

 ∂u n  ∂u n ∆t ∂ 2 u h2 ∂ 3 u
εnj = +c + 2
(xj , t + θt ∆t) + c (xj + θh, tn )
∂t j ∂x j 2 ∂t 6 ∂x3
(3.13)
∆t ∂ 2 u h2 ∂ 3 u
= (x j , t + θ t ∆t) + c (xj + θh, tn ).
2 ∂t2 6 ∂x3
Ainsi nous obtenons

εnj = O ∆t + h2 .

78 C HAPITRE 3 : Approximation par différences finies

Pour l’analyse de convergence, il sera utile d’introduire une version « continue en

espace » de l’erreur de troncature, à savoir, en reprenant les notations de la section
précédente,
n Uhn+1 − Uhn
εh = + Ah Uhn .
∆t
où Uh = Πh u(·, t ). Reprenant l’expression (3.13), on déduit une estimation L∞ de εnh
n n

(avec une constante numérique C convenable)

 ∂2u ∂3u 
kεnh kL∞ ≤ C ∆t sup (·, t) + h2 (·, tn ) . (3.14)
t∈[tn ,tn+1 ] ∂t2 L∞ ∂x3 L∞

Nous pouvons également obtenir une estimation équivalente de εnh pour la norme L2 ,
c’est à dire qu’il existe une constante C telle que
 ∂2u ∂3u 
kεnh kL2 ≤ C ∆t sup (·, t) + h2 (·, tn ) . (3.15)
t∈[tn ,tn+1 ] ∂t2 L2 ∂x3 L2

Le lecteur se convaincra aisément que les estimations (3.14) et (3.15) se généralisent

à toute norme Lp (la constante C ci-dessous dépendant de p)
 ∂2u ∂3u 
kεnh kLp ≤ C ∆t sup (·, t) + h2 (·, tn ) .
t∈[tn ,tn+1 ] ∂t2 Lp ∂x3 Lp

3.2.2 Propagation numérique et condition nécessaire de convergence

La formule (3.4) ( ou (3.8) ) démontre que le support de la solution numérique se
propage à vitesse finie au sens où :

supp u0h ⊂ [a, b] =⇒ supp unh ⊂ [a − nh, b + nh].

On peut même dire que la propagation du support de la solution numérique est donc
donnée par
h
Vnum (∆, h) = .
∆t
Si le support de la solution numérique se propage moins vite que que le support de
la solution exacte, la solution numérique n’aura aucune chance de s’approcher de la
solution exacte puisqu’il va exister une partie du support de la solution exacte sur la-
quelle elle sera toujours nulle (voir figure 3.3-(a)).

Une condition nécessaire de convergence du schéma est donc que la vitesse numé-
rique doit être supérieure à la vitesse de propagation de la solution exacte (voir figure
3.3-(b)) :
h c∆t
≥c ⇔ α= ≤ 1. (3.16)
∆t h
3.2 Le schéma explicite centré pour l’équation d’advection 79

Supp u(·, tn ) Supp u(·, tn )

n n
t t
Supp unh Supp unh

a b a b
Supp u0 Supp u0
c∆t c∆t
(a) Cas où Vnum < c (⇔ α = > 1) (b) Cas où Vnum ≥ c (⇔ α = ≤ 1)
h h

F IGURE 3.3 : Condition nécessaire de convergence : comparaison du support de u(·, t)

(zone grise) et du support de unh (zone hachurée).

R EMARQUE 3.2.1 (R EMARQUE SUR LA CONDITION NÉCESSAIRE DE CONVERGENCE )

La condition (3.16) a une interprétation géométrique simple. En effet, pour l’équa-
tion de transport, on sait que la connaissance de u à l’instant tn dans l’intervalle
[xj−1 , xj+1 ] permet de calculer la solution dans l’intervalle

xj−1 − ctn ≤ x − ct ≤ xj+1 − ctn .

La condition (3.16) signifie donc que le point (xj , tn+1 ) doit se trouver dans la bande
grisée de la Figure 3.4.

tn+1
h/c
∆t
tn
h
xj−1 xj xj+1
h
F IGURE 3.4 : Condition (3.16) : ∆t ≤
c

3.2.3 Analyse de la stabilité d’un schéma : principes généraux

Avant d’aborder la notion de stabilité du schéma numérique, une notion intuitive-
ment simple mais délicate à définir mathématiquement, il est utile de préciser ce
que nous entendrons par convergence du schéma. Intuitivement, la convergence si-
gnifie que, dans une certaine norme, l’erreur entre la solution exacte et la solution
approchée tend vers 0 quand ∆t et h tendent vers 0. C’est le sens plus précis du « tend
vers 0 » qu’il nous faut définir.
80 C HAPITRE 3 : Approximation par différences finies

On supposera dans toute la suite que nous nous intéressons à l’approximation de

la solution dans un intervalle de temps fixe [0, T ], T > 0. Nous nous limiterons à la
convergence L∞ en temps à valeur Lp en espace.

D ÉFINITION 3.2.2
On dit qu’un schéma d’approximation de (3.3) produisant une suite unh ∈ Lp
converge dans L∞ (0, T ; Lp) si et seulement si, pour toute donnée initiale u0 dans
Lp ,
lim sup ku(·, tn ) − unh kLp = 0.
(h,∆t)→0 tn ≤T

Selon le principe qu’une condition nécessaire pour qu’une suite converge est qu’elle
soit bornée, nous introduisons le principe de stabilité Lp (qui devrait s’appeler plutôt
stabilité L∞ (0, T ; Lp )).

D ÉFINITION 3.2.3
On dit qu’un schéma numérique est Lp -stable si et seulement si, pour tout temps
T > 0, il existe une constante C(T ) indépendante de h et ∆t, telle que pour toute
donnée initiale u0 dans Lp ,

∀ tn ≤ T, kunh kLp ≤ C(T ) ku0 kLp . (3.17)

R EMARQUE 3.2.2
Comme on le verra, le résultat de stabilité (3.17) (ou encore le résultat de conver-
gence) peut n’être satisfait que sous une condition de stabilité de la forme

(h, ∆t) ∈ S ⊂ R∗+ × R∗+ ,

où S définit le domaine de stabilité de la méthode. Évidemment tout ceci n’a de sens

que si l’ensemble S admet (0, 0) comme point d’accumulation.

Dans ce qui suit, nous allons nous intéresser aux schémas numériques à un pas pro-
duisant une suite unh satisfaisant une relation de récurrence du type

un+1
h = Sh (∆t) unh , Sh (∆t) ∈ L(Lp ). (3.18)

Le schéma explicite (3.6) entre donc dans cette catégorie.

R EMARQUE 3.2.3
Un schéma consistant de la forme (3.18) vérifie, au moins formellement

lim Sh (∆t) = I .
∆t→0
3.2 Le schéma explicite centré pour l’équation d’advection 81

R EMARQUE 3.2.4
Démontrer un résultat de stabilité Lp pour un schéma à un pas du type (3.18) corres-
pond à l’obtention d’une estimation uniforme en (h, ∆t) de
kSh (∆t)n vkLp
kSh (∆t)n kL(Lp ) = sup .
v∈Lp kvkLp

En effet, la définition équivaut à

∀ tn ≤ T, kSh (∆t)n kL(Lp ) ≤ C(T ). (3.19)

Une condition suffisante de stabilité, très utile dans la pratique, est donnée par le

T HÉORÈME 3.2.2
Une condition suffisante de stabilité Lp du schéma (3.18) est qu’il existe une
constante ν ≥ 0 telle que

∀ (∆t, h) ∈ S, kSh (∆t)kL(Lp ) ≤ 1 + ν ∆t. (3.20)

Auquel cas on a
n
kSh (∆t)n kL(Lp ) ≤ eν t ,

ce qui entraine en particulier (3.17) avec C(T ) = eν T .

D ÉMONSTRATION : Par définition de kSh (∆t)kL(Lp ) , il est clair que

kSh (∆t)n kL(Lp ) ≤ (1 + ν∆t)n .

Pour conclure, il suffit d’utiliser le résultat suivant : pour tout x > 0 et tout entier n
x
n
1+ ≤ ex .
n
En effet, pour n = 1, l’inégalité 1 + x ≤ ex traduit simplement le fait que le graphe de la
fonction x 7→ ex est au dessus de sa tangente au point d’abscisse x = 0. Plus généralement, la
fonction x 7→ xn étant croissante, il vient
x x x
n
1 + ≤ e n =⇒ 1+ ≤ ex .
n n


Notons que si la condition (3.20) est réalisée avec ν = 0, c’est à dire si

∀ (∆t, h) ∈ S, kSh (∆t)kL(Lp ) ≤ 1,

82 C HAPITRE 3 : Approximation par différences finies

alors on a
kunh kLp ≤ ku0h kLp ,
qui est l’équivalent discret de

ku(·, t)kLp ≤ ku0 kLp .

On parlera alors de stabilité uniforme.

Signalons ainsi que le théorème 3.2.2 admet une réciproque, dont la démonstration,
qui nécessite un peu de connaissances sur la théorie spectrale, dépasse le cadre de ce
cours.

T HÉORÈME 3.2.3
Pour qu’un schéma numérique du type (3.18) soit Lp stable, il faut qu’il existe une
constante ν ≥ 0 telle que

∀ (∆t, h) ∈ S, kSh (∆t)kL(Lp ) ≤ 1 + ν ∆t. (3.21)

3.2.4 Analyse de stabilité L2 du schéma explicite centré : la méthode

de Fourier-Von Neumann
Le schéma explicite (3.6) entre dans la classe (3.18) avec

Sh (∆t) = I − ∆t Ah . (3.22)

Pour étudier kSh (∆t)kL(L2 ) , nous utilisons la transformation de Fourier et introdui-

sons l’opérateur
Sbh (∆t) = F Sh (∆t) F −1 ∈ L(L2 ). (3.23)
On montre alors qu’il existe Sbh (k, ∆t) ∈ L∞ (R, C) tel que


∀k ∈ R, ûn+1 b n b n
(3.24)
h (k) = Sh (∆t)ûh (k) = Sh (k, ∆t) ûh (k),

où on a posé ûnh = F unh .

Grâce au théorème de Plancherel, nous avons

kSbh (∆t)kL(L2 ) = kSh (∆t)kL(L2 ) .

Le théorème 3.2.4 va nous donner une autre expression pour kSh (∆t)kL(L2 ) qui sera
beaucoup utilisée en pratique. Pour cela, il est utile d’introduire quelques définitions
et énoncer (sans preuve) quelques résultats de base (que nous invitons le lecteur à
démontrer à titre d’exercice) :
3.2 Le schéma explicite centré pour l’équation d’advection 83

• Étant donné bb(k) ∈ L∞ (R, C), on appelle opérateur de multiplication par bb,
b linéaire et continu de Lp (R; C) dans lui-même défini par :
l’opérateur B,


∀ v ∈ Lp (R, C), b (k) = bb(k) v(k),
Bv p. p. k ∈ R.

C’est un exercice très simple de montrer que :

∀ p ∈ [1, +∞], b L(Lp ) = kbbkL∞

kBk (3.25)

• Dans le cas p = 2, on appelle opérateur de symbole bb, l’opérateur

b F −1
B=F B ∈ L(L2 ).

Du théorème de Plancherel et de (3.25), nous déduisons aisément que

kBkL(L2 ) = kbbkL∞ (3.26)

D’après les relations (3.23-3.24), Sh (∆t) est un opérateur à symbole, on en déduit

donc d’après (3.26)

T HÉORÈME 3.2.4

∀ ∆t, kSh (∆t)kL(L2 ) = sup |Sbh (k, ∆t)|

k∈R

Il nous reste donc à étudier la norme L∞ de k 7→ Sbh (k, ∆t). Pour ce faire, il nous faut
caractériser la transformée de Fourier de Ah . Celui ci est également un cas particulier
d’opérateur à symbole :

L EMME 3.2.5
On a l’identité
∀ u ∈ L2 , d
A b b(k),
h u (k) = Ah (k) u ∀ k ∈ R, (3.27)
bh (k) donné par
et l’opérateur Ah a pour symbole A

bh (k) = ı c sin kh .
A
h

Ah
u Ah u

F F
bh
A
b
u bh u
A b
84 C HAPITRE 3 : Approximation par différences finies

D ÉMONSTRATION : La propriété (1.10) signifie que l’opérateur τa (défini par τa u(x) = u(x + a))
n’est autre que l’opérateur de symbole eiak . Comme
τh − τ−h
Ah = c ,
2h
on en déduit que Ah est l’opérateur de symbole
ikh − e−ikh
bh (k) = c e
A = ıc
sin kh
2h h
ce qui achève la démonstration. 

R EMARQUE 3.2.5
On peut remarquer que


bh (k) = ı c|k| 1 + O |k|2 h2 .
A

On retrouve alors sous une autre forme, le fait que l’opérateur Ah dont le symbole est
bh (k) approche l’opérateur c ∂ dont le symbole est ı ck , et ce au second ordre en h.
A ∂x

De (3.22), nous déduisons facilement le

C OROLLAIRE 3.2.6
L’opérateur Sh (∆t) est l’opérateur de symbole

c∆t
Sbh (k, ∆t) = 1 − ı α sin kh, avec α = .
h
Par conséquent, d’après le théorème 3.2.4

kSh (∆t)kL(L2 ) = sup |Sbh (k, ∆t)| = (1 + α2 )1/2 .

k∈R

Finalement, nous avons tous les éléments pour établir le théorème suivant

T HÉORÈME 3.2.7
Le schéma explicite centré (3.4) est inconditionnellement instable dans L2 .

Au final, le schéma explicite centré qui semblait à première vue naturel est donc, mal-
gré sa consistance, instable. La solution numérique calculée grâce à ce schéma ne
peut a priori pas converger vers la solution du problème continu (3.3).

La figure 3.5 illustre l’instabilité du schéma explicite centré. Pour plusieurs pas de
3.2 Le schéma explicite centré pour l’équation d’advection 85

temps, nous représentons en rouge la solution exacte (de l’équation d’advection (3.3)
pour une vitesse c particulière) et en bleu la solution numérique calculée en utilisant
le schéma explicite centré (3.4).
R EMARQUE 3.2.6 (R EMARQUE SUR LA MÉTHODE DE F OURIER -V ON N EUMANN )
Les calculs précèdents apparaissent tout naturellement quand on cherche la solution
explicite de (3.6) en appliquant la transformation de Fourier. Nous procédons donc
ici de la même façon que pour l’équation continue (3.3) à la section 1.2.2.

Appliquons la transformation de Fourier à (3.6). Si u

bnh désigne la transformée de Fou-
rier de unh , en utilisant (3.27), nous obtenons :
bn+1
u bnh (k)
h (k) − u sin kh n
+ ıc bh (k) = 0,
u
∆t h
c’est à dire
bn+1
u b bn (k) (3.28)
h (k) = Sh (k, ∆t) u
Cette formule est à comparer à la formule (1.11) pour la solution exacte : le symbole
Sbh (k, ∆t) apparaît donc comme le coefficient d’amplification du schéma.

Par suite
bnh (k) = Sbh (k, ∆t)n u
u b0(k).
Par le théorème de Plancherel
Z +∞ n
kunh k2L2 = 2π Sbh (k, ∆t) u0 (k)|2 dk
|b
−∞

et la recherche de l’estimation (3.17) (pour p = 2) mène ainsi naturellement au résul-

tat du théorème 3.2.7.

Retenons ici une recette de calcul pour la mise en oeuvre de la méthode de Fourier-
Von Neumann lors de l’analyse de stabilité L2 d’un schéma numérique approchant
une équation aux dérivées partielles d’évolution en une dimension d’espace :
(i) On cherche une solution de la forme
unj = unh (k) eikxj , j ∈ Z, n ≥ 0, (unh (k) ∈ C).

(ii) On détermine une relation de récurrence satisfaite par la suite unh (k),
(iii) On étudie le comportement de la famille de suites (unh (k), k ∈ R).
Remarquons alors que :
• cette méthode n’est applicable que si l’on s’intéresse à l’approximation d’une
équation linéaire à coefficients constants posée sur R, sur un maillage régulier,
• dans le cas des schémas numériques à un pas de temps, la relation de récur-
rence est nécessairement de la forme (3.28) et l’analyse de stabilité se réduit à
l’étude du coefficient d’amplification numérique.
86 C HAPITRE 3 : Approximation par différences finies

F IGURE 3.5 : Schéma explicite centré (inconditionnellement instable).

3.3 Le schéma de Lax-Friedrichs pour l’équation d’advection 87

3.3 Le schéma de Lax-Friedrichs pour l’équation d’ad-

vection
À la différence du schéma explicite centré, le schéma de Lax-Friedrichs cherche à
symétriser le traitement sur la variable temporelle tout en conservant le caractère
explicite du schéma, le caractère centré. On utilise donc :
∂u
• l’approximation de suivante
∂t
n+1
unj+1 + unj−1
∂u u j −
(xj , tn ) ≃ 2 ,
∂t ∆t
∂u
• une approximation centrée, d’ordre 2, pour , écrite à l’instant tn (voir Section
∂x
3.1)
∂u unj+1 − unj−1
(xj , tn ) ≃ Dh0 u (xj , tn ) = .
∂x 2h
L’équation du schéma de Lax-Friedrichs au point (xj , tn ) s’écrit donc
2un+1
j − unj+1 − unj−1 unj+1 − unj−1
+c = 0. (3.29)
2∆t 2h
Il s’agit d’un schéma à un pas (ou deux points) en temps et à deux points en espace.
Ce schéma est bien explicite. Il permet de calculer explicitement la solution à l’instant
n + 1 à partir de la solution à l’instant n
1−α n 1+α n
un+1
j = uj+1 − uj−1, j ∈ Z,
2 2
où nous rappelons que
c∆t
α= .
h
Evidemment, pour que (3.29) définisse entièrement la solution approchée, il suffit de
connaître les u0j , j ∈ Z, ce pour quoi on va, comme indiqué dans (3.5), utiliser la
condition initiale.
R EMARQUE 3.3.1 (I NTERPRÉTATION GÉOMÉTRIQUE DU SCHÉMA DE L AX F RIEDRICHS )
On rappelle que la solution exacte est constante sur les droites caractéristiques qui
sont ici plus précisément les droites de pente c. On a donc
u(xj , tn+1 ) = u(x∗ , tn ), avec x∗ = xj − c∆t.
On remarque ainsi aisément que (voir Figure 3.6) :
1+α 1−α
x∗ = xj−1 + xj+1 .
2 2
Le schéma de Lax Friedrichs consiste donc à considérer une interpolation linéaire sur
le maillage de la solution continue.
88 C HAPITRE 3 : Approximation par différences finies

tn+1
dx
dt
=c

tn
x∗
xj−1 xj xj+1

h − c∆t h + c∆t
= 1−α
2
· 2h = 1+α
2
· 2h

F IGURE 3.6 : Interprétation géométrique du schéma de Lax Friedrichs

3.3.1 Consistance et ordre du schéma

Nous commençons par l’étude de l’erreur de troncature qui nous conduit à l’analyse
de la consistance. Contrairement au schéma explicite centré, l’analyse est ici un peu
plus complexe.

Respectant la définition 3.2.1, l’erreur de troncature du schéma (3.29) est donnée par
U n +U n
n
Ujn+1 − j+1 2 j−1 n
Uj+1 n
− Uj−1
εj = +c , (3.30)
∆t 2h
où Ujn est définie par Ujn = u(xj , tn ), u étant une solution régulière de (3.3).

On peut alors démontrer

T HÉORÈME 3.3.1
Le schéma (3.29) est d’ordre 1 en temps et en espace.

D ÉMONSTRATION : On commence par remarquer que

n + Un
Uj+1 n − 2U n + U n
Uj+1
j−1 j j−1
= Ujn + ,
2 2
en utilisant la définition 3.2.1 de l’erreur de troncature pour le schéma explicite centré que
nous notons (εnj )EC , il vient facilement que l’erreur de troncature pour le schéma de Lax-
Friedrichs εnj est donnée par :
n n n
h2 Uj+1 − 2Uj + Uj−1
εnj = (εnj )EC + .
2∆t h2
Si on fixe α = c∆t
h et en utilisant l’opérateur d’approximation de la dérivée seconde donnée
dans (3.2), on obtient
ch h ∂ 2 u n
εnj = (εnj )EC + + O(h2 )]. (3.31)
2α ∂x2 j
3.3 Le schéma de Lax-Friedrichs pour l’équation d’advection 89

En utilisant le calcul de l’erreur de consistance pour le schéma explicite centré (voir la preuve
du Théorème 3.2.1)

(εnj )EC = O ∆t + h2 ,
on obtient finalement

εnj = O ∆t + h .


Pour l’analyse de convergence, nous utilisons la version « continue en espace » de

l’erreur de troncature, à savoir, en reprenant les notations de la section 3.2,
u(x+h,tn )+u(x−h,tn )
u(x, tn+1 ) − u(x + h, tn ) − u(x − h, tn )
εnh (x) = 2
+c .
∆t 2h
À partir des expressions (3.13) et (3.31), on déduit une estimation L∞ de εnh (avec une
constante numérique C dépendant de α)
 ∂2u ∂2u 
n
kεh kL∞ ≤ C(α) ∆t sup (·, t) + h sup (·, t) .
t∈[tn ,tn+1 ] ∂t2 L∞ t∈[tn ,tn+1 ] ∂x2 L∞

Nous pouvons également obtenir une estimation équivalente de εnh pour la norme Lp ,
c’est à dire qu’il existe une constante C(α) telle que
 ∂2u ∂2u 
n
kεh kLp ≤ C(α) ∆t sup (·, t) + h sup (·, t) ,
t∈[tn ,tn+1 ] ∂t2 Lp t∈[tn ,tn+1 ] ∂x2 Lp

ce qui pour le problème de Cauchy se ramène au

T HÉORÈME 3.3.2
Pour le schéma de Lax-Friedrichs (3.29), l’erreur de troncature est telle que pour
tout p, il existe une constante C(α) dépendant de α telle que

∂2u
kεnh kLp ≤ C(α) h (·, t) . (3.32)
∂x2 Lp

3.3.2 Vitesse de propagation numérique

La formule (3.29) démontre que le support de la solution numérique se propage à
vitesse finie au sens où
supp u0h ⊂ [a, b] =⇒ supp unh ⊂ [a − nh, b + nh]
La vitesse de propagation de la solution numérique est donc donnée par
h
Vnum (∆, h) = .
∆t
90 C HAPITRE 3 : Approximation par différences finies

Comme pour le schéma explicite centré, une condition nécessaire de convergence

du schéma est donc que la vitesse numérique doit être supérieure à la vitesse de pro-
pagation de la solution exacte (voir Figure 3.3) :
h c∆t
≥c ⇔ α= ≤ 1. (3.33)
∆t h

3.3.3 Analyse de stabilité

Le schéma de Lax Friedrichs (3.29) produit une suite unh satisfaisant la relation de
récurrence
un+1
h = Sh (∆t) unh (3.34)
avec
  1−α 1+α
Sh (∆t) u (x) = u(x + h) + u(x − h).
2 2
En appliquant la méthode de Fourier-Von Neumann développée à la section 3.2.3 ou
l’algorithme proposé en Remarque 3.2.6 sur la stabilité L2 , on montre que la relation
précédente mène au coefficient d’amplification
1 − α ıkh 1 + α −ıkh
Sbh (k, ∆t) = e + e ,
2 2
soit (voir Figure 3.7 pour la représentation dans le plan complexe de k 7→ Sbh (k, ∆t))

Sbh (k, ∆t) = cos(kh) − ıα sin(kh).

On rappelle qu’un schéma donné par une relation de récurrence d’ordre 1 est L2 -
stable si et seulement si il existe une constante ν ≥ 0 telle que

∀ (∆t, h) ∈ S, kSh (∆t)kL(L2 ) = sup |Sbh (k, ∆t)| ≤ 1 + ν ∆t.

k∈R

On en déduit alors le résultat de stabilité L2 pour le schéma de Lax-Friedrichs.

T HÉORÈME 3.3.3
Le schéma de Lax-Friedrichs (3.29) est stable (et même uniformément stable) pour
la norme L2 si et seulement si
c∆t
α= ≤ 1. (3.35)
h
Cette condition est appelée condition CFL (Courant-Friedrichs-Lewy).

D ÉMONSTRATION : Il suffit de remarquer que

sup |Sbh (k, ∆t)| = max(1, α).

k∈R


3.3 Le schéma de Lax-Friedrichs pour l’équation d’advection 91

Im z Im z
k 7→ Sbh (k, ∆t) k 7→ Sbh (k, ∆t)

1 1
Re z Re z

(a) Cas où α > 1 : le schéma est instable (b) Cas où α ≤ 1 : le schéma est stable

F IGURE 3.7 : Coefficient d’amplification dans le plan complexe en fonction de α =

c∆t
h

La stabilité L2 est donc obtenue par le calcul du coefficient d’amplification en utili-

sant la méthode de Fourier-Von Neumann. Plus généralement, pour la stabilité Lp , il
faut utiliser des arguments différents. Pour le schéma de Lax-Friedrichs,on utilise un
argument de convexité.

T HÉORÈME 3.3.4
Le schéma de Lax-Friedrichs (3.29) est stable (et même uniformément stable) pour
la norme Lp si et seulement si
c∆t
α= ≤1
h

D ÉMONSTRATION : Supposons que α ≤ 1. On rappelle que pour tout u

  1−α 1+α
Sh (∆t) u (x) = u(x + h) + u(x − h),
2 2
avec
1−α 1+α 1−α 1+α
≥ 0, ≥ 0, + = 1.
2 2 2 2
Par convexité de x 7→ |x|p pour p ≥ 1, on obtient pour toute fonction u ∈ Lp , d’une part que
Sh (∆t) u ∈ Lp et d’autre part que
Z Z Z
  1−α 1+α
| Sh (∆t) u (x)|p dx ≤ |u(x + h)|p dx + |u(x − h)|p dx
R 2 R 2 R

ce qui entraîne
kSh (∆t) ukLp ≤ kukLp .
En utilisant, le théorème 3.2.2, on en déduit que le schéma est stable pour la norme Lp . 
92 C HAPITRE 3 : Approximation par différences finies

La condition CFL (3.35) correspond exactement à la condition nécessaire de conver-

gence (3.33) obtenue en comparant la vitesse numérique et la vitesse de propagation
de la solution continue. Elle dit que le pas de temps doit être inférieur à h/c. Cette
condition est quelque peu contraignante : même si on se place à la limite de la condi-
tion de stabilité, lorsque qu’on divise le pas d’espace par 2, le pas de temps doit être
divisé par 2. Si on raffine le maillage en espace par 10, il faut raffiner le maillage en
temps par 10 !

L’intérêt du schéma explicite est que les calculs à effectuer à chaque pas de temps
sont très simples et peu coûteux. La contre-partie est que la condition de stabilité
contraint à choisir un pas de temps aussi petit que le pas en espace ce qui peut en-
gendrer un gros surcoût de calcul.

3.3.4 Convergence du schéma

Passons maintenant à la convergence du schéma (3.29) dans L∞ (0, T ; Lp ) en mon-
trant d’après la définition 3.2.3 que pour toute donnée initiale dans Lp

lim sup ku(·, tn ) − unh kLp = 0,

(h,δt)→0 tn ≤T

où u est la solution continue de (3.3) et unh est la suite produite par le schéma (3.29).

Nous introduisons donc l’erreur commise sur la solution

enh = u(·, tn ) − unh ,

et obtenons la convergence du schéma comme une conséquence de sa consistance

et de sa stabilité.

La première étape consiste donc à relier l’erreur enh à l’erreur de troncature εnh . C’est
particulièrement simple ici grâce à la linéarité des équations. En effet, des égalités
(3.34) et (3.30), nous déduisons que

u(·, tn+1) = Sh (∆t) u(·, tn ) + ∆t εnh ,

un+1
h = Sh (∆t) unh ,

soit par différence

en+1
h = Sh (∆t) enh + ∆t εnh . (3.36)
Ainsi, l’erreur enh apparaît comme la solution d’un problème d’évolution discret en-
tretenu par un terme source qui n’est autre que l’erreur de troncature.

Plus généralement, nous allons en fait montrer que si une suite enh vérifie la relation
(3.36) où l’opérateur Sh (∆t) satisfait l’estimation de stabilité Lp (3.19), alors on peut
3.3 Le schéma de Lax-Friedrichs pour l’équation d’advection 93

estimer enh en fonction de e0h et des εkh uniformément en ∆t, tant que tn ≤ T .

En effet, nous observons que

e1h = Sh (∆t) e0h + ∆t ε0h ,

 
e2h = Sh (∆t) e1h + ∆t ε1h = Sh (∆t) 2
e0h 1 0
+ ∆t εh + Sh (∆t) eh ,

et, de proche en proche,

X
n−1 
enh = Sh (∆t) n
e0h + ∆t Sh (∆t)k εhn−1−k .
k=0

En passant aux normes, on obtient alors

X
n−1 
kenh kLp = kSh (∆t) n
kL(Lp ) ke0h kLp + ∆t kSh (∆t) k
kL(Lp ) kεhn−1−k kLp .
k=0

Si on suppose l’hypothèse de stabilité (3.21) satisfaite, on a

kSh (∆t)kL(Lp ) ≤ 1 + ν ∆t.
Or on sait que
q
∀tq ≤ T, kSh (∆t)q kL(Lp ) ≤ eνt ≤ eνT .
Nous avons donc,
 
n
∀ t ≤ T, kenh kLp ≤ e νT 0 n k
keh kLp + t sup kεh kLp ,
k≤n−1

ce qui mène à l’estimation d’erreur

 
sup kenh kLp ≤ e νT
ke0h kLp + T sup kεnh kLp .
tn ≤T tn ≤T

Bien entendu, dans le cas d’un schéma uniformément stable (ν = 0), l’estimation
devient  
sup kenh kLp ≤ ke0h kLp + T sup kεnh kLp . (3.37)
tn ≤T tn ≤T

Nous pouvons appliquer ce résultat général pour établir des estimations d’erreur Lp
pour le schéma de Lax Friedrichs (3.29).

T HÉORÈME 3.3.5
c∆t
Soient u la solution de (3.3) et unh la solution de (3.29, 3.7). Si ≤ 1 et si, pour
h
p ∈ [1, +∞], on a l’estimation d’erreur L (0, T ; L ) :
∞ p

d2 u 0
sup ku(·, tn ) − unh kLp ≤ h C(α) T .
tn ≤T dx2 Lp
94 C HAPITRE 3 : Approximation par différences finies

D ÉMONSTRATION : D’après le théorème 3.3.4, sous la condition CFL (3.35), nous savons que
le schéma est uniformément stable dans Lp . Nous pouvons donc appliquer la formule (3.37)
dans laquelle nous injectons l’estimation (3.32) de l’erreur de troncature.

Le terme de (3.37) concernant eh0 est nul si la condition initiale imposée pour la solution nu-
mérique est (choix (i) de (3.7))

∀j, u0j = u(xj , t = 0) = u0 (xj ).

Si la condition initiale imposée est autre ( par exemple le choix (ii) de (3.7)), pour estimer le
terme en eh0 nous utilisons le résultat classique suivant (sur l’interpolation linéaire - exercice)

d2 u0
ke0h kLp = kΠh u0 − u0 kLp ≤ Ch k kLp .
dx2


D’après la Section 3.3.2 et le théorème, la condition CFL (3.35) est donc une condi-
tion nécessaire et suffisante de convergence du schéma de Lax-Friedrichs. Les figures
3.8 et 3.9 illustrent l’importance de cette condition dans la stabilité et donc la conver-
gence du schéma : pour plusieurs pas de temps, nous représentons en rouge la solu-
tion exacte (de l’équation d’advection (3.3) pour une vitesse c particulière) et en bleu
la solution numérique calculée en utilisant le schéma de Lax Friedrichs (3.29).

De manière générale, nous venons d’établir dans cette section, le théorème de Lax.

T HÉORÈME 3.3.6 (LAX)

Soit u une solution suffisamment régulière d’un système différentiel et unj la solution
numérique discrète obtenue par un schéma aux différences finies avec la donnée
initiale u0j = u0 (xj ). On suppose le schéma linéaire, consistant et, pour un certain
domaine de stabilité S admettant (0, 0) comme point d’accumulation, stable pour
une norme k k, alors le schéma est convergent au sens où

∀(h, ∆t) ∈ S, ∀T > 0, lim sup kenh k = 0,

∆t,h→0 tn ≤T

avec enh = unh − u(tn , ·) l’erreur commise sur la solution.

De plus, l’ordre de convergence est donné par l’erreur de troncature : c’est à dire que
si le schéma est d’ordre p en temps et d’ordre q en espace alors

∀n, enh = O(∆tp + hk ).

3.3 Le schéma de Lax-Friedrichs pour l’équation d’advection 95

F IGURE 3.8 : Schéma de Lax-Friedrichs avec c∆t = 0.3h

96 C HAPITRE 3 : Approximation par différences finies

F IGURE 3.9 : Schéma de Lax-Friedrichs avec c∆t = 1.01h

3.4 Le schéma saute-mouton pour l’équation des ondes 97

3.4 Le schéma saute-mouton pour l’équation des ondes

Nous considérons ici l’équation des ondes sur R à vitesse constante c
 2 2
 ∂ u 2∂ u

 − c = 0, pour t ≥ 0, x ∈ R

 ∂t2 ∂x2


u(t = 0, x) = u0 (x), pour x ∈ R (3.38)





 ∂u (t = 0, x) = u1 (x),

pour x ∈ R
∂t
Nous discrétisons toujours l’espace avec un pas h et le temps avec ∆t > 0 et cher-
chons toujours à calculer l’inconnue discrète un = (unj )j∈Z à chaque pas de temps.

3.4.1 Présentation du schéma

Le schéma saute-mouton est un schéma explicite centré :
un+1
j − 2unj + ujn−1 un − 2unj + unj−1
2 j+1
∀n, ∀j, − c = 0. (3.39)
(∆t)2 h2
La première remarque fondamentale sur le plan pratique, est que, son caractère ex-
plicite permet une mise en oeuvre peut couteuse en résolvant simplement :

un+1
j = 2unj − ujn−1 + α2 (unj+1 − 2unj + unj−1),

où le nombre sans dimension :

c∆t
α= ,
h
mesure le rapport entre la distance parcourue par l’onde pendant le temps ∆t (soit
c∆t) et la taille de la maille en espace h. Nous verrons un peu plus loin que, comme
pour l’equation d’advection, cette quantité est liée à la condition CFL des schémas
introduits.

Bien entendu pour pouvoir déterminer complètement la solution discrète, il faut la

connaître aux deux premiers instants. C’est pour la détermination de la solution aux
deux premiers instants que vont intervenir les conditions initiales. Les difficultés qui
en découlent sont caractéristiques de l’équation des ondes qui est d’ordre 2 en espace
et en temps.

3.4.2 Le schéma de démarrage : choix des conditions initiales ap-

prochées
Pour le problème continu, les données initiales sont définies par les valeurs de la
solution à t = 0, ainsi que celles de sa dérivée en temps au même instant.
98 C HAPITRE 3 : Approximation par différences finies

Pour le problème discrétisé, pour déterminer u0h et u1h on doit approcher u(xj , 0) et
u(xj , ∆t).

Pour u0h , il est naturel de choisir

∀j, u0j = u0 (xj ). (3.40)

Pour approcher u(xj , ∆t), on utilise simplement un développement de Taylor au pre-

mier ordre
∂u
u(xj , ∆t) = u(xj , 0) + ∆t (xj , 0) + O(∆t2 ),
∂t
ce qui mène à l’approximation

∀j, u1j = u0 (xj ) + ∆t u1 (xj ).

Nous verrons que la précision d’un tel schéma de démarrage est insuffisante. Pour en
avoir l’intuition, il suffit de remarquer que ce choix correspond à une approximation
d’ordre 1 en temps de la condition initiale
∂u
(x, 0) = u1 (x).
∂t
Ainsi cette dernière équation se réécrit
u1j − u0j
∀j, = u1 (xj ),
∆t
où le terme de gauche peut être vue comme une approximation décentrée (donc
d’ordre 1) de la dérivée en temps de la solution à l’instant t = 0.

Pour avoir un schéma de démarrage d’ordre 2 en temps (plus cohérent avec l’ordre
du schéma (3.39) qui, comme on le verra plus loin, est égal à 2), il suffit de pousser le
développement de Taylor un cran plus loin
∆t2 ∂ 2 u
u(xj , ∆t) = u0 (xj ) + ∆t u1 (xj ) + (xj , 0) + O(∆t3 ),
2 ∂t2
et, comme on ne dispose pas de la dérivée seconde de la solution en tant que donnée
initiale, d’utiliser l’équation des ondes écrite à l’instant t = 0 pour en déduire que
∂2u 2 0
2∂ u
(x, 0) = c (x),
∂t2 ∂x2
ce qui mène au schéma de démarrage

c2 ∆t2 u0j+1 − 2u0j + u0j−1

∀j, u1j = u0j 1
+ ∆t u (xj ) + . (3.41)
2 h2
Le schéma numérique est donc complètement défini par (3.39), (3.40) et (3.41).
3.4 Le schéma saute-mouton pour l’équation des ondes 99

3.4.3 Ordre du schéma : consistance et erreur de troncature

De la même façon que dans la définition 3.2.1, l’erreur de troncature du schéma (3.39)
est donnée par
Ujn+1 − 2Ujn + Ujn−1 U n − 2 Ujn + Uj−1
2 j+1
n
εnj
= − c .
∆t2 h2
où Ujn est définie par Ujn = u(xj , tn ), u étant une solution régulière de (3.38).

On peut alors montrer que

T HÉORÈME 3.4.1
Le schéma (3.39) est consistant à l’ordre 2 en temps et en espace :

εnj = O(∆t2 + h2 ).

D ÉMONSTRATION : Il suffit d’utiliser l’opérateur d’approximation de la dérivée seconde don-

née dans (3.2). 

3.4.4 Analyse de stabilité L2 par techniques énergétiques

Comme nous l’avons vu pour d’autres types d’équation, la stabilité d’un schéma (dans
une certaine norme) est un préalable incontournable à sa convergence (dans cette
même norme). Nous allons mener cette analyse dans le cadre le plus adapté à l’équa-
tion des ondes qui est le cadre L2 du fait de la périodicité des solutions de l’équation
des ondes. Ayant déjà étudié la méthode de Fourier-Von Neumann pour l’équation
d’advection, nous allons ici utiliser une autre méthode fondée sur une propriété fon-
damentale de l’équation des ondes : la conservation de l’énergie.

Un résultat de conservation d’énergie discrète

Il est assez naturel de se demander si le schéma utilisé satisfait une relation de conser-
vation équivalente à la conservation d’énergie pour la solution du problème continu.
En effet, nous avions vu au théorème 1.3.3 que l’équation des ondes vérifie pour tout
t > 0,
Z Z
1 ∂u 2 c2 ∂u 2
E(t) = E(0) avec E(t) = (t, x) dx + (t, x) dx.
2 R ∂t 2 R ∂x
1
Introduisons alors l’énergie discrète (à l’instant tn+ 2 = n + 12 ∆t) :
1 X  uj − unj 2 1
n+1
1
E n+ 2 = + ah (un+1 n
h , uh ), n ≥ 1, (3.42)
2 j∈Z ∆t 2
100 C HAPITRE 3 : Approximation par différences finies

avec X  uj+1 − uj  vj+1 − vj 

ah (u, v) = c2 .
j∈Z
h h

et où unh ,
définie sur (0, 1) est la version « continue en espace » de la solution discrète,
en reprenant les notations de la section 3.2.

En remarquant que ah (u, v) est un équivalent discret de

Z
2 ∂u ∂v
c dx,
R ∂x ∂x

on voit que la quantité (3.42) peut être vu comme une approximation de l’énergie
1
E(tn+ 2 ) où l’on a remplacé le produit scalaire L2 (R) sur la variable d’espace par le
produit scalaire ℓ2 sur le vecteur (unj )j :
X
kunh k2ℓ2 = unh , unh ℓ2
= (unj )2 .
j∈Z

Autrement dit
n+ 21 1 un+1 − unh 1
Eh = k h kℓ2 + ah (un+1 n
h , uh ).
2 ∆t 2

T HÉORÈME 3.4.2
Si la suite (unj )j est solution du schéma saute-mouton (3.39), l’énergie discrète
n+ 21
Eh se conserve au cours du temps :
n+ 12 n− 21
Eh = Eh , ∀n ≥ 1.

P REUVE : Nous laissons les détails de la démonstration au lecteur pour nous concentrer sur
les principes de la démonstration.

Rappelons que la conservation de l’énergie continue E(t) s’obtient en multipliant l’équation

des ondes par ∂u
∂t et en intégrant le résultat en espace. Nous procédons de manière équiva-
lente au niveau discret en multipliant par

un+1
h − unh
.
2∆t
Or on observe l’apparition d’identités remarquables
* + !
un+1
h − 2u n + un−1 un+1 − un
h h h h 1 un+1
h − unh unh − uhn−1
, = − ,
∆t2 2∆t 2∆t ∆t ℓ2 ∆t ℓ2
ℓ2
* +
unh − 2unh + unh un+1 − u n
1
2
, h h
= (ah (un+1 n n n−1
h , uh ) − ah (uh , uh )).
h 2∆t 2
∆t
ℓ
3.4 Le schéma saute-mouton pour l’équation des ondes 101

Ceci conduit bien en additionnant ces deux relations à

1 n+1/2 n−1/2
(E − Eh ) = 0,
2∆t h
c’est à dire le résultat annoncé. 

Principe de la stabilité à partir de la conservation d’énergie

Une fois cette conservation établie il nous faut maintenant l’utiliser pour assurer un
résultat de stabilité de type
kunh kL2 (R) ≤ ku0h kL2 (R) .

Le problème majeur dans l’étude discrète est que l’énergie discrète n’est pas néces-
sairement une quantité positive dû au terme
ah (un+1 n
h , uh ).

Ainsi même si l’énergie se conserve on ne peut pas être sûr qu’on conserve une norme
sur unh . Toutefois si le pas de temps est suffisamment petit (nous voyons apparaître ici
une condition de type CFL), un+1 h sera proche de unh et ah (un+1
h , uh ) sera donc presque
n

un carré. Dans ce cas, on pourra obtenir une forme « d’équivalence de norme » qui
permettra de relier la stabilité de la norme L2 (R) à celle de l’énergie.

On obtient alors le théorème suivant que nous admettrons

T HÉORÈME 3.4.3 (Admis)
Soit unh la solution discrète. Si la solution du problème est assez régulière alors
c∆t
α= <1 =⇒ ∀T > 0 lim sup kenh kL2 (R) = 0
h ∆t,h→0 tn ≤T

avec enh = unh − u(tn , ·) l’erreur commise sur la solution. Plus précisément, on montre
que si on pose l’énergie discrète de l’erreur
n+1/2 1 en+1 − enh 1
Eh (enh ) = h
+ ah (en+1 n
h , eh )
2 ∆t ℓ2 2
alors
n+1/2
Eh (enh ) > 0
et q q √ n
n+1/2 n 2 X
Eh (eh ) ≤ E 1/2 (e0 ) + √ ∆t kεkh kℓ2
2 1−α 2
k=1
Pn
où ∆t k=1 kε
k
kℓ2 est l’erreur de troncature du schéma sommée à chaque pas de temps.
102 C HAPITRE 3 : Approximation par différences finies

Nous n’entrons pas dans les détails de la démonstration car notre objectif était sim-
plement de vous montrer que l’écriture d’un schéma numérique, où les identités
remarquables permettent d’écrire une énergie discrète, sont un outil extrêmement
puissant pour l’obtention d’un schéma stable.
Chapitre 4

Discrétisation des équations

hyperboliques non linéaires 1D

4.1 Présentation du problème

Ce chapitre est consacré à l’approximation numérique par différences finies du pro-
blème de Cauchy :


 Trouver u(x, t) : R × R+ → R tel que
 ∂u ∂
+ f (u) = 0 , x ∈ R, t > 0, (4.1)

 ∂t ∂x

u(x, 0) = u0 (x), x ∈ R,

u2
où f est une fonction strictement convexe (par exemple, f (u) = pour l’équation
2
de Bürgers).

On se donne un pas de discrétisation en espace h > 0 et un pas de discrétisation

en temps ∆t > 0. On pose :
(
xj = jh, j ∈ Z,
tn = n∆t, n ∈ N,

et on se propose de calculer une approximation unj de u(xj , tn )1 où u est la solution de

(4.1).

Z xj + h
1 1 2
ou plutôt une approximation de u(x, tn )dx si l’on considère une approche par volumes
h xj − h
2
finis (cf. section 4.2.1)

103
104 C HAPITRE 4 : Discrétisation des équations hyperboliques non linéaires 1D

Dans toute la suite, nous nous intéressons à des schémas explicites à un pas de temps
et à trois pas d’espace qui s’écrivent

un+1
j = H(unj−1, unj , unj+1) , (4.2)

où H est une fonction régulière. Les schémas implicites sont difficiles à mettre en
oeuvre pour les problèmes non-linéaires et nous n’en parlerons pas ici.

4.2 Propriétés des schémas

4.2.1 Schémas sous forme conservative
Nous avons vu que la forme conservative :

∂u ∂
+ f (u) = 0 (4.3)
∂t ∂x
et la forme non conservative :
∂u ∂u
+ f ′ (u) =0
∂t ∂x
des lois de conservation ne sont équivalentes que pour les solutions u régulières du
problème et que la notion de solution faible repose sur la forme conservative de
l’équation. On suit le même principe pour le calcul numérique d’une solution faible ;
on utilisera un schéma « sous forme conservative ».

Le principe de base des schémas conservatifs Zest la méthode des volumes finis qui
1 xj+ 21
vise à trouver une approximation de uj (t) = u(x, t) dx à tous les temps t = tn
h x 1
j− 2

h
où xj± 1 = xj ± . Pour cela, on intègre tout d’abord l’équation (4.3) entre xj− 1 et xj+ 1
2 2 2 2

(et on divise par h), Z

d 1 xj+ 21 ∂
uj (t) + f (u)dx = 0
dt h x 1 ∂x
j− 2

ou de manière équivalente

d f (u(xj+ 1 , t)) − f (u(xj− 1 , t))

uj (t) + 2 2
= 0. (4.4)
dt h
Le schéma numérique conservatif est une approximation de l’équation (4.4) : on ap-
proche la quantité f (u(xj+ 1 , t)) (resp. f (u(xj− 1 , t))) en utilisant les ul (t). Un choix na-
2 2
turel consiste à prendre f (u(xj+ 1 , t)) ≈ g(uj (t), uj+1(t)) et f (u(xj− 1 , t)) ≈ g(uj−1(t), uj (t))
2 2
où g est le flux numérique qu’il faut alors définir. Après une discrétisation en temps
4.2 Propriétés des schémas 105

explicite, on obtient le schéma (4.5).

D ÉFINITION 4.2.1 (S CHÉMA CONSERVATIF )

On dit que le schéma (4.2) peut se mettre sous forme conservative s’il existe une
fonction g régulière telle que l’on ait :
∆t  n n
un+1
j = unj − g(uj , uj+1) − g(unj−1, unj ) (4.5)
h
La fonction g est appelée le flux numérique du schéma (elle est définie à une
constante additive près).

On utilisera souvent la notation abrégée suivante

 
un+1
j = u n
j − α g n
j+ 1
− g n
j− 1
.
2 2

où : 
 α ∆t
= ,
h
.
 gn = g unj , unj+1 .
j+ 21

L’importance de la forme conservative pour un schéma sera illustrée par le théorème

de Lax-Wendroff (théorème 4.2.2) énoncé plus loin. Mais avant tout, il nous faut pré-
ciser sous quelle condition, portant sur le flux numérique, un schéma sous forme
conservative est consistant avec l’équation (4.3).

4.2.2 Consistance
Les notions d’erreur de troncature, de consistance et d’ordre sont définies comme
dans le cas linéaire. L’erreur de troncature du schéma (4.5) s’écrit
1
1 n 
εnj = ū n+1 n
− ūj + n
ḡ 1 − ḡj− 1 ,
∆t j h j+ 2 2

où 
 ūnj = ū(xj , tn ).
 ḡ n 1 = ḡ(ūnj , ūnj+1),
j+ 2

et ū est la solution régulière de l’équation (4.3).

Remarquons en particulier que l’erreur de troncature se calcule à l’aide de solutions

régulières de l’équation (4.3). Ces notions n’auront donc pas de sens au voisinage des
chocs.
106 C HAPITRE 4 : Discrétisation des équations hyperboliques non linéaires 1D

L EMME 4.2.1 (C ONDITION DE CONSISTANCE )

Un schéma mis sous la forme conservative (4.5) est consistant avec l’équation (4.3)
si et seulement si :
g(u, u) = f (u) + cte ∀u ∈ R. (4.6)
Il est alors d’ordre 1 au moins.

D ÉMONSTRATION : Soit un schéma de la forme :

 
un+1
j = u n
j − α g n
j+ 1 − g n
j− 1 .
2 2

Calculons son erreur de troncature

1  n+1  1 
εnj = ūj − ūnj + n
ḡj+ 1 − ḡ
n
1
j− 2
,
∆t h 2

où l’on a utilisé les notations définies précédemment. Or on a


 ∂ ū
 ūn+1
j − ūnj = ∆t (xj , tn ) + O(∆t2 ),
∂t

 ūn − ūn ∂ ū
j+m j = mh (xj , tn ) + O(h2 ) avec m = ±1.
∂x
∂g ∂g
D’où, si et désignent les 2 dérivées partielles de g :
∂u ∂v


∂g n n
∂ ū

 g ūnj , ūnj+1 = g ūnj , ūnj + ū , ū h (xj , tn ) + O(h2 ),
∂v j j ∂x


 g ūn , ūn = g ūn , ūn −
∂g
∂ ū
j−1 j j j ūnj , ūnj h (xj , tn ) + O(h2 ).
∂u ∂x
∆t
On en déduit enfin, pour α = fixé :
h
 
n ∂ ū n ∂g n n ∂g n n ∂ ū
εj = (xj , t ) + (ū , ū ) + (ū , ū ) (xj , tn ) + O(h).
∂t ∂v j j ∂u j j ∂x

Le schéma est donc consistant (et au moins d’ordre 1) si et seulement si

∂g ∂g
(u, u) + (u, u) = f ′ (u) ∀u ∈ R,
∂u ∂v
d’où l’égalité (4.6). 

Insistons sur le fait qu’il existe des schémas consistants et non conservatifs.

4.2.3 Linéarisation du schéma - Stabilité L2

La stabilité asymptotique d’un schéma non-linéaire se définit de la manière suivante :
4.2 Propriétés des schémas 107

D ÉFINITION 4.2.2 (S TABILITÉ ASYMPTOTIQUE )

Le schéma (4.2) est asymptotiquement stable pour la norme k.k si pour tout couple
de données initiales (u0 , v 0 ), le couple de solutions (un , v n ) vérifie

∀n ∈ N, kun − v n k ≤ C u0 − v 0

où la constante C est indépendante de u0 , v 0 , n, ∆t et h.

Contrairement au cas linéaire, la stabilité L2 d’un schéma non-linéaire ne peut plus

s’étudier de manière simple à l’aide de la transformée de Fourier. Cependant, on peut
obtenir un critère heuristique de stabilité par linéarisation.

Considérons un schéma sous la forme conservative (4.5) et soient un et v n les solu-

tions correspondant aux données initiales u0 et v 0 .
( n+1

uj = unj − α g(unj, unj+1) − g(unj−1, unj ) ,

(4.7)
vjn+1 = vjn − α g(vjn , vj+1
n n
) − g(vj−1 , vjn ) .
On pose alors
wjn = unj − vjn .
Si wjn est supposé très petit devant unj , on peut calculer la différence entre les deux
égalités (4.7) et développer le résultat par rapport aux puissances de w. Si on ne
conserve que les termes du premier ordre, on obtient :
 
n+1 n ∂g n n ∂g n
wj = wj − α (u , u ) − (u , u ) wjn
n
(4.8)
∂u j j+1 ∂v j−1 j

∂g n n n ∂g n n n
+ (u , u )w − (u , u )w .
∂v j j+1 j+1 ∂u j−1 j j−1
Relativement à w, il s’agit d’un schéma linéaire à coefficients variables. On dit qu’il
s’agit du schéma linéarisé associé au schéma (4.5).

L’étude de la stabilité du schéma (4.8) ne peut pas se faire comme au chapitre pré-
cédent par analyse de Fourier, puisque les coefficients du schéma sont ici variables.
Pour pallier cette difficulté, il nous faut faire une nouvelle hypothèse : nous suppo-
sons que unj varie lentement en fonction de j. On peut alors considérer que, au voisi-
nage d’un point xk , le schéma (4.8) « ressemble » beaucoup au schéma suivant :
  
n+1 n ∂g ∂g n ∂g n ∂g n
wj = wj − α − wj + w − w (4.9)
∂u ∂v ∂v j+1 ∂u j−1
∂g ∂g

où toutes les dérivées et sont calculées en unk , unk . Ce schéma étant à coeffi-
∂u ∂v
cients constants, sa stabilité s’analyse de manière classique.
108 C HAPITRE 4 : Discrétisation des équations hyperboliques non linéaires 1D

D ÉFINITION 4.2.3 (S TABILITÉ LINÉAIRE )

On dit que le schéma (4.5) est linéairement stable L2 si le schéma (4.9) est stable L2
pour tout k ∈ Z.

R EMARQUE 4.2.1 (U NE STABILITÉ NON GARANTIE )

Le critère de stabilité linéaire que nous venons d’expliciter n’a qu’une valeur indica-
tive. En effet, les approximations que nous avons faites ne sont pas justifiées au voi-
sinage des chocs ni même dans les zones où la dérivée devient trop importante. En
résumé, les schémas linéairement instables sont toujours à proscrire. En revanche, il
est seulement plausible, mais nullement garanti, qu’un schéma linéairement stable
ait quelques propriétés de stabilité.

R EMARQUE 4.2.2 (U NE STABILITÉ LINÉAIRE QUI PEUT ÊTRE CONDITIONNELLE )

Certains schémas sont linéairement stables sous une condition portant sur le rap-
port des pas ∆t et h. Tout comme dans le cas linéaire, cette condition CFL peut s’in-
terpréter géométriquement. En effet, nous avons vu au cours précédent que, pour
l’équation
∂u ∂
+ f (u) = 0,
∂t ∂x
la connaissance de u à l’instant tn dans l’intervalle [xj−1 , xj+1 ] détermine la solution
dans l’ensemble du cône :

xj−1 + A(t − tn ) ≤ x ≤ xj+1 − A(t − tn ) (4.10)

 
n
où A = sup |a(ξ)|; |ξ| ≤ sup|u(x, t )| . La condition (4.10) signifie que le point (xj , tn+1 )
x∈R
doit se trouver dans ce cône.

4.2.4 Vers la convergence : le théorème de Lax-Wendroff

Nous allons maintenant démontrer un résultat important concernant les schémas
sous forme conservative. Il exprime le fait que si un schéma sous forme conserva-
tive consistant avec l’équation (4.3) converge, c’est nécessairement vers une solution
faible du problème.

On suppose que la donnée initiale u0 du problème (4.1) a été discrétisée comme suit :
Z
0 1 xj+ 12 0
uj = u (x) dx ∀j ∈ Z,
h x 1
j− 2

où :  
1
xj+ 1 = j+ h ∀j ∈ Z.
2 2
4.2 Propriétés des schémas 109

On construit alors, pour ∆t et h donnés, une suite un = unj j∈Z
par un schéma de la
forme
n+1 n ∆t  n n

uj = uj − g 1 − gj− 1 , (4.11)
h j+ 2 2

consistant avec l’équation (4.3), c’est-à-dire

g(u, u) = f (u) , ∀u ∈ R.

Notons qu’on a omis la constante pour simplifier la démonstration. On note u∆ la

fonction définie par :

u∆ (x, t) = unj pour xj− 1 ≤ x < xj+ 1 et tn ≤ t < tn+1 .

2 2

A chaque couple (∆t, h) correspond une fonction u∆ . On a alors le théorème suivant.

T HÉORÈME 4.2.2 (L AX -W ENDROFF )

Supposons qu’il existe une suite (∆tk , hk )k∈N , une suite de fonctions (u∆k )k∈N asso-
ciées, et une fonction u telles que, pour tout domaine K borné de R × [0, +∞[:

(i) u ∈ L∞ (K),
(ii) u∆k −→ u dans L∞ (K) quand k −→ +∞.

Alors u est une solution faible du problème (4.1).

D ÉMONSTRATION : Soit ϕ ∈ D(R2 ). On doit montrer que

Z +∞ Z +∞   Z +∞
∂ϕ ∂ϕ
u + f (u) dx dt + u0 (x)ϕ(x, 0) dx = 0 . (4.12)
x=−∞ t=0 ∂t ∂x x=−∞

Posons ϕnj = ϕ(xj , tn ), pour j ∈ Z, n ∈ N. En multipliant (4.11) par ϕnj et en sommant, on

obtient (on omet l’indice k pour alléger les notations) :
X  X 
h un+1
j − u n
j ϕn
j + ∆t g n
j+ 1 − g n
j− 1 ϕnj = 0.
2 2
j,n j,n

Par « sommation par parties », on obtient :

X   X X

h un+1
j
n n+1
ϕj − ϕj −h u0j ϕ0j + ∆t n
gj+ n n
1 ϕj − ϕj+1 = 0. (4.13)
2
j,n j

Si on note ϕ∆ et g∆ les fonctions définies par :


 ϕ∆ (x, t) = ϕnj avec xj− 1 ≤ x < xj+ 1 et tn ≤ t < tn+1 ,
2 2

 g∆ (x, t) n
= gj+ 1 avec xj ≤ x < xj+1 et tn ≤ t < tn+1 ,
2
110 C HAPITRE 4 : Discrétisation des équations hyperboliques non linéaires 1D

l’égalité (4.13) s’écrit aussi :

Z +∞ Z +∞ Z +∞
ϕ∆ (x, t − ∆t) − ϕ∆ (x, t)
u∆ (x, t) dx dt − u∆ (x, 0)ϕ∆ (x, 0)dx
x=−∞ t=∆t ∆t x=−∞
Z Z
+∞ +∞ ϕ∆ (x − h2 , t) − ϕ∆ (x + h2 , t)
+ g∆ (x, t) dx dt = 0 .
x=−∞ t=0 h

On vérifie alors facilement que la première intégrale converge vers :

Z +∞ Z +∞
∂ϕ
− u(x, t) (x, t)dx dt
x=−∞ t=0 ∂t

et la seconde vers : Z +∞
u0 (x)ϕ(x, 0)dx.
x=−∞

On remarque ensuite que g∆ vérifie par définition :

    
h h
g∆ (x, t) = g u∆ x − , t , u∆ x + , t ∀x, ∀t.
2 2

Donc g∆ converge uniformément sur tout compact vers la fonction g(u(x, t), u(x, t)) qui, d’après
la consistance du schéma est identiquement égale à f (u(x, t)). La troisième intégrale converge
donc vers : Z +∞ Z +∞
∂ϕ
− f (u(x, t)) (x, t)dx dt,
x=−∞ t=0 ∂x
d’où le résultat (4.12). 

Ce théorème nous assure que les schémas sous forme conservative sont adaptés au
calcul de solutions faibles. Notons qu’un schéma qui ne peut pas s’écrire sous forme
conservative peut très bien converger vers une fonction qui n’est pas solution faible
du problème.

4.2.5 Schémas entropiques

Un schéma sous forme conservative convergent ne converge pas toujours vers la
solution entropique du problème (voir exemple le schéma de Murman-Roe section
4.3.1).

Nous allons maintenant énoncer un critère (malheureusement difficile à vérifier en

pratique) qui assure la convergence (si cette dernière a lieu) vers la solution entro-
pique.
4.2 Propriétés des schémas 111

On rappelle tout d’abord que la solution entropique est caractérisée par le fait qu’elle
satisfait la condition d’entropie

∂ ∂
U(u) + F (u) ≤ 0 (4.14)
∂t ∂x

pour toute entropie (U, F ), c’est-à-dire pour tout couple de fonctions C 1 , U étant stric-
tement convexe et F étant définie par :

F ′ (v) = U ′ (v)f ′ (v) ∀v ∈ R.

Nous sommes en mesure de définir la notion de consistance d’un schéma avec une
condition entropique.

D ÉFINITION 4.2.4 (S CHÉMA CONSISTANT AVEC UNE CONDITION D ’ ENTROPIE )

Le schéma (4.5) est dit consistant avec la condition d’entropie (4.14) si et seule-
ment si il existe une fonction G(u, v) telle que :

(i) G est consistante avec le flux d’entropie F ,

G(u, u) = F (u) , ∀u ∈ R ;


(ii) Pour toute solution unj de (4.5), on a
 
Ujn+1 ≤ Ujn − α Gnj+ 1 − Gnj− 1
2 2

où l’on a posé :

 Ujn = U(unj ),

∀j ∈ Z, ∀n ∈ N .
 Gn 1 = G unj , unj+1 ,
j+ 2

La fonction G est appelée le flux d’entropie numérique.

On peut alors établir, par une démarche tout à fait analogue à celle que nous avons
suivie pour démontrer le théorème de Lax-Wendroff, le théorème suivant.

T HÉORÈME 4.2.3
Sous les hypothèses du théorème 4.2.2, si de plus le schéma considéré est
consistant avec toute condition d’entropie, alors la limite u est l’unique solution
entropique du problème (4.1).
112 C HAPITRE 4 : Discrétisation des équations hyperboliques non linéaires 1D

Ceci conduit finalement à définir les schémas entropiques.

D ÉFINITION 4.2.5 (S CHÉMA ENTROPIQUE )

On dit que le schéma (4.5) est entropique s’il est consistant avec toute condition
d’entropie.

Nous allons étudier une classe particulière de schémas entropiques : les schémas mo-
notones.

4.2.6 Schémas monotones

Condition de monotonie

On a vu dans le cours précédent que l’opérateur qui à u0 associe la solution entro-

pique du problème de Cauchy (4.1) est monotone. Autrement dit, si u et v sont les
solutions pour les données initiales u0 et v 0 , on a :

u0(x) ≥ v 0 (x) p.p. x ∈ R ⇒ u(x, t) ≥ v(x, t) p.p. x ∈ R, t > 0.

Il est naturel de s’intéresser aux schémas possédant la même propriété.

D ÉFINITION 4.2.6 (S CHÉMA MONOTONE )

Un schéma est dit monotone si
 n 
uj ≥ vjn , ∀j ∈ Z ⇒ un+1
j ≥ vjn+1 , ∀j ∈ Z .

Nous disposons d’un critère simple pour voir si un schéma est monotone ou pas.

L EMME 4.2.4 (C RITÈRE DE MONOTONIE )

Le schéma

un+1
j = H unj−1, unj , unj+1 (4.15)
est monotone si et seulement si H est une fonction croissante de chacun de ses
arguments.

D ÉMONSTRATION : La démonstration de ce lemme est immédiate. 

En pratique, la fonction H dépend des pas de discrétisation ∆t et h et la monotonie

∆t
n’est réalisée que sous une condition sur α = (condition de type CFL).
h
4.2 Propriétés des schémas 113

Propriétés des schémas monotones

Nous pouvons déjà énoncer quelques propriétés immédiates des schémas mono-
tones.

L EMME 4.2.5 (P ROPRIÉTÉS DES SCHÉMAS MONOTONES )

Soit un schéma monotone de la forme (4.15) et tel que H(u, u, u) = u pour tout
u ∈ R. Alors :

(i) il préserve la monotonie i.e.

(
∀j, unj ≤ unj+1 ⇒ ∀j, un+1
j ≤ un+1
j+1 ,

∀j, unj ≥ unj+1, ∀j ⇒ ∀j, un+1

j ≥ un+1
j+1 ,

(ii) il vérifie le principe du maximum suivant :

∀k ∈ Z inf unj ≤ un+1

k ≤ sup unj ,
j∈Z j∈Z

(iii) il est borné en norme l∞ :

kun+1 kℓ∞ ≤ kun kℓ∞

où kun kℓ∞ = sup|unj |.

j∈Z

D ÉMONSTRATION :

(i) il suffit de poser vjn = unj+1 (puis vjn = unj−1 ) et d’appliquer la définition.

(ii) De même, il suffit de poser : vjn = inf unk (puis vjn = sup unk ) ∀j et d’appliquer la défini-
k∈Z k∈Z
tion en utilisant l’identité H(u, u, u) = u.

Enfin, (iii) est une conséquence directe de (ii). 

En particulier, si la condition initiale est monotone, la solution approchée calculée à

l’aide d’un schéma monotone ne comportera jamais d’oscillations.

Une propriété essentielle des schémas monotones est donnée par le théorème sui-
vant que nous admettrons :
114 C HAPITRE 4 : Discrétisation des équations hyperboliques non linéaires 1D

T HÉORÈME 4.2.6 (U N RÉSULTAT DE CONVERGENCE )

Soit un schéma sous forme conservative consistant avec (4.3) et monotone. Alors,
de toute suite hk tendant vers 0, on peut extraire une sous-suite telle que la suite
de fonctions associées (u∆k ) converge presque partout vers une solution faible u
du problème (4.1) (les notations sont celles du théorème 4.2.2).

Autrement dit, les schémas monotones sont, en un certain sens, convergents. De

plus, la limite est bien la solution entropique du problème.
On a en effet le théorème suivant.

T HÉORÈME 4.2.7 (L IEN SCHÉMA ENTROPIQUE - SCHÉMA MONOTONE )

Un schéma sous forme conservative consistant avec (4.3) et monotone est entro-
pique.

D ÉMONSTRATION : On considère un schéma de la forme suivante :

un+1
j = H(unj−1 , unj , unj+1 )

où H(u, v, w) = v − α{g(v, w) − g(u, v)}. On doit montrer qu’il est consistant avec toute condi-
tion d’entropie, mais on admettra qu’il suffit de considérer les entropies de Kruzkov qui sont
de la forme : (
U (u) = |u − k|,
k ∈ R.
F (u) = signe(u − k)(f (u) − f (k)),
D’après la définition 4.2.4, il nous faut trouver une fonction G(u, v) telle que l’on ait :

G(u, u) = signe (u − k)(f (u) − f (k)) ∀u ∈ R (4.16)

et


|un+1
j − k| ≤ |unj − k| − α G(unj , unj+1 ) − G(unj−1 , unj ) (4.17)
pour toute solution (unj ) du schéma. Posons alors :

G(u, v) = g(u ∨ k, v ∨ k) − g(u ∧ k, v ∧ k)

où l’on a noté :
a ∨ b = max(a, b) et a ∧ b = min(a, b).
On a alors, puisque le schéma est consistant avec (4.3) :

G(u, u) = g(u ∨ k, u ∨ k) − g(u ∧ k, u ∧ k)

= f (u ∨ k) − f (u ∧ k)
= signe (u − k)(f (u) − f (k)) .

L’identité (4.16) est donc satisfaite.

4.2 Propriétés des schémas 115

Posons alors : vjn = unj ∨ k. Comme le schéma est monotone :

  

 vjn ≥ unj ∀j ⇒ H vj−1 n , vn , vn n+1
j j+1 ≥ uj
 

 vjn ≥ k ∀j ⇒ H vj−1 n , vn , vn
j j+1 ≥ H(k, k, k) = k

d’où finalement :
H(unj−1 ∨ k, unj ∨ k, unj+1 ∨ k) ≥ un+1
j ∨k. (4.18)
On démontrerait de même l’inégalité suivante :

H(unj−1 ∧ k, unj ∧ k, unj+1 ∧ k) ≤ un+1

j ∧ k. (4.19)

Mais comme |u − k| = u ∨ k − u ∧ k, on déduit de (4.18) et (4.19) que :

|un+1
j − k| ≤ H(unj−1 ∨ k, unj ∨ k, unj+1 ∨ k) − H(unj−1 ∧ k, unj ∧ k, unj+1 ∧ k),

et on vérifie aisément que cette dernière inégalité coïncide avec (4.17). 

Malheureusement, la portée des schémas monotones est limitée par leur ordre :

T HÉORÈME 4.2.8 (O RDRE DES SCHÉMAS MONOTONES )

Un schéma sous forme conservative consistant avec (4.3) et monotone est d’ordre
1 (exactement).

D ÉMONSTRATION : On considère comme précédemment un schéma de la forme :

un+1
j = H(unj−1 , unj , unj+1 ),

où H(u, v, w) = v − α{g(v, w) − g(u, v)}.

1) Soit ū une solution régulière de (4.3). En poursuivant le calcul amorcé lors de la démons-
tration du lemme 4.2.1, on vérifie aisément que l’erreur de troncature s’écrit
 n    
n ∆t ∂ 2 ū h2 ∂ ∂g ∂g ∂ ū n
εj = + α (ū, ū) − (ū, ū) + O(h2 ) .
2 ∂t2 j 2∆t ∂x ∂v ∂u ∂x j

En dérivant (4.3) nous avons

       
∂ 2 ū ∂ ∂ ∂ ∂ ū ∂ ∂ ∂ 2 ∂ ū
= − f (ū) = − c(ū) = c(ū) f (ū) = c(ū) .
∂t2 ∂t ∂x ∂x ∂t ∂x ∂x ∂x ∂x

Ainsi, l’erreur de troncature s’écrit finalement :

  
n ∆t ∂ ∂ ū n
εj = β(ū, α) + O(h2 )
2 ∂x ∂x j

où  
2 1 ∂g ∂g
β(ū, α) = c(ū) − (ū, ū) − (ū, ū) .
α ∂u ∂v
116 C HAPITRE 4 : Discrétisation des équations hyperboliques non linéaires 1D

2) Le schéma étant monotone, on a :

∂H ∂g
∀u ∈ R, (u, u) = α (u, u) ≥ 0
∂u ∂u  
∂H ∂g ∂g
∀u ∈ R, (u, u) = 1 − α (u, u) − (u, u) ≥ 0
∂v ∂u ∂v
∂H ∂g
∀u ∈ R, (u, u) = −α (u, u) ≥ 0.
∂w ∂v
On a donc finalement  ∂g

 (u, u) ≥ 0

 ∂u


∂g
∀u ∈ R, (u, u) ≤ 0 (4.20)

 ∂v



 ∂g (u, u) − ∂g (u, u) ≤ 1
∂u ∂v α
3) Le schéma étant consistant, on a :

∀u ∈ R, g(u, u) = f (u) + cste,

d’où
∂g ∂g
∀u ∈ R, (u, u) + (u, u) = c(u). (4.21)
∂u ∂v
De (4.20) et (4.21), on déduit finalement :
 2  
2 ∂g ∂g 1 ∂g ∂g
c(u) ≤ (u, u) − (u, u) ≤ (u, u) − (u, u) .
∂u ∂v α ∂u ∂v

Donc β(ū, α) est toujours négatif ou nul.

Pour que le schéma soit d’ordre supérieur ou égal à 2, il faut que β(u, α) soit nul pour tout u.
D’après ce qui précède, il faut pour cela que l’on ait :
 2  2
∂g ∂g ∂g ∂g
∀u ∈ R, (u, u) − (u, u) = (u, u) + (u, u) ,
∂u ∂v ∂u ∂v

et
∂g ∂g 1
(u, u) − (u, u) = .
∂u ∂v α
Ceci n’est vrai que si :
∂g 1 ∂g
(u, u) = et (u, u) = 0 ,
∂u α ∂v
où
∂g 1 ∂g
(u, u) = − et (u, u) = 0 .
∂v α ∂u
Plaçons nous par exemple dans le premier cas. D’après (4.21), on a alors :

1
∀u ∈ R, a(u) = .
α
4.3 Exemples 117

1
Il s’agit donc du cas linéaire de l’équation de transport avec la vitesse · Le schéma s’écrit :
α
un+1
j = unj−1 ,

et il est exact car il coïncide avec la méthode des caractéristiques.

Ainsi nous avons démontré qu’un schéma monotone est d’ordre 1 exactement, sauf si f (u) =
∆t
cu et c = 1 auquel cas il est d’ordre ∞.
h


4.3 Exemples
Nous allons illustrer les notions présentées précédemment en passant en revue les
principaux exemples de schémas pour les équations hyperboliques non linéaires en
dimension 1.

4.3.1 Quelques schémas d’ordre 1

Le schéma de Lax-Friedrichs
Le schéma de Lax dans le cas linéaire s’étend de manière naturelle au cas non-linéaire
sous la forme :
1 n
α
un+1
j = uj−1 + unj+1 − f (unj+1) − f (unj−1) .
2 2
Le schéma de Lax-Friedrichs possède les propriétés suivantes.

• Conservatif. On vérifie aisément qu’il peut se mettre sous la forme conservative

avec pour flux numérique :
 
1 1
g(u, v) = f (u) + f (v) + (u − v) .
2 α

• Consistant. Le schéma de Lax-Friedrichs est consistant d’ordre 1 au moins. En

effet,  
1 1
g(u, u) = f (u) + f (u) − (u − u) = f (u).
2 α
On verra ensuite qu’il est monotone donc qu’il est exactement d’ordre 1.

• Linéairement stable. Pour le schéma de Lax-Friedrichs, (4.8) s’écrit :

1 n  α
 n 
wjn+1 = n
wj−1 + wj+1 − c unk wj+1 n
− wj−1
2 2
118 C HAPITRE 4 : Discrétisation des équations hyperboliques non linéaires 1D

où c désigne la dérivée de f . C’est exactement

le schéma de Lax-Friedrichs li-
néaire avec une vitesse égale à c uk . On sait que ce schéma est stable L2 sous
n

la condition CFL :

|c unk |∆t
≤ 1. (4.22)
h
Le schéma non-linéaire de Lax-Friedrichs est donc linéairement stable L2 si
(4.22) est vérifiée pour tout k. Comme cette condition dépend de n, le pas de
temps doit être ajusté à chaque itération. Autrement dit, on doit avoir :

∆t(n)
sup|c unk | ≤ 1 ∀n ∈ N . (4.23)
k∈Z h

• Monotone. On a
1 α
H(u, v, w) = (u + w) − (f (w) − f (u)).
2 2
D’où :
∂H 1 α
(u, v, w) = + c(u),
∂u 2 2
∂H
(u, v, w) = 0,
∂v
∂H 1 α
(u, v, w) = − c(w).
∂w 2 2
Par conséquent, le schéma de Lax-Friedrichs est monotone sous la condition
CFL (4.23).

Le schéma d’Engquist-Osher
Il est défini comme suit :
( Z unj Z unj+1 )
α
un+1
j = unj − f (unj+1) − f (unj−1) + |c(ξ)|dξ − |c(ξ)| dξ
2 un
j−1 u n
j

où l’on a noté a la dérivée de f.

On retrouve le schéma décentré à gauche dans les zones où f est croissante et à droite
lorsque f est décroissante.

Le schéma d’Engquist-Osher possède les propriétés suivantes.

• Conservatif. Il peut s’écrire sous forme conservative à condition de poser :

 Z v 
1
g(u, v) = f (u) + f (v) − |c(ξ)|dξ .
2 u
4.3 Exemples 119

• Consistant. On vérifie aisément que f (u, u) = g(u). Le schéma d’Engquist-Osher

est donc consistant et d’ordre 1 au moins. On va voir qu’il est monotone donc
qu’il est exactement d’ordre 1.
• Linéairement stable. Avec une démarche similaire à celle effectuée pour le schéma
de Lax-Friedrichs, on montre qu’il est linéairement stable L2 sous la condition
CFL (4.23).
• Monotone. On a
 Z v Z w 
α
H(u, v, w) = v − f (w) − f (u) + |c(ξ)|dξ − |c(ξ)|dξ
2 u v

D’où :
∂H α α
(u, v, w) = c(u) + |c(u)|
∂u 2 2
∂H
(u, v, w) = 1 − α|c(v)|
∂v
∂H α α
(u, v, w) = − c(w) + |c(w)| .
∂w 2 2
Par conséquent, le schéma d’Engquist-Osher est monotone sous la condition
CFL (4.23).

Le schéma de Murman-Roe
Il est naturel de vouloir étendre au cas non-linéaire le schéma décentré. Mais il y a
une difficulté puisque nous avons vu que le schéma décentré à gauche est stable sous
une condition CFL si la vitesse c est positive. Si celle-ci est négative, il est impératif de
décentrer à droite. Autrement dit, le schéma non-linéaire suivant :

un+1
j = unj − α(f (unj ) − f (unj−1))

est linéairement instable dans toute zone où c(u) est négatif.

Pour remédier à cette difficulté, on peut par exemple considérer le schéma suivant,
appelé schéma de Murman-Roe



un+1
j = unj − α g unj , unj+1 − g unj−1, unj

où
1
g(u, v) = (f (u) + f (v) + |c(u, v)|(u − v)),
2
avec 
 f (u) − f (v) si u 6= v
c(u, v) = u−v .
 ′
f (u) si u = v.
120 C HAPITRE 4 : Discrétisation des équations hyperboliques non linéaires 1D

Tout comme le schéma précédent, ce schéma coïncide avec le schéma décentré lorsque
f est monotone.

Le schéma de Murman-Roe possède les propriétés suivantes.

• Conservatif. On remarque que le schéma de Murman-Roe est sous forme conser-

vative en prenant pour le flux numérique la fonction g définie précédemment.

• Consistant. Il est consistant car g(u, u) = f (u) pour tout u ∈ R.

• Linéairement stable dans certains cas. Il n’est pas différentiable aux points
(u, v) pour lesquels c(u, v) s’annule. On ne peut donc pas linéariser le schéma
au voisinage de ces points. On peut seulement affirmer que, sous la condition
CFL (4.23), le schéma de Murman-Roe est linéairement stable L2 dans les zones
où f est monotone.

• Non entropique et donc non monotone. Supposons enfin que l’on veuille ré-
soudre numériquement à l’aide de ce schéma le problème (4.1) avec la donnée
initiale : (
ug si x ≤ 0,
u0 (x) =
ud si x > 0,
où ug et ud sont tels que : f (ug ) = f (ud ). Pour la donnée initiale discrète suivante
(
ug si j ≤ 0,
u0j =
ud si j ≥ 1,

le schéma de Murman-Roe fournit une solution discrète stationnaire, i.e. :

unj = u0j ∀n ∈ N.

En effet, on a :

g(ug , ug ) = g(ud, ud ) = g(ug , ud ) = g(ud, ug ) = f (ug ) = f (ud ).

Ainsi, le schéma de Murman-Roe permet de calculer une solution du problème

donnée par
u(x, t) = u0 (x) ∀t ≥ 0.
C’est une solution faible car la relation de Rankine-Hugoniot est satisfaite.
Si ug ≥ ud , il s’agit bien de la solution entropique. En revanche, si ug < ud , on
sait que la solution entropique correspond à une onde de détente et non à ce
choc stationnaire. Le schéma de Murman-Roe ne permet pas de « sélectionner »
la solution entropique.
4.3 Exemples 121

4.3.2 Un exemple de schéma d’ordre 2 : le schéma de Lax-Wendroff

Nous allons construire un schéma d’ordre 2.

Un développement limité en temps nous donne (avec les notations habituelles) :

∂ ū (∆t)2 ∂ 2 ū
ūn+1
j = ūnj n
+ ∆t (xj , t ) + 2
(xj , tn ) + O(∆t3 ). (4.24)
∂t 2 ∂t
Mais de plus, on déduit de l’équation (4.3) les identités suivantes :
     
∂ 2 ū ∂ ∂ ∂ ∂ ū ∂ ∂
= − f (ū) = − c(ū) = c(ū) f (ū) .
∂t2 ∂t ∂x ∂x ∂t ∂x ∂x
L’équation (4.24) devient
 n   n
ūn+1
j − ūnj ∂ ∆t ∂ ∂
=− f (ū) + c(ū) f (ū) + O(∆t2 ),
∆t ∂x j 2 ∂x ∂x j

où la notation [ . ]nj signifie que l’on calcule la quantité entre les crochets au point xj
et au temps tn .

Nous allons ensuite discrétiser les termes entre crochets. Le premier terme se traite
de manière habituelle (discrétisation ordre 2)
 n
∂ f (ūnj+1) − f (ūnj−1)
(f (ū)) = + O(h2 ).
∂x j 2h

Pour le second, on applique deux fois de suite un schéma centré afin d’obtenir une
approximation d’ordre 2 :
  n ( n  n )
∂ ∂ 1 ∂ ∂
c(ū) f (ū) = c(ū) f (ū) − c(ū) f (ū) + O(h2 ),
∂x ∂x j 2h ∂x 1
j+ 2 ∂x 1
j− 2
 
1 n
f (ūnj+1) − f (ūnj ) n
f (ūnj ) − f (ūnj−1)
= cj+ 1 − cj− 1 + O(h2 ),
2h 2 h 2 h
 ūn + ūn 
j j+1
où l’on a posé : cnj+ 1 = c .
2 2
On obtient finalement le schéma du second ordre suivant
α n
α2  n 
un+1 = u n
j − f − f n
+ c 1 (f n
− f n
) − c n
(f
j− 21 j
n
− f n
) ,
j
2 j+1 j−1
2 j+ 2 j+1 j j−1

avec  n
 fj
 = f (unj )
 n n 
 n
u j + u j+1
 cj+ 1 = c
2 2
122 C HAPITRE 4 : Discrétisation des équations hyperboliques non linéaires 1D

• Conservatif. Le schéma de Lax-Wendroff peut s’écrire sous forme conservative

et son flux numérique est donné par
   
1 u+v
g(u, v) = f (u) + f (v) − αc (f (v) − f (u)) .
2 2

• Consistant. On vérifie facilement que g(u, u) = f (u).

• Linéairement stable. Sa version linéarisée correspond au schéma linéarisé de

Lax-Wendroff. Il est donc linéairement stable L2 sous la condition CFL (4.23).

• Non entropique et donc non monotone. Enfin, on peut montrer qu’il n’est pas
entropique en suivant la même démarche que pour le schéma de Murman-Roe.

4.3.3 Un exemple fondamental - Le schéma de Godounov

La démarche de construction de ce schéma est originale : elle repose sur la résolution
exacte de problèmes de Riemann locaux. La construction de un+1 à partir de un se fait
en deux étapes :

- Première étape : On note un (x) la fonction définie par

un (x) = unj si xj− 1 ≤ x < xj+ 1 , (4.25)

2 2

 
1
où xj+ 1 = j + h, ∀j ∈ Z. On calcule alors v n solution entropique du problème
2 2
suivant : 
 ∂v + ∂ f (v) = 0, x ∈ R, t ∈ [tn , tn+1 ]
∂t ∂x (4.26)

v(x, tn ) = un (x), x∈R

On montrera (cf. lemme 4.3.1) que ce problème peut être résolu de façon exacte si la
condition suivante est vérifiée :

∆t(n) 1
sup|a (unk )| ≤ · (4.27)
k∈Z ∆x 2

- Deuxième étape : On définit un+1 comme projection de v n sur les fonctions constantes
par morceaux :
Z
n+1 1 xj+ 12 n
uj = v (x) dx, ∀j ∈ Z. (4.28)
h x 1
j− 2
4.3 Exemples 123

L EMME 4.3.1
Sous la condition (4.27), la solution du problème (4.26) est donnée par :
 
n
x − xj+ 1
v (x, t) = wR 2
; unj , unj+1 si xj < x < xj+1 , (4.29)
t − tn
x 
où wR ; ug , ud est la solution du problème de Riemann pour les états (ug , ud ).
t

D ÉMONSTRATION : Le problème (4.26) consiste en une succession de problèmes de Riemann

de la forme

 ∂v ∂

 + f (v) = 0, x ∈ [xj , xj+1 ] , t ∈ [tn , tn+1 ],
 ∂t ∂x
(Pj ) ( n
 uj si x < xj+ 1 ,

 v (x, t n) = 2
 un si x > x 1 .
j+1 j+ 2

La solution de ce problème pour x ∈ R est donnée par :

x − x 
j+ 21
v(x, t) = wR ; unj , unj+1 · (4.30)
t− tn

Il s’agit d’une onde de détente si unj ≤ unj+1 et d’un choc dans le cas contraire. La « juxta-
position » de solutions de la forme (4.30) définit bien la solution de (4.26) si les différents
problèmes n’ont pas « interagi » avant l’instant tn+1 . Autrement dit, il faut que la solution v de
chacun des problèmes (Pj ) vérifie :

(
v(xj , t) = unj ,
n n+1
∀t ∈ [t , t [,
v(xj+1 , t) = un+1
j ,

Ceci s’écrit encore :

  
 h n n+1 n

 wR − ; uj , uj , = uj ,
 2τ 
∀τ ∈ [0, ∆t[,

 h n n+1
 wR ;u ,u , = un+1
2τ j j j

et ceci est automatiquement vérifié sous la condition :

h



≥ max |c unj |, |c unj+1 | .
2∆t
124 C HAPITRE 4 : Discrétisation des équations hyperboliques non linéaires 1D

choc détente
tn+1

uj-1n ujn n
uj+1

tn
xj-1 xj- xj xj+ xj+1

Autrement dit, la vitesse de propagation de la solution du problème (Pj ) étant :


 
Vj = max |c unj |, |c un+1
j | ,

l’onde de choc ou l’onde de détente n’auront pas atteint les bornes de l’intervalle xj et xj+1
avant l’instant :
h
Tj = tn + Vj .
2
Il suffit donc de choisir t n+1 tel que : t n+1 ≤ Tj , ∀j ∈ Z . 

• Conservatif. Nous allons maintenant vérifier que le schéma de Godounov dé-

crit précédemment peut se mettre sous forme conservative.

L EMME 4.3.2
(i) L’application ξ → f (wR (ξ; u, v)) est continue en 0, quels que soient u et v.

(ii) Le schéma de Godounov peut s’écrire sous forme conservative et son flux
numérique est donné par :

g(u, v) = f (wR (0; u, v)) . (4.31)

D ÉMONSTRATION :
f (u) − f (v)
(i) On sait que l’application ξ → wR (ξ; u, v) est continue partout sauf en ξ =
u−v
si u > v. Elle est donc continue en 0 sauf si : u > v et f (u) = f (v). Or, on a dans ce cas :
wR (0− ; u, v) = u et wR (0+ ; u, v) = v. D’où : f (wR (0− ; u, v)) = f (u) = f (w) = f (wR (0+ ; u, v)).

(ii) Intégrons l’équation (4.26) sur la maille

h i  
xj− 1 , xj+ 1 × tn , tn+1 .
2 2
4.3 Exemples 125

On a : Z Z  
xj + 21 tn+1
∂v ∂
0= + f (v) dx dt.
xj − 21 tn ∂t ∂x
A priori, cette écriture n’a de sens que si la maille considérée n’est pas traversée par une ligne
de choc. Dans le cas contraire, il faudrait intégrer de part et d’autre du choc. Cependant,
comme v satisfait la relation de Rankine-Hugoniot, on vérifie aisément que l’intégrale sur
la ligne de choc issue de l’intégration par parties est nulle. Tout se passe donc comme si v
était régulière !

On a alors
Z xj+ 1 Z tn+1  
2 n+1

n
0= v(x, t ) − v(x, t ) dx + f (v(xj+ 1 , t)) − f (v(xj− 1 , t)) dt.
2 2
xj− 1 tn
2

D’où, en vertu de (4.28), (4.29) et du point (i) :

 

0 = h un+1
j − u n n n n n
j + ∆t f (wR (0; uj , uj+1 )) − f (wR (0; uj−1 , uj )) .

C’est ce qu’il fallait démontrer. 

Remarquons que le schéma de Godounov coïncide avec le schéma décentré si f est

monotone. En effet, si par exemple


c unj ≥ 0 , ∀j ∈ Z,

alors
wR (0; unj , unj+1) = unj , ∀j ∈ Z
et le schéma de Godounov s’écrit


un+1
j = unj − α f (unj ) − f (unj−1) .

• Consistant d’ordre 1. En utilisant (4.31) on remarque que g(u, u) = f (u) et donc

que le schéma de Godounov est consistant d’ordre 1 au moins. Comme il est
monotone (cf. ci-après), il est d’ordre 1 exactement.

• Monotone et donc entropique. Le schéma de Godounov est monotone. En ef-

fet, la première étape qui consiste en une résolution exacte d’un problème de
Cauchy est monotone (voir cours précédent). La seconde étape l’est triviale-
ment.