12.3 Approximation par ´ el´ ements finis de l’´ equation de la chaleur
461
c’est-` a-dire,
Nh
j=1
. u j (t)
1
0
ϕ j ϕ i dx +
Nh
j=1
u j (t)a(ϕ j , ϕ i ) = F (ϕ i ),
i= 1, . . . , N h .
Avec les mˆ emes notations qu’en (12.9), on a
M
. u(t) + A fe u(t) = f fe (t),
(12.15)
o` u A fe = (a(ϕ j , ϕ i )), f fe (t) = (F (ϕ i )) et M = (m ij ) = (
1
0
ϕ j ϕ i dx) pour
i, j = 1, . . . , N h . La matrice M s’appelle matrice de masse . Comme elle est
inversible, on peut ´ ecrire le syst` eme d’EDO (12.15) sous forme normale :
. u(t) = −M
−1 A fe u(t) + M
−1 f fe (t).
(12.16)
Pour r´ esoudre (12.16) de mani` ere approch´ ee, on peut `
a nouveau appliquer le
θ-sch´ ema :
M
u
k+1
− u
k
∆t
+ A fe
θu
k+1 + (1 − θ)u
k
= θf
k+1
fe
+ (1 − θ)f
k
fe .
(12.17)
Comme d’habitude, l’exposant k indique que la quantit´ e est consid´ er´ ee au
temps t
k . De la mˆ eme mani` ere qu’avec les diff´ erences finies, on retrouve respectivement pour θ = 0, 1 et 1/2 les sch´ emas d’Euler explicite, implicite et le
sch´ ema de Crank-Nicolson. Seul ce dernier est du second ordre en ∆t.
Pour tout k, (12.17) est un syst` eme lin´ eaire associ´ e ` a la matrice
K =
1
∆t
M + θA fe .
Comme M et A fe sont sym´ etriques d´ efinies positives, la matrice K l’est aussi.
On peut donc effectuer, `
a t = 0, sa d´ ecomposition de Cholesky K = H
T H, o` u
H est une matrice triangulaire sup´ erieure (voir Section 3.4.2). On doit alors
r´ esoudre ` a chaque pas de temps les deux syst` emes lin´ eaires triangulaires de
taille de N h
⎧
⎨
⎩
H
T y =
1
∆t
M − (1 − θ)A fe
u
k + θf
k+1
fe
+ (1 − θ)f
k
fe ,
Hu
k+1 = y,
ce qui repr´ esente un coˆ ut de calcul de l’ordre N
2
h /2 op´ erations par syst` eme.
Quand θ = 0, une technique appel´ ee condensation de la masse (mass-lumping
en anglais), consistant `
a approcher M par une matrice diagonale inversible
M,
permet de d´ ecoupler le syst` eme d’´ equations (12.17). Dans le cas d’´ el´ ements
finis affines par morceaux, on obtient la matrice
M en appliquant la formule
composite du trap` eze aux noeuds {x i } pour ´ evaluer les int´ egrales
1
0
ϕ j ϕ i dx.
On obtient alors ˜
m ij = hδ ij , i, j = 1, . . . , N h (voir Exercice 2).
461
c’est-` a-dire,
Nh
j=1
. u j (t)
1
0
ϕ j ϕ i dx +
Nh
j=1
u j (t)a(ϕ j , ϕ i ) = F (ϕ i ),
i= 1, . . . , N h .
Avec les mˆ emes notations qu’en (12.9), on a
M
. u(t) + A fe u(t) = f fe (t),
(12.15)
o` u A fe = (a(ϕ j , ϕ i )), f fe (t) = (F (ϕ i )) et M = (m ij ) = (
1
0
ϕ j ϕ i dx) pour
i, j = 1, . . . , N h . La matrice M s’appelle matrice de masse . Comme elle est
inversible, on peut ´ ecrire le syst` eme d’EDO (12.15) sous forme normale :
. u(t) = −M
−1 A fe u(t) + M
−1 f fe (t).
(12.16)
Pour r´ esoudre (12.16) de mani` ere approch´ ee, on peut `
a nouveau appliquer le
θ-sch´ ema :
M
u
k+1
− u
k
∆t
+ A fe
θu
k+1 + (1 − θ)u
k
= θf
k+1
fe
+ (1 − θ)f
k
fe .
(12.17)
Comme d’habitude, l’exposant k indique que la quantit´ e est consid´ er´ ee au
temps t
k . De la mˆ eme mani` ere qu’avec les diff´ erences finies, on retrouve respectivement pour θ = 0, 1 et 1/2 les sch´ emas d’Euler explicite, implicite et le
sch´ ema de Crank-Nicolson. Seul ce dernier est du second ordre en ∆t.
Pour tout k, (12.17) est un syst` eme lin´ eaire associ´ e ` a la matrice
K =
1
∆t
M + θA fe .
Comme M et A fe sont sym´ etriques d´ efinies positives, la matrice K l’est aussi.
On peut donc effectuer, `
a t = 0, sa d´ ecomposition de Cholesky K = H
T H, o` u
H est une matrice triangulaire sup´ erieure (voir Section 3.4.2). On doit alors
r´ esoudre ` a chaque pas de temps les deux syst` emes lin´ eaires triangulaires de
taille de N h
⎧
⎨
⎩
H
T y =
1
∆t
M − (1 − θ)A fe
u
k + θf
k+1
fe
+ (1 − θ)f
k
fe ,
Hu
k+1 = y,
ce qui repr´ esente un coˆ ut de calcul de l’ordre N
2
h /2 op´ erations par syst` eme.
Quand θ = 0, une technique appel´ ee condensation de la masse (mass-lumping
en anglais), consistant `
a approcher M par une matrice diagonale inversible
M,
permet de d´ ecoupler le syst` eme d’´ equations (12.17). Dans le cas d’´ el´ ements
finis affines par morceaux, on obtient la matrice
M en appliquant la formule
composite du trap` eze aux noeuds {x i } pour ´ evaluer les int´ egrales
1
0
ϕ j ϕ i dx.
On obtient alors ˜
m ij = hδ ij , i, j = 1, . . . , N h (voir Exercice 2).