ÉLÉMENTS DE MODÉLISATION DÉTERMINISTE (M1)

SYLVIE BENZONI

L’objectif est de se familiariser avec quelques modèles classiques, aux domaines d’application
variés, et de voir quelques méthodes simples permettant leur simulation sur ordinateur. Cette partie
concerne exclusivement des modèles déterministes, de type équations différentielles (ordinaires puis
aux dérivées partielles).

Programme prévisionnel (dans le désordre):

1) Équations différentielles ordinaires. Équations d’Euler–Lagrange en mécanique, exemple du
pendule simple. Équation de van der Pol en électricité. Équation logistique, équations de Lotka-
Volterra en dynamique des populations. Comparaison avec les modèles discrets. Schéma d’Euler.
Portraits de phase.
2) Équations aux dérivées partielles. - Équation de Poisson en électrostatique. Équations des
fluides parfaits incompressibles et irrotationnels (potentiel, fonction de courant). - Équation de
transport. Équation des cordes vibrantes. Équation de continuité en mécanique des fluides, ap-
plication au trafic routier. - Équation de la chaleur. Équations de réaction-diffusion en biologie.
Simulation par différences finies. Conditions aux limites de Neumann et de Dirichlet.

1. Équations différentielles ordinaires

Ceci n’étant pas un cours d’équations différentielles, nous aborderons exclusivement des exemples
d’équations différentielles. Par ailleurs, les pré-requis en provenance d’autres disciplines seront
limités à leur strict minimum. Ce chapitre est divisé en trois parties, concernant trois domaines
d’application distincts (bien que l’on puisse leur trouver des points communs):
– la dynamique des populations, qui sera l’occasion d’étudier les propriétés qualitatives d’équa-
tions scalaires (modèles de Malthus et de Verhulst) et de systèmes de deux équations (modèles
prédateurs-proies de type Lotka-Volterra); on en profitera également pour introduire le schéma
d’Euler;
– la mécanique du point avec l’équation de l’oscillateur harmonique et celle du pendule simple,
dont on étudiera en détail le portrait de phase, puis la mécanique du solide avec le mouvement
d’une toupie,
– l’électricité avec les circuits RLC modélisés par l’équation de van der Pol, admettant une
solution périodique « non-triviale ».
Pour la dynamique des populations, strictement aucun pré-requis n’est nécessaire. Pour la méca-
nique du point il suffira de connaı̂tre la loi de Newton. Nous verrons en outre quelques éléments
concernant le formalisme lagrangien permettant la mise en équation de problèmes de mécanique
du solide, sur l’exemple de la toupie. Pour l’électricité, nous serons amenés à admettre les lois
classiques (Faraday) reliant différences de potentiel et intensités.
1.1. Dynamique des populations. Une référence particulièrement recommandée pour ce para-
graphe est le livre de Murray [7] (au moins le premier chapitre).
Un modèle mathématique en dynamique des populations, qu’il concerne des humains, du gibier,
des bactéries ou des végétaux, vise à prédire l’évolution future du nombre d’individus. Il est
Date: 2 mars 2011.
1
généralement construit sur la base de la situation au moment présent, voire de l’évolution passée,
et d’arguments plus ou moins empiriques.
1.1.1. Malthus. Le modèle le plus simple qui soit est dû à Malthus et s’exprime sous forme d’une
équation différentielle linéaire du premier ordre. Considérons une population dont le nombre d’in-
dividus est N(t ) à l’instant t . En pratique, il est impossible de compter le nombre d’individus à
chaque instant t , il faut donc imaginer qu’on les compte avec une périodicité donnée: un jour,
un mois, un an, ..., cela dépend de l’échelle de temps qui nous intéresse. Disons pour fixer les
idées que l’intervalle de temps entre chaque décompte (recensement exhaustif) est ∆t . Pendant cet
intervalle de temps, des individus naissent, meurent, arrivent ou partent (selon les flux migratoires
ou les prédateurs dans la nature, ou encore l’intervention du biologiste en laboratoire). Faisons
pour l’instant abstraction des arrivées et des départs. On appelle taux d’accroissement naturel de
la population entre t et t + ∆t la variation du nombre d’individus dûe aux naissances et aux décès,
divisée par ∆t et par le nombre total d’individus à l’instant « initial » t : la dimension physique
du taux d’accroissement naturel r est donc l’inverse d’un temps. Par définition,
N(t + ∆t ) − N(t )
(1) r = .
N(t ) ∆t
Notons que ce taux d’« accroissement » peut parfaitement être négatif. C’est en fait la différence
entre le taux de natalité et le taux de mortalité. Des valeurs typiques de r sont par exemple
13h −9h = 4h par an en France, et 10h −16h = −6h par an en Ukraine (source: INED/Banque
Mondiale). En général, r dépend évidemment de t et de ∆t . Cependant, l’hypothèse de base dans
le modèle de Malthus consiste à dire que r est constant. Dans ce cas, et si on connaı̂t la valeur
de r , on peut voir la relation (1) comme un moyen de calculer N(t k ) à tout instant t k = t 0 + k ∆t
connaissant N(t 0 ). En effet, (1) équivaut à
(2) N(t + ∆t ) = (1 + r ∆t ) N(t ) ,
ce qui veut simplement dire que la suite des valeurs (N(t k ))k∈N est géométrique de raison 1 + r ∆t .
Ceci n’est toutefois qu’une version discrète du modèle de Malthus, qui s’obtient en faisant tendre
∆t vers 0 dans (1). En supposant que ce passage à la limite est justifié, on obtient ainsi l’équation
différentielle ordinaire, linéaire du premier ordre:
dN
(3) = r N.
dt
Pour que cette équation soit correcte, il faut prendre garde à bien exprimer le temps t et l’inverse
de r dans la même unité: seconde, heure, année, siècle, peu importe a priori, pourvu qu’il soit
raisonnable de supposer r constant sur l’échelle de temps considérée. On connaı̂t parfaitement les
solutions de (3): elles sont de nature exponentielle, et plus précisément, toute solution de (3) est
donnée par
N(t ) = N(0) er t .
Autrement dit, une population régie par l’équation de Malthus (3) s’accroı̂t exponentiellement vite
si r > 0, et s’éteint exponentiellement vite si r < 0. Si l’on observe les données de la population
humaine mondiale, celle-ci semble plutôt relever, qualitativement au moins, du premier cas !
Sensibilité par rapport aux données initiales. Il est crucial de remarquer que l’équation (3), aussi
simple soit-elle, est extrêmement sensible aux données initiales si r > 0. En effet, la forme exponen-
tielle des solutions montre qu’une erreur (de mesure ou d’arrondi) à l’instant initial est multipliée
par er T au bout d’un temps T. Étant données les valeurs prises par la fonction exponentielle (pour
mémoire, e ≈ 2, 7, e10 ≈ 22000), s’il est raisonnable d’utiliser (3) comme modèle mathématique sur
des temps de l’ordre de 1/r , cela n’a aucun sens (quantitativement) pour des temps de l’ordre
de 10/r . Avec par exemple r = 4h par an, si l’on veut que l’erreur soit au pire multipliée par
2
10, cela limite la portée du modèle à ln(10) × 1000/4 ≈ 575 ans. Notez que le problème est ana-
logue avec le modèle discret (2): une erreur sur la valeur initiale de la suite géométrique de raison
1 + r ∆t est multipliée par (1 + r ∆t )k à l’instant t k = t 0 + k ∆t . Avec la valeur précédente r = 4h par
an et ∆t = 1 an, pour que l’erreur soit au pire multipliée par 10, il ne faut pas aller au delà de
ln(10)/ ln(1 + 4/1000) ≈ 577 ans.
Comparaison entre modélisation discrète et continue. Le modèle (2) s’obtient en discrétisant (3)
par le schéma d’Euler explicite. Un pas de discrétisation ∆t = 1 an peut paraı̂tre grossier mais
cela dépend de la valeur de r et de l’intervalle de temps auquel on s’intéresse. En effet, pour des
données initiales identiques, l’erreur relative à l’instant t k = t 0 + k ∆t entre la solution du modèle
continu (3) et la solution du modèle discret (2) est
εk := 1 − (1 + r ∆t )k e−r k ∆t .
De façon générale, si l’on fixe un temps T que l’on subdivise en K intervalles de longueur ∆t ,
l’erreur εK tend vers zéro lorsque K tend vers +∞ (c’est-à-dire ∆t tend vers zéro): ceci résulte
précisément de la formule d’Euler :
lim (1 + x/K)K = eK .
K→+∞

Bien sûr, à ∆t fixé, εk croı̂t avec k . Par exemple en prenant à nouveau r = 4h par an avec ∆t = 1,
on trouve
ε1 ≈ 8.10−6 , ε10 ≈ 8.10−4 , ε100 ≈ 6.10−2 , ε1000 ≈ 0, 9 .
Autrement dit, la cohérence entre le modèle discret et le modèle continu dure plus longtemps que
leur robustesse par rapport aux données initiales...
Introduction de « termes sources ». Nous avons parlé au début de possibles migrations ou inter-
ventions de l’expérimentateur. Si l’on tient compte du flux net d’individus (les entrants moins les
sortants) par unité de temps, disons f , l’équation (3) devient
dN
(4) =rN+ f .
dt
On dit que f est un terme source. De façon générale, si l’on autorise à la fois r et f à dépendre
du temps t , la solution de (4) est donnée par la formule de Duhamel
Rt Z t Rt
(5) N(t ) = e t 0 r (s)ds N(t 0 ) + e τ r (s)ds f ( τ) d τ .
t0

En particulier lorsque r et f sont constants, cette formule se réduit à

N(t ) = N(t 0 ) er (t −t0 ) + ( er (t −t0 ) − 1) f /r ,
où l’on voit que N(t ) s’annule en temps fini si f < −r N(t 0 ). Ce n’est évidemment pas surprenant:
si l’on soustrait suffisamment d’individus, la population finit par s’éteindre.
1.1.2. Verhulst. Au 19ème siècle, Verhulst proposa une version modifiée du modèle de Malthus,
avec un terme « auto-limitant », c’est-à-dire limitant la croissance lorsque la population devient
trop grande. Ceci revient en fait à ajouter un terme source f dépendant de N. L’hypothèse de
Verhulst est que ce terme dépend de façon quadratique de N, de sorte que la nouvelle équation
s’écrive
dN
(6) = r N (1 − N/K) ,
dt
où K représente la capacité de l’environnement à assurer la survie de l’espèce: tant que N n’a pas
atteint la valeur K, la population s’accroı̂t; si elle dépasse K à un certain moment, elle décroı̂t.
En fait, une analyse du portrait de phase de cette équation montre qu’une population initiale en
3
dessous de la capacité K ne peut jamais la dépasser si elle régie par (6), aussi appelée équation
logistique. Ceci est confirmé par la résolution explicite de cette équation. On trouve en effet
er t N(t 0 )
N(t ) = K .
K − N(t 0 ) + er t N(t 0 )
Vers les modèles prédateurs-proie. Une version raffinée de l’équation logistique est l’équation de
Ludwig:

dN B N2
(7) = r N (1 − N/K) − 2 .
dt A + N2
Elle prend en compte non seulement une auto-régulation quadratique mais aussi un terme source
s’annulant avec N et asymptotiquement constant lorsque la population est très grande. Il cor-
respond par exemple au prélèvement d’individus par des prédateurs, qui est faible lorsque la
population de proies est petite (les prédateurs vont chercher leur nourriture ailleurs), et limité par
la quantité de nourriture que les prédateurs sont capables d’absorber.
Adimensionnalisation. Afin d’étudier efficacement un modèle, il est bon de s’affranchir du choix
des unités de mesure, de sorte que les notions de quantités « petites » ou « grandes » prennent un
sens absolu. Pour cela, on procède à ce que l’on appelle une adimensionnalisation. Ceci consiste à
choisir des valeurs de référence (dans une unité de mesure quelconque) pour toutes les dimensions
physiques apparaissant dans le modèle, puis à définir de nouvelles inconnues et paramètres sans
dimension à partir de ces valeurs de référence. Cette approche a pour effet « secondaire » de réduire
le nombre de paramètres au strict minimum, d’après le théorème π de Buckingham (reposant sur
des arguments d’algèbre linéaire, le théorème du rang notamment, voir [6, p. 5]). Sur l’exemple
de l’équation (7), les dimensions « physiques » sont seulement au nombre de deux: le temps et
le nombre d’individus (même si ce nombre n’a pas vraiment une dimension physique, on peut
l’exprimer en milliers, millions ou milliards par exemple). Choisissons donc un temps de référence
T et un nombre d’individus de référence M. Les paramètres K et A représentent des nombres
d’individus. On leur associe donc naturellement les paramètres sans dimension K e := K/M et Ae :=
A/M. Le taux d’accroissement naturel r a la dimension de l’inverse du temps, on définit donc
re := Tr . Quant au paramètre B, il est homogène à d N/d t , ce qui nous conduit à poser Be := TB/M.
Enfin, l’inconnue N est un nombre d’individus. On considère donc l’inconnue sans dimension
Ne := N/M, vue comme une fonction du temps adimensionné e t := t /T. On a ainsi

dN
e T dN
= .
dt
e M dt
En substituant toutes ces nouvelles quantités dans (7), on obtient exactement la même équation,
avec des tildes:
dN
e B
eNe2
= re N
e (1 − N/
e K)
e − .
t
de e2 + N
A e2

À première vue, il semble qu’on n’ait rien gagné en terme de nombre de paramètres! Cependant,
on peut en supprimer deux en choisissant « astucieusement » M = A et T = A/B, ce qui nous
ramène à A
e = 1 et B e = 1. Ainsi, en posant q := K/A et en « oubliant les tildes pour simplifier »
(phrase classique dans cette approche), on obtient l’équation adimensionnalisée à seulement deux
paramètres:

dN N2
(8) = r N (1 − N/q) − .
dt 1 + N2
4
1.1.3. Lotka–Volterra. Les « vrais » modèles prédateurs-proies s’intéressent au nombre d’individus
dans les deux catégories de populations. Disons que N(t ) et P(t ) représentent respectivement le
nombre de proies et le nombre de prédateurs à l’instant t . Le système proposé par Volterra suppose
que:
• les proies ont un certain taux d’accroissement naturel indépendant du temps,
• le prélèvement de proies par les prédateurs est proportionnel au produit NP,
• les prédateurs ont un taux d’« accroissement » naturel indépendant du temps (ce taux étant
négatif car en l’absence de proie la population de prédateurs est amenée à dépérir),
• la présence de proies ajoute un terme source proportionnel au produit NP à l’équation sur P.
On obtient un système de deux équations différentielles ordinaires
dN


 = a N − b NP,
dt

(9)
 dP = −d P + c N P ,


dt
doté de plusieurs propriétés remarquables. Avant tout, le lecteur pourra à titre d’exercice1 adi-
mensionner le système (9) et montrer qu’il suffit d’un (et non quatre) paramètre sans dimension,
le système adimensionné s’écrivant
dN


 = N (1 − P) ,
dt

(10)
 dP = −δ P (1 − N) .


dt
Comme première propriété on voit que les axes du plan, N = 0 (absence de proies) et P = 0 (absence
de prédateurs), sont invariants par le flot (aussi bien de de (10) que de de (9), évidemment): cela
signifie que si N(t 0 ) = 0 alors N(t ) = 0 quel que soit t (et de même pour P); du point de vue de la
modélisation, c’est la moindre des choses (il ne peut y avoir de « génération spontanée » de proies
ou de prédateurs); du point de vue mathématique, c’est une conséquence de la partie unicité dans
le théorème de Cauchy–Lipschitz , qui s’applique sans problème ici car les seconds membres de (9)
sont des fonctions polynomiales de (N, P). On montre même plus par cet argument d’unicité: les
demi-droites {(N, 0) ; N > 0} et {(0, P) ; P > 0} sont des trajectoires, c’est-à-dire des courbes décrites
par des solutions de (9), et par suite, les trajectoires ne pouvant s’intersecter sans coı̈ncider, le
quart de plan {(N, P) ; N > 0 , P > 0} est aussi invariant par le flot. Une autre propriété mathématique
de (9) et de sa version adimensionnée (10), est leur conservativité, autrement dit l’existence d’une
intégrale première, qui est par définition une quantité constante le long des trajectoires. Cette
intégrale première se calcule aisément car elle est à variables séparées, c’est-à-dire somme d’une
fonction de N et d’une fonction de P. En effet, supposons que (N, P) soit une solution non triviale
de (10), c’est-à-dire pour laquelle ni N ni P ne s’annule, et telle que ni N ni P ne prennent la valeur
1. Alors
1 dN 1 dP
+ = 0,
N(1 − P) dt δP(1 − N) dt
d’où après multiplication par δ(1 − P)(1 − N),
d¡ ¢
δ(ln N − N) + ln P − P = 0 .
dt
Ceci montre précisément que la fonction H : (N, P) 7→ δ(ln N−N) + ln P−P est une intégrale première
de (10) (notez que les restrictions N 6= 1, P 6= 1, sont inutiles: elles sont apparues artificiellement
dans la manière de faire de la calcul). Cependant, aussi intéressante soit-elle mathématiquement, la
conservativité du système de Lotka–Volterra révèle ses limites du point de vue de la modélisation:
1Indication: poser T = 1/a , N
e = Nc/d , P
e = Pb/a .
5
les courbes de niveau de H forment en effet des courbes fermées dans le plan de phase {(N, P) ; N ≥
0 , P ≥ 0}, correspondant à des solutions périodiques assez peu réalistes. Pour plus de détails et des
modèles plus réalistes, on renvoie au chapitre 3 de [7].
Remarque 1. Le lecteur plus intéressé par la finance que par la dynamique des populations pourra
transposer ce qui précède en changeant nombre d’individus en masse d’argent, taux d’accroissement
naturel en taux d’intérêt, soldes migratoires en soldes budgétaires (différence entre recettes et
dépense), etc.
1.2. Mécanique. La modélisation en mécanique relève de lois assez bien établies. Celle que nous
utiliserons principalement en mécanique du point est la loi de Newton. Si un point matériel de
masse m est repéré dans l’espace Rd par sa position x(t ) ∈ Rd à l’instant t , sa vitesse est v(t ) =
x 0 (t ) ∈ Rd , son accélération est a(t ) = x 00 (t ) ∈ Rd , et s’il est soumis à une force F ∈ Rd , la loi de
Newton s’exprime par la célèbre relation F = ma , qui revient à une équation différentielle du
second ordre
d2 x
m = F.
dt 2
Lorsque la force F est constante, cette équation différentielle est triviale: elle s’intègre à vue pour
donner des solutions polynomiales d’ordre 2 en t . Ceci explique notamment les trajectoires pa-
raboliques des (petits) objets lancés en l’air, essentiellement soumis à la seule force de gravité
F = m~g , où ~
g est l’accélération de la pesanteur : à la surface de la Terre, ~ g est un vecteur dirigé
verticalement vers le bas (ces notions de « verticalité » et de « bas » dépendant évidemment du
point où l’on se situe, mais à une échelle spatiale très inférieure au rayon de la Terre, les variations
de ~
g sont négligeables), et de norme g ≈ 9.8m.s −2 . Lorsque la force F dépend de la position x (ce
qui est en principe le cas avec la force de gravité, comme on vient de le rappeler), la loi de Newton
devient une « vraie » équation différentielle du second ordre, en général non-linéaire:
d2 x F(x)
(11) 2
= .
dt m
1.2.1. Oscillateur harmonique. L’exemple le plus simple est celui d’un point matériel situé à l’ex-
trémité d’un ressort dont l’autre extrémité est fixée. Dans ce cas, (si le ressort est suffisamment
!x

E 0 x

Fig. 1. Ressort, modèle d’oscillateur harmonique.

rigide, ou contraint, pour avoir des mouvements rectilignes) d = 1. Si à chaque instant t , x(t ) est
l’abscisse de l’extrémité mobile dans un repère d’origine la position au repos (ressort ni étiré ni
comprimé), une expression couramment admise pour la force est alors
F(x) = − k x ,
où le paramètre k > 0 est la raideur du ressort (plus k est élevé plus la force de rappel est intense).
Pour une telle fonction F (linéaire), l’équation (11) devient, en posant pour simplifier κ := k/m ,
d2 x
(12) + κx = 0 ,
dt 2
appelée équation de l’oscillateur harmonique. On connaı̂t ses solutions explicitement. Plus préci-
sément, le problème de Cauchy pour (12) et une donnée initiale (x(0), x 0 (0)) = (x 0 , v 0 ) admet comme
solution unique la fonction
p p p
x : t 7→ x 0 cos( κt ) + (v 0 / κ) sin( κt ) .
6
La forme sinusoı̈dale de cette solution explique évidemment le terme d’oscillateur.

1.2.2. Pendule simple. Voyons maintenant un exemple où la force F dépend non-linéairement de
x . C’est encore un cas d’équation scalaire (d = 1) mais dans lequel x(t ) représente un angle, en
l’occurrence l’angle par rapport à la position verticale d’une tige fixée à son extrémité haute et
contrainte à se déplacer dans un plan. En fait, on suppose cette tige de masse négligeable sauf à
son extrémité libre M, où l’on met une masse m , dont on va voir que la valeur n’intervient pas
dans l’équation ! En revanche, la longueur ` de la tige intervient. L’équation du mouvement de ce
O

x ex

m
g

Fig. 2. Pendule simple (oscillant dans un plan).

dispositif, appelé pendule simple, s’écrit en effet

d2 x g
(13) + sin x = 0 .
dt 2 `
Attention, nous avons fait ici un petit glissement de notations par rapport au cadre général: les
« vraies » positions, vitesses et accélérations de M sont des vecteurs de R2 (et non les scalaires
x(t ), x 0 (t ), x 00 (t )); ceci devrait être clair dans le calcul qui suit.
On peut obtenir (13) par la loi de Newton (moyennant quelques hypothèses implicites, cf [1]),
bien que ce ne soit pas le plus facile. En effet, dans le bilan des forces exercées sur l’extrémité
libre du pendule, il faut tenir compte du fait que ce point est lié à la tige et non libre de se
déplacer dans l’espace. Par suite, la loi de Newton se réduit ici à l’égalité entre les projections
orthogonales à la tige de l’accélération du point et de l’accélération de la pesanteur. Voyons ce
que cela donne. Dans un repère orthonormé direct Oz y où O est le point de fixation du pendule,
z est la direction verticale dirigée vers le bas, l’extrémité libre M du pendule a pour coordonnées
(` cos x, ` sin x), de sorte que sa vitesse a pour coordonnées (−`x 0 sin x, `x 0 cos x), et son accélération
(−`x 00 sin x − `x 0 cos, `x 00 cos x − `x 0 sin x). Quant à l’accélération de la pesanteur, ses coordonnées
sont évidemment (g , 0). En égalant les produits scalaires des accélérations par le vecteur ~ e x de
coordonnées (− sin x, cos x) on obtient l’équation attendue
`x 00 = − g sin x .
On peut aussi obtenir cette équation de manière plus rapide grâce à la conservation de l’énergie
totale. En effet, l’énergie cinétique instantanée du point M est
1 −−→ 1
Ec := mkOM0 k2 = m `2 x 02 ,
2 2
et son énergie potentielle (nous reviendrons sur cette notion plus loin) est (à une constante près)
Ep := − m g z = − m ` cos x .
7
L’équation (13) se déduit donc de
d
(14) (Ec + Ep ) = 0
dt
(lorsque la masse m est non nulle!) aux instants où x 0 6= 0. Notez qu’« inversement » , (13) implique
(14) avec Ec et Ep définies comme ci-dessus. Autrement dit, l’énergie totale est une intégrale
première de l’équation du mouvement. Pour être plus précis, il est préférable d’écrire (13) comme
un système dans le plan de phase {(x, v = x 0 ) ∈ R2 }:
dx


 = v,
dt

(15)
 dv = − g sin x ,


dt `
et d’exprimer l’énergie totale en fonction de (x, v):
1
Ec + E p = m `2 v 2 − m ` cos x .
2
Cette fonction de (x, v) est une intégrale première de (15). C’est bien sûr un modèle idéalisé: en
pratique, le mouvement du pendule finira toujours par s’amortir. Un modèle plus réaliste devrait
tenir compte des frottements dans le bilan des forces, c’est-à-dire en fait d’une force F dépendant
aussi de la vitesse, ce qui détruirait la conservation de l’énergie (et en général aussi l’indépendance
de l’équation par rapport à la masse). Si l’on reste sur le modèle de pendule simple idéal, l’intégrale
première Ec +Ep permet de tracer facilement le portrait de phase de (15). On y observe notamment:
• des points d’équilibre: le point bas (x, v) = (0, 0), stable, et le point haut (x, v) = (π, 0), instable,
• des cycles correspondant à des oscillations périodiques autour du point d’équilibre bas (restant
strictement en dessous du point haut),
• des orbites « hétéroclines » (en fait homoclines si l’on identifie les points (π, 0) et (−π, 0) du plan
de phases), correspondant à un mouvement « atteignant » en temps infini le point d’équilibre
haut,
• des trajectoires correspondant à des tours complets effectués de façon périodique autour du
point d’équilibre haut.

1.2.3. Saut à ski. Intéressons nous à un skieur dévalant un pente en « ligne droite », c’est-à-dire en
suivant une courbe (pré-determinée) située dans un plan vertical (par exemple un tremplin de saut).
Supposons grossièrement que l’on peut assimiler ce skieur à un point matériel de masse m . Comme
pour le pendule, on va voir que si l’on néglige les frottements la masse m n’intervient pas dans
l’équation du mouvement, équation que l’on obtient très facilement en utilisant la conservation de
l’énergie totale. En effet, si l’on suppose pour fixer les idées que le skieur part avec une vitesse
nulle d’un point O, si l’on exprime les coordonnées de sa position M un repère orthonormé direct
Oz y où z est la direction verticale dirigée vers le bas, son énergie totale initiale est nulle (cette
valeur étant liée au choix arbitraire de la constante dans l’énergie potentielle) et à chaque instant
t elle vaut
1
Ec (t ) + Ep (t ) , avec Ec (t ) := m (y 0 (t )2 + z 0 (t )2 ) , Ep (t ) := − m g z(t ) .
2
Il faut de plus se souvenir que la trajectoire n’est pas libre mais donnée par une courbe (plane).
Disons qu’elle soit d’équation z = ζ(y) (notez que ce paramétrage par y autorise le tremplin à
remonter). La conservation de l’énergie totale fournit alors comme équation du mouvement:
(1 + ζ0 (y(t ))2 ) y 0 (t )2 − 2 g ζ(y(t )) = 0 ,
8
soit encore, en dérivant et en supposant que y 0 ne s’annule pas (il est assez rare qu’un sauteur
s’immobilise sur le tremplin!),
d2 y ³ dy ´2 ¶ ζ0 (y)
µ
00
(16) = g − ζ (y) .
dt 2 dt 1 + ζ0 (y)2
On a là un exemple d’équation différentielle du second ordre bien plus compliqué à première vue
que ceux rencontrés jusqu’à présent (notamment parce que la dérivée première de y intervient).
Cependant on sait déjà, étant donnée la façon dont on l’a obtenue, que cette équation ou plus
exactement le système équivalent dans le plan de phase

dy
 dt = w ,




dw ¢ ζ0 (y)
= g − ζ00 (y) w 2
 ¡

 .
dt 1 + ζ0 (y)2


admet
1
Ec + E p = m (1 + ζ0 (y)2 ) w 2 − m g ζ(y)
2
comme intégrale première. En fait, l’équation (16) est ce qu’on appelle une équation d’Euler–
Lagrange pour le lagrangien
1
L : (y, w) 7→ L(y, w) := m (1 + ζ0 (y)2 ) w 2 + m g ζ(y)
2
(notez le changement de signe par rapport à l’énergie totale), c’est-à-dire que (16) s’écrit de façon
abstraite à l’aide de L:
d ∂L ³ dy ´ ∂L ³ dy ´
µ ¶
(17) y, = y, .
dt ∂w dt ∂y dt
(La vérification est laissée au lecteur.) On a gardé la masse m dans le lagrangien par souci d’homo-
généité physique (un lagrangien ayant en général la dimension d’une énergie) mais elle n’intervient
pas ici dans l’équation du mouvement. Bien entendu, il s’agit comme pour le pendule simple d’un
modèle idéalisé: tout amateur de ski sait bien que plus on est lourd plus on va vite (si les sauteurs
sont maigres, c’est pour mieux « voler » après le tremplin; pour aller plus vite sur le tremplin il
vaudrait mieux qu’ils soient lourds comme les descendeurs).
Remarque 2. Un problème célèbre2 de calcul des variations est de déterminer la forme optimale
du tremplin, c’est-à-dire la fonction ζ minimisant le temps de parcours entre le point de départ O
et le point de décollage, ces deux points étant fixés. Or, sachant d’après ce qui précède que
s
dy 2 g ζ(y)
= ,
dt 1 + ζ0 (y)2
si a est l’abscisse du point de décollage, ce temps de parcours est donné par
s
Z a 1 + ζ0 (y)2
T= dy .
0 2 g ζ(y)
Une condition nécessaire de minimisation de T sous les conditions ζ(0) = 0, ζ(a) = h (ordonnée du
point de décollage) s’exprime aussi au moyen d’une équation d’Euler–Lagrange, en l’occurrence
celle associée au lagrangien s
1 + ξ2
P : (ζ, ξ) 7→ P(ζ, ξ) := ,
2g ζ
2si ce n’est le premier: il est connu sous le nom de problème de la brachistochrone.
9
qui s’écrit
d ∂P ³ dζ ´ ∂P ³ dζ ´
µ ¶
ζ, = ζ, .
dy ∂ξ dy ∂ζ dy
On montre que la courbe optimale est une cycloı̈de. On voit déjà, sans aucun calcul compliqué,
que cela ne peut pas être une droite. En effet, notons que
∂P ³ ´ ξ
ζ, ξ = P(ζ, ξ) ,
∂ξ 1 + ξ2
et donc, si l’on avait une solution linéaire ζ : y 7→ ξ0 y de l’équation d’Euler–Lagrange associée à
P, on aurait
ξ20 ∂P ∂P
(ξ y, ξ0 ) =
2 ∂ζ 0
(ξ0 y, ξ0 ) .
1 + ξ0 ∂ζ

1.2.4. Planètes. Supposons qu’une « petite » planète (ou satellite) de masse m soit dans le champ
de gravité d’une « grosse » planète (ou astre, comme le soleil) de masse M. Les adjectifs entre
guillemets sont surtout là pour justifier qu’on identifie la petite planète à un point. Si x ∈ R3
désigne la position de cette petite planète dans un référentiel lié à la grosse planète, la loi de la
gravitation universelle détermine la force exercée par la grosse planète sur la petite:
GMm x
F(x) = − .
kxk2 kxk
où G est la constante d’attraction universelle (G = 6, 672.10−11 m3 .s−2 .kg−1 ). Notez que l’expression
F = m~g utilisée plus haut pour un objet dans l’attraction de la Terre revient à négliger les variations
de x par rapport au rayon terrestre R et à poser g = GM/R2 où M est la masse de la Terre
(l’application numérique permettant de calculer la valeur de g est laissée au lecteur). Si l’on
conserve l’expression exacte de F, la loi de Newton (11) donne comme équation du mouvement
d2 x GM x
2
=− ,
dt kxk2 kxk
que l’on peut aussi écrire
d2 x 1 GMm
2
= ∇V(x) , V(x) := − .
dt m kxk
On dit que la force F dérive du potentiel V. C’est ainsi que l’on définit l’énergie potentielle: lorsque
la force F dérive d’un potentiel V, la loi de Newton s’écrit dans l’espace des phases {(x, v) ∈ Rd ×Rd }
dx


 = v,
dt


 dv = 1 ∇V(x) ,


dt m
et ce système d’équations différentielles admet l’énergie totale
1
m kvk2 + V(x)
2
comme intégrale première; le premier terme est l’énergie cinétique et le second est l’énergie po-
tentielle, définie à une constante près. On ne reconnaı̂t pas immédiatement dans V l’expression de
l’énergie potentielle utilisée plus haut à propos du pendule, Ep = m g h pour un point matériel de
masse m et d’altitude h . En fait, cette dernière résulte de l’approximation suivante:
GMm GMm GM
V=− ≈− + m 2 h, h ¿ R.
R+h R R
10
1.2.5. Toupie. Au paragraphe précédent nous avons assimilé les planètes à des points matériels.
Voyons maintenant un « vrai » problème de mécanique du solide, à savoir le mouvement d’une
toupie. On va obtenir les équations du mouvement comme équations d’Euler–Lagrange: il « suffit »
pour cela de savoir exprimer l’énergie cinétique et l’énergie potentielle en fonction des données
cinématiques (position et vitesse).
On considère plus précisément un solide de révolution (c’est-à-dire invariant par rotation autour
de son axe), homogène (c’est-à-dire de masse uniformément répartie), avec un point fixe apparte-
nant à l’axe de révolution (la pointe de la toupie). Ce solide S a un centre de gravité G, défini de
façon unique par Ñ
−−→
GM dλ(M) = 0 ,
S
où λ désigne la mesure de Lebesgue divisée par le volume de S (de sorte que λ(S) = 1). S’il est de
masse m , son énergie potentielle dûe à la gravité est (à une constante près, avec l’approximation
mentionnée ci-dessus)
Ñ
−−→ −−→
Ep = − m~
g · OM dλ(M) = − m ~
g · OG .
S
Pour décrire la position du solide S, il suffit de connaı̂tre la position de G, que l’on peut repérer
grâce à deux angles, appelés précession et nutation et que l’on note ici respectivement ψ et θ,
ainsi que l’angle de rotation propre 3 du solide autour de son axe. Avec les notations de la figure
→
−
3, θ est l’angle entre les vecteurs Z (vertical orienté vers le haut) et → −z (axe du solide), et ψ est
→
− →
− −z )/k→
− →
l’angle entre les vecteurs X (l’un des axes horizontaux) et → −
u := ( Z ∧ → Z ∧ −z k. Ces trois angles
sont appelés angles d’Euler . On voit que l’énergie potentielle s’exprime à l’aide de θ seulement: si

Z
z
v

y
G
# u
"
x
O Y

!
u

Fig. 3. Toupie et angles d’Euler.

−−→
` = kOGk, Ep = m g ` cos θ. Quant à l’énergie cinétique, elle est donnée de façon générale par
1
Ñ
Ec = m k~v (M)k2 dλ(M) .
S 2

Pour la calculer explicitement, il faut savoir qu’en tout point M du solide, la vitesse ~
v (M) s’exprime
→
−
au moyen du vecteur vitesse de rotation instantanée Ω par
→
− −−→
~
v (M) = Ω ∧ OM .
11
Ce vecteur vitesse de rotation instantanée s’exprime simplement à l’aide des dérivées par rapport
au temps des angles d’Euler:
→
−
u + ψ̇~
Ω = 3̇~z + θ̇ ~ Z.
Autrement dit,
1 →
− −−→ 1 →
− → −
Ñ
Ec = m k Ω ∧ OMk2 dλ(M) =: JO ( Ω , · Ω )
S 2 2
La forme quadratique JO est diagonalisable dans une base orthonormée appelée repère principal
d’inertie. Comme S est un solide de révolution, (~ x ,~ z) est un repère principal d’inertie, et les
y ,~
valeurs propres associées, appelées moments principaux d’inertie, sont des nombres A, B = A et C,
positifs ou nuls puisque JO est positive. De plus, si on note ~ v := ~
z ∧~
u , (~
u ,~ z) est aussi un repère
v ,~
principal d’inertie. Par suite,
1 →
− 1 →
− 1 →
−
Ec = A kP~u ( Ω )k2 + A kP~v ( Ω )k2 + C kP~z ( Ω )k2 ,
2 2 2
où P~u , P~v et P~z désignent les projections orthogonales sur les droites dirigées respectivement par
~
u~ v et ~ z . Or par définition de θ on a
→
−
Z = sin θ~
v + cos θ~
z,
donc →
− →
− →
−
P~u ( Ω ) = θ̇ ~
u, P~v ( Ω ) = ψ̇ sin θ ~
v, z,
P~z ( Ω ) = (3̇ + ψ̇ cos θ)~
et finalement,
1 1
A (θ̇2 + ψ̇2 sin2 θ) + C (3̇ + ψ̇ cos θ)2 .
Ec =
2 2
En résumé, le lagrangien L = Ec − Ep associé au mouvement de cette toupie s’écrit
1 1
L(θ, θ̇, ψ̇, 3̇) = A (θ̇2 + ψ̇2 sin2 θ) + C (3̇ + ψ̇ cos θ)2 − m g ` cos θ ,
2 2
et les équations d’Euler–Lagrange
d ∂L ¡ ∂L ¡ ¢ d ∂L ¡ d ∂L ¡
µ ¶ µ ¶ µ ¶
¢ ¢ ¢
θ, θ̇, ψ̇, 3̇ = θ, θ̇, ψ̇, 3̇ , θ, θ̇, ψ̇, 3̇ = 0 , θ, θ̇, ψ̇, 3̇ = 0
dt ∂θ̇ ∂θ dt ∂ψ̇ dt ∂3̇
donnent
³¡ ´
A ψ̇ cos θ − C (3̇ + ψ̇ cos θ) ψ̇ + m g ` sin θ ,
¢
A θ̈ =
A ψ̇ sin2 θ + C (3̇ + ψ̇ cos θ) cos θ = cste .
3̇ + ψ̇ cos θ = cste .
On a ainsi directement deux intégrales premières (les deux dernières lignes, qui sont dûes à ce que
ψ et 3 sont des variables cycliques, c’est-à-dire que L n’en dépend pas, bien qu’il dépende de ψ̇ et
3̇). Par ailleurs, l’énergie totale Ec + Ep est aussi une intégrale première, c’est-à-dire que
1 1
A (θ̇2 + ψ̇2 sin2 θ) + C (3̇ + ψ̇ cos θ)2 + m g ` cos θ = cste .
2 2
Les équations du mouvement se réduisent ainsi à une équation scalaire sur θ, de la forme:
β − γ cos θ 2 2mg `
µ ¶
2
θ̇ + + cos θ = α ,
sin θ A
où α, β, γ sont des constantes du mouvement (dépendant des conditions initiales, ainsi que de
A et C). En posant χ = cos θ, on se ramène à une équation de la forme χ̇2 = f (χ) où f est une
fonction polynômiale de degré 3 ayant trois zéros réels. Cette p équation, pourtant scalaire, admet
une solution périodique: ceci est dû au fait que la fonction χ 7→ f (χ) n’est
p pas Lipschitziennep aux
points où f s’annule, de sorte qu’il n’y a pas unicité des solutions de χ̇ = f (χ) (ni de χ̇ = − f (χ)).
12
On montre ainsi par une analyse purement qualitative l’existence de solutions périodiques dans
l’espace physique.
Remarque 3. La vitesse et l’accélération d’un point mobile P sur un solide lui-même en mouve-
ment (un humain sur la terre, une petite bête sur la toupie, etc.) sont données par les « lois de
composition » suivantes:
~
v (P) = ~
v e (P) + ~
v r (P) ,
~
a (P) = ~
a e (P) + ~
a r (P) + ~
a c (P) ,
où ~v e (P) et ~
a e (P) sont la vitesse et l’accélération d’ entraı̂nement, c’est-à-dire, à chaque instant
t , la vitesse et l’accélération du point lié au solide se trouvant au point P(t ) (ce point change au
cours du temps), ~ v r (P) et ~
a r (P) sont la vitesse et l’accélération d’ relative, c’est-à-dire la vitesse et
l’accélération dans un référentiel lié au solide, et enfin ~ a c (P) est l’accélération de Coriolis. Cette
dernière est définie comme:
→−
~
a c (P) = 2 Ω ∧ ~
v r (P) .

1.3. Électricité. Abandonnons maintenant le vaste domaine de la mécanique pour étudier les
circuits électriques, qui servent aussi à modéliser des phénomènes biologiques (battements du
coeur par exemple). L’état d’un circuit électrique composé de résistances, bobines et condensateurs,
peut être décrit par l’intensité I et la différence de potentiel U dans chacun de ces composants.
Les différentes lois de l’électricité conduisent à un système d’équations différentielles pour ces
quantités. Considérons par exemple un circuit fermé comprenant un composant de chaque sorte,
dans l’ordre résistance-bobine-condensateur (voir la figure 4), la bobine ayant pour inductance L,

L
C

Fig. 4. Circuit RLC.

le condensateur ayant pour capacité C, et le comportement de la résistance étant régi une relation
algébrique (loi d’Ohm généralisée) :
UR = F(IR ) .
Les équations régissant l’évolution de ce circuit sont :
dIL

 dt = UL ,
L




 C dUC = IC ,



dt
assorties des relations de compatibilité
UC = U L + UR , IR = IL = − IC .
13
(L’orientation du circuit étant arbitraire, les équations sont fort heureusement invariantes par
changement d’orientation.) En éliminant les autres inconnues, on se ramène à un système de deux
équations pour x = IL et y = UC par exemple :
dx

 L dt = y − F(x) ,




 C dy = − x .



dt
Remarque 4. Comme expliqué plus haut, on peut adimensionner ce système. Considérons une
intensité de référence I, un potentiel de référence U et un temps de référence T, et posons
x y t F(x)
xe := , ye := , e t := , F(e e x ) := .
I U T U
On voit alors qu’une fonction t 7→ (x(t ), y(t )) est solution du système de départ si et seulement si
la fonction e
t 7→ (e
x (e t )) est solution de
t ), ye(e
x
de TU ¡
 ¢

 = ye − F(e
e x) ,
 de t LI


 d ye = − IT xe ,



de t CU
où il reste deux paramètres sans dimension, p
à savoir (TU)/(LI)
p et (IT)/(CU). Notons qu’ils valent
tous deux 1 si l’on choisit en particulier T = LC et U/I = L/C, de sorte que le système adimen-
sionné n’a alors plus aucun paramètre.

2. Équations aux dérivées partielles

Dans ce chapitre, nous allons nous intéresser à des problèmes où la fonction inconnue dépend à
la fois du temps et de la position dans l’espace, problèmes que l’on peut modéliser par une équation
aux dérivées partielles d’évolution.
2.1. Équation de continuité. Comme point de départ nous allons considérer une équation « uni-
verselle », au sens où elle ne résulte d’aucune approximation et intervient dans différents contextes
physiques. Il s’agit d’une loi de conservation appelée équation de continuité. Cette loi fait interve-
nir une quantité k et son flux ~ q ∈ Rd pour un phénomène en dimension d d’espace. Par exemple:
• dans un modèle simple de trafic routier appelé modèle de Lighthill-Whitham-Richards, d = 1,
k(x, t ) représente la concentration de véhicules, c’est-à-dire le nombre de véhicules par unité
de longueur au point x ∈ R de la route (vue de suffisamment loin) à l’instant t et q(x, t ) est le
débit de véhicules par unité de temps,
• en mécanique des fluides, d = 3, k(x, t ) est la masse volumique du fluide au point x ∈ R3 à
l’instant t et ~ q (x, t ) = (q 1 , q 2 , q 3 )(x, t ) ∈ R3 avec q j (x, t ) le débit massique (masse par unité de
temps et unité d’aire); si l’on s’intéresse au transport de polluants par exemple, k(x, t ) peut
être la concentration volumique d’une espèce chimique (nombre de moles par unité de volume)
et ~q son débit (nombre de moles par unité de temps et unité d’aire)
• en électrodynamique, d = 1, 2 ou 3, k(x, t ) est la densité de charge électrique qui s’exprime en
Coulomb (C) par unité de longueur, d’aire ou de volume, selon que d = 1, 2 ou 3, et ~ q (x, t ) est
la densité de courant électrique, dont les composantes s’expriment en Ampère (A) si d = 1, en
Ampère par unité de longueur si d = 2, en Ampère par unité d’aire si d = 3, en supposant dans
tous les cas que l’unité de temps est la seconde (sachant qu’un Ampère vaut exactement un
Coulomb par seconde); on peut aussi utiliser une autre unité de densité de charge électrique,
à savoir le Faraday, qui est la charge électrique d’une mole de charges élémentaires (1F ≈
14
 96485C), auquel cas il faut remplacer les Ampère dans la densité de courant électrique par des
Faraday par unité de temps; on peut également définir k(x, t ) comme le nombre de charges
élémentaires (sachant qu’une charge élémentaire vaut environ 1, 602×10−19 C) et ~
q (x, t ) comme
le débit de charges élémentaires.
• en thermique, d = 1, 2 ou 3, k(x, t ) est la température du matériau au point x ∈ Rd à l’instant
t et q(x, t ) est le flux de chaleur .
Quel que soit le contexte physique, l’équation de continuité est valable en l’absence de source
extérieure et s’écrit:
d
∂t k + div~ div~ ∂x j q j .
X
(18) q = 0, q :=
j =1

Voyons comment on l’obtient en dimension 1 d’espace pour simplifier. Pour fixer les idées, disons
que k(x, t ) est un nombre de particules (atomes, molécules, charges électriques, véhicules, etc.) par
unité de longueur et que q(x, t ) est le débit de ces particules, c’est-à-dire un nombre de particules
par unité de temps. L’équation de continuité
(19) ∂t k + ∂x q = 0
se déduit du bilan de conservation suivant. Soient a , b , s , T tels que a ≤ b , s ≤ T. Le nombre de
particules situées dans l’intervalle [a, b] à l’instant T est
Z b
k(x, T) dx .
a
En l’absence de sources extérieures (qui viendraient ajouter ou retirer des particules), il est égal
au nombre particules situées dans l’intervalle [a, b] à l’instant s , à savoir
Z b
k(x, s) dx
a
retranché du nombre de particules sortant par le point b entre s et T, à savoir
Z T
q(b, t ) dt ,
s
auquel on ajoute le nombre de particules entant par le point a entre s et T, à savoir
Z T
q(a, t ) dt .
s
Autrement dit, on a le bilan de conservation:
Z b Z b Z T Z T
k(x, T) dx = k(x, s) dx − q(b, t ) dt + q(a, t ) dt ,
a a s s
que l’on peut aussi écrire I
k dx − q dt = 0 ,
∂R
où R = [a, b] × [s, T]. Or d’après la formule de Green–Riemann
I Z
k dx − q dt = − (∂t k + ∂x q) dx dt .
∂R R
Donc Z
(∂t k + ∂x q) dx dt = 0 , quel que soit le rectangle R .
R
Ceci montre que l’équation de continuité (19) doit être satisfaite presque partout.
L’obtention de l’équation de continuité (18) en dimension d d’espace est tout à fait analogue,
en remplaçant le rectangle R par un pavé de Rd × R.
15
2.2. Convection versus diffusion. Cependant l’équation (18) n’est pas encore un modèle ma-
thématique à proprement parler, car elle porte en fait sur d + 1 inconnues: les fonctions k et q j ,
j ∈ {1, . . . , d }. C’est le lien entre k et ~
q , appelé loi de fermeture, qui définit un modèle. On distingue
en fait deux catégories de tels modèles: ceux décrivant un phénomène de convection, correspondant
au transport de particules suivant un champ de vitesse ~ v , auquel cas
(20) ~
q = k~
v,
et ceux décrivant un phénomène de diffusion, dans lesquels
(21) ~
q = − λ∇k ,
où λ > 0 est le coefficient de diffusion. Cette dernière relation signifie que la quantité k a tendance
à s’« homogénéiser », son flux étant orienté en direction des zones où k est moins élevée. Lorsque
k est une température, (21) est la loi de Fourier et λ la diffusivité thermique du matériau. Lorsque
k est la concentration d’une espèce chimique, (21) est la loi de Fick , inspirée de la loi de Fourier.
Quant à la relation (20), elle n’est pas encore une loi de fermeture: on a simplement troqué
l’inconnue ~ q pour ~ v . L’équation qui en résulte,
∂t k + div(k ~
v) = 0,
peut être couplée à
– la donnée pure et simple de ~ v , par exemple lorsque k est la concentration d’un polluant
transporté par un fluide de vitesse prédéterminée ~
v,
~
– une loi empirique exprimant v en fonction de k , ce qui est le cas par exemple dans le modèle
de Lighthill-Whitham-Richards pour le trafic routier: v = v max (1 − k/kmax ),
– une autre loi de conservation, comme celle de la quantité de mouvement
v ) + div(k ~
∂t (k~ v ⊗~
v + p I) = 0 ,
où p est la pression du fluide.
Dans le cas d’un transport (20) par ~
v fonction de k , l’équation de continuité (18) s’écrit sous
forme d’une loi de conservation non-linéaire
(22) ∂t k + div(k ~
v (k)) = 0 ,
tandis que lorsque ~
q vérifie (21), (18) devient ce qu’on appelle l’équation de la chaleur
(23) ∂t k = λ∆k .
Il se trouve que (22) et (23) ont des propriétés mathématiques très différentes (hormis le fait que
l’une soit linéaire et l’autre non): la première propage l’information à vitesse finie (en ce sens
elle relève des EDPs dites hyperboliques) et pas l’équation de la chaleur; cette dernière régularise
instantanément la condition initiale (si k est solution de (23) avec k(0, ·) ∈ L∞ (Rd ) alors pour tout
t > 0, k(0, ·) ∈ C ∞ (Rd ): ceci se montre facilement grâce à la transformation de Fourier); au contraire,
une solution de (22) initialement de classe C ∞ devient « génériquement » discontinue en temps
fini lorsque l’équation est vraiment non-linéaire, c’est-à-dire lorsque ~ v (k) dépend effectivement
de k (et lorsque ~ v (k) est constant, la régularité de la donnée initiale est propagée sans aucune
amélioration). Nous allons préciser les propriétés de (22) sur un exemple.

2.3. Trafic routier. Considérons le modèle de Lighthill-Whitham-Richards:

∂t k + ∂x (k v max (1 − k/k max )) = 0 .
Pour simplifier les calculs à venir, commençons par l’adimensionnaliser en posant
t := t /T , xe := t /X , ke := k/k max .
e
16
Si l’on choisit de plus X et T tels que v max = X/T, il ne reste plus aucun paramètre dans le modèle
adimensionné (où l’on omet bien sûr les tildes):
(24) ∂t k + ∂x (k (1 − k)) = 0 .

(Attention, c’est assez pratique de ne pas avoir de paramètre, mais on perd de vue la nature
physique des variables, si bien que les erreurs de calcul éventuelles sont plus difficiles à repérer.)
On se fixe comme objectif de pouvoir calculer « à la main » la solution de (24) correspondant à
l’écoulement d’une file d’attente, c’est-à-dire avec une donnée initiale de la forme
½
0 si x < −` ou x > 0,
(25) k 0 (x) =
1 si − ` ≤ x ≤ 0,

où ` est la longueur de la file d’attente. (Rappelons que 1 est la concentration adimensionnée
maximale.) Nous aurons besoin pour cela de trois ingrédients:
• la méthode des caractéristiques, permettant de calculer la solution dans les zones où elle est
régulière,
• les ondes de détente, permettant de calculer la solution autour à partir d’un point où la donnée
initiale a une discontinuité décroissante (comme c’est le cas pour k0 donnée par (25) en x = 0),
• les ondes de choc, permettant de définir des solutions discontinues.

2.3.1. Méthode des caractéristiques. Une solution régulière de (24) doit vérifier l’équation dite
quasi-linéaire
(26) ∂t k + (1 − 2k) ∂x k = 0 .

En guise d’introduction à la méthode des caractéristiques, commençons par considérer l’équation

linéarisée autour d’une constante k . C’est l’équation obtenue en première approximation pour
une solution de la forme k = k + ε kl + o(ε). En ne retenant que les termes d’ordre ε on arrive à
l’équation de transport linéaire sur kl :
∂t k l + (1 − 2k) ∂x k l = 0 .

Celle-ci se résout facilement car elle signifie

d
k l (y + (1 − 2k)t , t ) = 0 , ∀y ∈ R ,
dt
d’où
k l (y + (1 − 2k)t , t ) = k 0 (y) , ∀y ∈ R .
Autrement dit, kl est constante sur les droites d’équation x = y + (1 − 2k)t , appelées caractéris-
tiques, ou encore, les valeurs de la donnée initiale k0 se propagent à la vitesse 1 − 2k (en variables
physiques, cette vitesse vaut v max (1 − 2k/kmax )). Il est tout à fait remarquable que cette propriété
(de propagation sur des droites de la donnée initiale) s’étende aux solutions (régulières) de l’équa-
tion non-linéaire (26). Supposons en effet que (x, t ) 7→ k(x, t ) soit une telle solution. Soit alors une
courbe, appelée courbe caractéristique, paramétrée par t 7→ χ(t ), où χ est solution de l’équation
différentielle ordinaire
dχ
= (1 − 2k(χ, t )) .
dt
Alors par dérivation de fonctions composées,
d
k(χ(t ), t ) = ∂t (χ(t ), t ) + (1 − 2k(χ(t ), t )) ∂x k(χ(t ), t ) = 0 .
dt
17
Ceci montre que k(χ(t ), t ) est constant et, par contrecoup, que la courbe caractéristique est une
droite. Autrement dit, une solution régulière du problème de Cauchy
∂t k + (1 − 2k) ∂x k = 0 ,
½

k(x, 0) = k 0 (x)
est telle que
k(y + (1 − 2k 0 (y))t , t ) = k 0 (y) , ∀y ∈ R .
On remarque que si k 0 ∈ Cb1 (R) et si t ∈ [0, T] avec 2T maxR k 00
< 1, la fonction y 7→ y + (1 − 2k 0 (y))t
est une bijection strictement croissante de R dans R. Donc pour tout x ∈ R il existe un unique
y ∈ R (le « pied » de la caractéristique passant par x ) tel que x = y +(1−2k 0 (y))t , ce qui détermine
sans ambiguı̈té k(x, t ) = k0 (y). Notons que la condition 2T maxR k00 < 1 est satisfaite pour tout
T ≥ 0 lorsque k 0 est décroissante: dans ce cas, la méthode que l’on vient de décrire fournit une
une unique solution régulière globale. En revanche, s’il existe un point où k00 prend une valeur
strictement positive, le temps
1
T∗ :=
2 maxR k 00
est un temps d’explosion, à partir duquel la méthode des caractéristiques n’est plus valable: on
vérifie en calculant ∂x k par le théorème des fonctions implicites que
lim ∂x k(y + (1 − 2k 0 (y))t , t ) = +∞
t %T∗

si y est un point où k00 atteint son maximum.

En pratique, on peut toujours appliquer la méthode des caractéristiques localement lorsque k0
est C 1 par morceaux.

2.3.2. Ondes de détente. Lorsque k0 a une discontinuité décroissante en un point y 0 , la méthode

des caractéristiques ne permet pas de calculer complètement la solution car elle ne fournit aucune
valeur dans le secteur compris entre les (demi-)droites d’équation x = y 0 + (1 − 2k0 (y 0 −))t et x =
y 0 + (1 − 2k 0 (y 0 +))t . On va combler ce « vide » par une onde de détente. L’idée pour la construire
est basée sur l’observation suivante. On suppose pour simplifier que y 0 = 0 (ce qui est le cas pour
la donnée initiale (25)), et on note k± = k0 (0±), en supposant k− > k+ . Le problème de Cauchy
 ∂t k + (1 −
½ 2k) ∂x k = 0 ,


(27) k − si x < 0,
 k(x, 0) =
k + si x > 0,
(que l’on appelle problème de Riemann), est invariant par la dilatation (x, t ) 7→ (αx, αt ) quel que
soit α > 0. Par suite, à supposer qu’il admette une solution unique (ce qui est le cas), celle-
ci vérifie nécessairement k(x, t ) = k(αx, αt ) quel que soit α > 0. En particulier pour α = 1/t , on
voit que k(x, t ) = k(x/t , 1). Autrement dit, cette solution ne dépend que de x/t : on dit qu’elle
est auto-similaire. Sachant cela, on peut facilement la calculer. Cherchons en effet k sous la
forme k(x, t ) = κ(x/t ) pour t > 0, où la fonction ξ 7→ κ(ξ) est dérivable. Par dérivation de fonctions
composées, on obtient
(−x/t + 1 − 2κ(x/t ))κ0 (x/t ) = 0 , ∀x ∈ R , ∀t > 0 .
On en déduit
1 − 2κ(ξ) = ξ , ∀ξ ∈ R ; κ0 (ξ) 6= 0 .
Par ailleurs, pour que k satisfasse la condition initiale on doit avoir
lim κ(ξ) = k ± .
ξ→±∞
18
La solution à ce problème en κ est donnée par:

 k− si ξ < 1 − 2k− ,
κ(ξ) = (1 − ξ)/2 si 1 − 2k− ≤ ξ ≤ 1 − 2k+ ,
k+ si x > 1 − 2k+ ,


ce qui donne, pour tout t > 0,


 k− si x < (1 − 2k− )t ,
k(x, t ) = (t − x)/2t si (1 − 2k− )t ≤ x ≤ (1 − 2k+ )t ,
k+ si x > (1 − 2k+ )t .


Cette solution est précisément ce que l’on appelle une onde de détente. (Le fait qu’elle soit affine
dans le secteur vient de la forme particulière de la loi q = k(1 − k), polynomiale du second degré.
En général, une onde de détente est déterminée par une équation implicite.)
2.3.3. Ondes de choc. Lorsque k− < k+ , la résolution du problème de Cauchy (27) est très dif-
férente de ce qui précède. On peut commencer par se convaincre facilement que la méthode des
caractéristiques est vouée à l’échec: en effet, les droites caractéristiques issues des points à gauche
de la discontinuité initiale ont une vitesse 1 − 2k− strictement supérieure à la vitesse 1 − 2k+ des
droites caractéristiques issues des points à droite, ce qui fait que les premières s’intersectent avec les
secondes. Pour pallier ce problème, on va chercher une solution auto-similaire discontinue appelée
onde de choc. Autrement dit, on cherche une vitesse σ telle que la fonction
½
k − si x < σt ,
(28) (x, t ) 7→ k(x, t ) =
k + si x > σt ,
soit une solution de (24). Attention, pour définir les ondes de choc on ne peut pas utiliser la forme
quasilinéaire (26) de l’équation, car elle n’a tout simplement pas de sens pour une fonction k
discontinue en espace (cela reviendrait à multiplier une fonction de Heaviside par une masse de
Dirac). En revanche, on peut donner un sens (faible) à (24) pour toute fonction k admettant une
trace intégrable sur les droites verticales et horizontales du plan (x, t ): on l’a même déjà rencontrée
plus haut, il s’agit de la formule intégrale
I
k dx − q dt = 0 , ∀R = [a, b] × [s, T] .
∂R
Si l’on suppose pour fixer les idées que la vitesse cherchée σ est positive, choisissons a = 0, b = σ,
s = 0, T = 1. D’après la formule intégrale ci-dessus, pour que la fonction discontinue (28) soit
solution (faible) de (24) il faut que
(k + − k − ) σ = q + − q − , q ± := k ± (1 − k ± ) ,
c’est-à-dire encore
σ = 1 − k− − k+ .
Inversement, on vérifie que cette relation, appelée condition de Rankine–Hugoniot, suffit pour que
(28) soit solution (faible) de (24). (Noter que lorsque k+ tend vers k− , la vitesse σ du choc tend
vers la vitesse caractéristique 1−2k− .) Plus généralement, on montre qu’une fonction de classe C 1
par morceaux, discontinue à travers une seule courbe Γ, appelée courbe de choc paramétrée par
t 7→ γ(t ) est solution (faible) de (24) si et seulement si elle vérifie
• l’équation (24) (ou (26) de façon équivalente) en tout point (x, t ) tel que x 6= γ(t ) et
• la condition de Rankine–Hugoniot à travers Γ, à savoir:
dγ
(29) = 1 − k(γ(t )−, t ) − k(γ(t )+, t ) .
dt
Attention, contrairement aux courbes caractéristiques, les courbes de choc n’ont aucune raison
d’être des droites en général.
19
2.3.4. Résolution du problème de la file d’attente. On considère le problème de Cauchy (24)-(25),
que l’on propose de résoudre explicitement à titre d’exercice. Ce problème modélise par exemple
l’écoulement d’une file d’attente arrêtée à un feu de circulation situé au point x = 0, qui passe au
vert à l’instant t = 0.
(1) Montrer, notamment à l’aide de la méthode des caractéristiques, qu’il existe un temps T0 tel
que, pour tout t ∈]0, T0 [, la solution à l’instant t est composée d’une zone où la concentration
est nulle, séparée par un choc stationnaire (c’est-à-dire de vitesse nulle) d’une zone où est
la concentration est maximale, elle-même séparée par une onde de déténte d’une autre zone
où la concentration est nulle.
(2) Calculer avec ce modèle le temps auquel le dernier véhicule de la file d’attente démarre.
(3) Calculer la solution au delà du temps T0 , en montrant notamment grâce à la condition
de Rankine–Hugoniot que la courbe de choc est incurvée par l’onde de détente. (Cette
condition fournit une équation différentielle ordinaire qu’il faudra résoudre pour calculer
la position du choc à chaque instant.)
(4) Calculer le temps T1 auquel le dernier véhicule de la file d’attente franchit le feu.
(5) Exprimer T0 et T1 dans les variables physiques. Faire une application numérique avec par
exemple v max = 50 km/h et ` = 100m.
2.4. Mélange des genres. On peut bien sûr considérer un modèle prenant en compte à la fois
de la convection et de la diffusion en posant
(30) ~
q = k~
v − λ ∇k
dans (18), auquel cas on obtient une équation dite de convection-diffusion:
∂t k + div(k ~
v ) = λ∆k .
Ceci correspond par exemple à la propagation de la chaleur à la fois par diffusion et par convection
dans un fluide en mouvement. On peut également prendre en compte un terme source, dépendant
éventuellement de k , auquel cas on parle plutôt de terme de réaction (chimique par exemple). Une
équation de la forme
∂t k = λ∆k + f (k)
est appelée équation de réaction-diffusion. De nombreux modèles en biologie s’écrivent à l’aide
d’équation de réaction-diffusion (voir [7]).

2.5. Équation des ondes. Nous allons maintenant nous intéresser à une autre équation qui, si
elle n’est pas universelle, apparaı̂t comme un modèle satisfaisant dans divers contextes physiques.
Elle partage avec les équations de transport par convection la propriété de propager l’information
à vitesse finie. Il en est ainsi de l’équation des ondes en acoustique pour la propagation du son (à
une vitesse dépendant des variables thermodynamiques mais disons de l’ordre de 340 m/s dans
l’air ambiant) et en électromagnétisme pour la propagation de la lumière (à la vitesse vertigineuse
de 300 000 km/s). A l’origine, elle a été obtenue par d’Alembert au XVIIIème siècle pour les cordes
vibrantes. Nous verrons aussi deux autres exemples.
2.5.1. Cordes vibrantes. Dans le mouvement d’une corde en tension, les paramètres physiques mis
en jeu sont la densité linéaire ρ0 (masse par unité de longueur) de la corde, reliée à la masse
volumique µ0 et à la section σ0 par
ρ0 = σ0 µ0 ,
et T0 la tension initiale de la corde, nombre strictement positif ayant la dimension physique d’une
force.
20
 Notons ~ u (x, t ) ∈ R3 le déplacement transversal de la corde à l’instant t , par rapport à une position
de référence x~e1 ∈ R3 , x ∈ R, et supposons le déplacement longitudinal négligeable. Autrement dit,
le point situé en x~e1 dans la position de référence se retrouve à l’instant t en w ~ (x, t ) = x~e1 + ~
u (x, t )
avec ~u (x, t ) ⊥ ~e1 . Si l’on isole mentalement un morceau de corde situé au repos dans un intervalle
[a, b], chacune de ses extrémités est soumise à une force colinéaire à ∂x w =~e1 + ∂x u . C’est ainsi
que l’on définit la la tension T(x, t ) > 0 de la corde en w ~ (x, t ), de sorte que la résultante des forces
appliquées au morceau de corde est
~F = T(b, t ) ∂x w
~ (b, t ) − T(a, t ) ∂x w
~ (a, t ) .
Rb
Or la somme des masses × accélérations le long du morceau de corde est a ~ (x, t ) dx . En
ρ0 ∂2t w
intégrant la loi de Newton le long de ce morceau on obtient donc
Z b
~ (b, t ) − T(a, t ) ∂x w
T(b, t ) ∂x w ~ (a, t ) = ρ0 ∂2t w
~ (x, t ) dx .
a

En projetant cette égalité vectorielle sur la direction ~e1 on en déduit

T(b, t ) − T(a, t ) = 0 ,
quels que soient a et b , ce qui montre que T(x, t ) ne dépend en fait pas de x . Pour la suite, on
suppose que la tension T ne dépend pas non plus de t , c’est-à-dire que T(x, t ) = T0 quels que soient
x et t . La projection sur ~e⊥
1 de l’égalité vectorielle ci-dessus donne alors
Z b
u (b, t ) − ∂x ~
T0 (∂x ~ u (a, t )) = ρ0 ∂2t ~
u (x, t ) dx ,
a
ou encore
Z b
u (x, t ) − ρ0 ∂2t ~
(T0 ∂2x ~ u (x, t )) dx .
a
Ceci étant vrai quels que soient a et b , on en déduit que ~
u vérifie l’équation des ondes mono-
dimensionnelle
(31) ∂2t u − c 2 ∂2x u = 0 ,
avec la vitesse s s
T0 T0
c := = .
ρ0 σ0 µ0
Remarque 5. On voit ici l’importance de la positivité de T0 (si T0 était négatif l’équation
T0 ∂2x ~
u − ρ0 ∂2t ~
u = 0 serait de nature très différente de l’équation des ondes: le problème de Cauchy
correspondant serait mal posé). Ceci n’est pas surprenant car une corde qui n’est pas en tension
s’affaisse et ne peut pas pas vibrer.
Remarque 6. On verra plus loin que si la corde est de longueur L, la vitesse c est reliée à la
fréquence temporelle ω des ondes « élémentaires » (aussi appelées harmoniques) par c = 2Lω.
Autrement dit,
s
T0
ω= 2
.
4L σ0 µ0
Du point de vue acoustique, plus ω est grand plus le son produit est perçu comme aigü. On voit
ainsi sur l’expression de ω qu’une corde produira un son d’autant plus aigü que la tension sera
grande, ou que sa section sera petite, ou encore que sa longueur sera petite. Ceci est conforme à
l’expérience avec les instruments de musique à corde tels que violon, guitare ou piano.
21
2.5.2. Barres élastiques. Contrairement à ce qui se passe pour une corde, les déplacements dans
une barre élastique rigide sont essentiellement longitudinaux, c’est-à-dire qu’un point situé en x~e1
dans la position de référence se retrouve après compression ou étirement en w ~ (x, t ) = (x + u(x, t ))~e1 .
On définit encore T(x, t ) la tension de la barre en w ~ (x, t ), mais cette fois elle n’a pas de signe
défini (la barre pouvant être indifféremment en compression ou en étirement). Une loi de l’élasticité
affirme que pour faire varier de δl un morceau de longueur l 0 il faut une variation de tension δT
proportionnelle à δl /l 0 . Quantitativement, on définit E0 le module d’Young du matériau tel que
δl
δT = E0 σ0 ,
l0

où σ0 désigne à nouveau la section (de la barre). Par définition, E0 est un nombre positif homogène
à une pression (c’est-à-dire une force par unité d’aire). En appliquant cette loi à un morceau
[x, x + δx], qui devient [x + u(x, t ), x + δx + u(x + δx, t )] à l’instant t , on obtient

u(x + δx, t ) − u(x, t )

T(x, t ) − T0 (x) = E0 σ0 ,
δx
d’où à la limite lorsque δx tend vers 0 :

T(x, t ) = T0 (x) + E0 σ0 ∂x u .

Pour la suite, on supposera T0 indépendant de x . D’après la loi de Newton appliquée au morceau

de barre compris entre les abscisses a et b , on a
Z b
T(b, t ) − T(a, t ) = ρ0 ∂2t u(x, t ) dx ,
a

ce que l’on peut encore écrire grâce à l’expression de T en fonction de ∂x u ,

Z b
(E0 σ0 ∂2x u − ρ0 ∂2t u(x, t )) dx = 0 .
a

Ceci étant vrai quels que soient a et b , on en déduit que u satisfait l’équation des ondes (31) avec
la vitesse
s
E0
c= .
µ0
Si l’on s’intéresse à la densité ρ le long de la barre, on voit assez facilement qu’elle est donnée par

ρ(x, t ) = ρ0 ( 1 − ∂x u ) .

En effet, pour chaque morceau de longueur initiale l 0 on a

δρ δl
ρ l = ρ0 l 0 d’où + = 0.
ρ0 l 0
En appliquant cette relation au morceau [x, x + δx] et en faisant tendre δx vers 0, on en déduit
ρ − ρ0
+ ∂x u = 0 .
ρ0
Par suite, en supposant la densité initiale ρ0 homogène, on voit par dérivation que ρ satisfait la
même équation des ondes que u (et T).
22
2.5.3. Tuyaux sonores. Un peu d’intuition physique (en pensant à nos tympans par exemple)
montre que la tension T exercée par un fluide est reliée à la pression p par p = − T/σ0 . D’où,
δT δp 1 δv
= =− ,
T0 p0 χ0 v 0
où χ0 est le coefficient de compressibilité (sans dimension), et v le volume. Or dans un tube de
section constante,
δv δl
= .
v0 l0
Donc on a une loi analogue à celle de l’élasticité, avec
p0
E0 = .
χ0
En particulier, pour un gaz parfait adiabatique,

p v γ = cte ,

d’où χ0 = 1/γ et E0 = γ p 0 . On trouve comme vitesse de propagation

γ p0
r
c= ,
µ0
ce qui l’expression bien connue de la vitesse du son.
Application numérique. Dans l’air, assimilé à un gaz di-atomique, on a approximativement
γ = 7/5 (ceci s’obtient en raisonnant sur le nombre n de degrés de liberté des molécules; de façon
générale, γ = (5 + n)/(3 + n)). La loi des gaz parfaits
µRT
p= , R = 8, 3144 J.K−1 .mol−1 , M = 28, 8.10−3 [Link]−1 ,
M
permet de calculer
s
RT
c= γ ≈ 332 m.s −1
M
à une température de 273 K.

2.5.4. Équation des ondes multi-dimensionnelle. L’analogue de (31) en dimension d d’espace

s’écrit

(32) ∂2t u − c 2 ∆u = 0 ,

Par exemple pour d = 2, cette équation peut modéliser la vibration d’une peau (au lieu d’une
corde, ce qui revient pour un musicien à délaisser le violon, la guitare ou le piano pour une
timbale), la vibration d’une plaque (au lieu d’une barre: on troque le triangle pour une cymbale),
la propagation de vagues à la surface d’un lac. Quant à l’équation des ondes en dimension d = 3,
c’est l’équation fondamentale de l’acoustique (avec u la vitesse (macroscopique) du fluide ambiant
ou la perturbation de pression) et de l’électromagnétisme (avec u le champ électromagnétique).
L’équation des ondes (et ses variantes, non-linéaires notamment) peut faire l’objet d’un cours
entier de master. Nous allons nous contenter de quelques propriétés de base, de formules de réso-
lution explicite et d’éléments pour son approximation numérique.
23
2.5.5. Propriétés de base de l’équation des ondes. L’équation (32) est conservative, au sens où
l’énergie totale
1 1
Z Z
2
E := E c + E p , avec E c := (∂t u) dx , E p := c 2 k∇uk2 dx ,
2 R 2 R

est conservée au cours du temps. En effet, supposons que u ∈ C ([0, T]; L2 (Rd )) ∩ C ([0, T]; H2 (Rd ))
2

soit une solution. Alors, en multipliant l’équation par ∂t u , on obtient après intégration par parties
du second morceau :
d 1 d 1
Z Z
(∂t u)2 dx + c 2 k∇uk2 dx = 0 .
dt 2 R dt 2 R
Par ailleurs, l’équation (32) propage l’information à vitesse finie c , comme on va le voir sur les
formules explicites.
2.5.6. Formules de résolution explicite.
Forme générale des solutions en dimension 1. Plusieurs approches sont possibles. Nous
allons en voir trois. La première consiste à remarquer la décomposition de l’opérateur des ondes
(appelé aussi d’Alembertien) en
∂2t − c 2 ∂2x = ( ∂t − c ∂x ) ( ∂t + c ∂x ) .
Par suite, pour résoudre
∂2t u − c 2 ∂2x u = 0
il suffit de résoudre successivement
( ∂t − c ∂x ) v = 0 ,
( ∂t + c ∂x ) u = v .
Or la première équation est une simple équation de transport, dont la solution ne dépend que de
(x + c t ). Pour s’en convaincre, il suffit de vérifier que
d
v(y − c t , t ) = ∂t v − c ∂x v = 0 ,
dt
et donc v(y − c t , t ) = v(y, 0) quel que soit y . En renversant les notations, on a v(x, t ) = v(x + c t , 0)
quel que soit x . Notons plus simplement h(y) = v(y, 0). Il faut ensuite résoudre l’équation de
transport avec terme source :
( ∂t + c ∂x ) u(x, t ) = h(x + c t ) .
La solution générale de l’équation homogène est u(x, t ) = g (x − c t ) (par le même argument que
ci-dessus). De plus, en notant f une primitive de h/2c , on a une solution particulière « évidente »
u(x, t ) = f (x + c t ). Finalement, par linéarité, la solution recherchée est de la forme :
u(x, t ) = f (x + c t ) + g (x − c t ) .
Une seconde méthode consiste à faire le changement de variables
(x, t ) 7→ (y, z) := ( x + c t , x − c t ) .
On a
∂x = ∂ y + ∂z , ∂t = c ( ∂ y − ∂z ) ,
d’où
∂2t − c 2 ∂2x = c 2 ( ∂ y − ∂z )2 − ( ∂ y + ∂z )2 = − 2 c 2 ∂2y z .
¡ ¢

On doit donc résoudre le problème

∂2y z ue = 0
dans les nouvelles variables, dont la solution est évidemment de la forme
u(y,
e z) = f (y) + g (z) .
24
Problème de Cauchy sur R. Soit à résoudre le problème de Cauchy :
∂t u − c 2 ∂2x u = 0
 2





u(x, 0) = φ(x) ,




∂t u(x, 0) = ψ(x) .


Il s’agit d’exprimer les fonctions f et g apparaissant dans la forme générale de la solution en

fonction de φ et ψ. Les deux équations à satisfaire sont
½
φ(x) = f (x) + g (x) ,
ψ(x) = c f 0 (x) − c g 0 (x) .
En dérivant la première, on en déduit
2c f 0 (x) = c φ0 (x) + ψ(x) ,
½

2c g 0 (x) = c φ0 (x) − ψ(x) ,

d’où Z x
 1 1
 f (x) = φ(x) +
 ψ(y) dy + cste ,
2 2c Z0x
1 1
 g (x) = φ(x) − ψ(y) dy + cste’ ,

2 2c 0
Finalement, en revenant à l’équation non dérivée φ(x) = f (x) + g (x), on en déduit
Z x +c t
1 1
u(x, t ) = f (x + c t ) + g (x − c t ) = ( φ(x + c t ) + φ(x − c t ) ) + ψ(y) dy .
2 2c x − c t
Cette formule est appelée formule de d’Alembert. Réciproquement, si φ est au moins deux fois
dérivable, et ψ au moins une fois dérivable, la formule de D’Alembert fournit la solution cherchée.
(En fait, l’équation des ondes étant linéaire, elle admet des solutions « faibles » pour des données
initiales φ et ψ moins régulières.)
On peut également résoudre le problème de Cauchy par une troisième méthode, qui a en outre le
mérite d’autoriser un terme source (correspondant à un « forçage » : on peut penser par exemple
à la gravité). Considérons l’équation des ondes avec terme source :
∂2t u − c 2 ∂2x u = f ,
et cherchons la valeur en (x 0 , t 0 ) de la solution du problème de Cauchy associé aux données initiales
u(x, 0) = φ(x), ∂t u(x, 0) = ψ(x). On appelle « cône » de dépendance (qui est en fait ici un triangle)
de (x 0 , t 0 ) l’ensemble
∆ := { (x, t ) ; x 0 − c (t 0 − t ) ≤ x ≤ x 0 + c (t 0 − t ) } .
Par la formule de Green on a
t

t0

∆
x
x0
25
 Z Z
∂2t u c 2 ∂2x u ∂t u dx + c 2 ∂x u dt .
¡ ¢ ¡ ¢
− =−
∆ ∂∆
On décompose bien sûr cette intégrale en trois morceaux. L’un vaut simplement
Z x 0 +ct
ψ(y) dy .
x 0 −ct

D’autre part, le long du segment

{ (x, t ) ; x = x 0 − c (t − t 0 )} ,
on a
∂t u dx + c 2 ∂x u dt = − c (∂t u dt + ∂x u dx ) .
Donc l’intégrale correspondante vaut
− c ( u(x 0 , t 0 ) − φ(x 0 + c t 0 ) ) .
De la même façon, l’intégrale restante vaut
− c ( u(x 0 , t 0 ) − φ(x 0 − c t 0 ) ) .
En additionnant, on obtient donc
Z Z x 0 +ct
∂2t u c 2 ∂2x u
¡ ¢
− dx dt = ψ(y) d y − 2 c u(x 0 , t 0 ) + c ( φ(x 0 + c t0 ) + φ(x 0 − c t0 ) ) ,
∆ x 0 −ct

d’où la formule
1 1
Z x0 + c t0 1
Z
u(x 0 , t 0 ) = ( φ(x 0 + c t 0 ) + φ(x 0 − c t 0 ) ) + ψ(y) dy + f dx dt
2 2c x0 − c t0 2c ∆
(dont la formule de d’Alembert est évidemment un cas particulier, avec f ≡ 0).
Problème de Dirichlet sur R+ . La résolution du problème de Dirichlet homogène


 ∂2t u − c 2 ∂2x u = 0 , x > 0 ,





 u(x, 0) = φ(x) , x > 0,



∂t u(x, 0) = ψ(x) x > 0,










 u(0, t ) = 0 , t > 0,
lorsque φ et ψ satisfont les conditions de compatibilité φ(0) = 0 et ψ(0) = 0, se ramène à celle
du problème de Cauchy par la méthode dite des images. En effet, si u est solution, alors (x, t ) 7→
u(−x, t ) est solution du problème symétrique, posé sur R− , avec φ et ψ prolongées en des fonctions
impaires : φ(x) = −φ(−x) et ψ(x) = −ψ(−x). Donc la fonction u définie comme la superposition
des deux solutions est une solution sur R tout entier, sauf peut-être en 0. Réciproquement, soit
1 1
Z x +c t
u(x, t ) = ( φ(x + c t ) + φ(x − c t ) ) + ψ(y) dy .
2 2c x −c t
Cette fonction satisfait l’équation des ondes partout où elle est deux fois dérivable. C’est le cas
si φ et ψ sont respectivement deux fois et une fois dérivable sur R+ , sauf peut-être sur le cône
caractéristique
{ (x, t ) ; x ± c t = 0} .
(Si φ satisfait la condition de compatibilité supplémentaire φ0 (0) = 0, il n’y pas de problème.) De
plus, lorsque φ et ψ sont impaires, il en est de même de u . Par conséquent, c’est bien une solution
du problème de Dirichlet homogène, sauf peut-être sur la demi-droite { x = c t }. Lorsque x > c t ,
26
les fonctions φ et ψ sont évaluées en des points positifs. Lorsque x < c t , on peut vouloir exprimer
u(x, t ) à l’aide des fonctions φ et ψ originales, définies sur R+ . On trouve ainsi la formule
1 1
Z c t +x
u(x, t ) = ( φ(c t + x) − φ(c t − x) ) + ψ(y) dy .
2 2c c t −x

Problème de Dirichlet sur un intervalle borné [0, L]. La « méthode des images » se généralise
en prolongeant φ et ψ en fonctions impaires et 2L-périodiques, comme sur le dessin ci-dessous, ce

−L 0 L 2L 3L

qui permet de calculer u à nouveau par la formule de d’Alembert :

1 1
Z x +c t
u(x, t ) = ( φ(x + c t ) + φ(x − c t ) ) + ψ(y) dy .
2 2c x −c t

Cela fournit bien une solution de l’équation des ondes, en dehors de

{ (x, t ) ; x ± c t ∈ L Z } .
Cette solution est impaire et 2L-périodique, donc satisfait les conditions de Dirichlet homogènes
en x = 0 et x = L.
Une autre approche consiste à utiliser les séries de Fourier. En effet, les fonctions prolongées φ
et ψ ont des séries de Fourier en sin (puisqu’elles sont impaires et s’annulent en ±L):
X X
φ(x) = a n sin(πnx/L) , ψ(x) = b n sin(πnx/L) ,
n∈N∗ n∈N∗

où les sommes sont à comprendre au sens de la convergence dans L2pér ([−L, L]), ou en un sens
plus fort si φ et ψ sont assez régulières (la série de Fourier d’une fonction de classe C 2 converge
uniformément). On cherche alors u sous la forme
X
u(x, t ) = c n (t ) sin(πnx/L) ,
n∈N∗

ce qui donne la famille d’équations différentielles ordinaires

c n00 + (c πn/L)2 c n = 0 ,
dont les solutions sont de la forme
c n (t ) = αn cos(πc t n/L) + βn sin(πc t n/L) .
Afin de satisfaire les données initiales, on doit avoir c n (0) = an , c n0 (0) = bn , ce qui donne
L bn
αn = a n , βn = .
πc n
On constate en particulier (comme évoqué plus haut) que la fréquence3 temporelle de u est égale
à ω = c/(2L). De plus, si cette dernière est la fréquence de la « première harmonique » c 1 , on voit
que c n a même comme fréquence le multiple n ω.

3inverse de la période
27
Problème de Cauchy sur R3 . Soit à résoudre le problème de Cauchy :
∂t u − c 2 ∆ u = 0 , x ∈ R3 , t > 0 ,
 2





u(x, 0) = φ(x) , x ∈ R3 ,




∂t u(x, 0) = ψ(x) , x ∈ R3 .


On va utiliser la méthode des moyennes sphériques, qui permet de se ramener à un problème en

dimension 1, que l’on sait résoudre d’après de qui précède. Pour cela on associe à u la fonction de
(r, t ) ∈ R+ × R+
1
Z
ū(r, t ) := u,
4πr 2 Sr
où Sr est la sphère de rayon r :
Sr := { x ∈ R3 ; kxk = r } .
Pour obtenir l’équation satisfaite par u , on applique la formule de Green sur la boule
Br := { x ∈ R3 ; kxk ≤ r } ,
∂u
Z Z
∆u = ,
Br Sr ∂n
où n est le vecteur normal unitaire sur Sr = ∂ Br , sortant de Br . Pour que u soit solution de
l’équation des ondes, il faut donc que
∂u
Z Z
∂2t u = c 2 .
Br Sr ∂n
En utilisant les coordonnées sphériques habituelles, représentées ci-dessous

θ 1
0
0
1

r
ϕ

cette égalité s’écrit de façon équivalente :

Z r Z 2πZ π Z 2πZ π
∂2t u 2
ρ sin θ dθ d3 dρ = c r 2 2
∂r u sin θ dθ d3 .
0 0 0 0 0
Or, par définition, Z 2πZ π
1
ū = u sin θ dθ d3 .
4π 0 0
On déduit donc de l’équation précédente :
Z r
ρ2 ∂2t ū dρ = c 2 r 2 ∂r ū .
0
28
D’où, en dérivant une fois par rapport à r :
r 2 ∂2t ū = c 2 ∂r ( r 2 ∂r ū ) ,
c’est-à-dire
c2
∂2t ū − c 2 ∂2r ū − 2 ∂r ū = 0 .
r
(Ceci n’est rien d’autre que l’équation des ondes axisymétrique.) En posant v := r ū , on se ramène
à une équation des ondes monodimensionnelle ordinaire :
∂2t v − c 2 ∂2r v = 0 .
Puisque v doit s’annuler en r = 0, cela revient à chercher la solution du problème de Dirichlet dans
R+ , avec v(r, 0) = r φ̄(r ), ∂t v(r, 0) = r ψ̄(r ) et les notations évidentes
1 1
Z Z
φ̄(r, t ) := φ, ψ̄(r, t ) := ψ.
4πr 2 Sr 4πr 2 Sr

On a une formule explicite pour v :

1 1
Z c t +r
v(r, t ) = ( (c t + r ) φ̄(c t + r ) − (c t − r ) φ̄(c t − r ) ) + y ψ̄(y) dy , 0 ≤ r ≤ ct,
2 2c c t −r
que l’on peut réécrire
µ
1
Z c t +r 1
¶ Z c t +r
v(r, t ) = ∂t y φ̄(y) dy + y ψ̄(y) dy .
2c c t −r 2c c t −r
On récupère ensuite u(0, t ) = ū(0, t ) par dérivation, puisque ū(0, t ) = ∂r v(0, t ). Comme
µ
1
Z c t +r ¶
1
Z
∂r y ψ̄(y) dy |r =0 = t ψ̄(c t ) = ψ,
2c c t −r 4 π c2 t Sc t

avec la même formule pour φ, on en déduit

µ ¶
1 1
Z Z
u(0, t ) = ∂t φ + ψ.
4 π c2 t Sc t 4 π c2 t Sc t

Bien sûr, on a la même formule lorsqu’on translate 0 en x 0 . Finalement, on a obtenu la formule

générale, dite de Kirchhoff :
µ ¶
1 1
Z Z
u(x 0 , t ) = ∂t φ + ψ.
4 π c2 t kx−x 0 k=ct 4 π c2 t kx−x 0 k=ct

(Il est possible de montrer cette formule en utilisant la transformation de Fourier.)

On observe sur la formule de Kirchhoff que la solution au point (x 0 , t ) ne dépend des données
initiales qu’au voisinage de la sphère { x ; kx − x 0 k = c t }. Autrement dit, si les données initiales sont
concentrées dans { x ; kx − x 0 k < R }, la solution au point (x 0 , t ) est nulle pour tout t > R/c . C’est ce
qu’on appelle le principe de Huyghens . Il est faux en dimension 2 comme on va le voir.
Problème de Cauchy sur R2 . On va se servir de la formule de Kirchhoff en dimension 3 pour
obtenir une formule en dimension 2. Il suffit en effet de résoudre le problème étendu à R3 , avec φ
et ψ indépendantes de la troisièmeRvariable x 3 =: z . Pour une fonction ψ dépendant seulement de
(x 1 , x 2 ) =: (x, y), l’intégrale double Sr ψ se décompose en
Z Z
ψ=2 ψ,
Sr Sr ∩{z > 0}

et la demi-sphère Sr ∩ {z > 0} est paramétrée par (x, y) :

q
2 2 2
Sr ∩ {z > 0} = { (x, y, z) ; x + y < r et z = r 2 − x2 − y 2 } ,
29
l’élément de surface étant
r
p dxdy .
r2 − x2 − y 2
D’où
ψ(x, y)
Z Z
ψ = 2r p dx dy .
Sr x 2 +y 2 ≤r r 2 − x2 − y 2
On en déduit la formule explicite en dimension 2 :
à !
1 φ(x, y)
Z
u(x 0 , y 0 , t 0 ) = ∂t p dx dy
2 π c (x−x0 )2 +(y−y 0 )2 ≤c 2 t 2 c 2 t 2 − x 2 − y 2
1 ψ(x, y)
Z
+ p dx dy .
2 π c (x−x0 )2 +(y−y 0 )2 ≤c 2 t 2 c 2 t 2 − x 2 − y 2
On remarque sur cette formule que, même avec des données initiales concentrées dans { x ; kx − x 0 k < R },
la solution au point (x 0 , t ) ne « s’éteint jamais ».
2.5.7. Approximation numérique. Les formules de résolution explicites que l’on vient de voir ne
sont pas toujours commodes à mettre en œuvre numériquement, d’autant qu’elles portent sur le
problème de Cauchy posé dans tout l’espace, alors qu’en pratique on doit résoudre le problème
avec des conditions aux limites. Même si on a vu comment traiter certains problèmes aux limites
par la méthode des images, cela nécessite un calcul d’intégrale en chaque point de résolution, ce
qui n’est pas façon la plus efficace de résoudre le problème. Nous allons construire et étudier une
méthode de résolution approchée, très simple à implémenter et qui donne de bons résultats.
Schéma de base. Une méthode élémentaire de discrétisation des EDP est celle dite des différences
finies, consistant à remplacer les dérivées exactes par des dérivées approchées discrètes au moyen
de la formule de Taylor. En ce qui concerne l’équation des ondes, on doit approcher des dérivées
secondes, ce que l’on peut faire en remarquant que pour toute fonction f de classe C 2 , la formule
de Taylor donne
f (y + h) − 2 f (y) + f (y − h)
f 00 (y) = + O (h) .
h2
On se limitera dans cette présentation à la dimension 1.
Étant donnés un pas d’espace ∆x et un pas de temps ∆t , on cherche à calculer u nj , supposé
approcher la solution exacte de l’équation des ondes
∂2t u − c 2 ∂2x u = 0
au point (x j , t n ), avec x j +1 − x j = ∆x et t n+1 − t n = ∆t . En appliquant la formule précédente dans
la direction du temps et dans la direction de l’espace, et en substituant les approximations ainsi
obtenues dans l’équation des ondes on obtient le schéma centré :
u n+1
j
− 2 u nj + u n−1
j
u nj+1 − 2 u nj + u nj−1
2
−c = 0.
∆t 2 ∆x 2
L’ordre en temps, respectivement en espace, d’un schéma est l’ordre en ∆t , respectivement en ∆x ,
de l’erreur de troncature, obtenue en appliquant le schéma à la solution exacte :
u(x j , t n+1 ) − 2 u(x j , t n ) + u(x j , t n−1 ) u(x j +1 , t n ) − 2 u(x j , t n ) + u(x j −1 , t n )
e nj := −c 2
.
∆t 2 ∆x 2
Dans le cas présent, on a
e nj = O (∆t ) + O (∆x) ,
le schéma est donc d’ordre 1 : c’est le minimum que l’on puisse faire pour espérer approcher
effectivement la solution exacte !
30
 Une seconde propriété importante est la stabilité d’un schéma. Une condition nécessaire de
stabilité dans `2 s’obtient facilement grâce à la transformée de Fourier discrète.
Définition 1. La transformée de Fourier discrète d’une suite u = (u n )n∈Z ∈ `2 (Z) est l’unique
fonction Fd (u) = ue ∈ L2 (T) dont les coefficients de Fourier sont les u n . Plus précisément, ue est la
limite dans L2 (T) des sommes partielles de la série
u n e− 2 i π n ω .
X
n∈Z

D’après la formule de Parseval, on a :

Z 1
2
e ω)|2 dω ,
X
|u n | = |u(
n∈Z 0

ce qui signifie que la transformation de Fourier discrète :

Fd : `2 (Z) → L2 (T)
u = (u n )n∈Z 7→ ue
est une isométrie. On observe que pour l’opérateur de « shift » :
T : u = (u j ) j ∈Z 7→ U = (U j := u j +1 ) j ∈Z , e (ζ) = e2 i πζ u(
on a (Fd ◦ T (u)) (ζ) = U e ζ) .
Le schéma centré se réécrivant
∆t 2
u n+1
=j +s ( u nj+1
+ 2(1 − u nj−1 )
− s := c s ) u nj
, u n−1
j , 2
∆x 2
(noter que s est un nombre sans dimension) ou encore, en posant v nj := u n−1
j
:
v n+1
   n 
vj

j 0 1
= ,
    
 
n
u n+1
j
−1 2 ( 1 − s ) + s (T + T )
−1
uj

par transformation de Fourier discrète (en espace seulement !), on obtient

vbn+1 (ζ) vbn (ζ)
   
  = H(ζ)  ,
n+1 n
ub (ζ) ub (ζ)
où H(ζ) est ce qu’on appelle la matrice d’amplification du schéma. Elle vaut simplement
   
0 1 0 1
H(ζ) =  = .
2 i πζ − 2 i πζ 2 ( 1 − s ( 1 − cos 2 πζ ) )
−1 2(1 − s ) + s (e +e ) −1
L’itération ci-dessus se résout explicitement en
vbn (ζ) vb0 (ζ)
   
  = H(ζ)n  .
n 0
ub (ζ) ub (ζ)
On pourra donc contrôler la norme L2 de la solution à l’itération n en fonction de la norme L2
initiale si et seulement si la norme de la norme de la matrice H(ζ)n est uniformément bornée
(en n et ζ). Par la propriété d’isométrie de la transformation de Fourier discrète, ceci fournit
une condition nécessaire et suffisante de stabilité `2 du schéma. Dans l’état, elle n’est cependant
pas très exploitable. On n’a pas précisé la norme matricielle à utiliser. En fait, il suffit qu’il en
existe une telle que H(ζ)n soit bornée. Une condition nécessaire pour cela, dite condition de von
31
Neumann, est que le rayon spectral de H(ζ) soit inférieur ou égal à 1. Or les valeurs propres de
H(ζ) sont les racines de
λ2 − 2 b(ζ) λ + 1 = 0 , b(ζ) := 1 − s ( 1 − cos 2 πζ ) .
Le produit de ces racines vaut donc 1. Si jamais elles sont réelles, l’une sera nécessairement de
valeur absolue supérieure à 1, et la condition de von Neumann sera violée. Il faut donc que le
discriminant de cette équation soit négatif (ou nul), c’est-à-dire que |b(ζ)| ≤ 1 quel que soit ζ. Ceci
revient à demander s ≤ 1. Réécrite en fonction de ∆t et ∆x , cette condition est
c ∆t ≤ ∆x .
Ce type de condition, majorant le pas de temps en fonction du pas d’espace, est appelée condition
de Courant-Friedrichs-Lewy (CFL en abrégé). Par extension, le nombre sans dimension c ∆t /∆x
est souvent appelé nombre CFL. On constate en pratique que si cette condition n’est pas satisfaite,
le schéma « diverge » (au sens où les valeurs numériques obtenues en guise de solution approchée
deviennent démesurément grandes) au bout de quelques itérations.

32
 Index
énergie condition de, 32
cinétique, 7, 10 nombre de, 32
potentielle, 7, 10 courbe de choc, 19
totale, 7 cycle, 8
équation cyclique
aux dérivées partielles d’évolution, 14 variable, 12
de réaction-diffusion, 20 cycloı̈de, 10
de transport linéaire, 17
quasi-linéaire, 17 d’Alembert
conservative, 24 formule de , 25
de continuité, 14 d’Alembertien, 24
de la chaleur, 16 débit, 14
de Ludwig, 4 densité
de Malthus, 2 de charge électrique, 14
de Verhulst, 3 de courant électrique, 14
linéarisée, 17 différences finies, 30
logistique, 4 diffusion, 16
équation des ondes axisymétrique, 29 diffusivité thermique, 16
équation différentielle discrétisation
du second ordre, 6 pas de, 3
linéaire du premier ordre, 2 schéma de, 3
linéaire du premier ordre avec terme source, 3 discret
équation scalaire, 7 modèle, 2
Duhamel
équation de transport, 24 formule de, 3
accélération, 6 entraı̂nement
accélération de la pesanteur, 6 vitesse et accélération, 13
adimensionnalisation, 4 erreur de troncature, 30
Ampère, 14 Euler
angles d’Euler, 11 angles d’, 11
auto-similaire, 18, 19 Euler–Lagrange
équation d’, 9
brachistochrone, 9
Faraday, 14
cône caractéristique, 26
file d’attente, 17
cône de dépendance, 25
flux, 14
calcul des variations, 9
flux de chaleur, 15
Cauchy
force, 6
problème de, 6, 18
forme optimale, 9
centre de gravité, 11
formule d’Euler, 3
champ
Fourier
de gravité, 10
série de, 27
de vitesse, 16
cinématique, 11 Green
coefficient de diffusion, 16 formule de , 25
conservation Green–Riemann
de l’énergie totale, 7 formule de, 15
conservativité, 5
convection-diffusion, 20 harmonique, 21
convention, 16 Huygens
corde vibrante, 20 principe de, 29
Coriolis hyperbolique
accélération de, 13 EDP, 16
Coulomb, 14
Courant-Friedrichs-Lewy inertie
33
 moment principal d’, 12 d’accroissement naturel, 2
repère principal, 12 de mortalité, 2
intégrale première, 5, 12 de natalité, 2
invariance par le flot, 5 temps d’explosion, 18
tension, 21
Kirchhoff terme source, 3, 20
formule de , 29 théorème
π de Buckingham, 4
Lighthill-Whitham-Richards
de Cauchy–Lipschitz, 5
modèle de, 16
trajectoire, 5
loi
trajectoire parabolique, 6
de conservation, 14
transformation de Fourier discrète, 31
de conservation non-linéaire, 16
de la gravitation universelle, 10 variables séparées, 5
de fermeture, 16 vitesse, 6
de Fick, 16 von Neumann
de Fourier, 16 condition de , 32
de Newton, 6

méthode
des caractéristiques, 17
méthode des images, 26
matrice d’amplification, 31
moyenne sphérique, 28

Newton
loi de, 6
nutation, 11

onde de choc, 17, 19

onde de détente, 17, 18
orbite hétérocline, 8
ordre, 30
oscillateur harmonique, 6

pendule simple, 7
phases
espace des, 10
plan de phase, 6, 8
points d’équilibre, 8
portrait de phase, 3, 8
position, 6
potentiel, 10
précession, 11

réaction, 20
Rankine–Hugoniot
condition de, 19
relative
vitesse et accélération, 13
Riemann
problème de, 18
rotation propre, 11

schéma d’Euler explicite, 3

solution périodique, 6, 13
stabilité, 31

taux
34
 Références
[1] S. S. Antman. The simple pendulum is not so simple. SIAM Review, vol. 40, n° 4, pp. 927–930, 1998.
[2] S. Benzoni-Gavage. Calcul différentiel et équations différentielles. Cours et exercices corrigés. Dunod, Collection
Sciences Sup, 2010.
[3] M. Braun. Differential equations and their applications, volume 11 of Texts in Applied Mathematics. Springer-
Verlag, New York, fourth edition, 1993. An introduction to applied mathematics.
[4] C. Chicone. Ordinary differential equations with applications, volume 34 of Texts in Applied Mathematics.
Springer-Verlag, New York, 1999.
[5] J.-P. Demailly. Analyse numérique et équations différentielles. Collection Grenoble Sciences. Presses Universi-
taires de Grenoble, Grenoble, 1996.
[6] J. David Logan. Applied mathematics. Wiley-Interscience [John Wiley & Sons], Hoboken, NJ, third edition,
2006.
[7] J. D. Murray. Mathematical biology. I, volume 17 of Interdisciplinary Applied Mathematics. Springer-Verlag,
New York, third edition, 2002. An introduction.
[8] L. Schwartz. Méthodes mathématiques pour les sciences physiques. Hermann, Paris, 1961.

35