Université Mohammed Premier

Ecole Nationale des Sciences Appliquées

Oujda

Modélisation Mathématique
Filière Génie Civil

Première année

2022-2023

Imad EL MAHI

ADRESSE : COMPLEXE UNIVERSITAIRE, HAY EL QODS, B.P 669, 60000 OUJDA, MAROC
Table des matières

1 Equations de la physique - Classication des EDP 3

1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Généralités sur les EDP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Classication des EDP linéaires du 2ème ordre : . . . . . . . . . . . . . . . . . . . 6
1.4 Quelques modèles d'équations hyperboliques et paraboliques . . . . . . . . . . . . 8
1.4.1 Equations hyperboliques . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4.2 Equations paraboliques . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.5 Méthode des caractéristiques pour les EDP du 1er ordre . . . . . . . . . . . . . . 10

2 Méthode des diérences nies 17

2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2 Principe de la méthode des diérences nies - Ordre de précision . . . . . . . . . 17
2.2.1 Principe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2.2 Ordre de la méthode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.3 Diérences nies pour une équation diérentielle du 1er ordre . . . . . . . . . . . 18
2.4 Diérences nies pour un problème elliptique 1D . . . . . . . . . . . . . . . . . . 22
2.4.1 Modèle mathématique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.4.2 Approximation aux diérences nies . . . . . . . . . . . . . . . . . . . . . 22
2.5 Diérences nies pour un problème parabolique 1D (Equation de la chaleur) . . . 24
2.5.1 Modèle mathématique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.5.2 Approximation diérences nies . . . . . . . . . . . . . . . . . . . . . . . . 25
2.6 Consistance, stabilité et convergence des schémas aux diérences nies . . . . . . 27
2.6.1 Erreur de troncature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.6.2 Consistance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.6.3 Stabilité d'un schéma numérique . . . . . . . . . . . . . . . . . . . . . . . 29
2.6.4 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.6.5 Théorème de Lax . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.6.6 Stabilité au sens de Fourier - Von Neumann . . . . . . . . . . . . . . . . . 29

2
Chapitre 1

Equations de la physique -
Classication des EDP
1.1 Introduction
La modélisation consiste à remplacer un système complexe en un objet simple qui reproduit
les comportements principaux du système à étudier. Ce modèle peut être un modèle réduit, une
expérience au laboratoire, un modèle numérique, etc.
On s'intéresse ici à la modélisation numérique qui consiste à remplacer le problème physique,
chimique, biologique, économique ou industriel par des équations mathématiques. Ces équations
s'écrivent souvent sous forme d'EDP (Équations aux dérivées partielles) et sont généralement
non linéaires et posées en deux ou trois dimensions d'espaces avec des conditions aux limites
parfois compliquées. Leurs solutions analytiques (exactes) sont souvent très diciles à calculer
voire impossible. On fait donc recours aux méthodes numériques pour approcher les solutions de
ces équations.

La modélisation numérique se décompose en 5 étapes :

1. Établissement d'un modèle mathématique qui décrit le problème physique ou industriel

que l'on veut résoudre.

2. Discrétisation des équations, ce qui permet d'avoir un modèle mathématique posé en un

nombre ni de n÷uds obtenus par une discrétisation (maillage) du domaine du calcul.

3. Méthodes numériques pour la résolution du problème discret. Ce point est crucial car la
validité du modèle dépend étroitement de la précision des méthodes utilisées. A ce niveau,
souvent une analyse d'erreurs et de convergence de la méthode est requise.

4. Mise en ÷uvre informatique. Ceci consiste à reproduire les algorithmes développés par un
langage de programmation (Fortran, Matlab, Python, ...).

5. Analyse, interprétation des résultats et évaluation de la précision du modèle et des mé-

thodes utilisées. Ceci se fait en général en comparant avec les études expérimentales exis-
tantes, ou par des études d'analyse mathématique d'ordre et convergence des méthodes
utilisées.

1.2 Généralités sur les EDP

Dénition 1.2.1 (Dénition d'une EDP) Soit u = u(x1 , x2 , . . . , xn ) une fonction de plu-
sieurs variables indépendantes. Par exemple, x1 peut représenter la variables temps t et x2 , x3 ,
x4 peuvent représenter les variables d'espace x, y et z .

3
Une équation aux dérivées partielles (ou en abrégé EDP) pour la fonction u est une relation
liant la fonction u, les variables x1 , x2 , . . . , xn et les dérivées partielles de u par rapport à ces
variables. Elle s'écrit en général sous la forme :
∂u ∂u ∂ 2 u ∂ 2 u
F (x1 , x2 , . . . , xn , u, ,..., , , , . . .) = 0 (1.1)
∂x1 ∂xn ∂x21 ∂x1 ∂x2

Remarque 1:
Les EDP interviennent dans beaucoup d'applications en physique, biologie, économie, ingénierie,
génie civil, etc :

ˆ En électromagnétisme → Équations de Maxwell

ˆ En mécanique des uides → Équations de Navier-Stokes
ˆ En mécanique quantique → Équations de Shrödinger
ˆ En nance → Équations de Black & Scholes
ˆ En génie civil → Équations d'équilibre, équations de propagation, équations de l'hydrau-
lique.

Exemple : Equation de diusion (1D)

∂u ∂ 2 u
− 2 =0 (E)
∂t ∂x
1
u1 (t, x) = t + x2 est une solution de (E) dénie sur R2 .
2
1 − x2 ∗
u2 (t, x) = √ e 4t est aussi solution de (E) sur R+ × R (c'est-à-dire pour t>0 et x ∈ R) dite
4πt
solution fondamentale ou noyau de la chaleur.
1
On peut représenter par exemple la solution u2 en fonction de x pour t= .
4

et on a
+∞ +∞
1 √
Z Z
1 1 2
u2 ( , x) dx = √ e−x dx = √ π = 1
−∞ 4 −∞ π π

Remarque 2:
Il est clair que pour avoir une solution unique de l'EDP (E), il faut rajouter à l'équation une
condition initiale et des conditions aux limites.

4
Exercice 1 Montrer que la solution de s'exprime, quand le domaine est non borné (inni),
(E)
comme le produit de convolution de la condition initiale par le noyau de la chaleur :
Z +∞
u(t, x) = u(0, x) ? u2 (t, x) = u(0, x − y) u2 (t, y) dy
−∞

Dénition 1.2.2 (Ordre d'une EDP) On appelle ordre d'une EDP l'ordre le plus élevé des
dérivées que contient l'EDP.

Dénition 1.2.3 (EDP linéaire - EDP non linéaire) Une EDP est dite linéaire quand elle
l'est par rapport à u et par rapport à toutes ses dérivées partielles.
En fait, les EDP linéaires ne contiennent aucun produit de u avec elle même ou avec ses dérivées
partielles alors que les EDP non linéaires peuvent en avoir.
 Par exemple, une EDP linéaire du 1er ordre pour la fonction u(x, y) a la forme générale :

∂u ∂u
A +B + Cu = D (1.2)
∂x ∂y
où A, B , C et D peuvent être des fonctions de x et y .

 La forme générale d'une EDP linéaire du 2ème ordre pour u(x, y) est :

∂2u ∂2u ∂2u ∂u ∂u

A2
+ B + C 2
+D +E + Fu = G (1.3)
∂x ∂x∂y ∂y ∂x ∂y
où A, B , C , D, E , F et G peuvent dépendre de x et y .

Dénition 1.2.4 (EDP homogène - EDP non homogène) Une EDP est dite homogène lors-
qu'elle ne contient que des termes faisant intervenir la fonction u et ses dérivées partielles.
 Par exemple, l'EDP linéaire du 1er ordre (1.2) serait homogène si D = 0. Celle du 2ème
ordre (1.3) est homogène si G = 0.

Exemples :
∂u ∂u
+ =0 EDP du 1er ordre, linéaire et homogène (à coecients constants).
∂x ∂y

∂2u ∂2u
+ =0 (Eq. d'équilibre) est une EDP du 2ème ordre, linéaire et homogène.
∂x2 ∂y 2

∂u ∂u
+ + sin(u) = 0 1er ordre, non linéaire et homogène.
∂x ∂y

∂u ∂u
x + y2 + 4xyu = 1 1er ordre, linéaire et non homogène.
∂x ∂y

∂2u ∂2u
− =0 (Eq. de propagation des ondes) 2ème ordre, linéaire et homogène.
∂t2 ∂x2
∂u ∂u ∂2u
+u =ν 2 2ème ordre, non linéaire et homogène.
∂t ∂x ∂x
∂u 3 ∂u ∂ 2 u
( ) + · =0 2ème ordre, non linéaire et homogène.
∂x ∂x ∂y 2

5
1.3 Classication des EDP linéaires du 2ème ordre :
On considère une EDP linéaire du 2ème ordre sous forme générale en (2D ) :

∂2u ∂2u ∂2u ∂u ∂u

A 2 +B + C 2 +D +E + Fu = G (1.4)
∂x ∂x∂y ∂y ∂x ∂y
| {z }
(P)

où A, B , C , D , E , F et G peuvent dépendre des variables x et y.

La partie (P) s'appelle partie principale de l'EDP (1.4). En fait, c'est elle qui donne la nature
de l'EDP.

On note par ∆ = B 2 − 4AC le discriminant associé à l'EDP (1.4).

Classication :
On dira que l'EDP (1.4) est

 hyperbolique si ∆ > 0,
 parabolique si ∆ = 0,
 elliptique si ∆ < 0.

Remarques :
1. Les termes hyperbolique, parabolique, elliptique proviennent du fait que si l'on sup-
pose que A, B , C , D , E et F constants et que l'on cherche une solution ϕ(x, y) de
l'équation homogène associée à (1.4) (G = 0) de la forme :

ϕ(x, y) = exp(r1 x+r2 y)

On obtient l'équation algébrique pour r1 , r2 , dite équation caractéristique, suivante :

Ar12 + Br1 r2 + Cr22 + Dr1 + Er2 + F = 0

qui est de type :

AX 2 + BXY + CY 2 + DX + EY + F = 0 (A)

et qui représente l'équation d'une conique.

On rappelle que la conique d'équation (A) représente

une hyperbole si ∆ > 0, une parabole si ∆ = 0, une ellipse si ∆ < 0.

2. Si les coecients A, B et C dépendent des variables (x, y) alors la dénition du type de

l'EDP (1.4) devient locale, c'est-à-dire que l'EDP (1.4) sera

 hyperbolique au point (x0 , y0 ) si ∆(x0 , y0 ) > 0,

 parabolique au point (x0 , y0 ) si ∆(x0 , y0 ) = 0,
 elliptique au point (x0 , y0 ) si ∆(x0 , y0 ) < 0.
avec ∆(x0 , y0 ) = B 2 (x0 , y0 ) − 4A(x0 , y0 )C(x0 , y0 )
C'est-à-dire que l'EDP peut être hyperbolique dans certaines régions, parabolique ou
elliptique dans d'autres régions.

6
Exemples :
1. Equation des ondes (1D) :
∂2u 2
2∂ u
− c =0 où c représente la vitesse du son.
∂t2 ∂x2
Variables −→ (t, x)
Coecients : A = 1 ; B = 0 ; C = −c2 ; D = E = F = G = 0
Le discriminant de l'EDP : ∆ = B 2 − 4AC = 4c2 > 0
L'équation des ondes est donc une EDP hyperbolique.

2. Equation de diusion (1D) :

∂u ∂2u
−ν 2 =0 (ν >0 est le coecient de diusion)
∂t ∂x
Variables −→ (t, x)
Coecients : A = 0; B = 0; C = −ν ; D = 1; E=F =G=0
∆= B2 − 4AC = 0 =⇒ L'équation de diusion est parabolique.

3. Equation de Laplace :
∂2u ∂2u
+ 2 =0
∂x2 ∂y

Variables −→ (x, y)
Coecients : A = 1; B = 0; C = 1; D = E = F = G = 0
∆ = −4 < 0 =⇒ L'équation de Laplace est elliptique.

4. Equation de Tricomi :

∂2u ∂2u
y − 2 =0
∂x2 ∂y

Variables −→ (x, y)
Coecients : A = y ; B = 0 ; C = −1 ; D = E = F = G = 0
Le discriminant de l'EDP : ∆ = 4y
Si y>0 alors l'EDP est hyperbolique.

Si y=0 alors l'EDP est parabolique.

Si y<0 alors l'EDP est elliptique.

L'équation de Tricomi peut modéliser par exemple les écoulements transsoniques non
visqueux. On voit que c'est une EDP mixte : hyperbolique sur le demi plan supérieur
y > 0, parabolique sur l'axe horizontal y = 0 et elliptique sur le demi plan inférieur y < 0.

7
1.4 Quelques modèles d'équations hyperboliques et paraboliques
1.4.1 Equations hyperboliques
Les équations hyperboliques modélisent les phénomènes de propagation d'ondes sans aucune
dissipation, c'est-à-dire que l'énergie produite par ces ondes ne s'atténue pas au cours du temps.
On parle de convection pure. On les trouve dans divers problèmes de la physique.
Dans le cas linéaire, ces équations peuvent modéliser la propagation du son dans un milieu
homogène. En électromagnétisme, les équations de Maxwell sont hyperboliques et linéaires.
Dans le cas non linéaire, ces équations proviennent de lois de conservation de certaines quantités.
C'est le cas par exemple des équations d'Euler qui expriment la conservation de la masse, de la
quantité de mouvement et de l'énergie pour un uide parfait (non visqueux). En hydraulique,
les équations de Saint-Venant modélisent les écoulements de l'eau peu profonde (océans, rivières,
écoulements suite aux rupture de barrages, etc).

Exemples types de problèmes hyperboliques

1. Equation des ondes :

∂2u 2
2∂ u
−

c =0 (1.5)
∂t2 ∂x2
u(t = 0, x) = u0 (x) donnée (condition initiale)


Cette équation modélise la propagation d'une onde sonore à vitesse c dans milieu homo-
gène. On l'appelle aussi équation de la corde vibrante.
2
Si l'on prend c = 1 et comme condition initiale la fonction u0 (x) dénie par u0 (x) = e−x ,

alors à t = 1 sec et à t = 5sec, on obtient en résolvant l'équation (1.5) :

8
 2. Equation de convection :

 ∂u ∂u
+c =0
∂t ∂x (1.6)
 u(t = 0, x) = u0 (x) donnée

C'est l'équation de transport de la quantité u sans dissipation avec la vitesse c.

On peut montrer que la solution de cette équation hyperbolique s'écrit :

u(t, x) = u0 (x − ct) ←− condition initiale u0 (x) transportée avec la vitesse c.

c'est-à-dire que si l'on prend la condition initiale u0 (x) comme suit :

alors la simulation avec une vitesse c = 2m/sec, donne la solution suivante à t = 3 sec :

1.4.2 Equations paraboliques

Ces équations modélisent la diusion d'une quantité u dans un milieu. Cela peut être la
diusion de la température dans une barre métallique ou dans un milieu bidimensionnel, ou la
diusion d'une quantité u d'un polluant dans une rivière.

Exemple :

 ∂u ∂2u
−ν 2 =0 (1.7)
 ∂t ∂x
u(t = 0, x) = u0 (x)

9
 Cette équation permet de simuler la diusion de la quantité u(t, x) où ν est le coecient de
diusion. Il est clair que plus ν est grand plus la quantité u(t, x) diuse plus vite.
2
Si l'on prend à t = 0, la condition initiale u(0, x) = u0 (x) = e−x , alors à t > 0, on aura :

1.5 Méthode des caractéristiques pour les EDP du 1er ordre

On considère l'EDP linéaire du 1er ordre, sous sa forme générale :

∂u ∂u
P +Q +Cu=D (1.8)
∂x ∂y
où P , Q, C et D sont des fonctions de x et y.
(1.8) s'écrit aussi :
∂u ∂u
P +Q =R
∂x ∂y
avec R=D−Cu

Le principe de la méthode des caractéristiques pour la résolution de l'EDP (1.8) et de chercher

des courbes dites courbes caractéristiques le long desquelles l'EDP (1.8) se ramène à une simple
équation diérentielle ordinaire (EDO).

s
)
(s)
s ),y
(x(
Courbe caractéristique

Si on arrive à résoudre l'EDO sur la courbe caractéristique Γ, alors on pourra répéter la procédure
en partant de la courbe voisine, etc et donc obtenir la solution u(x, y) dans tout le domaine.

On va illustrer le principe sur l'exemple suivant :

10
Equation de convection

 ∂u ∂u
+c =0 pour x ∈ R ; t > 0
∂t ∂x
 u(0, x) = u0 (x) ←− condition initiale

Nous allons chercher une courbe caractéristique Γ ((t(s), x(s)), s étant le paramètre qui décrit la
courbe, le long de laquelle l'EDP devient un système d'EDO.

Dérivons u le long de la courbe Γ :

du ∂u dt ∂u dx
= +
ds ∂t ds ∂x ds
Or d'après l'EDP
∂u ∂u
= −c
∂t ∂x
On a alors :  
du ∂u dx dt
= −c
ds ∂x ds ds

dx dt
On voit que si on impose : −c =0 c'est-à-dire dx = c dt
ds ds
dx
ou encore =c
dt
on obtient :
du
=0
ds
Les courbes caractéristiques sont donc les droites d'équation :

dx
=c
dt
et sur ces courbes caractéristiques la solution vérie :

du = 0

c'est-à-dire que u reste constante sur chaque caractéristique (on parle d'invariants de Riemann).

Le système d'EDO à résoudre est donc :


 dx
=c qui donne la courbe caractéristique Γ
dt
 du =0 qui donne la solution sur cette courbe caractéristique

Courbes caractéristiques
dx
= c ⇔ x = ct + ξ (ξ ∈ R)
dt

11
Solution
Sur chaque courbe caractéristique (Γ) : x − ct = ξ , on a :

du = 0 =⇒ u(t, x) = cte = f (ξ) ←− c-à-d u ne dépend que de ξ

↑
f ∈ C 1 (R)
Soit alors
u(t, x) = f (x − ct)

Cette solution doit être retrouvée aussi à t = 0.

or à t = 0, on a :

u(0, x) = u0 (x) =⇒ f (x) = u0 (x)

c-à-d f ≡ u0
On obtient nalement :
u(t, x) = u0 (x − ct)

Retour au cas général

∂u ∂u
P +Q =R
∂x ∂y
Cherchons une courbe caractéristique Γ = Γ(x(s), y(s)) le long de laquelle l'EDP se réduit à une
équation diérentielle.

Dérivons u le long de la courbe Γ :

du ∂u dx ∂u dy
= +
ds ∂x ds ∂y ds
Soit
du ∂u dx ∂u dy
P = P +P
ds ∂x ds ∂y ds
∂u dx ∂u dy
= (R − Q ) +P
∂y ds ∂y ds
∂u dy dx dx
= (P −Q )+R
∂y ds ds ds

12
On remarque alors que dès lors qu'on impose

dy dx
P −Q =0
ds ds
on aura
du dx
P =R
ds ds
Autrement dit, le long des courbes caractéristiques :

P dy = Q dx (1.9)

la solution vérie :
P du = R dx (1.10)

Comme moyen de retenir ces relations, il sut d'écrire :

dx dy du
= = (1.11)
P Q R

Attention : Lorsque P , Q, ou R sont nuls, il vaut mieux utiliser (1.9) et (1.10)).

(

Remarques :
1. La relation P dy = Q dx est appelée direction caractéristique. La pente locale de la ca-
ractéristique est :
dy Q
=
dx P
2. Dans le cas où R = 0, la méthode des caractéristiques donne : du = 0 le long de la
caractéristique, c'est-à-dire que u est constance le long de la caractéristiques (on appelle
cela un invariant de Riemann).

3. La méthode des caractéristiques est particulièrement adaptée aux problèmes hyperbo-

liques (régis par un transport convectif d'une quantité).

Exercice 2 Utiliser la méthode des caractéristiques pour résoudre les EDP suivantes :
∂u ∂u
+ =0 (1)
∂x ∂y
∂u ∂u
+y =0 (2)
∂x ∂y
∂u ∂u
x +2 − 2u = 0 (3)
∂x ∂y

∂u ∂u
+ =0 (1)
∂x ∂y

Cherchons une courbe caractéristique Γ = Γ (x(s), y(s)) le long de laquelle l'EDP se réduit à une
équation diérentielle. Dérivons le long de la courbe Γ

13
 du ∂u dx ∂u dy
= +
ds ∂x ds ∂y ds
∂u dx ∂u dy
= − +
∂y ds ∂y ds
 
∂u dy dx
= −
∂y ds ds
On voit donc que dès lors qu'on pose

dy dx
− =0
ds ds
c'est-à-dire le long de la courbe caractéristique

dy = dx

ou encore
y − x = cte = ξ
la solution u vérie :
du
=0
ds
c'est-à-dire u est constante le long de chaque courbe Γ (invariant de Riemann).

Autrement dit, u ne dépend que de ξ d'où :

La solution générale de l'EDP (1) est : u(x, y) = f (ξ) = f (y − x) où f ∈ C 1 (R)

Plus simple pour la résolution

On peut écrire directement le système d'EDO (P = 1 ; Q = 1 ; R = 0)
dx dy

 1 = 1
 
 y = x + ξ ←− courbes caractérisiquess


dx dy du
= = ⇔ ⇔
1 1 0
 du = dx du = 0 ←− solution
 

 sur ces courbes
0 1

14
 ∂u ∂u
+y =0 (2)
∂x ∂y

Manière 1
On dérive le long d'une courbe caractéristique Γ :

du ∂u dx ∂u dy
= +
ds ∂x ds ∂y ds
∂u dx ∂u dy
= (−y ) +
∂y ds ∂y ds
∂u dy dx
= ( −y )
∂y ds ds
Le long de la courbe caractéristique

dy dx
−y c-à-d dy = y dx
ds ds
La solution vérie
du
=0 c-à-d du = 0
ds
Les courbes caractéristiques sont

dy
dy = y dx ⇔ = dx ln |y| = x + K
y

y = ξex
ou encore
ye−x = ξ
Le long de ces courbes, la solution u est constante c'est-à-dire :

La solution générale de l'EDP (2) est : u(x, y) = f (ξ) = f (ye−x ) où f ∈ C 1 (R)

Les courbes caractéristiques sont :

15
Plus simple
Écrire le système d'EDP sur les caractéristiques (P = 1 ; Q = y ; R = 0)

dx dy
=

 
 y dx = dy ←− caractéristiques

dx dy du  1
 y
= = ⇔ ⇔
1 y 0
du dx du = 0 ←−
 

 solution

 =
0 1

∂u ∂u
x +2 − 2u = 0 (3)
∂x ∂y

On dérive le long d'une caractéristique :

du ∂u dx ∂u dy
= +
ds ∂x ds ∂y ds
∂u dx 1 ∂u dy
= + (2u − x )
∂x ds 2  ∂x ds
∂u dx x dy dy
= − +u
∂x ds 2 ds ds

Les courbes caractéristiques ont pour équation :

dx x dy dx 1
− =0 c-à-d = dy
ds 2 ds x 2
1
ln |x| = y + K ou x = ξey/2
2

La solution vérie le long de ces caractéristiques :

du dy du
=u c-à-d = dy
ds ds u
ln |u| = y + K soit u = Cey

16
Chapitre 2

Méthode des diérences nies

2.1 Introduction
On distingue parmi les méthodes numériques pour la résolution des équations aux dérivées
partielles trois grandes familles de méthodes :
ˆ La méthode des diérences nies qui s'utilise surtout sur des domaines réguliers (rectangles
en 2D et hexaèdres en 3D).
ˆ La méthode des volumes nis qui peut être utilisée sur des géométries complexes pour des
EDP représentant la conservation de certaines quantités physiques (Exemple : Équations
de Navier-Stokes).
ˆ La méthode des éléments nis qui peut être aussi appliquée sur des géométries complexes
et qui est largement utilisée pour la mécanique des uides et surtout en mécanique des
structures.

2.2 Principe de la méthode des diérences nies - Ordre de pré-

cision
2.2.1 Principe
On considère une EDP modélisant un problème de la physique ou autre. La méthode des
diérences nies appliquée à l'EDP consiste à approximer les dérivées partielles dans l'EDP
en utilisant des développements limités de Taylor. En procédant ainsi, les dérivées partielles de
l'EDP sont donc remplacées par leurs approximations diérences nies constituées par les valeurs
de l'inconnu en un certain nombre de points du maillage. On aboutit généralement à un schéma
itératif déni par des équations algébriques ou à des systèmes linéaires.

2.2.2 Ordre de la méthode

On se place en 3D. Soit u(x, y, z, t) une fonction désignant l'inconnue d'une EDP provenant
de la physique.

On a par dénition :

∂u u(x + ∆x, y, z, t) − u(x, y, z, t)

(x, y, z, t) = lim ←− Cette relation est exacte.
∂x ∆x→0 ∆x
∂u
Donnons maintenant une approximation de (x, y, z, t).
∂x

17
Faisons un DL de u(x + ∆x, y, z, t) au voisinage de (x, y, z, t) avec ∆x → 0.

∂u (∆x)2 ∂ 2 u
u(x + ∆x, y, z, t) = u(x, y, z, t) + ∆x (x, y, z, t) + (x, y, z, t)
∂x 2! ∂x2
(∆x)n ∂ n u
(x, y, z, t) + O (∆x)n+1 )


+......... + n
n! ∂x
Si on tronque ce DL à l'ordre 1 en ∆x, on aura :

∂u
u(x + ∆x, y, z, t) = u(x, y, z, t) + ∆x (x, y, z, t) + O((∆x)2 )
∂x
On obtient alors :

∂u u(x + ∆x, y, z, t) − u(x, y, z, t)

(x, y, z, t) = + O(∆x)
∂x ∆x
∂u
Nous venons en fait de donner une une approximation diérences nies de qui est :
∂x
∂u u(x + ∆x, y, z, t) − u(x, y, z, t)
(x, y, z, t) '
∂x ∆x
avec une précision d'ordre 1 (puissance de ∆x dans l'erreur de troncature).

O(∆x) est appelé Erreur de troncature.

Et on a la dénition :

Dénition 2.2.1 On appelle ordre d'une méthode, la puissance de ∆x avec laquelle l'erreur de
troncature tend vers 0.

2.3 Diérences nies pour une équation diérentielle du 1er ordre

Nous allons appliquer la méthode des diérences nies pour la résolution numérique de l'équa-
tion diérentielle suivante :

u0 (x) + u(x) = x − 1

pour x ∈ ]0, 2]
(2.1)
u(0) = 1

La méthode de diérences nies se compose en 3 étapes :

Etape 1 : Discrétisation du domaine : Maillage

On subdivise l'intervalle [0, 2] en N points xi (i = 1, . . . , N ) régulièrement espacés.

a = x1 x2 x3 xi−1 xi xi+1 xN = b
b−a
On note par ∆x = le pas du maillage (ici a = 0 et b = 2).
N −1
Les n÷uds du maillage sont :

x1 = a ; x2 = a + ∆x ; . . . ; xi = a + (i − 1)∆x

18
Le problème continue (2.1) qui consiste à chercher la fonction u(x) sur l'intervalle [a, b] de di-
mension inni se ramène alors à la recherche de N valeurs discrètes ui en chaque n÷ud xi du
maillage.

On dit qu'on a remplacé le problème continu par un problème discret.

Le problème discret s'écrit alors :


 Chercheru(xi ) pour i = 2, . . . , N vériant
u0 (xi ) + u(xi ) = xi − 1 (2.2)
u(x1 ) = 1


Etape 2 : Construction du schéma numérique

Cette étape consiste à approximer les opérateurs dérivées u0 (xi ) en chaque n÷ud xi du
maillage.

Notons par ui = u(xi ) la valeur discrète de u(x) au n÷ud xi et u0i = u0 (xi ) la dérivée au n÷ud
xi .

Schéma d'ordre 1 pour l'approximation de u0i

On a déjà vu que :
u(x + ∆x) − u(x)
u0 (x) = + O(∆x)
∆x
Au n÷ud xi , on aura :
ui+1 − ui
u0i = + O(∆x)
∆x
Une approximation de l'ordre 1 de u0i est donc :

ui+1 − ui
u0i '
∆x
qui est un schéma aux diérences nies décentré "aval".

On peut aussi construire un schéma d'ordre 1 décentré "amont" en écrivant :

ui − ui−1
u0i '
∆x

Schéma d'ordre supérieur pour l'approximation de u0i

Les deux choix d'approximation d'ordre 1 pour u0i ne sont pas uniques. On peut construire
des schémas d'ordre supérieur en manipulant les DL au voisinage de xi .

Ecrivons les DL de ui+1 et ui−1 au voisinage de xi :

(∆x)2 00
ui+1 = u(xi + ∆x) = ui + ∆x u0i + ui + O (∆x)3 )


(A)
2!
(∆x)2 00
ui−1 = u(xi − ∆x) = ui − ∆x u0i + ui + O (∆x)3 )


(B)
2!
En faisant (A)-(B), on obtient :

ui+1 − ui−1 = 2∆x u0i + O (∆x)3

19
c'est-à-dire :
ui+1 − ui−1
u0i = + O (∆x)2


2∆x
Autrement dit, une approximation d'ordre 2 "centrée" de u0i est :

ui+1 − ui−1
u0i '
2∆x

Remarque 1:
1. En général, pour obtenir des ordres encore supérieurs, il faut utiliser plusieurs n÷uds
voisins de xi . Le nombre de n÷uds nécessaires pour l'écriture du schéma numérique est
appelé "stencil".

2. An de pouvoir conclure sur l'ordre d'un schéma, il est conseillé de pousser les dévelop-
pements limités susamment loin.

3. En théorie, la précision d'un schéma d'ordre supérieur est plus élevée. Mais attention, ceci
n'est vrai que dans la limite ∆x, ∆t → 0. D'ailleurs, parfois un schéma d'ordre 1, peut
s'avérer être plus précis qu'un schéma d'ordre supérieur à ∆x, ∆t nis.

Par exemple, nous avons l'inégalité :

x ≤ 10 x2 ⇐⇒ x ≥ 0.1

C'est-à-dire que pour x ≥ 0.1, la précision à l'ordre 1 est plus élevée.

Exercice 3 Ecrire un schéma aux diérences nies d'ordre 3 par l'approximation de la dérivée
u0i au n÷ud xi (Indication : utiliser les valeurs de u aux n÷uds xi , xi−1 , xi+1 , et xi+2 ).

Réponse
−ui+2 + 6ui+1 − 3ui − 2ui−1
u0i = + O (∆x)3


6 ∆x
on peut aussi utiliser les valeurs de u aux n÷uds xi−2 , xi−1 , xi , xi+1 , et on obtient :

−ui−2 + 6ui−1 − 3ui − 2ui+1

u0i = + O (∆x)3


−6 ∆x

Détail corrigé
On fait le DL de ui+1 , ui+2 , ui−1 au voisinage de xi à l'ordre 3 :

ui+1 = u(xi+1 ) = u(xi + ∆x)

(∆x)2 00 (∆x)3 (3)
= ui + ∆x u0i + ui + O (∆x)4


ui +
2 3!
ui+2 = u(xi+2 ) = u(xi + 2∆x)
(2∆x)2 00 (2∆x)3 (3)
= ui + 2∆x u0i + ui + O (∆x)4


ui +
2 3!
ui−1 = u(xi−1 ) = u(xi − ∆x)
(∆x)2 00 (∆x)3 (3)
= ui − ∆x u0i + ui + O (∆x)4


ui −
2 3!

20
On fait une combinaison linéaire de ui+1 , ui+2 et ui−1 :

(∆x)2
αui+1 + βui+2 + γui−1 = (α + β + γ)ui + ∆x(α + 2β − γ)u0i + (α + 4β + γ)u00i (2.3)
2
(∆x)3 (3)
(α + 8β − γ)ui + O (∆x)4


+
3!
(3)
Pour obtenir u0i en fonction de ui , ui−1 , ui+1 et ui+2 , on annule on les coecients de u00i et ui .
On impose alors que :

 α + 4β + γ = 0
α + 8β − γ = 0 ⇒ 2α + 12β = 0 ⇒ α = −6β
α + 2β − γ 6= 0


On a un système de 2 équations à 3 inconnues α, β et γ. On peu choiri par exemple β = 1, et

on obtient α = −6β = −6 et γ = −α − 4β = 2.
L'égalité (2.3) s'écrit alors :

−6ui+1 + ui+2 + 2ui−1 = −3ui − 6∆x u0i + O (∆x)4

Ce qui donne :  
6ui+1 − ui+2 − 2ui−1 − 3ui
u0i + O (∆x)3


=
6 ∆x

→ Revenons maintenant à notre équation diérentielle du 1er ordre.

Nous avons vu que le problème discret s'écrit :


 Trouver ui pour i = 2, . . . , N tel que :
u0i + ui = xi − 1
u1 = 1


Si on utilise le schéma aux diérences nies décentré amont :

ui − ui−1
u0i '
∆x
alors le problème discret devient :

 ui pour i = 2, . . . , N
Trouver tel que :
ui − ui−1

+ ui = xi − 1

 ∆x
u1 = 1
c'est-à-dire : 
1
ui = (ui−1 + ∆x(xi − 1)) pour i = 2, . . . , N

1 + ∆x
 u =1
1

Connaissant u1 , on peut calculer facilement u2 , puis à partir de u2 on peut calculer u3 et ainsi

de suite. Ce schéma est itératif et est très simple à mettre en ÷uvre.

Exercice 4 Programmer cette méthode aux diérences nies sous Matlab, puis comparer la so-
lution numérique avec la solution exacte.
La solution exacte de l'équation diérentielle étant :
u(x) = 3e−x + x − 2

21
2.4 Diérences nies pour un problème elliptique 1D
2.4.1 Modèle mathématique
On considère une barre métallique de longueur 1m chauée par une source de chaleur f et
dont les deux extrémités sont plongées dans un milieu où la température est nulle (la glace par
exemple).

On note par u(x) la température en un point x.

A l'état stationnaire, u est la solution de l'EDP suivante :

−u00 (x) = f (x)


pour x ∈]0, 1[
(2.4)
u(0) = u(1) = 0 ←− conditions aux limites du type Dirichlet.

f étant une fonction donnée supposée continue sur [0, 1].

2.4.2 Approximation aux diérences nies

Etape 1 : Maillage et problème discret

x1 = 0 x2 x3 xi−1 xi xi+1 xN = 1

On subdivise le domaine de calcul [0, 1] en N points xi espacés de ∆x.

On a :

1
xi = xi−1 + ∆x = (i − 1)∆x avec ∆x =
N −1
Le problème continue (2.4) est remplacé par le problème discret suivant :


 Chercher u(xi ) pour i = 2, . . . , N − 1 vériant
−u00 (xi ) = f (xi ) (2.5)
u(x1 ) = u(xN ) = 0


Etape 2 : Schéma numérique

Donnons une approximation d'ordre 2 de u00 (xi ) :

Ecrivons les DL de ui+1 et ui−1 au voisinage de xi

(∆x)2 00 (∆x)3 (3)

ui+1 = ui + ∆x u0i + ui + O (∆x)4


ui + (A)
2! 3!
(∆x)2 00 (∆x)3 (3)
ui−1 = ui − ∆x u0i + ui + O (∆x)4


ui − (B)
2! 3!
00
= 2ui + (∆x) ui + O (∆x)4
2


(A)+(B) =⇒ ui−1 + ui+1

soit :
ui−1 − 2ui + ui+1
u00i = + O (∆x)2


(∆x) 2

22
c'est-à-dire qu'une approximation à l'ordre 2 de u00i est :

ui−1 − 2ui + ui+1

u00i '
(∆x)2

Avec ce choix d'approximation de u00i , on peut approcher le problème (2.5) par le problème discret
suivant : 

 u2 , u3 , . . . , un−1 ∈ R /
Trouver
 ui−1 − 2ui + ui+1
− = fi pour i = 2, . . . , N − 1 (2.6)

 (∆x)2
u1 = uN = 0


Etape 3 : Passage au problème matriciel

Maintenant, il sut d'écrire le problème discret (2.6) sous forme matricielle, l'inconnue étant
le vecteur U = (u2 , u3 , . . . , uN −1 ).
Mais avant, il faut tenir compte des conditions aux limites u1 = uN = 0.

En i = 2, on a :
u1 − 2u2 + u3
− = f2
(∆x)2

Puisque u1 = 0 alors  
1
− (−2u2 + u3 ) = f2
(∆x)2

En i = N − 1, on a :
1
− (uN −2 − 2uN −1 + uN ) = fN −1
(∆x)2

Or uN = 0, alors  
1
− (uN −2 − 2uN −1 ) = fN −1
(∆x)2

Et pour 3 ≤ i ≤ N − 2, on a :

 
1
− (ui−1 − 2ui + ui+1 ) = fi
(∆x)2

Le système linéaire s'écrit alors :


1

 − (−2u2 + u3 ) = f2
(∆x)2




 1
− (ui−1 − 2ui + ui+1 ) = fi pour 3≤i≤N −2

 (∆x)2

 1
 − (uN −2 − 2uN −1 ) = fN −1


(∆x)2

Sous forme matricielle

23
 −2 1 0 ··· ··· ···
    
0 u2 f2
 1 −2 1 0 ··· ··· 0   u3  f3
    
0
 1 −2 1 0 ··· 
0  
  u4 
 f4
1  .
 .. . .. . .. . .. . .. .. . 
.
.
.
  .. 
− . .
 .  =  . 
    
(∆x)2 
 .  .   . 
.. .. .. ..
 .. . . . . 0  .   . 
  .   . 
 0 ··· ··· 0 1 −2 1  uN −2  fN −2 
0 ··· ··· ··· 0 1 −2 uN −1 fN −1
Matrice tridiagonale
| {z }

⇐⇒
AU = b

Pour résoudre ce système linéaire on peut utiliser :

ˆ une méthode directe : Elimination de Gauss ou Cholesky.

ˆ une méthode itérative : Jacobi, Gauss-Seidel ou relaxation.

ˆ On peut aussi utiliser, comme la matrice est tridiagonale, l'algorithme de Thomas basé
sur la factorisation LU et qui est moins coûteux en nombre d'opérations (voir plus tard).

Exercice 5 On prend f (x) = ex − x + 1 , de sorte que la solution exacte est

x3 x2 2
u(x) = −ex + − + (e − )x + 1
6 2 3
Ecrire un programme Matlab qui calcule la solution par le schéma aux diérences nies et trace
la solution exacte et numérique sur une même gure.

2.5 Diérences nies pour un problème parabolique 1D (Equa-

tion de la chaleur)
2.5.1 Modèle mathématique
Considérons le problème de conduction de la chaleur dans une barre métallique de longueur
1 m, dont les extrémités sont soumises à des températures Tg et Td . On suppose qu'à l'instant
initial, la température le long de la barre vaut f (x).

Notons par u(x, t) la température au point x de la barre à l'instant t. u vérie l'équation de la

chaleur suivante avec conditions initiales (CI) et conditions aux limites (CL) :


∂u ∂2u
−D 2 =0 ∀ x ∈ ]0, 1[ ; ∀ t ≥ 0



∂t ∂x (2.7)

 u(x, 0) = f (x) ← donnée ∀ x ∈]0, 1[
u(0, t) = Tg ; u(1, t) = Td ∀t ≥ 0


D étant le coecient de diusion.

24
2.5.2 Approximation diérences nies
Etape 1 : Maillage (Discrétisation spatiale et temporelle)

a = x1 x2 x3 xi−1 xi xi+1 xN = b

ˆ [0, 1] en N points (x1 , x2 , . . . , xN ) espacés de ∆x.

On discrétise l'intervalle
1
∆x = , et xi = xi−1 + ∆x = (i − 1)∆x
N −1
ˆ Le temps t est aussi discrétisé en des intervalles ∆t (t0 , t1 , . . . , tn , . . .)
avec t
n+1 = tn + ∆t ∆t étant le pas du temps.
n
On note par ui la température au n÷ud xi à l'instant t = n∆t.
n

Etape 2 : Diérences nies (Approximation des dérivées)

Dérivée spatiale
∂2u
Une approximation à l'ordre 2 centrée de au n÷ud xi s'écrit :
∂x2
∂2u ui−1 − 2ui + ui+1
'
∂x2 i (∆x)2
Donc le problème approché s'écrit au point xi :

∂u ui−1 − 2ui + ui+1

= D
∂t i (∆x)2
du
⇔ = f (u)
dt i

Dérivée temporelle :
Pour la discrétisation temporelle, il y a deux manières de procéder : soit utiliser un schéma
explicite, soit un schéma implicite.
∂u
Tout d'abord, une approximation à l'ordre 1 en temps de s'écrit :
∂t i

∂u un+1
i − uni
'
∂t i ∆t

Schéma d'Euler explicite

Dans ce cas, les termes de l'opérateur spatiale sont évalués à l'instant tn . On a donc :

∂2u uni−1 − 2uni + uni+1

'
∂x2 i (∆x)2

Le problème discret s'écrit alors :

 n+1
u − uni un − 2uni + uni+1
 i − D i−1 =0 pour i = 2, . . . , N − 1


∆t (∆x)2
(2.8)
u0 = f (xi ) pour i = 2, . . . , N − 1
 in


u1 = Tg ; unN = Td ∀n ∈ N

25
On obtient alors le schéma :

 
D∆t n
un+1 = uni + (u − 2uni + uni+1 ) ∀i = 2, . . . , N − 1 (2.9)
i
(∆x)2 i−1

Comme on connait u0i i = 1, . . . , N ), alors on peut calculer u1i à l'instant t1 = ∆t par le

( pour
2
procédé itératif (2.9), puis ui et ainsi de suite.
n+1 n+1 = (n + 1)∆t se déduisent directe-
En fait le schéma est dit explicite car les ui à l'instant t
n n
ment des ui à l'instant t = n∆t qui sont connues. Il n'y a aucun système à résoudre.

Remarque 1:
Le schéma (2.9) est très simple à mettre en oeuvre, cependant on montrera par la suite que pour
être stable, il nécessite une restriction sur le pas du temps ∆t. Ce dernier doit être choisi de telle
sorte que la condition suivante soit vériée :

2D
∆t ≤1
(∆x)2

Schéma d'Euler implicite

Dans ce cas, l'opérateur spatial est évalué à l'instant tn+1 . Et on a :

∂2u un+1
i−1 − 2un+1
i + un+1
i+1
'
∂x2 i (∆x)2

Le problème discret s'écrit alors :


 un+1
i − uni un+1 n+1
i−1 − 2ui + un+1
i+1
−D =0 pour i = 2, . . . , N − 1


(∆x)2

∆t (2.10)
u0 = f (xi ) pour i = 2, . . . , N − 1
 uin = T ; un = T



1 g N d ∀n ∈ N

Soit

−λun+1 n+1
− λun+1 n
 
i−1 + (1 + 2λ)ui i+1 = ui pour i = 2, . . . , N − 1 (2.11)

D∆t
avec λ= .
(∆x)2

Le schéma aux diérences nies implicite (2.11) peut être écrit sous la forme matricielle suivante :

  
 un+1 
−λ ··· ··· ··· un2 + λTg

1 + 2λ 0 0 2
 n+1  
 −λ 1 + 2λ −λ 0 ··· ··· 0  u un3 
  3n+1  
 
un4

 0
 −λ 1 + 2λ −λ 0 · · · 0   u4  
 
 


 .. .. .. .. .. . . . 
.   ..   .
  
 . . . . . . .

 
  ..  = 
   
 . . . . . .

 .. .. .. .. .. .
0  .   .

    
. . .  .   .

.. .. ..   ..   .

 0 0 .

    
n+1 un
 0 ··· ··· 0 −λ 1 + 2λ −λ

uN −2 
    N −2 
0 ··· ··· ··· 0 −λ 1 + 2λ n
uN −1 + λTd
un+1 N −1

26
On voit donc qu'à chaque instant, c'est-à-dire qu'à chaque passage de tn à tn+1 , ce schéma
nécessite la résolution d'un système linéaire dont la matrice est tridiagonale.

Néanmoins, ce schéma possède l'avantage d'être inconditionnellement stable c'est-à-dire stable

quelque soit le choix du pas du temps ∆t (comme on le verra par la suite). L'un des avantages
des schémas implicites est qu'ils se permettent l'utilisation de grands pas de temps et donc sont
bien adaptés pour la simulation des problèmes longs en temps physique.

Remarque 2:
Pour la résolution du système linéaire avec matrice tridiagonale, on peut utiliser l'algorithme de
Thomas, qui eectue la factorisation LU de Gauss et qui nécessite un coût de calcul de 8N − 7
opérations (3(N − 1) opérations pour la factorisation et 5N − 4 opérations pour la substitution).

Algorithme de Thomas
 
a1 c1 0 · · · · · · 0
 b2 a2 c2 0 · · · 0 
 
 .. .. .. 
A=0
 . . . 

.
 ..

cn−1 
0 · · · · · · 0 bn an

On décompose : A = LU avec :

   
1 0 ··· 0 α1 c1 · · · 0
. .
. .
  
β2 .  .
   
.. . .. .. . 
L=0 . ; U =0 . 
 
. . . . .
   
. .
 ..  ..
 
 cn−1 
0 ··· 0 βn 1 0 ··· αn

On aura alors :

α1 = a1
b2 b2
β2 α1 = b2 =⇒ β2 = =
α1 a1
b2
β2 c1 + α2 = a2 =⇒ α2 = a2 − β2 c1 = a2 − c1
a1
.
.
.

2.6 Consistance, stabilité et convergence des schémas aux dié-

rences nies
2.6.1 Erreur de troncature
Dénition 2.6.1 On appelle erreur de troncature d'un schéma numérique la quantité, notée ET ,
obtenue par la diérence entre l'équation exacte et l'équation approchée par le schéma numérique.

27
Autrement dit, si l'équation s'écrit :
L(u) = 0
où L est un opérateur diérentiel, et l'équation approchée s'écrit :

Lh (u) = 0

alors
ET = L(u) − Lh (u)

Exemple 2.6.1 Dans le cas de l'équation de la chaleur :

∂u ∂2u
−D 2 =0 ⇐⇒ L(u) = 0 (2.12)
∂t ∂x
avec le schéma numérique :
un+1 − uni un − 2uni + uni+1
i
− D i−1 =0 ⇐⇒ Lh (u) = 0 (2.13)
∆t (∆x)2

l'erreur de troncature s'écrit :

∂2u un+1 − uni un − 2uni + uni+1
   
∂u i
− D i−1 = O(∆t) + O (∆x)2


ET = −D 2 −
∂t ∂x ∆t (∆x)2

On voit que l'erreur de troncature en temps est en O(∆t) et celle en espace est en O (∆x)2 .

Nous avons donc un schéma d'ordre 1 en temps et d'ordre 2 en espace.

2.6.2 Consistance
En fait, pour une même équation diérentielle ou aux dérivées partielles, il existe une multi-
tude de schémas numériques permettant de résoudre cette équation. Il faut donc assurer que le
schéma représente bien l'équation exacte dans la limite où le pas d'espace ∆x et la pas de temps
∆t tendent vers zero, et partout dans le domaine considéré. On parle de consistance du schéma.

Dénition 2.6.2 Un schéma aux diérences nies est dit consistant si et seulement si :
lim ET = 0
∆x→0
∆t→0

Exemple 2.6.2 Le schéma (2.13) pour l'équation de la chaleur (2.12) est consistant. En eet :

lim ET = lim O(∆t) + O (∆x)2

∆t→0 ∆t→0
∆x→0 ∆x→0
or
|O(∆t)| ≤ λ1 |∆t| −→ 0 =⇒ lim O(∆t) = 0
∆t→0 ∆t→0
et
O( ∆x)2 ≤ λ2 (∆x)2 −→ 0 =⇒ lim O (∆x)2 = 0



∆x→0 ∆x→0

28
2.6.3 Stabilité d'un schéma numérique
Dénition 2.6.3 Un schéma numérique est dit stable si la solution numérique calculée par ce
schéma reste bornée dans le temps et dans l'espace.
Autrement dit, les erreurs produites par le schéma (Erreur de troncature, erreurs d'arrondi, ...)
ne doivent pas croître lors de la procédure numérique.

La stabilité du schéma numérique est parfois attachée aux valeurs du pas du temps ∆t et du pas
d'espace ∆x.
En fait, trois cas sont possibles pour un schéma numérique :

ˆ Il peut être inconditionnellement stable, c'est-à-dire stable quelque soit le choix de ∆t et

∆x (généralement c'est le cas des schémas implicites).

ˆ Conditionnellement stable (cas de schémas explicites).

ˆ Il peut être instable.

2.6.4 Convergence
Dénition 2.6.4 Un schéma numérique sera dit convergent pour une norme donné (notée k.k),
si pour toute solution initiale u0 et pour tout n ∈ N, on a :
lim kunh − unex k = 0
∆x→0
∆t→0

où unex est la solution exacte au temps tn et unh la solution numérique.

2.6.5 Théorème de Lax

En général, il n'est pas facile de montrer la convergence d'une méthode numérique en calculant
lim kunh − unex k du fait qu'on ne connait pas toujours unex . La consistance et la stabilité d'un
∆x→0
∆t→0
schéma sont en général beaucoup plus facile à étudier que sa convergence.

On utilise le théorème de Lax (Richtmyer et Norton 1967) suivant :

Théorème de Lax
Pour un problème linéaire aux valeurs initiales, la solution numérique d'un schéma itératif en
temps aux diérences nies converge vers la solution exacte si le schéma est consistant et stable.

2.6.6 Stabilité au sens de Fourier - Von Neumann

Dénition 2.6.5 Un schéma numérique est dit L2 -stable si pour tout temps T > 0, il existe une
constante C(T ) indépendante de ∆x et ∆t telle que :
Pour toute donnée initiale u0 ∈ L2 ; et ∀n ∈ N avec n∆t ≤ T , on a :
kunh kL2 ≤ C(T )ku0 kL2

où unh est la solution approchée à l'instant tn = n∆t.

→ Pour l'étude de la stabilité au sens L2 , on utilise la méthode de Fourier - Von Neumann.

La méthode consiste à chercher des solutions de l'EDP sous la forme :

29
 unj = C n eiξxj avec xj = j∆x et ξ est le nombre d'onde.

puis injecter cette solution dans le schéma numérique, et étudier le comportement de la suite de
fonction (C n )n∈N .
On obtient une relation de récurrence sur les C n. Et on a :
n
si la suite de fonctions (C )n∈N est bornée pour la norme L2 alors le schéma numérique est stable,
sinon il est instable.

Exemple 2.6.3 Considérons l'équation de la chaleur :

∂u ∂2u
−D 2 =0
∂t ∂x
et soit le schéma aux diérences nies explicite suivant :
un+1
j − unj unj−1 − 2unj + unj+1
−D =0
∆t (∆x)2

qui peut s'écrire aussi :

∆tD
un+1
j = unj + λ(unj−1 − 2unj + unj+1 ) avec λ=
(∆x)2

Étudions la stabilité au sons de Fourier - Von Neumann :

on pose :
unj = C n eiξj∆x
on injecte cette solution dans le schéma, on aura :

C n+1 eiξj∆x = C n eiξj∆x + λ(C n eiξ(j−1)∆x − 2C n eiξj∆x + C n eiξ(j+1)∆x )

C n+1 = [1 − 2λ + λ(e−iξ∆x + eiξ∆x )] C n

= [1 − 2λ + 2λ cos(ξ∆x)] C n
= [1 − 2λ(1 − cos(ξ∆x))] C n or 1 − cos(2x) = 2 sin2 (x)
ξ∆x
C n+1 = [1 − 4λ sin2 ( )] C n
2
posons :

ξ∆x
A = 1 − 4λ sin2 ( ) A est appelé facteur d'amplication.
2
on a alors :
C n+1 = AC n
c'est-à-dire :
C n = AC n−1 = A2 C n−2 = . . . = An C 0
En démarrant par une condition initiale u0j = C 0 eiξj∆x bornée
0
c'est-à-dire telle que |C | ≤ M , pour que la solution unj à l'instant tn = n∆t soit bornée ∀n ∈ N,
il faut que |A| ≤ 1 ; ∀ξ ∈ R

30
 c'est-à-dire :

ξ∆x
−1 ≤ 1 − 4λ sin2 ( )≤1 ∀ξ ∈ R
2
or on a toujours

ξ∆x ξ∆x
1 − 4λ sin2 ( )≤1 car 4λ sin2 ( ≥0
2 2
Il reste :

ξ∆x
−1 ≤ 1 − 4λ sin2 ( )
2
ξ∆x
⇔ −2 ≤ −4λ sin2 ( )
2
ξ∆x
⇔ 2λ sin2 ( ) ≤ 1 ∀ξ ∈ R
2
ξ∆x
⇔ sup(2λ sin2 ( )) ≤ 1
ξ∈R 2
2∆tD
⇔ 2λ ≤ 1 ⇐⇒ ≤1
(∆x)2

Le schéma proposé est donc stable si :

∆tD 1
≤
(∆x)2 2

Exemple 2.6.4 Considérons cette fois-ci un schéma implicite pour l'équation de la chaleur :
un+1
j − unj un+1 n+1
j−1 − 2uj + un+1
j+1
−D =0
∆t (∆x)2

Regardons ce qu'il en est par la stabilité L2 .

Le schéma s'écrit :
∆tD
un+1
j = unj + λ(un+1 n+1
j−1 − 2uj + un+1
j+1 ) avec λ=
(∆x)2
on prend :
unj = C n eiξj∆x
Le schéma s'écrit :
C n+1 eiξj∆x = C n eiξj∆x + λ(C n+1 eiξ(j−1)∆x − 2C n+1 eiξj∆x + C n+1 eiξ(j+1)∆x )

(1 + 2λ) C n+1 − λ(e−iξ∆x + eiξ∆x ) C n+1 = C n

1
C n+1 = Cn
(1 + 2λ) − λ(e−iξ∆x + eiξ∆x )
1
C n+1 = Cn
ξ∆x
1 + 4λ sin2 ( )
2
= A Cn

31
or
1
∀ξ ∈ R ; |A| = ≤1
ξ∆x2
1 + 4λ sin ( )
2
Il s'en suit que le schéma implicite est stable dans L2 (R).
On dit qu'il est inconditionnellement stable (c'est-à-dire quelques soient les paramètres ∆t et
∆x).

Exemple 2.6.5 On considère l'équation de transport avec la CI suivante :


 ∂u ∂u
+c = 0 ∀x ∈ R ; ∀t > 0
∂t ∂x
 u(x, 0) = u (x) ∀x ∈ R
0

On suppose que la vitesse c>0.

Étudions la stabilité L2 du schéma décentré amont suivant :
un+1
j − unj unj − unj−1
+c =0
∆t ∆x
qui s'écrit aussi :
c∆t
un+1
j = unj − λ(unj − unj−1 ) avec λ=
∆x
on prend :
unj = C n eiξj∆x
Injectons dans le schéma :
C n+1 eiξj∆x = C n eiξj∆x − λ(C n eiξj∆x − C n eiξ(j−1)∆x )

| − λ +{zλe
−iξ∆x
C n+1 = (1 }) C
n

A
On doit avoir :
|A| ≤ 1 ∀ξ ∈ R

|A|2 = |1 − λ + λ cos(ξ∆x) − iλ sin(ξ∆x)|2

= (1 − λ + λ cos(ξ∆x))2 + (λ sin(ξ∆x))2
= 1 − 2λ + 2λ2 + 2λ cos(ξ∆x) − 2λ2 cos(ξ∆x)

|A|2 ≤ 1 ⇔ 1 + 2λ2 (1 − cos(ξ∆x)) − 2λ(1 − cos(ξ∆x)) ≤ 1 ∀ξ ∈ R

⇔ 2λ(1 − cos(ξ∆x))(λ − 1) ≤ 0

or
∀ξ ∈ R ; 2λ(1 − cos(ξ∆x)) ≥ 0 car λ>0
On doit donc avoir :
λ−1≤0
c'est-à-dire :
c∆t
λ= ≤1
∆x

32
Exercice 6 On reprend l'équation d'advection précédente (avec c > 0).

 ∂u ∂u
+c = 0 ∀x ∈ R ; ∀t > 0
∂t ∂x
 u(x, 0) = u0 (x) ∀x ∈ R
Étudier l'ordre, la consistance et la stabilité de Fourier - Von Neumann des schémas suivants :

Schéma 1 (décentré aval) :

un+1
j − unj unj+1 − unj
+c =0
∆t ∆x
Schéma 2 (centré) :
un+1
j − unj unj+1 − unj−1
+c =0
∆t 2∆x
Schéma 3 (Leap-frog ou saute-mouton) :

un+1
j − ujn−1 unj+1 − unj−1
+c =0
2∆t 2∆x

Schéma 1 (décentré aval)

un+1
j − unj unj+1 − unj
+c =0
∆t ∆x

Ordre du schéma et consistance

n
∂u
ˆ un+1
j = u(xj , tn + ∆t) = unj + ∆t + O((∆t)2 )
∂t j

∂u n un+1
j − unj
=⇒ = + O(∆t)
∂t j ∆t
n
∂u
ˆ unj+1 = u(xj + ∆x, tn ) = unj + ∆x + O((∆x)2 )
∂x j

∂u n unj+1 − unj
=⇒ = + O(∆x)
∂x j ∆x

On déduit l'erreur de troncature :

∂u ∂u ujn+1 − unj unj+1 − unj

ET = ( +c )−( +c ) = O(∆t) + O(∆x)
∂t ∂x ∆t ∆x
Le schéma (1) est donc d'ordre 1 en temps et en espace.

Il est consistant. En eet,

|O(∆t)| ≤ λ1 |∆t| −→ 0 et |O(∆x)| ≤ λ2 |∆x| −→ 0

∆t→0 ∆x→0

c'est-à-dire
lim ET = 0
∆t→0
∆x→0

33
Stabilité
Le schéma (1) s'écrit :

c∆t
un+1
j = unj − λ(unj+1 − unj ) avec λ=
∆x
on pose
unj = C n eiξj∆x
on injecte cette solution dans le schéma :

C n+1 eiξj∆x = C n eiξj∆x − λ(C n eiξ(j+1)∆x − C n eiξj∆x )

C n+1 = [1
|+λ−
iξ∆x
{zλe }] C
n

Le schéma (1) est stable si et seulement si

|A| ≤ 1 ∀ξ ∈ R

|A|2 ≤ 1 ⇔ |(1 + λ) − λ(cos(ξ∆x) + i sin(ξ∆x))|2 ≤ 1

⇔ ((1 + λ) − λ cos(ξ∆x))2 + (λ sin(ξ∆x))2 ≤ 1
⇔ (1 + λ)2 − 2λ(1 + λ) cos(ξ∆x) + λ2 cos2 (ξ∆x) + λ2 sin2 (ξ∆x) ≤ 1
⇔ 1 + 2λ + λ2 − 2λ(1 + λ) cos(ξ∆x) + λ2 ≤ 1
⇔ 2λ(1 + λ)(1 − cos(ξ∆x))∀ξ∈R ≤ 0
c∆t
Impossible car λ= ≥0 et 1 − cos(ξ∆x) ≥ 0
∆x
Le schéma (1) est donc instable.

Schéma 2 (centré)

un+1
j − unj unj+1 − unj−1
+c =0
∆t 2∆x

Ordre et consistance

∂u n un+1
j − unj
= + O(∆t)
∂t j ∆t

n n
∂u (∆x)2 ∂ 2 u
unj+1 = unj + ∆x + O (∆x)3


+ (2.14)
∂x j 2 ∂x2 j
n n
∂u (∆x)2 ∂ 2 u
unj−1 = unj − ∆x + O (∆x)3


+ (2.15)
∂x j 2 ∂x2 j

∂u n unj+1 − unj−1
+ O (∆x)2


(2.14)-(2.15) =⇒ =
∂x j 2∆x

= O(∆t) + O (∆x)2


ET
|ET | ≤ |O(∆t)| + |O(∆x)2 | ≤ λ1 ∆t + λ2 (∆x)2 −→ 0
∆t/∆x→0

Le schéma (2) est donc consistant, d'ordre 1 en temps et 2 en espace.

34
Stabilité
Le schéma s'écrit :

c∆t
un+1
j = unj − λ(unj+1 − unj−1 ) avec λ=
2∆x
On pose :
unj = C n eiξj∆x
On l'injecte dans le schéma

h i
C n+1 eiξj∆x = iξj∆x
e − λ(eiξ(j+1)∆x − eiξ(j−1)∆x ) C n

h i
C n+1 = 1 − λ(eiξ∆x − e−iξ∆x ) C n

C n+1 = [1 − 2i sin(ξ∆x)] C n
| {z }
A

Or
|A|2 = |1 − 2iλ sin(ξ∆x)|2
= 1 + 4λ2 sin2 (ξ∆x) ≥ 1 ∀ξ ∈ R
Donc le schéma (2) est instable.

Schéma 3 (Leap-frog)

un+1
j − ujn−1 unj+1 − unj−1
+c =0
2∆t 2∆x

n n
∂u (∆t)2 ∂ 2 u
un+1 unj + O (∆t)3


j = + ∆t + (2.16)
∂t j 2 ∂t2 j
n n
∂u (∆t)2 ∂ 2 u
un−1 = unj − ∆t + O (∆t)3


j + (2.17)
∂t j 2 ∂t2 j

∂u n un+1
j − ujn−1
+ O (∆t)2


(2.16)-(2.17) =⇒ =
∂t j 2∆t

De même on a :
∂u n unj+1 − unj−1
+ O (∆x)2


=
∂x j 2∆x

= O((∆t)2 ) + O (∆x)2


ET

=⇒ Schéma consistant d'ordre 2 en temps et en espace.

35
Stabilité
Le schéma s'écrit :
c∆t
un+1
j = un−1
j − λ(unj+1 − unj−1 ) avec λ=
∆x
On pose :
unj = C n eiξj∆x
on l'injecte dans le schéma

C n+1 eiξj∆x = C n−1 eiξj∆x − λ(C n eiξ(j+1)∆x − C n eiξ(j−1)∆x )

C n+1 = C n−1 − 2iλ sin(ξ∆x) C n

C n+1 + 2iλ sin(ξ∆x) C n − C n−1 = 0

L'équation caractéristique de cette suite récurrente est :

r2 + 2iλ sin(ξ∆x) r − 1 = 0 (2.18)

Le terme général de la suite (C n )n∈N s'écrit alors :

C n = αr1n + βr2n
où r1 et r2 sont les racines de l'équation caractéristique (2.18)
( p
r1 = −iλ sin(ξ∆x) − 1 − λ2 sin2 (ξ∆x)
p
r2 = −iλ sin(ξ∆x) + 1 − λ2 sin2 (ξ∆x)
Pour avoir la stabilité c'est-à-dire la suite (C n ) bornée ∀n ∈ N, il faut et il sut que les racines
de l'équation caractéristique (2.18) soient toutes les deux de module inférieur ou égale à 1.

Si λ2 sin2 (ξ∆x) ≤ 1 alors

1 − λ2 sin2 (ξ∆x) ≥ 0
Donc :
|r1 |2 = |r2 |2 = (−λ sin(ξ∆x))2 + (1 − λ2 sin2 (ξ∆x)) = 1
Le schéma est donc stable.

Si λ2 sin2 (ξ∆x) > 1 alors

q
r1 = −iλ sin(ξ∆x) − i λ2 sin2 (ξ∆x) − 1
q
r2 = −iλ sin(ξ∆x) + i λ2 sin2 (ξ∆x) − 1
Dans ce cas, les deux racines sont imaginaires purs. De plus, leur produit r1 · r2 = −1. Donc
l'une des deux racines a un module > 1.
Le schéma est donc instable.

La stabilité est donc vériée si et seulement si ∀ξ ∈ R

λ2 sin2 (ξ∆x) ≤ 1

⇔ sup(λ2 sin2 (ξ∆x)) ≤ 1

ξ∈R

c∆t
⇔ λ2 ≤ 1 ⇔ λ= ≤1
∆x

36