CHAPITRE 1 METHODE DES DIFFERENCES FINIES

CHAPITRE 1

1 - Classifications des équations aux dérivées partielles (EDP) :

1.1 - Définitions :
Une équation aux dérivées partielles (E.D.P) est une relation faisant intervenir
les variables indépendantes x1,……,xn, la fonction u et ses dérivées partielles. Par
exemple, si u est une fonction de deux variables, une E.D.P peut s'écrire par la
relation :

f (x, y, ux , uy , uxx, uyy, uxy, ux2 y , uxy2 , ux3 , uy3 ,......)  0

* On appelle ordre de l'E.D.P l'ordre le plus élevé des dérivées partielles

intervenant dans l'E.D.P, par exemple :

ux2 y  3uxx  x u yy  ux  u  c  0 est d'ordre 3

(uxx  u yy )2  4( uxy )2  c  0 est d'ordre 2

* L'E.D.P est dite linéaire si f est linéaire par rapport à ses arguments u et ses
dérivées partielles, et si les coefficients qui les
lient ne dépendent que de (x,y); sinon elle est non linéaire. Par exemple, l'E.D.P
du second ordre :

al uxx + a2 uxy + a3 uyy + a4 ux + a5 uy + a6 u + a7 = 0 I-1

est linéaire si les ai ne dépendent que de (x,y).

Adapted by A. Boukhari 1
El-Oued University 2013.
CHAPITRE 1 METHODE DES DIFFERENCES FINIES

1.2 - Classification mathématique des E.D.P linéaires du 2nd ordre (cas de

deux variables indépendantes) :
De très nombreux phénomènes physiques se traduisent par des E.D.P linéaires
du second ordre du type (I.1) qui peuvent s'écrire sous la forme :
a uxx + 2 b uxy + c uyy = d I-2
où a, b et c ne dépendent que de (x,y) et d est une fonction linéaire de (x, y, u,
ux,uy )

Il y a trois types d'équations aux dérivées partielles représentés par

l’équation (I.2). Lorsque la quantité  (= b2 - 4 ac) < 0 l'équation (I.2) est dite du
type elliptique. Lorsque   0 elle est dite du type parabolique. Lorsque   0 , elle
est dite du type hyperbolique. Cette appellation est faite par analogie avec
l'équation générale du second ordre en géométrie analytique :
a x2 + 2 b xy + c y2 = d I-3

Ainsi, selon le signe du discriminant b2 - 4 ac, nous obtenons différentes formes

géométriques :

Adapted by A. Boukhari 2
El-Oued University 2013.
CHAPITRE 1 METHODE DES DIFFERENCES FINIES

1.3 - Classification mathématique dans le cas général (n variables

indépendantes) :

Si u est une fonction de n variables indépendantes, les E.D.P. linéaires du

second ordre sont du type

* Si tous les ai sont non nuls et de même signe, l'E.D.P est de type elliptique.
* Si tous les ai sont non nuls et sont; à une exception près, de même signe, l'E.D.P.
est de type hyperbolique.
* Si un seul des ai est nul (noté ai0) et tous les autres de même signe et si bi0 est non
nul l'E.D.P. est de type parabolique.

Les fonctions ai et bi étant dépendantes des variables (xl ,…, xn), la

classification est évidemment fonction du point (xl ,…, xn) considéré. Une E.D.P.
peut donc être de différents types suivant les points considérés : on dit qu'elle est de
type mixte.

Exemples :
Soient U(x,y) une fonction de deux variables et V(x,y,t) une fonction de
trois variables.

Adapted by A. Boukhari 3
El-Oued University 2013.
CHAPITRE 1 METHODE DES DIFFERENCES FINIES

 2U  2U
 2 0 est une E.D.P elliptique
x 2 y
 2V  2V  2V est une E.D.P hyperbolique
 2  2
t 2 x y

V  2V  2V
 2  2 est une E.D.P parabolique
t x y

 2U  2U
x 2  0 est une E.D.P Elliptique pour x0
x y 2 hyperbolique pour x 0
parabolique pour x0

2 - Classification physique des E.D.P :

De nombreux phénomènes physiques se rangent dans l'une des classes

suivantes :

 Les problèmes d'équilibre étudient l'état stationnaire d'un phénomène

(champ, chaleur, etc...) dans un domaine borné ou non. Ils sont gouvernés

par des E.D.P elliptiques.

 Les problèmes de valeurs propres sont en général des extensions des

problèmes d'équilibre dans lesquels les valeurs critiques de certains
paramètres doivent être déterminées. C'est le cas par exemple de la
résonance des circuits électriques.
Les problèmes d'évolution étudient l'évolution avec le temps d'un phénomène
(champ, chaleur, vibration, etc...) à partir d'un état initial donné. Ils sont gouvernés
par des E.D.P hyperboliques ou des E.D.P paraboliques.

Adapted by A. Boukhari 4
El-Oued University 2013.
CHAPITRE 1 METHODE DES DIFFERENCES FINIES

Exemples :

a) Equation de la chaleur :
La conduction thermique à l'intérieur d'un domaine D bidimensionnel
provoque un changement de la température T (t,x,y), qui est régi, en l'absence de
source de chaleur par l'E.D.P :

où k,  , c sont respectivement la conductivité thermique, la masse

volumique et la chaleur spécifique du solide constituant le domaine D.
Lorsque k dépend seulement de la position (x,y), l'E.D.P est linéaire; si k
dépend de la température T, l'E.D.P est non linéaire.
Dans la majorité des cas rencontrés, on considère k comme constante et alors
l'équation précédente peut être réécrite sous la forme :

k
 est le coefficient de diffusion thermique
c
 2T  2T désigne le Laplacien de T
T  2  2
x y

Adapted by A. Boukhari 5
El-Oued University 2013.
CHAPITRE 1 METHODE DES DIFFERENCES FINIES

b) Equation des ondes :

Considérons une membrane mince occupant au repos un domaine D
bidimensionnel de frontière F. Si on suppose :
 que la membrane est homogène et parfaitement élastique
 qu'elle est soumise à une tension T constante
 que les mouvements latéraux (dans le plan de la membrane) sont
négligeables
 que les déplacements verticaux (perpendiculaires au plan de la
membrane) z (t,x,y) sont de petite amplitude.

Alors les déplacements z(t , x , y) sont gouvernés par l'E.D.P

où   masse volumique de la membrane.

f = force extérieure à laquelle elle est soumise, par exemple la
gravitation.
Ici aussi, les conditions à la frontière F peuvent être de différents types
L'équation des ondes intervient aussi dans de nombreux phénomènes autres que le
mouvement vertical d'une membrane horizontale : courant électrique dans un
circuit, ondes acoustiques ou champs électriques dans les des d'ondes, etc...

Adapted by A. Boukhari 6
El-Oued University 2013.
CHAPITRE 1 METHODE DES DIFFERENCES FINIES

c) Equations de Laplace et de Poisson :

Plusieurs phénomènes sont régis par des E.D.P de la forme suivante :

I-4

Appelée équation de Laplace. Par exemple, en régime stationnaire de la conduction

de la chaleur, la température T(x , y) est régie par l'E.D.P

I-5

où k est la conductivité thermique et f une fonction donnant l'apport des

sources / puits de chaleur extérieurs.
Lorsque k ne dépend pas de la position (x , y), l'équation (I-5) devient :

Il s'agit de l'équation de Poisson ( si A=0).

L'équation de Poisson est rencontrée aussi en électrostatique sous la forme

 est le potentiel électrique et  la densité de charge.

La même équation est utilisée aussi en élasticité pour calculer la fonction 

dont on peut avoir l'angle de torsion d'un cylindre soumis à des forces de rotation
opposées.
d) Equation du téléphone :
La variation du voltage V le long d'un fil de transmission peut être prédite par
l'E.D.P

où R, L désignent la résistance et l'inductance par unité de longueur. c et G

désignent la capacité et la conductivité par unité de longueur.
Adapted by A. Boukhari 7
El-Oued University 2013.
CHAPITRE 1 METHODE DES DIFFERENCES FINIES

3 - Types de conditions aux limites :

D'après les exemples donnés ci-dessus nous pouvons remarquer que les
conditions aux limites sont de trois types :

1. Lorsqu'on précise la valeur de la solution u du problème sur la frontière F du

domaine D

la condition aux limites est dite du type Dirichlet.

2. Lorsqu'on impose le flux de la solution u à travers la frontière :

la condition aux limites est alors dite du type Neumann.

3. Si on impose une relation entre la valeur de u sur la frontière et le flux à

travers celle-ci :

la condition aux limites est dite de type mixte ou de 3ème type.

Adapted by A. Boukhari 8
El-Oued University 2013.
CHAPITRE 1 METHODE DES DIFFERENCES FINIES

4 - Discrétisation du domaine :
Il est bien connu qu’un problème aux dérivées partielles nécessite la donnée :
 d'un domaine D
 d'une équation aux dérivées partielles (E.D.P)
 de conditions aux limites
 de conditions initiales (pour les problèmes d'évolution)
Pour obtenir une approximation numérique de la solution de ce problème
nous devons approcher les dérivées partielles de l'E.D.P en chaque noeud du
domaine discrétisé (maillage) en utilisant les valeurs de la variable dépendante en
ce noeud et aux noeuds avoisinants.
Les calculs par différences finies sont effectués suivant un maillage obtenu
par un double réseau de parallèles aux axes et régulièrement espacés. L'intersection
de 2 droites du maillage définit un noeud M de cordonnées (xM , yM). Si les
parallèles à l'axe x sont espacées de Δ x  h et les parallèles à l'axe y de Δ y  k , le
noeud a comme coordonnées x  i Δ x  i h , et y  j Δ y  j k ou d'une manière
M M
condensée (i , j).
Ainsi la fonction U(x,y) prend au point M(xM , yM) la valeur U( i Δ x , j Δ y )

=U (i h , j k) ou U i , j ou U ij .

Fig.I-1 : Maillage

Adapted by A. Boukhari 9
El-Oued University 2013.
CHAPITRE 1 METHODE DES DIFFERENCES FINIES

5 - Approximation des dérivées :

Soit U(x,y) une fonction de deux variables indépendantes que nous
supposerons suffisamment différentiable. Si nous écrivons son développement en
séries de Taylor en un point ( x+h , y+k) , nous avons :

avec le résidu Rn donné par :

ou encore

Cette equation signifie qu'il existe un nombre positif constant M tel que

Le noeud ( i Δ x , j Δ y ) est entouré par les noeuds avoisinants montrés sur la

figure I-1. En développant en séries de Taylor pour U i 1, j , U i 1, j , U i  2, j et U i  2, j autour
de la valeur centrale U i , j , nous obtenons :

I-6

I-7

I-8

I-9

Ici etc... et toutes les dérivées sont évaluées au

noeud (i , j).
Adapted by A. Boukhari 10
El-Oued University 2013.
CHAPITRE 1 METHODE DES DIFFERENCES FINIES
5.1 - Approximation à l'ordre 1 en h de la dérivée première :
En négligeant les termes d'ordre 2 et plus dans les équations (I-6) et (I-7)
nous obtenons :

Nous avons donc approché Ux par des différences finies d'ordre 1

progressives (ou à droite) et régressives (ou à gauche) respectivement.
5.2 - Approximation à l'ordre 2 en h de la dérivée première par
les différences centrées :
En soustrayant (I-6) de (I-7) et en négligeant les termes d'ordre 4 et plus nous
obtenons :

C'est l'approximation de Ux par des différences finies centrées d'ordre 2.

5.3 - Approximation à l'ordre 1 en h de la dérivée seconde :
En multipliant l'équation (I-7) par (-2), en ajoutant le résultat à l'équation (I-9)
et en négligeant les termes d'ordre 3 et plus, nous obtenons :

C'est l'approximation de la dérivée seconde à l'ordre 1 en h par les différences

à droite. De même, en utilisant les équations (I-6) et (I-8), nous obtenons
l'approximation à l'ordre 1 en h de Uxx par les différences à gauche :

Adapted by A. Boukhari 11
El-Oued University 2013.
CHAPITRE 1 METHODE DES DIFFERENCES FINIES

5.4 - Approximation à l'ordre 2 en h de la dérivée seconde :

En ajoutant les équations (I-6) et (I-7) et en négligeant les termes d'ordre 4 et
plus, nous obtenons :

C'est l'approximation de la dérivée seconde à l'ordre 2 en h par les différences

centrées.
5.5 - Approximation à l'ordre 2 en h de la dérivée première par
les différences à droite ou à gauche :

De l'équation (I-7) nous pouvons tirer :

En substituant Uxx |i,j par sa approximation à l'ordre 1 en h par les différences

à droite, nous obtenons :

C'est l'approximation à l'ordre 2 en h de la dérivée première Ux par les

différences à droite.
De même en utilisant les équations (I-6) et l’approximation de Uxx |i,j à l'ordre
1 en h par les différences à gauche, nous obtenons l'approximation à l'ordre 2 en h
de Ux par les différences à gauche :

Adapted by A. Boukhari 12
El-Oued University 2013.
CHAPITRE 1 METHODE DES DIFFERENCES FINIES

5.6 - Approximation de Uy et Uyy :

Les approximations de Uy et Uyy sont obtenues de la même façon que celles

de Ux et Uxx et sont données ci-dessous :

Adapted by A. Boukhari 13
El-Oued University 2013.
CHAPITRE 1 METHODE DES DIFFERENCES FINIES

5.7 - Approximation de la dérivée croisée Uxy :

En utilisant l'approximation de Ux par des différences finies centrées d'ordre 2, nous

obtenons :

Finalement :

5.8 - Approximation à l'ordre n en h des dérivées :

En prenant un nombre de points avoisinants de plus en plus grand, on peut

obtenir un nombre illimité d'autres approximations pour chaque dérivée.
Cependant les formes ci-dessus sont les plus compactes et seront donc utilisées
dans la suite du cours.

Adapted by A. Boukhari 14
El-Oued University 2013.
 CHAPITRE 2

MÉTHODES DE RESOLUTION DES E.D.P.

PARABOLIQUES

Il existe essentiellement trois types de méthodes de base pour résoudre les

E.D.P. paraboliques : méthodes explicites, implicites et du type Crank-Nickolson
Nous exposerons chacune d'elles en prenant l'exemple d'une équation linéaire à une
variable d'espace x

II-1
II-2
II-3
II-4

où T est le temps où on veut la solution.

L'équation (II-1) décrit par exemple, en coordonnées adimensionnelles, la
conduction transitoire de la chaleur dans une barre isolée avec une distribution de
température au temps t = 0 et ayant les extrémités maintenues à des températures
qui peuvent être fonction du temps.
Dans les équations II-2 - II-4, f(x) est la condition initiale et g0(t) et g1(t) sont
les conditions aux limites qui sont du type Dirichlet.

Adapted by A. Boukhari 15
El-Oued University 2013.
CHAPITRE 2 METHODE DES DIFFERENCES FINIES

1 - Méthode explicite de résolution :

Afin de résoudre l'équation (II-1) par les différences finies, nous diviserons
1
l'intervalle [0 , 1] en M intervalles Δx  h  par les points
M
1
x0 = 0, x1,…., xM = 1, et l'intervalle [0 , T] en N intervalles Δt  k  par les
N
points t0 = 0, t1,…., tN où M et N sont des nombres entiers arbitrairement choisis.
Le maillage obtenu est représenté sur la figure II-1a.

[Link]-1a : Maillage dans le cas général ; [Link]-1b : La méthode explicite

n est l’indice du temps

Pour chaque noeud intérieur au domaine (c'est-à-dire les noeuds qui n'ont pas
i = 0 , i = M ou n = 0) les dérivées de l'équation (II-1) sont remplacées par les
approximations suggérées dans le chapitre 1. Ainsi, en désignant par Ui,n (ou U in ) la
valeur approchée de la fonction U au noeud (i , n) et en utilisant les différences à
droite pour la dérivée par rapport à t et les différences centrées pour la dérivée
spatiale, nous obtenons :
U in  1  U in U in 1  2 U in  U in 1 II-5

k h2

Adapted by A. Boukhari 16
El-Oued University 2013.
CHAPITRE 2 METHODE DES DIFFERENCES FINIES
ou, en définissant   t /(x ) 2  k / h 2 , alors :

U in 1   U in1  (1  2 ) U ni   U in1 II-6

Dans la figure II-1b, les croix et les cercles indiquent les noeuds intervenant
dans la discrétisation de Ut et Uxx respectivement.
Si tous les U in sont connus à l'instant tn, l'équation (II-6) permet à U in 1 d'être
calculé directement (c'est-à-dire explicitement) à l'instant tn+1 pour 1  i  M  1 .
Pour les noeuds des frontières i = 0, i = M, nous avons d'après les conditions
aux limites (II-3) et (II-4)
U 0n 1  g 0 ( t n 1 )
II-7
U nM1  g1 ( t n 1 )
Puisque les conditions initiales de v sont imposées au temps t = 0 par
l'équation (II-2) :
U 0i  f ( x i )
Les valeurs de U peuvent évidemment être obtenues dans tous les noeuds par
l’application répétée de (II-6) et (II-7). Nous devons calculer toutes les valeurs de U
à l’instant n avant d'avancer à l'instant n + 1.

Exemple II-1 :
Considérons le problème de la conduction de la chaleur des équations II-1- II-4
avec les conditions simples :

1
en choisissant Δx  0.2 et Δ t  0.01 on aura   et l'équation II-6 devient :
4
U in1  2 U in  U in1
U n 1
i  II-8
4

Adapted by A. Boukhari 17
El-Oued University 2013.
CHAPITRE 2 METHODE DES DIFFERENCES FINIES
Les valeurs calculées sont données dans la table suivante :

Une diffusion graduelle de la chaleur dans la barre est mise en évidence par
l'accroissement général de température.

1.1- Condition de convergence de la méthode explicite :

~
L'erreur de discrétisation locale w in est la différence de la valeur exacte U in au

noeud i,n et la valeur approchée U in obtenue par la méthode des différences finies :
~
w in  Uin  Uin II-9

Nous dirons que la schéma des différences (II-6) converge si l'erreur w in tend
vers 0 lorsque Δ x et Δ t tendent vers zéro.
Du développement en série de Taylor, en supposant que U possède un
nombre suffisant de dérivées partielles, nous avons :
k2
U n 1
i  U  k U t  U tt  O(k 3 )
n
i II-10
2!
h2 h3 h4 h5
U n
i 1  U  h Ux 
n
i U xx  U xxx  U xxxx  U xxxxx  O(h 6 ) II-11
2! 3! 4! 5!
Dans (II-10) et (II-11) les dérivées sont évaluées en ( i Δ x , n Δ t ).
En utilisant l'équation II-1, Ut = Uxx avec (II-10) et (II-11) nous aboutirons à :
U in 1   U in1  (1  2  ) U in   U in1  z in II-12

z in h h2
où  U tt  U xxxx  O(k 2 )  O(h 4 ) II-13
k 2 12

Adapted by A. Boukhari 18
El-Oued University 2013.
CHAPITRE 2 METHODE DES DIFFERENCES FINIES
En soustrayant (II-6) de (II-12) et en appliquant (II-9), nous obtenons :
w ni 1   w in1  (1  2  ) w in   w in1  z in II-14
Maintenant supposons que 0    1 / 2 ; les coefficients  et (1   ) sont
positifs et nous avons donc l'inégalité :
w in 1   w ni 1  (1  2 ) w in   w in1  z in
II-15
Posons w (max
n)
 max w in 1 i  M  1

et z (max
n)
 max z in

On a alors :
n 1)
w (max  w (max
n)
 z (max
n)
II-16
En particulier, puisque wmax(0) = 0
w max ( N )  N . Z où Z  max( z (max
n)
)
c'est à dire :
t ( x ) 2
w max ( N )  T U tt  U xxxx  O(t ) 2   O(x )4  II-17
2 12 max

Donc si   ] 0 , 1/2] , l'erreur de discrétisation est en O(h2) et la méthode converge

lorsque h  0.
Finalement, puisque Utt = Utxx = Uxxt = Uxxxx alors pour le choix spécial de   1 / 6
l'erreur de discrétisation est en O(h4) et la méthode converge encore plus vite.
L'inconvénient de la méthode explicite est donc qu'elle nécessite de choisir
t suffisamment petit (de telle façon à satisfaire la condition 0   1 / 2 ) sinon la
solution de l'équation II-1 devient instable. A cause de ce risque d'instabilité, il est
préférable d'utiliser l'une des méthodes suivantes.
2 - Méthode Implicite :
L'équation explicite (II-6) à été obtenue en écrivant les deux membres de (II-1)
à l'instant tn, où la solution est connue. Nous obtenons une équation implicite en
écrivant les deux membres de (II-1) à l'instant tn+1 où la solution n'est pas connue, ce
qui donne :

Adapted by A. Boukhari 19
El-Oued University 2013.
CHAPITRE 2 METHODE DES DIFFERENCES FINIES

U in  1  U in U in11  2 U in  1  U in11
 II-18
k h2
Ou encore :
 λ U in11  (1  2 λ) U in 1  λ U in11  U in II-19
Les noeuds intervenant dans la discrétisation de Ut et Uxx sont montrés sur la
figure 3.3 :

[Link]-2 : La méthode implicite

L'équation (II-19) écrite pour chaque point 1 i  M  1 résulte en un système

de M-1 équations simultanées à M+1 inconnues U 0n 1 ,....., U nM1 . Le nombre

d'inconnues est ensuite réduit à M-1 en incorporant les conditions aux limites :
U 0n 1  g 0 ( t n 1 )
U nM1  g1 ( t n 1 )
On obtient alors le système à matrice tridiagonale :
(1  2 λ) U 1n 1  λ U n2 1  U 1n  λ g 0 ( t i 1 ) pour i 1

 λ U in11  (1  2 λ) U in 1  λ U in11  U in pour 2  i  M  2

 λ U nM12  (1  2 λ) U nM11  U nM 1  λ g1 (t i 1 ) pour i  M  1

Ce système peut être résolu par une méthode quelconque de résolution des
systèmes algébriques linéaires. La méthode la plus adaptée pour ce type de systèmes
est en fait le double balayage de Choleski.
L'algorithme (II-19) permet de déterminer de proche en proche des
approximations de la solution du problème (II-1) aux divers instants tn, n = 1,…..,N.
2.1- Convergence de la méthode Implicite :
On peut montrer, en procédant comme dans le cas de la méthode explicite,
que la méthode implicite est inconditionnellement stable, c'est à dire qu'elle
converge qu'elle que soit la valeur de  .

Adapted by A. Boukhari 20
El-Oued University 2013.
CHAPITRE 2 METHODE DES DIFFERENCES FINIES
Exemple II-2 :
Déterminer, par la méthode implicite, la solution de l'équation (II-1) au temps
t  t avec les conditions suivantes :
II-20

II-21

II-22

utiliser Δx  0.1 et Δ t  0.01

Solution :

[Link]-3 : Le maillage de l’exemple II-2

La valeur de  est :

L'équation (II-19) s'écrit donc :

 U in11  3 U in 1  U in11  U ni II-23
Puisqu'il y a symétrie par rapport à i = 5, il suffit d'appliquer l'équation II-23 aux
noeuds (i , 1) avec i = 1,…..,5. Les valeurs de U1i avec i = 6,…..,9 sont obtenues par
symétrie.
Ainsi, pour i = 1 , l'équation (II-23) devient :
 U10  3 U11  U12  U10

et comme U10  0.2 (d'après la condition initiale II-20), et U10  0 (d'après la condition
aux limites II-22), alors :
3 U11  U12  0.2 II-24

Adapted by A. Boukhari 21
El-Oued University 2013.
CHAPITRE 2 METHODE DES DIFFERENCES FINIES
Pour i=2, 3, 4 on obtient respectivement après incorporation de la condition
initiale (II-20) :
 U11  3 U12  U13  0.4 II-25

 U 12  3 U13  U 14  0.6 II-26

 U 13  3 U14  U 15  0.8 II-27

Pour i = 5, on obtient :
 U14  3 U15  U16  1.0
mais, puisque U6 = U4 par symétrie, alors :
 2 U14  3 U15  1.0 II-28
Le système d'équations algébriques (II-24) - (II-28) est à matrice tridiagonale
et peut donc être résolu par la méthode du double balayage de Choleski. Les valeurs
obtenues sont :
U 11  0.1967 ; U12  0.3902 ; U 13  0.5740 ; U 14  0.7317 ; U15  0.8211
et par symétrie, nous avons :
U16  U14  0.7317 ; U17  U13  0.5740 ; U18  U12  0.3902 ; U19  U11  0.1967
3- La méthode de Crank-Nicholson :

[Link]-4 : La méthode de Crank-Nicholson

La méthode de Crank-Nicholson consiste à écrire les deux membres de

l’équation (II-1) au temps n + 1/2 (montré par un cercle sur la figure II-4). Le
premier membre est approché par les différences centrées et le second est exprimé
comme la demi-somme des approximations aux instants n et n + 1. Nous obtenons
alors : U n  1  U n 1 1
i i
  2x U in  1   2x U in
k 2 2
1  U i  1  2 U in  1  U in11 
n 1
1  U in 1  2 U in  U in 1 
    2 
   
2
Adapted by A. Boukhari
2 h h2 22
El-Oued University 2013.
CHAPITRE 2 METHODE DES DIFFERENCES FINIES
C’est - à - dire:
 λ U in11  2 (1  λ) U in 1  λ U in11  λ U in1  2 (1  λ) U ni  λ U in1 II-29
Le second membre de (II-29) étant connu, on remarque que (II-29) est en fait
une équation implicite, qui se traite exactement comme il à été indiqué à la section
précédente (concernant la méthode implicite).
L'avantage de ce schéma est que pour une valeur donnée de x , l'erreur de
troncature sur le terme t est nettement plus petite que dans les méthodes implicite
et explicite. Ce gain est obtenu au prix d'une complication mineure des calculs.
La méthode de Crank-Nicholson est inconditionnellement stable et converge
quelque soit la valeur du rapport  .
Exemple II-3 :
Utiliser la méthode de Crank-Nicholson pour résoudre le problème de
l'exemple II.2.
Solution :
En substituant pour  dans l'équation (II-29) nous obtenons :
 U in11  4 U in 1  U in11  U in1  U in1 II-30

En appliquant (II-30) pour i = 1,…..,5 et en tenant compte des conditions

initiales et aux limites, on aura :
 4 U 11  U 12  0.40

 U11  4 U12  U13  0.80

 U12  4 U13  U14  1.20

 U13  4 U14  U15  1.60

 2 U14  4 U15  1.60

La résolution de ce système donne alors :
U11  0.1989 ; U12  0.3956 ; U13  0.5834 ; U14  0.7381 ; U15  0.7691

Les valeurs de U 16 jusqu’à U19 sont obtenues par symétrie.

Adapted by A. Boukhari 23
El-Oued University 2013.
CHAPITRE 2 METHODE DES DIFFERENCES FINIES

4- Schémas explicites inconditionnellement stables :

Alors que la méthode explicite familière de la section 3.1 n'est stable qu’avec

la condition 0    1 2 , ils existent plusieurs autres méthodes explicites qui sont

libres de cette restriction. Des exemples de telles méthodes sont illustrées ci-dessous

pour la résolution de l'équation U t  U xx mais toutes peuvent être généralisées à plus

d'une dimension d'espace.
4.1- Méthode de Dufort - Frankel

[Link]-5 : Schéma de Dufort - Frankel

U in 1 - Uin 1 Uin1  U in 1  U ni 1  Uin1


2 t x 2
Cette méthode fait intervenir trois niveaux de temps comme montré sur la
figure (II-5).
4.2- Méthode de Saul'yev
Elle consiste à effectuer successivement deux types d'approximations. En
partant des U ni et en procédant dans le sens des x positives, nous obtenons les U in 1
en utilisant le schéma :
U in 1 - U in Uin11  Uin 1  U ni  Uin1

t x 2
Nous obtenons ensuite les U in  2 en procédant dans la direction des x négatives
et en utilisant le schéma suivant :
U in  2 - U in 1 Uin11  U in 1  Uin  2  Uin12

t x 2
Adapted by A. Boukhari 24
El-Oued University 2013.
CHAPITRE 2 METHODE DES DIFFERENCES FINIES

5- Application des méthodes explicite et implicite au cas général

d'une équation aux dérivées partielles linéaire :
Soit, à résoudre l'E.D.P
f (x, t)  2f (x, t) f (x, t)
 A(x)  B(x)  C(x)f (x, t)  D(x)
t x 2
x II-31
La discrétisation par la méthode explicite donne :
1 n 1 n A B
(f i  f i )  i2 (fin1  2fin  f in1 )  i (f in1  fin1 )  Ci fin  Di
t x 2x
Ou :
 A t B t   2Ai t   A t B t 
fin 1   i 2  i  fin1  1   Ci t  f in   i 2  i  fin1  Di t
 x 2 x   x   x 2x 
2
II-32
On peut montrer que la condition de stabilité de ce schéma est :
II-33

La méthode s'applique, sans aucune modification de principe, aux équations

dont les coefficients A, B, C, D dépendent non seulement de x, mais aussi de t et de
f, ainsi qu'aux fonctions de plusieurs variables d'espace.
L'application de la méthode implicite à l'équation (II-31) donne :
1 n 1 n A B
(f i  f i )  i2 (fin11  2fin 1  fin11 )  i (f in11  fin11 )  Ci fin 1  D i
t x 2x
ou encore :
 Ai Bi  n 1  1 2Ai  n 1  Ai Bi  n 1  f in 

 x 2x  f  
 t x  C  f   x  f     D 
2x 
 
 t
i 1 i i i 1 i

2 2
II-34
Ce schéma est universellement stable. Il peut facilement être généralisé au
cas où les coefficients A, B, C, D, dépendent de t. Aussi, on peut l'adapter au cas où
ces coefficients dépendent de f : il faut dans ce cas résoudre à chaque instant tk
l'équation (II-34) par itération : on se donne une fonction f0(x) qu'on injecte dans A,
B, C et D, la résolution de (II-34) donne alors une autre fonction f1(x), etc.... jusqu'à
la convergence.
Adapted by A. Boukhari 25
El-Oued University 2013.
CHAPITRE 4 METHODE DES DIFFERENCES FINIES

CHAPITRE 3
MÉTHODES DE RESOLUTION DES EDP
ELLIPTIQUES

1- Discrétisation à l'aide des différences finies :

La discrétisation d'une E.D.P elliptique à l'aide des différences finies comporte les étapes
suivantes :
1. On définit un maillage qui couvre le domaine et sa frontière.
2. On approche les dérivées partielles, à l'aide des différences finies, dans tous les noeuds
intérieurs au domaine : les équations algébriques obtenues contiennent des valeurs de la
fonction aux points appartenant à la frontière.
3. On utilise les conditions aux limites pour éliminer les valeurs aux points de la frontière.
Nous obtiendrons alors un système de N équations à N inconnues dont la résolution est faite
par l'une des méthodes données ci-dessous.
Exemple III.1 :
Déterminer la fonction f (x, y) dans le domaine rectangulaire 0  x  a , 0  y  b . f satisfait
à l'équation de Laplace avec les conditions aux limites suivantes :
f (x  0, y)  f 0 III-1

f (x  a, y)  f a III-2

f (c  x  a, y  0)  fa III-3

f y (x  c, y  0)  0 III-4

f y (x, y  b)  0 III-5

[Link]-1 : Domaine de définition de f

Adapted by A. Boukhari 26
El-Oued University 2013.
CHAPITRE 4 METHODE DES DIFFERENCES FINIES

Solution :
1. On définit un maillage (Fig. III-2) qui coïncide avec les frontières du domaine : on utilise
n + 1 pas sur x ( 0  i  n  1 ) de valeur  x = a/(n+1). Le pas  x est choisi comme sous
multiple de a et (a - c) de façon que x = c corresponde au pième pas sur x. On utilise aussi
m+1 pas sur y ( 0  k  m  1 ) de valeur  y=b/(m+1).

[Link]-2 : Maillage pour l’exemple III.1.

2. Dans chaque noeud intérieur (1  i  n et 1  k  m ) on approche l'équation f  0 à l'aide

des différences finies, ce qui donne :
1 1
(fi 1,k  2f i,k  fi 1,k )  2 (f i,k 1  2f i,k  f i,k 1 )  0 III-6
x 2
y
3. Les conditions aux limites sur les frontières i= 0, i = n+1, k= 0 et k = m+1 sont maintenant
incorporées dans les équations des noeuds situés sur les lignes i=1, i=n, k=1 et k=m
respectivement.
Noeuds (i, 1), i=1, n
Pour i =1 la condition aux limites imposée est f y  0 (équ. III-4) sur l'axe y =0 et f=f0 sur

l'axe x =0 (équ. III-1).

L'équation (III-6) avec ces conditions donne :

Adapted by A. Boukhari 27
El-Oued University 2013.
CHAPITRE 4 METHODE DES DIFFERENCES FINIES

III-7
Pour i= 2 à p-1 On obtient :

soit :
III-8

On continue ainsi à incorporer les conditions aux limites, ce qui donne : pour p  i  n  1
III-9

pour i=n
III-10
'
La même méthode est utilisée pour les noeuds (i, m), (1, k) et (n, k). Nous obtiendrons
ainsi un système de n équations linéaires pour chaque ligne k (1  k  m ) soit au total m x n = N
équations à N inconnues. Les inconnues sont les valeurs de f aux noeuds intérieurs du domaine.
NB : La condition aux limites f y  0 sur la frontière k=m est approchée en utilisant les

différences à gauche et non pas les différences à droite comme nous l'avons fait pour les
noeuds (i, k = 1).

2 - Méthode de résolution directe

Le système obtenu peut être résolu par la méthode d'élimination de Gauss comme illustrée
par l'exemple suivant .
Exemple III.2 :
Déterminer le champ des valeurs prises par la fonction T(x, y) dans le domaine et avec les
conditions aux limites suivantes :
 2T  2T
 0 dans un carré de côté C III-11
x 2 y2
T0 le long de x = 0 et y = 0 III-12
T  x3 le long de y = C III-13

Adapted by A. Boukhari 28
El-Oued University 2013.
CHAPITRE 4 METHODE DES DIFFERENCES FINIES
T  16y le long de x = C III-14
Ce problème correspond, par exemple, au cas concret suivant : détermine la température à
l'intérieur d'une plaque carrée (ou dans une barre infinie ou parfaitement isolée à chaque extrémité
de section carré) où 2 côtés sont à 0 °C et les 2 autres sont à températures imposées suivant les
lois T = x3 et T = 16y.
Solution :
En choisissant le même pas dans les deux directions h  x  y , la discrétisation pour
les noeuds à l'intérieur du maillage donne :

Soit :

III-15
L'application de ce schéma aux points à l'intérieur du domaine dans le cas d'un maillage
4x4 donne les équations :

Les quantités T0,0; T1,0; T2,0; T3,0; T4,0; T4,1; T4,2; T4,3; T4,4; T3,4; T2,4; T1,4; T0,4; T0,3; T0,2 et T0,1
sont connues d'après les conditions aux limites. On a donc un système d'équations linéaires: 9
équations à 9 inconnues. En posant Ul= T1,1, U2= T2,1, U3= T3,1, U4= T1,2, U5= T2,2, U6= T3,2,
U7=T1,3, U8= T2,3 et U9= T3,3 (voir figure III-3), on peut écrire :
Adapted by A. Boukhari 29
El-Oued University 2013.
CHAPITRE 4 METHODE DES DIFFERENCES FINIES

La résolution de ce système par la méthode d'élimination de Gauss donne [U1, U2, U3, U4, U5,
U6, U7, U8, U9] = [2,607; 5,964; 10,500; 4,464; 10,750; 20,036; 4,50; 12,536; 26,893].

[Link]-3 : Maillage pour

l’exemple III.2.

Le problème essentiel de la méthode directe est la place mémoire avec en conséquence un

temps élevé d'accès aux mémoires périphériques. Ceci devient dans la pratique prohibitif si le
nombre de variables indépendantes est supérieur à deux.
Par rapport à la méthode itérative exposée ci-dessous, la méthode directe possède les
avantages et inconvénients suivants :
 Avantages :
1. temps calcul en général plus petit pour une même précision
2. dans certains cas la méthode itérative peut ne pas converger.
 Inconvénients :
1. Occupation mémoire importante
2. risque d'erreur d'arrondi importante si certains pivots sont trop petits,
problème qu'on ne rencontre en pratique pas dans la méthode itérative.
3. la méthode directe n'est pas applicable aux équations non - linéaires.

Adapted by A. Boukhari 30
El-Oued University 2013.
CHAPITRE 4 METHODE DES DIFFERENCES FINIES

3 - Méthode itérative de Gauss - Seidel

Après multiplication par x 2 et regroupement des termes, l'équation (III-6) s'écrit :

Ou

III-16

Le procédé itératif est alors appliqué au système d'équations (III-16) comme suit :
(0)
i. on se donne une distribution fi,k avec i =1, n et k =1, m.
(1) (0)
ii. on calcule pour chaque couple (i, k) une nouvelle valeur fi,k à l'aide des valeurs f h,1
(1)
choisies ou f h,1 calculées. h et 1 étant les indices des noeuds avoisinants de (i, k).
(2) (3)
iii. on continue à calculer fi,k , fi,k etc .... de la même façon jusqu'à la convergence, c.-à-d. :

où j est le numéro d'itération, et  est le critère de convergence.

Dans le cas de l'équation de Laplace, la méthode de Gauss-Seidel converge toujours. Pour
d'autres E.D.P plus compliquées elle pourra diverger. Cette méthode est en général préférable aux
méthodes directes car elle est plus facile à programmer. Aussi, elle demande moins de place
mémoire et peut être appliquée à des matrices de grande taille et à des équations algébriques non-
linéaires.
La convergence peut être accélérée en utilisant un facteur de surrelaxation  , avec
1   2 .
( j) ( j) ( j1)
fi,k fi,k fi,k
Soit la valeur de fi,k à la jème itération, l'introduction du dans (III-16) donne :
( j1) ( j1) 
fi,k fi,k
au lieu d'utiliser pour l'itération suivante, on utilise donnée par :
III- 1 7

Exemple III.3 :
Résoudre le problème de l'exemple III.2 par la méthode de Gauss-Seidel avec surrelaxation.
Prendre un carré de côté c = 4 et imax = kmax = 8. Faire le calcul avec deux facteurs de surrelaxation
 =1.4 et  =1.45, Conclusion.
Adapted by A. Boukhari 31
El-Oued University 2013.
CHAPITRE 4 METHODE DES DIFFERENCES FINIES
Solution :
Soit R le résidu au noeud (i, k),

La combinaison des équations (III-15) et (III-17) donne :

III-18
où P est le nombre d'itérations.
En choisissant  =1.4 on obtient après 21 itérations les résultats de la table III.1 avec un
résidu maximal Rmax = 6.4144 x 10-3 . Lorsqu'on augmente le facteur de surrelaxation à 1.45 on
obtient de meilleures résultats (Rmax = 3.13187 x 10-3 seulement) en un nombre plus petit
d'itérations (table III.2).

[Link]-1 : Valeurs calculées lorsque  =1.4.

[Link]-2 : Valeurs calculées lorsque  =1.45.

Adapted by A. Boukhari 32
El-Oued University 2013.
CHAPITRE 4 METHODE DES DIFFERENCES FINIES
4- Méthode itérative du double balayage de Choleski
C'est une méthode très efficace de résolution des équations elliptiques. Ayant obtenu à
(p)
l'itération numéro p une estimation de fi,k on écrit l'équation (III-6) à l'itération p+1:

III-19

C'est à dire que nous calculons les valeurs de f sur la ligne k en utilisant les anciennes
valeurs de f sur les lignes k+1 et k-1 (valeurs obtenues à l'itération (p)). L'équation (III-19) écrite
pour tous les noeuds (i, k) (1  i  n ) de la ligne k donne un système de n équations algébriques à
n inconnues fi,k, f2,k, …fn,k . Le système est à matrice tridiagonale et est résolu par double balayage.
Les autres lignes sont traitées successivement de la même manière. On passe ensuite aux
colonnes i =1, ..., n : ici on calcule les valeurs de f sur la colonne i en utilisant les anciennes
valeurs de f sur les colonnes i-1 et i+1.
Le balayage des lignes et des colonnes constitue une itération. On effectue autant
d'itérations qu'il est nécessaire pour obtenir la convergence. Pour assurer une convergence plus
rapide, on doit balayer la première ligne et la première colonne, puis la 2ème ligne et la 2ème
colonne, etc…, et utiliser chaque fois les valeurs qui viennent d'être déterminées, au lieu des
anciennes valeurs.
5- Approximation des dérivées aux noeuds avoisinants une frontière
irrégulière
Lorsque la frontière rencontre le maillage rectangulaire en des points qui ne sont pas des
noeuds du maillage (voir figure III-4) les formules précédemment utilisées pour approcher les
dérivées aux noeuds à proximité du maillage ne sont plus valables car le pas varie de part et
d'autre de ces noeuds. Ce paragraphe concerne les approximations par les différences finies des
dérivées en un point tel que ‘0’, proche de la frontière F sur laquelle les valeurs de f sont
supposées connues. Ainsi, pour le cas simple d'un maillage carré de pas h le théorème de Taylor
donne :

 2f0
L'élimination de donne :
x 2

Adapted by A. Boukhari 33
El-Oued University 2013.
CHAPITRE 4 METHODE DES DIFFERENCES FINIES
f 0
de même, l'élimination de mène à :
x

Donc, l'équation de Laplace au point 0 s'écrit :

Figure III-4.

Exemple III.4 :
Résoudre l'équation de Laplace
U xx  U yy  0 dans D

U(x, y)  g(x, y) sur F
dans le domaine montré dans la figure III.5

Figure III-5.
Solution :
La discrétisation pour les noeuds 1, 3 et 4 est normale, c'est à dire différences centrées
pour Uxx et Uyy. Pour le noeud 2 , la discrétisation est normale pour Uxx mais spéciale pour Uyy.
On peut calculer 2 soit de la relation 2 h  y6  y 2 où y6 est obtenu de l'équation de la droite

‘57’, ou en utilisant le théorème de Thalès pour les triangles ‘267’ et ‘157’. On obtient dans les
deux cas 2 =1/2. La dérivée seconde Uyy au noeud 2 est donc donnée par :

Adapted by A. Boukhari 34
El-Oued University 2013.
CHAPITRE 4 METHODE DES DIFFERENCES FINIES
Nous obtenons le système :

EXERCICES RESOLUS
Exercice 1 :
Les contraintes de cisaillement et les déplacements des point dans un long cylindre
uniforme en torsion peuvent être calculés de la fonction des contraintes  qui est constante
autour du périmètre C d'une section droite dans le plan XOY et qui satisfait
 2  2
l’équation :   2  0 en chaque point P(x, y) de la section.
 X 2 Y 2
Le problème peut être écrit sous forme adimensionnelle en mettant x =X/L, y =Y/L et
   / L2 où L est une longueur représentative quelconque. Ce qui donne l'équation
adimensionnelle :
 2  2
 20
x 2 y 2
Pour un cylindre dont la section et les conditions aux limites sont montrés sur la figure ci-
dessous, calculer la valeur de  aux noeuds d'un maillage pour lequel x     y

Solution :
Puisqu'il y a symétrie par rapport à Ox, Oy et la droite OF, on calcule seulement les
valeurs de  aux noeuds a, b, c, d et e. L'approximation par les différences finies à l'équation
différentielle au noeud (i, j) est :

Adapted by A. Boukhari 35
El-Oued University 2013.
CHAPITRE 4 METHODE DES DIFFERENCES FINIES

Puisque dans notre cas h= k=1/2, on obtient :

L'application de cette équation en chacun des noeuds donne :

La solution de ces équations par l'élimination de Gauss est :

Exercice 2 :
Le courant dans un enroulement électrique donne naissance à une production interne de
chaleur qui est dissipée par rayonnement à travers les frontières. La distribution de température
stationnaire dans un enroulement ayant une section rectangulaire peut donc être approchée par la
solution de l'équation de la conduction stationnaire à travers un rectangle, avec production interne
de la chaleur et rayonnement de son périmètre vers le milieu qui l'entoure.
A cause du fait que la conductivité thermique K parallèlement à Ox diffère de sa valeur
2U 2U q
 K parallèlement à Oy, l'équation de la température U est :    ; où q est le
x 2
y 2
K
nombre d'unités de chaleur généré par unité de surface par seconde. Résoudre cette équation dans
q
le cas où :    et  16 et avec les conditions aux limites montrées sur la figure (III-6). La
K
section est un carré de côté 2. Prendre h =1/4 = k.
Solution :

Adapted by A. Boukhari 36
El-Oued University 2013. Figure III-6
CHAPITRE 4 METHODE DES DIFFERENCES FINIES

D'après les conditions aux limites, la solution est symétrique par rapport aux axes Ox et
Oy.(figure III-6). L'approximation de l'équation par les différences finies au noeud (i, j) est :
1 3
2
(Ui 1, j  2U i, j  U i 1, j )  2 (Ui, j 1  2Ui, j  Ui, j1 )  16
h h
Puisque h = k =1/4 alors :
U i 1, j  U i 1, j  3U i, j1  8Ui, j  3Ui, j1  1

En dénotant les noeuds par un seul indice comme montré sur la figure, nous obtenons :
Noeuds (1), (2), (3) et (4):

D'après la condition aux limites on a:

L'élimination de U-5 donne :

1
9 U1  2U 2  6U 5  1
2
Similairement, les points 2 , 3 et 4 donnent :
1
U1  9 U 2  U 3  6U 6  1
2
1
U 2  9 U 3  U 4  6U 7  1
2
1
U 3  9 U 4  6U8  1
2
Noeud (5) :
3U1  8U 5  2U 6  3U 9  1

Noeud (6):
3U 2  U 5  8U 6  U 7  3U10  1
Des équations similaires peuvent être facilement écrites pour les noeuds 7 à 20, donnant
vingt équations algébriques linéaires à 20 inconnues. Le système, sous forme matricielle est :
A U  b où b est le vecteur colonne dont les composants sont U1, U2,…,U20 ; b est un vecteur
colonne dont chaque composant est égal à -1, et A est une matrice qui peut être exprimée sous
forme partitionnée

Adapted by A. Boukhari 37
El-Oued University 2013.
CHAPITRE 4 METHODE DES DIFFERENCES FINIES

Où

La solution est :

Exercice 3 :
2U 2U
La fonction U satisfait l'équation :   32U  0 en chaque point intérieur au
x 2 y 2
carré x = ±1, y = ±1 et est sujette aux conditions aux limites :

Déterminer le système d'équations à résoudre dans le cas où x  1/ 4  y .

Solution :
En utilisant un maillage similaire à celui de l'exercice 2, et en prenant en considération la
symétrie par rapport à l'axe des y nous obtenons un système de 35 équations linéaires qu'on peut
écrire sous la forme matricielle : A U  b

Où

et A est une matrice dont la forme partitionnée est :

où

Adapted by A. Boukhari 38
El-Oued University 2013.
CHAPITRE 4 METHODE DES DIFFERENCES FINIES

et

N.B : Ce problème est un cas particulier, sous forme adimensionnelle de l’équation :

2U 2 U (U  U 0 )
 2  2H  0 qui détermine la température stationnaire aux points d'une plaque
x 2
y KD
mince et plate rayonnant la chaleur de sa surface dans un milieu à une température U0. D
représente son épaisseur, k sa conductivité thermique et H la conductance de la surface.
Exercice 4 :
Le mouvement lent, permanent d'un fluide visqueux dans un tube cylindrique, dont la
section est une surface S limitée par une courbe fermée C, peut être calculé de la fonction  qui
satisfait l'équation de Laplace dans tous les points de S et est égale à 1/2 r2 sur la courbe C, où
(r,  sont les coordonnées polaires d'un point dans le plan de S.
Calculer une solution par les différences finies de l'écoulement à travers le secteur
circulaire limité par les lignes  =0 et  =0,8 radians et le cercle r =1, aux noeuds du maillage
défini par: r=1/3 i,(i=1,2),  =0,2 j,(j=1,2,3).
Solution :
Le problème est symétrique par rapport à  = 0,4. Donc il y a quatre inconnues, 1 en

(1/3; 0,2),  2 en (1/3, 0,4),  3 en (2/3, 0,2) et  4 en (2/3, 0,4). Les valeurs à la frontière et la
forme polaire de l'équation de Laplace, donnent les équations :
1 7
521  25 2  1  3  1  0
2 18
1
501  52 2  1  4  0
2
3 1 1 1
1  14  3  6  4  2  0
4 2 4 72
3 1 1 5
 2  12  3  14  4   0
4 2 2 8
Leur solution est:

Exercice 5 :
Une lamelle semi-circulaire de rayon 2h et de conductivité thermique uniforme a son
diamètre maintenu à la température de 0° et sa circonférence à la température de 100°. Calculer

Adapted by A. Boukhari 39
El-Oued University 2013.
CHAPITRE 4 METHODE DES DIFFERENCES FINIES
par les différences finies les températures (dans le cas stationnaire) aux noeuds d'un maillage carré
de côté h.
Solution :

  3 1

 2  2
 0
x 2 A
y 2 A

(a)

 2  2
 0
x 2 B
y 2 B

(b)

De (a) et (b) nous obtenons :

A  70, 5 et B  60, 2
Bien que le maillage utilisé est grossier, les températures données sont correctes à presque
1 ° près. La solution analytique de ce problème par la méthode de séparation de variables et en
fonction des coordonnées polaires (r,  est :
2n 1
400  1  r 
   
 n 1 (2n  1)  2h 
sin(2n  1)

et donne : A  70, 5 et B  59 .

Adapted by A. Boukhari 40
El-Oued University 2013.
CHAPITRE 4 METHODE DES DIFFERENCES FINIES

CHAPITRE 4
MÉTHODES DE RESOLUTION DES EDP
HYPERBOLIQUES

1- Introduction :
Il existe deux méthodes de résolution des équations aux dérivées partielles hyperboliques :
1. la méthode des différences finies
2. la méthode des caractéristiques
La méthode des différences finies est appliquée pour les équations dont les solutions sont
continues. Lorsqu'une fonction ou l'une de ses dérivées est discontinue en un point de la frontière,
on sait que cette discontinuité se propage le long des caractéristiques passant par ce point. Dans ce
cas nous devons utiliser impérativement la méthode des caractéristiques car elle donne des
résultats plus précis que la méthode des différences finies.
2- Résolution par les méthodes des différences finies
Les méthodes des différences seront exposées pour l'équation simple (IV-1). Elles restent
entièrement valables pour des équations plus compliquées.
2.1- Méthode explicite :
Soit à résoudre l'équation des ondes :
 2f 1  2f
 IV-1
x 2 c2 t 2
f (x, t  0)  F(x) pour 0  x  a IV-2

 f
 t (x, t  0)  G(x) pour 0  x  a IV-3
où c est une constante.
On choisit un maillage de l'espace (x, t) correspondant à des pas x et t (figure IV-1). En
supposant f(x, t) connue jusqu'à l'instant k, la discrétisation par les différences centrées au noeud
(i, k) donne :

IV-4

IV-5
Adapted by A. Boukhari 41
El-Oued University 2013.
CHAPITRE 4 METHODE DES DIFFERENCES FINIES
en portant (IV-4) et (IV-5) dans (IV-1), on obtient :

ou encore, en regroupant les termes :

IV-6

Si f est connue pour tout x  ix jusqu'en t  kt , le second membre de (IV-6) est connu,
ce qui détermine fi, k+1. Le problème du démarrage reste à régler car en mettant k=0 dans (IV-6),
on voit clairement sur la figure (IV-1) qu'il est nécessaire de connaître fi,0 et fi,-1.
fi,0 est obtenue de la condition initiale :
fi,0=Fi IV-7
et fi,1 est déterminée par la donnée de la dérivée à t =0. En approchant cette dérivée par les
différences centrées, nous aurons :
 f  1
 t   2t (fi,1  fi,1 ) IV-8
 i,0
Les équations (IV-6) écrites en k =0 deviennent alors, en tenant compte des conditions
(IV-7) et (IV-8) :
IV-9

IV-10
L'élimination de fi,1 entre ces deux équations donne fi,-1 pour tout i. On pourra alors
calculer facilement fi,1 par l'équation (IV-9) : Il suffit de mettre i =1 à n.
Ayant calculé les valeurs aux noeuds (i, 1), nous passons ensuite aux noeuds (i, 2). Pour
cela, il faut appliquer l'équation (IV-6) en k=1 : on voit alors que fi,2 ne peut être déterminée que
pour i =2 à n-1. De même pour les noeuds (i, 3), fi,3 peut être déterminée seulement pour i=3 à n-2
et ainsi de suite, de sorte que le domaine où f peut être déterminée se rétrécit lorsque t croît (figure
IV-2).
La solution ne peut être déterminée complètement (c'est à dire pour tout t  0 et tout x tel
que 0  x  a ) que si l'on impose des conditions aux limites en x =0 et x =a par exemple :

En effet, avec ces conditions on peut appliquer (IV-6) pour toutes les valeurs de i (i =1 à n)
quelque soit k. On peut montrer que la solution de (IV-6) n'est pas stable si t  (x / c) . C'est la

Adapted by A. Boukhari 42
El-Oued University 2013.
CHAPITRE 4 METHODE DES DIFFERENCES FINIES
raison pour laquelle il faut employer les méthodes implicites si on désire calculer f en des pas
assez grands dans le temps.

Adapted by A. Boukhari 43
El-Oued University 2013.
CHAPITRE 4 METHODE DES DIFFERENCES FINIES

[Link]-1 [Link]-2

Exemple IV.1
soit l'E.D.P. Utt=Uxx avec les conditions aux limites U(0,t)=U(1,t)=0 pour t>0
et les conditions initiales U(x, 0)  0.5x(1  x) et Ut (x, 0)  0 pour 0  x  1 .

Résoudre ce problème par la méthode des différences finies.

Prendre    et x  h  0,1 .
Solution
Calculons d'abord les valeurs au temps t=0 (k=1) : d’après la condition initiale
U(x,0)=0,5x(1-x) avec h=0,1.
U 0,0  0; U1,0  0.5  0.1(1  0.1)  0.045; U 2,0  0.5  0.2(1  0.2)  0.080;
U 3,0  0.5  0.21  0.105; U4,0  0.120; U5,0  0.125; U 6,0  0.120;
U 7,0  0.105; U8,0  0.080; U9,0  0.045; U10,0  0.0;

Pour calculer les valeurs U0,1 et U10,1, on va utiliser les conditions aux limites du problème:
U0,1=U10,1=0. Les autres valeurs U1,1 (i =1 à 9) sont calculées à partir des formules (IV-9) et (IV-
10) en posant    et Gi=0:
Ui,1 = 0.5x(Ui-1,0 + Ui+1,0)
Ainsi U1,1 = 0.5 (U0,0 + U2,0) = 0.5 x 0.080 = 0.040
U2,1 = 0.5 (U1,0 + U3,0) = 0.5 x 0.150 = 0.075, U3,1 = 0.100, U4,1 = 0.115,
U5,1 = 0.120, U6,1 = 0.115, U7,1 = 0.100, U8,1 = 0.075, U9,1 = 0.040
Pour k>1, on utilise seulement les conditions aux limites du problème et la formule IV-6
avec    . Par exemple, lorsque k=2, l'équation IV-6 devient :
Ui,2 = Ui-1,1 + Ui+1,1 – Ui,0 (i =1, 2,…,9) , ainsi :
U1,2 = U0,1 + U2,1 – U1,0 = 0 + 0.075 – 0.045 = 0.030
U2,2 =U1,1 +U3,1 –U2,0 = 0.040 + 0.100 – 0.080 = 0,060
U3,2 = 0.085 ; U4,2 = 0.100 ; U5,2 = 0.105 ; U6,2 = 0.100 ; U7, 2 = 0.085 ; U8,2 = 0.060 ; U9,2 = 0.030
et ainsi de suite pour les autres valeurs de k.
Adapted by A. Boukhari 44
El-Oued University 2013.
CHAPITRE 4 METHODE DES DIFFERENCES FINIES
2.2 - Méthode Implicite :
La dérivée par rapport à t est approchée comme dans la méthode explicite (équation IV-4).
 2f 
 2  est cependant évaluée en prenant la demi-somme des valeurs de cette dérivée en k-1 et
 x i,k
k+l. on écrit donc (voir fig IV-3) :

 2f  1   2 f    2f  
 2    2   2   IV-11
 x i,k 2  x i,k 1  x i,k 1 
C'est à dire :
  2f 
2  i 1,k 1
 2f i,k 1  fi 1,k 1    fi 1,k 1  2fi,k 1  fi 1,k 1  
1
 2   f
 x i,k 2x

[Link]-3

En portant (IV-11) et (IV-8) dans (IV-1), il vient :


fi,k 1  2fi,k  fi,k 1   fi 1,k 1  2fi,k 1  fi 1,k 1    fi 1,k 1  2fi,k 1  fi 1,k 1  
2 

où   c 2 t 2 / x 2 On obtient en réarrangeant les termes :

  
 fi 1,k 1  (1  fi,k 1  fi 1,k 1  2f i,k  (fi 1,k 1  2fi,k 1  fi 1,k 1 ) IV-12
2 2 2
On obtient un système à matrice tridiagonale qui peut être résolu par double balayage, ou
par méthode itérative, à condition de connaître les conditions aux limites, par exemple en x =0 et
x =a : ceci n'était pas nécessaire dans la méthode explicite.
L'avantage principal de la méthode implicite est qu'elle est stable quelque soit t : on peut
donc utiliser des pas sur le temps plus importants que dans la méthode explicite, ce qui permet de
gagner du temps, au prix d'une complication mineure des calculs.

Adapted by A. Boukhari 45
El-Oued University 2013.
CHAPITRE 5 METHODE DES DIFFERENCES FINIES

CHAPITRE 5
PROPRIETES DES SCHEMAS NUMERIQUES
Nous avons montré dans les chapitres précédents que différents schémas formulés pour le
même problème ne se comportent pas de la même manière. En particulier la stabilité et la
précision de la solution ne sont pas les mêmes. Celles-ci sont seulement deux parmi plusieurs
propriétés par lesquelles les schémas numériques sont classifiés et caractérisés.
Dans ce chapitre, nous considérerons brièvement les propriétés les plus importantes des
formulations de différences finies : la consistance, la stabilité, la convergence et les propriétés de
transport et de conservation.
1- Consistance :
Considérons une fonction U du temps et d'une seule variable d'espace x. Notons l'équation
aux dérivées partielles régissant U par :
A(U)  0 V-1
Nous pouvons avoir une équation approchée de (V-1) en remplaçant les dérivées de (V-1)
par leur approximations, ce qui donne le schéma numérique :

A(V) 0 V-2
Par exemple, partant de l'E.D.P
U  2 U
A(U)   0 V-3
t x 2
Nous pouvons considérer le schéma numérique :

 V n 1  Vin Vin1  2Vin  Vin1

A(V)  i  0 V-4
s h2
où les Vin sont des approximations de U(a+ih, ns)= U in , i = 0,…,Nx ; n=0,…,Nt ;

ba T
h et s 
Nx Ns

Définition de la consistance :
Le schéma (V-2) est consistant avec l'équation (V-1) si on a
 n)  0
max (i,n ) A(U in )  A(U quand (h,s)  0
i

Exemple V.1 : Etudier la consistance du schéma (V-4) avec l'équation (V-3)

Solution :

Adapted by A. Boukhari 46
El-Oued University 2013.
CHAPITRE 5 METHODE DES DIFFERENCES FINIES

Exemple V.2 :
Etudier la consistance avec l équation (V-1) du schéma du Dufort-Frankel.

 V n 1  Vin 1 Vin1  Vin 1  Vin 1  Vin1

A(V)  i  V-5
2s h2
Solution : on a :
U s 2  2 U s3  3 U
U in 1  U ni  s    O(s 4 ) V-6
t 2 t 2
6 t 3

U h 2  2 U h 3  3 U
U in1  Uin  h    O(h 4 ) V-7
 x 2 x 2
6 x 3

Toutes les dérivées étant prises au point (a+ih , ns).

Uin 1  Uin 1 U
  O(s 2 )
2s t
Uin1  U in 1  Uin 1  Uin1  2 U s 2  2 U s4
   O(h 2
)  O( )
h2 x 2 h 2 t 2 h2

 s4 s2  2 U 
Donc on a: A(U)  A(U)  O(s )  O(h )  O( 2 )  2  2 
2 2

h h  t  a  ih,ns

On doit donc considérer deux cas :

s s2 s4
 Si reste constant quand (h, s)  0 on a  0 et  0 ; le schéma (V-5)est alors
h2 h2 h2
conditionnellement consistant avec l'équation (V-1).
s
 Si  C  Const quand (h, s)  0, on a alors
h
s4 2  U
2

A(U)  A(U)  O(s )  O(h )  O( 2 )  C
2 2

h t 2
et le schéma (V-5) n'est plus consistant avec l'équation (V-1) mais avec l'équation
U  2 U 2 U
 2  C2 2  0 V-8
t x t
Adapted by A. Boukhari 47
El-Oued University 2013.
CHAPITRE 5 METHODE DES DIFFERENCES FINIES
2- Définition de la Stabilité :
Une formulation est stable si l'erreur de discrétisation e n  U n  V n est bornée (c.-à-d.
n'augmente pas monotoniquement) quand n - l'indice du temps ou le nombre d'itérations-
augmente.
3- Analyse de la Stabilité :
Il existe plusieurs méthodes pour tester la stabilité des schémas numériques. Parmi celles-
ci, les méthodes de Von-Neumann (ou de Fourier) et la méthode matricielle sont les plus utilisées.
3.1- Méthode de Von-Neumann
Si au cours de l'application d'un schéma numérique, on remplace une valeur VjP

(approximation de U(a+jh , ps) par une valeur perturbée,

Wjp  Vjp   pj V-10

puis on continue les calculs avec Wjp , on obtiendra alors des valeurs Wjn différentes des Vjn pour

tous les indices (j, n) qui interviennent dans ces calculs. La déviation au point (j, n) est définie par
la relation : nj  Wjn  Vjn V-11

Dans la pratique, la perturbation  pj ; provient des erreurs d'arrondi ou de calculs-machine

et peut s'introduire en plusieurs noeuds du maillage. Il est facile pour certains types de problèmes
linéaires d'étudier la déviation qui correspond à une perturbation ayant la forme :

et qui est appliquée à tous les points d'une ligne du maillage pour un temps t =n s donné et pour
tous les x = a + k h ; k= 0,...,Nx.
Lorsque la perturbation est développable en série trigonométrique ou de Fourier, la
déviation totale est obtenue en faisant la somme des déviations de chacun des termes de la série.
Dans le cas particulier où les variables sont séparées c'est à dire lorsque :

V-12
Les valeurs de Vjn1 sont données par : V-13

La substitution de (V-12) et (V-13) dans la relation (V-2) du schéma numérique nous

donne :
V-14
Une condition suffisante de stabilité est que :   1, n, j,  V-15

Cette méthode est souvent appelée : Méthode de Fourier ou méthode de Von-Neumann.

Adapted by A. Boukhari 48
El-Oued University 2013.
CHAPITRE 5 METHODE DES DIFFERENCES FINIES
Exemple V.3 : Appliquer la méthode de Von-Neumann aux schémas explicite et implicite du
chapitre deux.
Solution : L'application du schéma explicite :

donne, en notant que t = ns et x = jh :

Pour que   1 , nous devons avoir pour tout  :

h 1
C'est-à-dire : sin 2  
2 2
h 1
sin 2  peut prendre sa valeur maximum 1, il suffit donc que   . retrouvant la condition déjà
2 2
rencontrée (Cf. Chapitre 2). L'application de la méthode de Von-Neumann au schéma implicite :

donne :

1
d'où  
h
1  4 sin 2 
2
Nous avons donc toujours   1 . Le schéma est donc inconditionnellement stable.

3.2- Méthode d'analyse matricielle :

La méthode de Von-Neumann ne tient pas compte des conditions aux limites. Pour faire ça
nous devons écrire les divers schémas numériques sous la forme :

V-16
où Vn+1 est le vecteur de composantes Vjn 1 , j =0,…Nx et où les matrices A et B tiennent compte

des conditions aux limites.

Si les valeurs propres de la matrice M = A-1 B sont toutes de module inférieur ou égale à 1,
la méthode définie par l'équation (V-2) est stable.
Adapted by A. Boukhari 49
El-Oued University 2013.
CHAPITRE 5 METHODE DES DIFFERENCES FINIES
Exemple V.4 : Etudier la stabilité du schéma explicite classique

V-17
pour lequel la matrice A, de dimensions (N-1) x (N-1), est

où r = k/h2 >0. Supposer que les valeurs aux frontières sont connues (problème de Dirichlet).

Solution :
La matrice M peut être réécrite :

où IN-1est la matrice unité d'ordre (N-1) et TN-1 une matrice (N-1) x (N-1) dont les valeurs propres
 s sont données par :

Les valeurs propres de M sont alors :

Donc les équations V-17 seront stables si maxs 1  4r sin 2 s  1

2N

L'inégalité de gauche donne :

Adapted by A. Boukhari 50
El-Oued University 2013.
CHAPITRE 5 METHODE DES DIFFERENCES FINIES

quand h  0, N   et sin 2 (N  1) 1
2N
Donc r  1 / 2
4- Convergence
Définition :
La méthode (V-2) est dite convergente si les valeurs Vin obtenues vérifient :

pour tester la convergence d'une méthode de différences finies on se base généralement sur le
théorème d'équivalence de Lax donné ci-dessous sans démonstration.
Théorème d'équivalence de Lax :
 (V) une méthode
Soient A(U) = 0 un problème d'E.D.P bien posé au sens de Cauchy, et A
d'approximation par les différences finies supposée consistante avec A(U). Alors la stabilité de A
est une condition nécessaire et suffisante de convergence.
5- Propriété transportive (ou de transport)
Une formulation est dite de posséder la propriété transportive si une perturbation est
"advectée" par la formulation seulement dans la direction de la vitesse.
Exemple V.5 :
considérons la formulation :

V-18
 
appliquée au problème schématisé ci-dessous (fig. V-1a) et dont l'E.D.P est : U  0.
t y
Physiquement, la perturbation peut seulement être convectée vers la droite par U.

Considérons pour c=1 la solution en n=1 et n=2 (fig. V-1b et V-1c). On peut remarquer
que le schéma mène à une propagation de perturbation vers l'amont et donc ne possède pas la
propriété transportive. Ceci est le cas de toutes les différences centrées pour la convection.

Adapted by A. Boukhari 51
El-Oued University 2013.
CHAPITRE 5 METHODE DES DIFFERENCES FINIES

Fig. V-1
6- Propriété de conservation :
Cet aspect est considéré en conjonction avec l'approche de solution par la méthode des
"volumes de contrôle". Ici nous allons considérer uniquement une formulation plus générale.
 Une approximation est dite conservative si elle préserve la relation :
V-19

Fig. V-2

L'équation (V-19) est une conséquence de l'application du théorème de divergence de

Gauss à l'équation transitoire de convection/diffusion :

   (     V-20
t
Alors que (V-20) exprime le principe de conservation à un niveau microscopique,

Adapted by A. Boukhari 52
El-Oued University 2013.
CHAPITRE 5 METHODE DES DIFFERENCES FINIES
l'équation (V-19) exprime ce principe pour un volume fini. L'équation (V-19) dit tout simplement

que le taux d'accumulation de  dans le volume V est égale au flux net de  par convection V et
par diffusion (  ) à travers la surface S.
Exemple V.6 :
La formulation
V-21
n'est pas conservative alors que le schéma
V-22

l'est.
La formulation (V-21) provient de
  j1   j1

y 2 y
Alors que (V-22) provient de l'approximation :
Vj1 j1  Vj1 j1
(V 1  (V 1 
j j 2
2 2

EXERCICES RESOLUS
Exercice 1:
U  2 U
L'équation   0 est approchée au point (i h, j k) par l'équation aux différences :
t x 2

Montrer que l'erreur de troncation locale en ce point est

Discuter la consistance de ce schéma avec l'E.D.P lorsque :

i) k= r h et ii) k= r h2 où r est une constante positive et  un paramètre variable.
Solution :
Le développement en séries de Taylor des termes Ui,j+1, Ui,j-1, Ui+1,j et Ui-1,j autour du point
(i h, j k) et la substitution dans :

mène a :

Adapted by A. Boukhari 53
El-Oued University 2013.
CHAPITRE 5 METHODE DES DIFFERENCES FINIES

d'où le résultat.

1erCas: k = r h
quand h  0 :

Quand   1/2 le troisième terme tend vers l'infini. Quand  =1/2 la valeur limite de Ti,j est :

Dans ce cas l'équation de différences finies est consistante avec l'équation hyperbolique :

U  2 U
Donc l'équation de différences finies est toujours inconsistante avec  quand k= r h.
t x 2
2ème Cas : k = r h2

Quand   1/2 le schéma de différences est consistant avec l'équation parabolique :

C'est seulement quand  =1/2 que le schéma de différence devient consistant avec
l'équation différentielle donnée. Ceci est alors le schéma bien connu de Dufort-Frankel, qui est
aussi stable pour toute valeur de r > 0.

Exercice 2 :
 2U 2 U
L'équation hyperbolique 2  2 est approchée en ( p h, q k) par le schéma de
t x
différence explicite :

Trouver la condition de stabilité de ce schéma.

Solution :

Adapted by A. Boukhari 54
El-Oued University 2013.
CHAPITRE 5 METHODE DES DIFFERENCES FINIES

Adapted by A. Boukhari 55
El-Oued University 2013.