Cours de la Méthode des Différences Finies

Méthode des différences finies

(Finite Difference Method)
Plan de cours :
I- Introduction à la méthode de différences finies (MDF)
1 - Définition de la méthode ;
2 – Discrétisation et les conditions aux limites;
3 – Points (Nœuds) réguliers ;
4 - Points (Nœuds) irréguliers ;
II – Passage d’une équation aux dérives partielles de type elliptique à une équation
différences finies :
1) – Cas régulier ; 2) - Cas irréguliers ;
III – Etude du système algébrique résultant :
IV – Application de la MDF aux équations paraboliques
V – Application de la MDF aux équations hyperboliques

I.1 – Définition de la méthode

L’idée de cette méthode numérique la plus universelle est de faire remplacer
approximativement l’équation aux dérivées partielles (équation différentielle) par un système
d’équations algébriques résultant d’un développement en séries de Taylor de la fonction à
approximer. La méthode des différences finies (MDF) est basée essentiellement sur la
discrétisation du domaine d’un dispositif (Composant, machine, etc.) à un ensemble de nœuds,
appelées maillage (grille ou réseau).
Avantages : Grande simplicité d'écriture et faible coût de calcul ;
Inconvénient : Limitation à des géométries simples.
I.2 – Discrétisation et les conditions aux limites : C’est une opération de transformation
de la géométrie originale du dispositif à un ensemble de nœuds (grille) ou le dispositif doit garder
la forme initiale (originale), ceci se fait par le traçage des lignes parallèles aux axes x et y et leur
intersection donne les nœuds du maillage. Nœuds (x,y) Frontière
Г
y y Lx
ym Frontière Г yn

Domaine
D Ly

hy
xn x hx xm x

Fig.1: Géométrie originale Fig.2: Géométrie discrétisée

1
 Cours de la Méthode des Différences Finies
Suivant la fig.2, on a :
Lx
Lx  (9  1).h x  (m  1).h x  h x  (1)
m 1
Ly
L y  (6  1).h y  (n  1).h y  h y  (2)
n 1
Les équations se résolvent généralement dans un domaine fermé, pour garantir l’unicité de
la solution, il convient de fournir des conditions aux limites de Dirichlet (donner la fonction), de
Neumann (fournir sa dérivée normale) ou mixtes (combinaisons
linéaires des deux). Point irrégulier
I.3 – Points réguliers : Ce sont des points (nœuds) qui se
trouvent à l’intérieur de la géométrie loin des frontières. Le calcul
numérique de la fonction dans ce genre de points, se fait par
l’application directe de l’équation aux différences finies.
I.4 – Points irréguliers : Ce sont des points (noeuds) qui se
trouvent proche de la frontière avec une distance égal à une
fraction du pas ‘hx ou hy’. Le calcul numérique de la fonction dans ce genre des noeuds demande
un traitement particulier.
II. Passage d’une équation aux dérivées partielles de type elliptique à une équation de
différences finies
II.1 Cas régulier
On suppose que notre problème est gouverné par le modèle mathématique suivant :
 2 ( x, y)  2 ( x, y)
 0 (3) [Modèle de Laplace]
x 2 y 2

Pour la transformation de l’équation aux dérivées partielles à une EDF, on utilise le

développement de Taylor. On arrête de développement à l’ordre 4 :

4 φ(x,y+hy)

3 1 2 φ(x-hx,y) φ(x+hx,y)
φ(x,y)
φ(x,y)
φ(x,y-hy)
5
( x, y) h 2x  2 ( x, y) h 3x  3 ( x, y) h 4x  4 ( x, y)
( x  h x , y)  ( x, y)  h x .  .  .  .  (h 5x ) (4)
x 2!  x
2
3!  x
3
4!  x
4

( x , y) h 2x  2 ( x , y) h 3x  3 ( x, y) h 4x  4 ( x , y)
( x  h x , y)  ( x, y)  h x .  .  .  .  (h 5x ) (5)
x 2!  x
2 3!  x
3 4!  x
4

Où : hx est suffisamment petit  h 4x  0  (h 5x )  0 .

En additionnant ces deux dernières relations (4 et 5) et en divisant par h 2x , on obtient:
 2 ( x, y)
.( x  h x , y)  ( x  h x , y)  2( x, y)
1
 (6)
 x
2
h 2x
2
 Cours de la Méthode des Différences Finies
On refait la même chose pour la variable y :
2 3 4
( x, y) h y  2 ( x, y) h y  3 ( x, y) h y  4 ( x, y)
( x, y  h y )  ( x, y)  h y .  .  .  .  (h 5y ) (7)
y 2!  y
2 3!  y
3 4!  y
4

2 3 4
( x, y) h y  2 ( x, y) h y  3 ( x, y) h y  4 ( x, y)
( x, y  h y )  ( x, y)  h y .  .  .  .  (h 5y ) (8)
y 2!  y
2 3!  y
3 4!  y
4

Où : hy est suffisamment petit  h 4y  0  (h5y )  0 .

En additionnant aussi ces deux dernières relations (7 et 8) et en divisant par h 2y , on obtient:

 2 ( x, y)

1

. ( x, y  h y )  ( x, y  h y )  2( x, y)  (9)
 y
2
h 2y

On remplace les équations (6) et (9) dans (3) :

1
.( x  h x , y)  ( x  h x , y)  2( x, y) 
1
 
. ( x, y  h y )  ( x, y  h y )  2( x, y)  0 
h 2x h 2y

h 2x .h 2y  1 
( x , y)  . .( x  h x , y)  ( x  h x , y) 
1
. ( x , y  h 
y )  ( x , y  h y 
)  (10.a)
2(h 2x  h 2y )  h 2x h 2y 

L’équation (10.a) représente une équation aux différences finies à cinq points.

Si en prenant les pas h x  h y  h , on aura :

( x , y)  .( x  h , y)  ( x  h , y)  ( x , y  h )  ( x , y  h )
1
(10.b)
4
II.2 Cas irrégulier
Le modèle de Laplace :
 2 ( x, y)  2 ( x, y)
 0 (3) φ(x, y + b.h)
x 2 y 2
b.h
h a.h
2
( x, y) (ah )  2( x, y) φ(x-h, y) φ(x+ah, y)
( x  ah, y)  ( x, y)  (ah ).  .  ((ah )3 ) (11)
x 2!  x
2 h
φ(x, y - h)
2
( x, y) h  2( x, y)
( x  h, y)  ( x, y)  h.  .  (h 3 ) (12)
x 2!  x
2

En faisant (11) + a x (12) :

a (a  1)h 2  2( x, y)
(x  a.h, y)  a.(x  h, y)  (a  1).( x, y)  . (13)
2  2x
On peut écrire l’équation (13) comme suit :
 2( x, y) 2.( x  ah, y)  a.( x  h, y)  (a  1).( x , y)
 (14)
 x 2
a (a  1)h 2

3
 Cours de la Méthode des Différences Finies
De la même manière, on obtient :
 2( x, y) 2.( x, y  bh )  b.( x, y  h )  (b  1).( x, y)
 (15)
2y b(b  1)h 2

On remplace les équations (14) et (15) dans (3) :

2.( x  ah, y)  a.( x  h, y)  (a  1).( x, y) 2.( x, y  bh )  b.( x, y  h )  (b  1).( x, y)
 0 
a (a  1)h 2 b(b  1)h 2

ab  1 1 1 1 
(x, y)  . .(x  ah, y)  .(x  h, y)  .(x, y  bh )  (x, y  h) (16)
(a  b)  a (a  1) a 1 b(b  1) b 1 
On a :
(x  ah, y)  A et (x, y  bh)  B
On remplace dans l’équation (16) :
ab  1 1 1 1 
(x, y)  . .A  .(x  h, y)  .B  (x, y  h) (17)
(a  b)  a (a  1) a 1 b(b  1) b 1 
Alors pour des raisons de programmation, l’équation (17) s’écrit :
ab  1 1 1 1 
(i, j)  ( x, y)  . . A  .(i  1, j)  . B  (i, j  1) (18)
(a  b)  a (a  1) a 1 b(b  1) b 1 
En repassant en notation indicielle, la relation (10.b) peut s’écrire sous la forme suivante :

φ(i+1, j)
φ(x, y+h)

φ(x, y) φ(i, j-1) φ(i, j) φ(i, j+1)

φ (x-h, y) φ(x-h, y)

φ(i-1, j)
φ (x, y-h)
Fig.4- Sous notation indicielle
Fig.3

(i, j)  .(i  1, j)  (i  1, j)  (i, j  1)  (i, j  1)

1
(19)
4
III. Etude du système algébrique résultant
Le système d’équations matriciel du problème aux différences finies est donné par :
[ K].[F] = [B]. (20)
Avec : [K] : La matrice dépendant du pas (h) ;
[F] : Le vecteur des inconnus aux nœuds ;
[B] : Le vecteur tenant compte des conditions aux limites.
Parmi les méthodes de résolution de l’équation (20) disponibles, ou trouve les méthodes
directes qui conduisent à une solution en un nombre d’opérations connus à l’avance et les
méthodes indirectes qui donne une solution par un processus itératif.
Pour les méthodes directes, les plus utilisées sont la méthode de Gauss ou la méthode de
Cramer.

4
 Cours de la Méthode des Différences Finies
Parmi les méthodes indirectes on trouve celle de Newton, de Jacobi, de Gauss
Seidel, …etc.
Méthode itérative (indirecte) : Imaginer si un nombre de nœuds très important, le nombre
des équations sera très important aussi. C’est pourquoi on introduira cette méthode.

(i, j)  .(i  1, j)  (i  1, j)  (i, j  1)  (i, j  1)

1
(19)
4
Dans le cas où (i, j) représente la solution exacte, l’équation (19) s’écrite :

.(i  1, j)  (i  1, j)  (i, j  1)  (i, j  1)  (i, j)  0

1
(21)
4

Si (i, j) représente une approximation  l’équation (21) devient le résultat d’un résidu
(Rés) :

Rés (i, j)  .(i  1, j)  (i  1, j)  (i, j  1)  (i, j  1)  (i, j)

1
(22)
4
La fonction de (i, j) deviendra à la nième itération :

 n (i, j)  (n 1) (i, j)  Rés (i, j) (23)

La convergence est assurée par :

 n (i, j)   (n 1) (i, j)   (24)

Young et Frankel ont proposé une méthode pour accélérer la convergence :

 n (i, j)  (n 1) (i, j)  w.Rés (i, j) (25)

Avec : 1  w  2 : Facteur de convergence.

IV. Application de la MDF aux équations paraboliques

On suppose que notre problème est gouverné par le modèle mathématique
monodimensionnel suivant :
  2
 k. 2 (26)
t x
On transforme l’équation (26) à une équation aux différences finies à quatre nœuds (points).
Explicité
(i, j  1)  (i, j) (i  1, j)  2.(i, j)  (i  1, j) j+1
 k. (27)
t (x ) 2 j
Avec : x i  (i  1).x ; t j  ( j  1).t ; i  1, 2, 3, ..., n ; j  1, 2, 3, ..., m .
i-1 i i+1
Dans la figure ci-contre: X change / à t ; O change / à x.
t
Si on pose : r  k. .
(x) 2
5
 Cours de la Méthode des Différences Finies
L’équation (27) devient :
(i, j  1)  (i, j)  r.(i  1, j)  2.(i, j)  (i  1, j)
(28)
 r.(i  1, j)  (1  2.r ).(i, j)  r.(i  1, j)

Soit sous forme matricielle :

 (2, j  1)  1  2.r r 0  0 0 0   (2, j)   d 
 (3, j  1)   r 1  2.r r  0 0  
0   (3, j)   0
    
 (4, j  1)   0 r 1  2.r  0 0 0   (4, j)  0
      
          .    r.  
 (n  3, j  1)   0 0 0  1  2.r r 0   (n  3, j)  0
      
(n  2, j  1)  0 0 0  r 1  2.r r  (n  2, j) 0
 (n  1, j  1)   0  1  2.r   (n  1, j)   
   0 0 0 r  f
Où :  0  (i,1 ) : Condition initiale;  d  (1, j ) ; f  (n, j ) : Conditions aux limites aux
extrémités ; n : Nombre de nœud sur l’axe d’espace x.
A chaque itération, le vecteur des inconnues se détermine directement.
1 1 1
Pour avoir une solution stable, il faut que r  , les valeurs r  ou r  ont été utilisées
2 4 6
pour minimiser les erreurs du calcul et diminuer les oscillations.
Si on remplace j → j +1 :
Implicite
 
2
(i  1, j  1)  2.(i, j  1)  (i  1, j  1)
 (29)
x 2
(x) 2 j+1

L’expression devient implicite, on aura : j

(i, j)  r.(i  1, j  1)  (i  1, j  1) i-1 i
(i, j  1)  (30) i+1
(1  2r )

L’expression (30) s’applique à tous les points sauf au premier et au dernier, ou elle doit être
modifiée pour refléter les C.L.
On peut réécrit l’expression (30) sous forme matricielle suivante:
1  2.r r 0  0 0 0   (2, j  1)   (2, j)   d 
  r 1  2.r r  0 0    
0   (3, j  1)   (3, j)   0
  
 0  r 1  2.r  0 0 0   (4, j  1)   (4, j)  0
      
              r.  
 0 0 0  1  2.r r 0   (n  3, j  1)   (n  3, j)  0
      
 0 0 0   r 1  2.r  r  (n  2, j  1) (n  2, j) 0
 0   r 1  2.r   (n  1, j  1)   (n  1, j)   
 0 0 0  f
A chaque itération, le vecteur des inconnues se détermine par la résolution d’un système
linéaire.

6
 Cours de la Méthode des Différences Finies
La molécule du calcul de (30) écrite pour tous les nœuds intérieur, elle nous donnera un
système tri diagonal de (n-2) équations et (n-2) inconnus qu’on doit résoudre à chaque instant t.
Avec ce schéma, on peut choisir r de façon comme on veut.
Imaginer si un nombre de nœuds très important, le nombre des équations sera très important
aussi. C’est pourquoi on introduira cette méthode itérative.
(i, j)  r.(i  1, j  1)  (i  1, j  1)
(i, j  1)  (30)
(1  2r )

Dans le cas où (i, j) représente la solution exacte, l’équation (30) s’écrite :

(i, j)  r.(i  1, j  1)  (i  1, j  1)
 (i, j  1)  0 (31)
(1  2r )

Si (i, j) représente une approximation  l’équation (31) devient le résultat d’un résidu
(Rés) :
(i, j)  r.(i  1, j  1)  (i  1, j  1)
Rés (i, j)   (i, j  1) (32)
(1  2r )

La fonction de (i, j) deviendra à la nième itération :

 n (i, j)  (n 1) (i, j)  Rés (i, j) (33)

La convergence est assurée par :

 n (i, j)   (n 1) (i, j)   (34)

Crank – Nicholson ont proposé un schéma inconditionnellement

j+1
stable (Δt on le prend comme on veut de même pour Δx et r), dans ce
j + 1/2
1
schéma la dérivé première prendre j  (point intermédiaire) :
2 j
 (i, j  1)  (i, j)
 , (35) i-1 i i+1
t t
et la dérivée 2 (spatial) peut être déterminer par la moyenne aux instants j et j + 1.

 2 1  (i  1, j)  2.(i, j)  (i  1, j) (i  1, j  1)  2.(i, j  1)  (i  1, j  1) 

 .   (36)
x 2 2  (x ) 2 (x ) 2 

En remplaçant (35) et (36) dans (26), on obtient :

2.(1  r).(i, j  1)  r.(i  1, j  1)  (i  1, j  1)  2.(1  r).(i, j)  r.(i  1, j)  (i  1, j) (37)

Après quelques simplifications, on aura :

r.(i  1, j  1)  (i  1, j  1)  2.(1  r ).(i, j)  r.(i  1, j)  (i  1, j)
(i, j  1)  (38)
2.(1  r )

7
 Cours de la Méthode des Différences Finies
Valable pour tous les nœuds intérieurs sauf pour le premier et le dernier, où il faut tenir le
compte des C.L.
V. Application de la MDF aux équations hyperboliques
On suppose que notre problème est gouverné par le modèle mathématique
monodimensionnel suivant :

 2  2
 2. (39)
t 2 x 2
Avec : μ : Vitesse de propagation de l’onde.
En remplaçant les dérivées par les approximations centrales, on aura :
(i, j  1)  2.(i, j)  (i, j  1) (i  1, j)  2.(i, j)  (i  1, j)
 2.
(t ) 2
(x) 2

Avec : x i  (i  1).x ; t j  ( j  1).t ; i  1, 2, 3, ..., n ; j  1, 2, 3, ..., m .

Après la simplification de l’éq.(39), on aura :
(i, j  1)  2.(1  r).(i, j)  r.(i  1, j)  (i  1, j)  (i, j  1) (40)
2
 t 
Avec : r   2 . .
 x 
Pour la stabilité du schéma, il faut que r  1 .
Un schéma populaire pour r = 1 :
(i, j  1)  (i  1, j)  (i  1, j)  (i, j  1) (41)

On peut utiliser la forme matricielle suivante pour r = 1 :

 (2, j  1)  0 1 0  0 0 0  (2, j)  1 0 0  0 0 0  (2, j  1)   d 

 (3, j  1)  1
   0 1  0 0 0  (3, j)  0 1 0  0 0 0  (3, j  1)   0 
 (4, j  1)  0 1 0  0 0 0  (4, j)  0 0 1  0 0 0  (4, j  1)   0 
         
           .           .    
 (n  3, j  1)  0 0 0  0 1 0  (n  3, j)  0 0 0  1 0 0  (n  3, j  1)   0 
         
(n  2, j  1) 0 0 0  1 0 1 (n  2, j) 0 0 0  0 1 0 (n  2, j  1)  0 
 (n  1, j  1)  0
   0 0  0 1 0  (n  1, j)  0 0 0  0 0 1  (n  1, j  1)   f 

Avec :  d  (1, j) et f  (n, j) : Conditions aux limites.

C’est toujours un schéma explicité.

Ces deux schémas nécessitent la connaissance de φ à temps j et j – 1 pour calculer φ à j + 1,
ce sont des schémas à deux pas. On doit donc trouver comment démarrer la résolution de
ces schémas car à t = Δt, on doit trouver φ à t = - Δt.

Si on donne :  0 (Condition de Neumann), on utilise l’approximation centrale :
t t 0

8
 Cours de la Méthode des Différences Finies
 (i,2)  (i,0)
  0  (i, 0)  (i, 2)
t 2.t
On remplace dans l’équation (40) :
.(i  1, 1)  (i  1, 1)
r
(i, 2)  (1  r ).(i, 1)  (42)
2
Si on prend r = 1, l’équation (42) devient :
(i,2)  .(i  1, 1)  (i  1, 1)
1
(43)
2
On peut utiliser la forme matricielle suivante pour r = 1 :

 (2, 2)  0 1 0  0 0 0  (2, 1)   d 
 (3, 2)  1 0 1  0 0  
0  (3, 1)   0
    
 (4, 2)  0 1 0  0 0 0  (4, 1)  0
  1    1  
    2 .        .    2 .  
 (n  3, 2)  0 0 0  0 1 0  (n  3, 1)  0
      
(n  2, 2) 0 0 0  1 0 1 (n  2, 1) 0
 (n  1, 2)  0    
0  (n  1, 1)   
   0 0 0 1  f

Les équations (42) et (43) ne sont valides que pour t = 0, pour obtenir les autres valeurs
φ(i, j) on doit utiliser l’équation (40) ou (41).
D’autres méthodes implicites aussi existent, mais elles résultent d’un grand nombre
d’équations à résoudre et leurs utilisations nécessitent des simplifications.

N.B. :

Première dérivée :
( x )  (x) 
.( x  h x )  ( x ) ou bien  .(x i 1 )  (x i )
1 1
Dérivée avec un pas avant :   
x hx  x  i h x
( x )  (x) 
.( x )  ( x  h x ) ou bien  .( x i )  ( x i 1 )
1 1
Dérivée avec un pas arrière :   
x hx  x  i h x
( x )  (x) 
.( x  h x )  ( x  h x ) ou bien  .(x i 1 )  (x i 1 )
1 1
Dérivée centrale :   
x 2.h x  x  i 2.h x
Deuxième dérivée :
Le schéma d’ordre deux dit « centré » pour approximer la dérivée seconde de φ :

 2 ( x )   2 ( x ) 
 2 .( x  h x )  ( x  h x )  2( x ) ou bien  2   2 .( x i 1 )  ( x i 1 )  2( x i )
1 1
 x
2
hx   x i h x
Il existe aussi une formulation « avant » et « arrière » pour la dérivée seconde, toute deux d’ordre 1 :

 2 ( x )   2 ( x ) 
 2 .( x  2.h x )  ( x )  2( x  h x ) ou bien  2   2 .( x i  2 )  ( x i )  2( x i 1 )
1 1
 x
2
hx   x i h x

9
 Cours de la Méthode des Différences Finies
 ( x )
2
  ( x) 
2
 2 .( x )  ( x  2.h x )  2( x  h x ) ou bien  2   2 .( x i )  ( x i 2 )  2( x i 1 )
1 1
et
 x
2
hx   x i h x
Dérivées croisées

 2 ( x, y)
Déterminons une approximation de la dérivée croisée de la fonction de 2 variables
xy
( x, y) . La discrétisation du domaine de calcul est bidimensionnelle et fait intervenir deux pas d’espace
supposés constants hx et hy dans les directions x et y.
Le principe est aussi basé sur les développements de Taylor :

( x, y) ( x, y) h 2x  2 ( x, y) h y  2 ( x, y)  2 ( x, y)
2

( x  h x , y  h y )  ( x, y)  h x  hy  .  .  h .h .
x y 2x 2y xy
x y
2! 2!
( x, y) (x, y) h 2x  2 ( x, y) h y  2 ( x, y)  2 ( x, y)
2

( x  h x , y  h y )  ( x, y)  h x hy  .  .  h .h .
x y 2x 2y xy
x y
2! 2!

( x, y) ( x, y) h 2x  2 ( x, y) h y  2 ( x, y)  2 ( x, y)
2

( x  h x , y  h y )  ( x, y)  h x hy  .  .  h .h .
x y 2x 2y xy
x y
2! 2!

( x, y) ( x, y) h 2x  2 ( x, y) h y  2 ( x, y)  2 ( x, y)
2

( x  h x , y  h y )  ( x, y)  h x hy  .  .  h .h .
x y 2x 2y xy
x y
2! 2!

En réalisant une combinaison des quatre équations précédentes ((1) + (2) - (3) - (4)), on obtient une
approximation de la dérivée croisée de 1er ordre:

 2 ( x , y) ( x  h x , y  h y )  ( x  h x , y  h y )  ( x  h x , y  h y )  ( x  h x , y  h y )

xy 4.h x .h y

On réécrit l’équation précédente en notation indicielle:

 2 ( x , y) (i  1, j  1)  (i  1, j  1)  (i  1, j  1)  (i  1, j  1)


xy ( i , j)
4.h x .h y

Exemple simple 1D avec conditions de Dirichlet

Considérons l’équation différentielle suivante :
 u" ( x )  f ( x ), x   0, 1 

 u (0)  a et u (1)  b
Avec : f : Fonction continue.
Le maillage est construit en introduisant n + 1 nœuds xi avec i = 0, 1, 2, …, n, régulièrement
espacés avec un pas hx. La quantité ui indiquera la valeur de la fonction u(x) au nœud xi.
L’équation à résoudre s’écrit, sous forme discrète en chaque nœud xi :

2u
 u" ( x ) i    f (x i )  f i
x 2 i

10
 Cours de la Méthode des Différences Finies
u2
u  2u i  u i 1
On peut écrire :  u"(x) i     i 1
x i
2
h 2x

2u i  u i 1  u i 1
On réécrit l’équation à résoudre sous la forme suivante :  fi
h 2x

Pour les différentes valeurs de i, on obtient :

 2u 1  u 2  a
  f1
 h 2x
 2u 2  u 3  u 1  f
 2
 h 2x
 

 2u n 1  b  u n  2  f n 1
 h 2x

On peut utiliser la formulation matricielle suivante :

 2 1 0  0 0 0  u1   f 1  a / h 2x 
 
 1 2  1
  0 0 0   u 2   f2 
 0 1 2  0 0 0  u 3   f3 
   2 
           h x   
0 0 0  2  1 0   u n 3   f n 3 
    
0 0 0   1 2  1 u n  2   f n 2 
0  
 0  1 2   u n 1   2
 0 0 f n 1  b / h x 
Après la résolution du système d’équations, on obtient les valeurs: u1 ; u2 ; u3 ;… ; un-2 ; un-1.
Exemple simple 1D avec conditions de Dirichlet et Neumann
Considérons l’équation différentielle suivante :

 u" ( x )  f ( x ), x   0, 1 

u (0)  a C.L de Dirichlet
u ' (1)  b
 C.L de Neumann

Les modifications du problème discrétisé par rapport au cas précédent sont les suivantes. Tout
d'abord, le nombre d'inconnues a changé. Il y a une inconnue au bord en x = 1. Le problème discret a
donc maintenant, sur la base du même maillage que précédemment, n inconnues ui pour i variant de 1 à n.

En appliquant la condition de N avec l’utilisation d’approximation d’ordre 1 (dérivée avant) :

u u  un
u ' (1)   n 1 b  u n 1  b. h x  u n
x x 1 hx
  
Dérivée "avant"

L’équation à résoudre aussi s’écrit, sous forme discrète en chaque nœud xi:

11
 Cours de la Méthode des Différences Finies
 u 2
 u" ( x ) i    f (x i )  f i
x 2 i

 2u u  2u i  u i1
On peut écrire :  u" ( x ) i   2   i1
x h2
i x
Dérivée centrale

2u i  u i 1  u i 1
On réécrit l’équation à résoudre sous la forme suivante :  fi
h 2x

Pour les différentes valeurs de i avec l’application des CL, on obtient :

 2u 1  u 2  a
 2
 f1
 h x

 2u 2  u 3  u 1  f
 h 2x
2

 

 2u n  2  u n 1  u n 3
  f n 2
 h 2x

 2u n 1  u n  u n  2
  f n 1
h 2x

 2u n  u n 1  u n 1  2u n  u n  b.h x   u n 1  u n  b.h x  u n 1  f
 h 2x h 2x h 2x
n

On peut utiliser la formulation matricielle suivante :

 2 1 0  0 0 0  u1   f1  a / h 2x 
 1 2  1     
  0 0 0  u 2   f2 
 0 1 2  0 0 0  u 3   f3 
    
            h 2x   
0 0 0  2  1 0  u n  2   f n 2 
    
0 0 0   1 2  1  u n 1   f n 1 
0  
 0  1 1  u n   f  b / h x 
 0 0  n
Après la résolution du système d’équations, on obtient les valeurs : u1 ; u2 ; u3 ;… ; un-1 ; un.
Exemple sur la discrétisation de l’équation de la chaleur 2D stationnaire
On considère le problème bidimensionnel stationnaire de la conduction de la chaleur dans un
 
domaine rectangulaire 0, L x x 0, L y . Le champ de température T(x, y) vérifie l’équation de Laplace :

  2 T ( x , y)  2 T ( x , y)
T( x, y)   0 
et ( x, y)  0, L x  x 0, L y 
 x 2 y 2
T(0, y)  Tg et T(L x , y)  Td 0  y  Ly ,

T( x, 0)  Tb et T( x, L y )  Th 0  x  Lx

Le domaine de calcul est discrétisé en n x m nœuds (xi, yj) (i = 1, 2, …, n et j = 1, 2, …, m). On
pose que les pas d’espace ∆x et ∆y. La température discrète au nœud (xi, yj) sera notée Tij = T(xi, yj).
12
 Cours de la Méthode des Différences Finies
Les deux termes de l’équation de Laplace peuvent s’écrire comme suit :
 2T Ti 1, j  2Ti , j  Ti 1, j

x 2 i, j
x 2

 2T Ti , j1  2Ti , j  Ti , j1


y 2 i, j
y 2

La formation discrétisée est alors pour i variant de 2 à n-1 et j variant 2 à m-1. On peut écrire
l’équation de Laplace approximée comme suit :
 
y 2 Ti 1, j  Ti 1, j   x 2 Ti, j1  Ti, j1   2 x 2  y 2 Ti, j  0

En posant A  x 2  y 2 et cas d’exemple n = m = 5, on obtient un système d’équations sous

forme matricielle suivant :
 2A y 2 0 x 2 0 0 0 0 0  T22  x 2 Tb  y 2 Tg 
 2    
 y  2A y 2 0 x 2 0 0 0 0  T32   x 2 Tb 
 0 y 2
 2A 0 0 x 2
0 0 
0 T42   x Tb  y Td 
2 2

 2     
 x 0 0  2A y 2 0 x 2 0 0  T23   y 2 Tg 
 0 x 2
0 y 2
 2A y 2
0 x 2
0   
T33    0 
     
 0 0 x 2 0 y 2  2A 0 0 x 2  T43   y 2 Td 
 0 0 0 x 2
0 0  2A y 2
0  T   x Th  y Tg 
2 2
   24   
 0 0 0 0 x 2 0 y 2  2A y 2  T34   x 2 Th 
     
x 2 y 2  2A  T44  x Th  y Td 
2 2
 0 0 0 0 0 0

Géométrie discrétisée
y
Th Th Th
∆x
∆y

Tg Td
T24 T34 T44
Ly Tg Td
T23 T33 T43
Tg Td
T22 T32 T42

Tb Tb Tb x

Lx

Dans le cas où les pas sont identiques ∆x = ∆y, la formulation devient pour i variant de 2 à
n-1 et j variant 2 à m-1. On peut réécrire l’équation de Laplace approximée comme suit :
Ti 1, j  Ti 1, j  Ti , j1  Ti , j1  4.Ti , j  0

On peut obtenir le système d’équations sous forme matricielle suivant pour les valeurs
de n = m = 5 et des températures : Tb =10 oC; Tg = 20 oC; Td = 40 oC; Th = 30 oC.

13
 Cours de la Méthode des Différences Finies
 4 1 0 1 0 0 0 0 0  T22  Tb  Tg 
 1 4 1 0 1 0 0 0   
0  T32   T 
  b 
0 1 4 0 0 1 0 0 0  T42  Tb  Td 
    
1 0 0 4 1 0 1 0 0  T23   Tg 
0 1 0 1 4 1 0 1 0  T33     0 
    
0 0 1 0 1 4 0 0 1  T43   Td 
0 0 0 1 0 0 4 1 0  T24  T  T 
    h g

 0 0 0 0 1 0 1  4 1   34 
T  Th 
0     
 0 0 0 0 1 0 1  4 T44  Th  Td 

 4 1 0 1 0 0 0 0 0  T22  30 T22  17,8571 

 1 4 1 0 1 0 0 0   
0  T32  10  T  19,4643 
    32   
0 1 4 0 0 1 0 0 0  T42  50 T42  25,0000
        
1 0 0 4 1 0 1 0 0  T23  20 T23   21,9643 
 0 1 0 1 4 1 0 1 0  T33     0   T33   25,0000
        
0 0 1 0 1 4 0 0 1  T43  40 T43  30,5357 
0 0 0 1 0 0 4 1  
0 T24  50 T  25,0000
      24   
0 0 0 0 1 0 1  4 1  T34  30 T34  28,0357
0   
1  4 T44  70 T   32,1429 
 0 0 0 0 1 0    44   

Exemple 1:
Trouver la distribution du potentiel V(x, y) entre deux armatures d’un condensateur. Où le
problème est gouverné par le modèle mathématique de Laplace suivant :
 2 V ( x , y)  2 V ( x , y)
 0
x 2 y 2
Où les conditions aux limites sont indiquées sur la géométrie du problème (fig.a) et on
prend les pas hx= hy= h et les valeurs des tensions appliquées suivantes : V1  24 V ; V2  12 V .
a) Calculer manuellement pour les valeurs n = 5 et m = 4 ;
b) En utilisant les expressions (22), (23) et (24) et l’organigramme ci-dessous, écrire un
programme sous MATLAB pour résoudre ce problème pour les valeurs n = m = 30 et
ε = 10-10 ;
c) Représenter graphiquement les lignes équipotentielles et les lignes de champ électrique à
l’intérieur.
Exemple 2 :
Nous appliquons les principes précédents à la détermination du potentiel à l’intérieur du
coaxial rectangulaire suivant (fig. b) :
Le conducteur extérieur étant au potentiel -10 et le conducteur intérieur au potentiel +10.
Le potentiel à l’intérieur du coaxial obéit à l’équation de Laplace.
a. Calculer manuellement pour les valeurs n = m = 7 ;

14
 Cours de la Méthode des Différences Finies
b. Ecrire un programme sous MATLAB pour résoudre ce problème pour les valeurs
n = m = 60 et ε = 10-10 ;
c. Représenter graphiquement les lignes équipotentielles et les lignes de champ électrique
à l’intérieur du coaxial.
Avec : A / C = B / D = 3.
∂V/∂n = 0

V1 V2

B D

C
V1 V2
A
Fig. a Fig. b
Réponse de l’exemple 1:
a). En utilisant les équations (6) et (9), l’équation à résoudre est :
V( x , y)  .V( x  h , y)  V( x  h, y)  V( x , y  h )  V( x, y  h )
1
4
Ou bien sous forme indicielle :

V(i, j)  .V(i  1, j)  V(i  1, j)  V(i, j  1)  V(i, j  1)

1
4
Cette équation est valable pour tous nœuds intérieurs, dans la 1 ère et la dernière colonnes, on
introduit les conditions de Dirichlet et dans la 1ère et la dernière lignes, on les corrige pour tenir en compte
les conditions aux limites de Neumann.
V V V(i, j  1)  V(i, j  1)
   0  V(i, j  1)  V(i, j  1) (Variation sur l’axe y)
n y 2.h

Dérivée centrale

Les équations ci-dessous sont valables pour la 1ère et la dernière lignes :


V(i, 1)  4 .V(i  1, 1)  V(i  1, 1)  2.V(i, 2)
1

 V(i, m)  .V(i  1, m)  V(i  1, m)  2.V(i, m  1)
1
 4
En utilisant les équations ci-dessus dans les différends nœuds, on obtient le système d’équations
suivant:


V21  4 .V1  V31  2.V22 ; V31  4 .V21  V41  2.V32 ; V41  4 .V31  V2  2.V42 ;
1 1 1

 V22  .V1  V32  V21  V23 ; V32  .V22  V42  V31  V33 ; V42  .V32  V2  V41  V43 ;
1 1 1
 4 4 4

 V23  .V1  V33  V22  V24 ; V33  .V23  V43  V32  V34 ; V43  .V33  V2  V42  V44 ;
1 1 1
 4 4 4
 V  1 .V  V  2.V ; V  1 .V  V  2.V ; V  1 .V  V  2.V .
 24 4 1 34 23 34
4
24 44 33 44
4
34 2 43

15
 Cours de la Méthode des Différences Finies
∂V/∂n = 0 V1 V24 V34 V44 V2

h
V1 V23 V33 V43 V2
V1 V2 h
V1 V22 V32 V42 V2

V1 V21 V31 V41 V2

V1 V2
→ y
x
V1 V2

Géométrie originale Géométrie discrétisée

On peut réécrire le système ci-dessus sous forme matricielle :

 4 1 0 2 0 0 0 0 0 0 0 0   V21   V1 
 1 4 1 0 2 0 0 0 0 0 0   
0   V31  0
  
0 1 4 0 0 2 0 0 0 0 0 0   V41  V2 
    
1 0 0 4 1 0 1 0 0 0 0 0  V22   V1 
0 1 0 1 4 1 0 1 0 0 0 0   V32  0
    
0 0 1 0 1 4 0 0 1 0 0 0  V42  V 
0      2 
0 0 1 0 0 4 1 0 1 0 0 V23 V
    1
0 0 0 0 1 0 1 4 1 0 1 0   V33  0
0 0 0 0 0 1 0 1 4 0 0  
1 V43  V 
    2
0 0 0 0 0 0 2 0 0 4 1 0  V24   V1 
    
0 0 0 0 0 0 0 2 0 1  4 1   V34  0
 0 0 0 0 0 0 0 0 2 0  
1  4 V44  V2 

Après la résolution du système d’équations obtenu, on obtient les valeurs suivantes :

V21  V22  V23  21 V ; V31  V32  V33  18 V ; V41  V42  V43  15 V

b). Programme sous environnement MATLAB :

%%%% Programme de calcul du champ électrique par la MDF
n=input ('entrer le nombre de colonnes');
m=input ('entrer le nombre de ligne');
p=input ('entrer le nombre iterations');
b=input ('entrer la valeur de la tention1');
a=input ('entrer la valeur de la tention2');
v=zeros(m,n); res=zeros(m,n);
%Conditions de Dirichlet
v(:,1)=b ;v(:,n)=a ;
for k=1:p
%Condition de Neumann aux lignes (première et dernière)
for i=2:n-1
res(1,i)=1/4*[2*v(2,i)+v(1,i-1)+v(1,i+1)]-v(1,i);
res(m,i)=1/4*[2*v(m-1,i)+v(m,i-1)+v(m,i+1)]-v(m,i);
end
%Calcul de residu des noeuds intérieurs
for i=2:n-1
16
 Cours de la Méthode des Différences Finies
for j=2:m-1;
res(j,i)=1/4*[v(j,i+1)+v(j,i-1)+v(j+1,i)+v(j-1,i)]-v(j,i);
end
end
% Test de convergence
if max(abs(res))<1e-10;
break
end
v=v+res;
end
v
k
% Représentation graphique
[chx,chy]=gradient(v)
chx=-chx
chy=-chy
set(gcf,'color',[1 1 1])
contour(v,n)
hold on
colormap(hsv(n))
quiver(chx,chy)
figure (2)
mesh(v)% Traçage 3D
C). 30

25
25
20
20
15
15
10
10 30
5 0 20
10
10
5 10 15 20 25 30 20
30 0

Avec : Lignes horizontales de champs ; Lignes équipotentielles (verticales) (Figure 1).

Exemple 03 :
t
Résoudre l’équation parabolique suivante : 0,2
  2
 k. 2 0,1
t x x
à l’aide des trois schémas étudies, où k = 1, Δx = 2 cm 0 2 4 6 8
et Δt = 0,1 s. C.L. → (0, t )  100 et (8, t )  50 et
C.I. → ( x,0)  0 pour 0< x <8. Faire jusqu’à t  0,2 s .
Exemple 04 : t
Résoudre l’équation hyperbolique suivante: 0,50

 2  2 0,25
 2.
t 2
x 2 x

où: μ = 1, r = 1 et Δt=0,25. C.L. → (0, t )  (1, t )  0 et 0 0,25 0,50 0,75 1

C.I. → ( x,0)  sin x . Faire jusqu’à t  0,50 s . Comparer la réponse avec la solution exacte
( x, t )  sin x. cos t .

17
 Cours de la Méthode des Différences Finies

Début

Entrer le nombre de lignes n

Entrer le nombre de colonnes m
Entrer le nombre d’itérations p
Entrer la valeur de la tension appliquée V

Initialisation de la matrice du potentiel (V)

& la valeur de ε

Introduire les
conditions de Dirichlet

k=1

k = k+1 Introduire les conditions de

Neumann

V=V+rés Calcul du résidu aux points

intérieurs du maillage

Non
Max[abs( rés)] < ε

Oui

Lire (V)

Tracer les équipotentielles

et les lignes du champs

Fin

Organigramme principal de résolution du problème (Exemple 1)

18
 Cours de la Méthode des Différences Finies

Solution de l’exemple 2 :
V =[ -10.0000 -10.0000 -10.0000 -10.0000 -10.0000 -10.0000 -10.0000
-10.0000 -5.8333 -1.6667 -0.8333 -1.6667 -5.8333 -10.0000
-10.0000 -1.6667 10.0000 10.0000 10.0000 -1.6667 -10.0000
-10.0000 -0.8333 10.0000 10.0000 10.0000 -0.8333 -10.0000
-10.0000 -1.6667 10.0000 10.0000 10.0000 -1.6667 -10.0000
-10.0000 -5.8333 -1.6667 -0.8333 -1.6667 -5.8333 -10.0000
-10.0000 -10.0000 -10.0000 -10.0000 -10.0000 -10.0000 -10.0000 ]
Solution de l’exemple 3 :
Schéma explicite :

Distance (x) en cm
φ(x, t)
0 2 4 6 8

0 100 0 0 0 50
Temps (t)
en s 0,1 100 2,5000 0 1,2500 50

0,2 100 4,8750 0,0938 2,4375 50

Schéma implicite :

Distance (x) en cm
φ(x, t)
0 2 4 6 8

0 100 0 0 0 50
Temps (t)
en s 0,1 100 2,3830 0,0851 1,1925 50

0,2 100 4,6563 0,2475 2,3321 50

Solution de l’exemple 4 :

Distance (x) en cm
φ(x, t)
0 0,25 0,50 0,75 1

0 0 0,7071 1,000 0,7071 0

Temps (t)
en s 0,25 0 0,5000 0,7071 0,5000 0

0,50 0 0 0 0 0

19