Université Mohammed Premier

Ecole Nationale des Sciences Appliquées

Oujda

Modélisation Mathématique

Filière Génie Civil

Première année

2019-2020

Imad EL MAHI
Table des matières

2 Méthode des diérences nies 3

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

2
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

3
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

4
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 :
 Chercher u(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

5
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 1 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!

6
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!
Pour obtenir u0i en fonction de ui , ui−1 , ui+1 et ui+2 , on annule on les coecients de u00i et u(3)
i .
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 :


u0 + ui = xi − 1
 i
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 :
Trouver ui pour i = 2, . . . , N 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 2 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

7
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 1 m 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 (A)


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


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

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


(∆x) 2

8
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 :
 Trouver u2 , u3 , . . . , un−1 ∈ R /


 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

9
 −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).

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
u(x, 0) =
∂x
f (x) ← donnée ∀ x ∈]0, 1[
(2.7)


u(0, t) = Tg ; u(1, t) = Td ∀t ≥ 0


D étant le coecient de diusion.

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

• On discrétise l'intervalle [0, 1] en N points (x1 , x2 , . . . , xN ) espacés de ∆x.

10
 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 tn+1 = tn + ∆t ∆t étant le pas du temps.
On note par ui la température au n÷ud xi à l'instant tn = 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 )
 in

 pour i = 2, . . . , N − 1
u1 = Tg ; unN = Td ∀n ∈ N

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

11
Comme on connait u0i ( pour i = 1, . . . , N ), alors on peut calculer u1i à l'instant t1 = ∆t par le
procédé itératif (2.9), puis u2i et ainsi de suite.
En fait le schéma est dit explicite car les un+1
i à l'instant tn+1 = (n + 1)∆t se déduisent directe-
ment des ui à l'instant t = n∆t qui sont connues. Il n'y a aucun système à résoudre.
n n

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 n+1
i−1 − 2ui + un+1
i+1
'
∂x2 i (∆x)2

Le problème discret s'écrit alors :

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

Soit

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

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
  un+1 un3
  
 −λ 1 + 2λ −λ 0 · · · · · · 0 3  
  n+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+1 n
uN −1 + λTd
u N −1

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.

12
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 
. . .
 
 . . . . . .

0
A= 
.

. . 
cn−1 
0 · · · · · · 0 bn an

On décompose : A = L U 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.
Autrement dit, si l'équation s'écrit :
L(u) = 0

13
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

14
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 :

15
 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 :
si la suite de fonctions (C n )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
c'est-à-dire telle que |C 0 | ≤ M , pour que la solution unj à l'instant tn = n∆t soit bornée ∀n ∈ N,
il faut que |A| ≤ 1 ; ∀ξ ∈ R

16
 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

17
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

18
Exercice 3 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

19
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 {zλe }] C
iξ∆x n
|+λ−
A

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
(2.14)-(2.15) =⇒ + O (∆x)2


=
∂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.

20
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 = eiξj∆x − λ(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 (2.16)


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


j +
∂t j 2 ∂t2 j

∂u n un+1 − ujn−1
(2.16)-(2.17) =⇒ j
+ O (∆t)2


=
∂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.

21
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

22