Chapitre 4

Eléments finis de Lagrange

La construction d’ une méthode d’éléments finis nécessite la donnée d’un maillage, de noeuds et d’un espace de
polynômes, qui doivent être choisis de manière cohérente. Les éléments finis de type Lagrange font intervenir
comme ”degrés de liberté” (c.à.d. les valeurs qui permettent de déterminer entièrement une fonction) les valeurs de
la fonction aux noeuds. Ils sont très largement utilisés dans les applications. Il existe d’autres familles d’éléments
finis, comme par exemple les éléments finis de type Hermite qui font également intervenir les valeurs des dérivées
directionnelles. Dans le cadre de ce cours, nous n’aborderons que les éléments finis de type Lagrange, et nous
renvoyons aux ouvrages cités en introduction pour d’autres éléments.

4.1 Espace d’approximation

4.1.1 Cohérence “locale”
4.1.1 Cohérence “locale”
Soit T un maillage de Ω, pour tout élément K de T , on note ΣK l’ensemble des noeuds de l’élément. On suppose
que chaque élément a Nℓ noeuds K : ΣK = {a1 , . . . , aNℓ }, qui ne sont pas forcément ses sommets. On note P un
espace de dimension finie constitué de polynômes, qui définit la méthode d’éléments finis choisie.

Définition 4.1 (Unisolvance, élément fini de Lagrange) Soit K un élément et ΣK = (ai )i=1,...,Nℓ un ensemble
de noeuds de K. Soit P un espace de polynômes de dimension finie. On dit que le triplet (K, ΣK , P ) est un élément
fini de Lagrange si ΣK est P -unisolvant, c’est à dire si pour tout (α1 , . . . , αNℓ ) ∈ IR Nℓ , il existe un unique élément
f ∈ P tel que f (ai ) = αi ∀i = 1 . . . Nℓ . Pour i = 1, . . . , Nℓ , on appelle degré de liberté la forme linéaire ζi
définie par ζi (p) = p(ai ), pour tout p ∈ P . La propriété d’unisolvance équivaut à dire que la famille (ζi )i=1,...,Nℓ
forme une base de P ′ (espace dual de P ).

La P -unisolvance revient à dire que toute fonction de P est entièrement déterminée par ses valeurs aux noeuds.

Exemple : l’élément fini de Lagrange P1 Prenons par exemple, en dimension 1, l’élément K = [a1 , a2 ], avec
ΣK = {a1 , a2 }, et P = P1 (ensemble des polynômes de degré inférieur ou égal à 1). Le triplet (K, ΣK , P ) est
unisolvant s’il existe une unique fonction f de P telle que :
,
f (a1 ) = α1
f (a2 ) = α2

Or toute fonction f de P s’exprime sous la forme f (x) = λx + µ et le système


 λa1 + µ = α1

λa2 + µ = α2

détermine λ et µ de manière unique.

De même si on considère le cas d = 2. On prend comme élément K un triangle et comme noeuds les trois sommets,
a1 , a2 , a3 du triangle. Soit P = P1 = {f : IR → IR; f (x) = λx1 + µx2 + ν} l’ensemble des fonctions affines.
Alors le triplet (K, ΣK , P ) est un élément fini de Lagrange car f ∈ P est entièrement déterminée par f (a1 ), f (a2 )
et f (a3 ).

136
4.1. ESPACE D’APPROXIMATION CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

a3 a3

a1 a2 a1 a2
φ1 φ3

a3

a1 a2
φ2

F IGURE 4.1 – fonctions de base locales pour l’élément fini de Lagrange P1 en dimension 2

f (x)
 x
a1 a2

F IGURE 4.2 – Interpolée P 1 sur [a1 , a2 ] (en trait pointillé) d’une fonction régulière (en trait continu)

Définition 4.2 (Fonctions de base locales) Si (K, ΣK , P ) est un élément fini de Lagrange, alors toute fonction f
de P peut s’écrire :
%Nℓ
f= f (ai )fi
i=1

avec fi ∈ P et fi (aj ) = δij . Les fonctions fi sont appelées fonctions de base locales.

Pour l’élément fini de Lagrange P1 en dimension 2 considéré plus haut, les fonctions de base locales sont décrites
sur la figure 4.1
Définition 4.3 (Interpolée) Soit (K, ΣK , P ) un élément fini de Lagrange, et soit v ∈ C(K, IR). L’interpolée de v
est la fonction Πv ∈ P définie par :
Nℓ
%
Πv = v(ai )fi
i=1

On montre sur la figure 4.2 un exemple d’interpolée pour l’élément fini de Lagrange P1 en dimension 1. L’étude
de "v − Πv" va nous permettre d’établir une majoration de l’erreur de consistance d(u, HN ).
Remarque 4.4 Pour que le triplet (K, ΣK , P ) soit un élément fini de Lagrange, ilfaut, mais il ne suffit pas, que
dim P = cardΣK . Par exemple si P = P1 et qu’on prend comme noeuds du triangle deux sommets et le milieu
de l’arête joignant les deux sommets, (voir figure 4.3), (K, ΣK , P ) n’est pas un élément fini de Lagrange.

Proposition 4.5 (Critère de détermination) Soit (K, Σ, P ) un triplet constitué d’un élément, d’un ensemble de
noeuds et d’un espace de polynômes, tel que :

dim P = cardΣ = Nℓ (4.1)

Analyse numérique des EDP, M1 137 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.1. ESPACE D’APPROXIMATION CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

1 2 3

F IGURE 4.3 – Exemple de triangle à trois noeuds qui n’est pas un élément fini de Lagrange)

Alors
si ∃!f ∈ P ; f = 0 sur Σ (4.2)
ou si
∀i ∈ {1 . . . Nℓ }∃fi ∈ P fi (aj ) = δij (4.3)
alors (K, Σ, P ) est un élément fini de Lagrange.

Démonstration : Soit :
φ : P → IR Nℓ

f !→ (f (ai ))ti=1,Nℓ .
 f !→ (f (ai ))i=1,Nℓ .

L’application φ est linéaire de P dans IR Nℓ , et, par hypothèse card sΣ = dim P . Donc φ est une application
linéaire continue de P dans IR Nℓ , avec dimP = dim(IR Nℓ ) = Nℓ . Si (K, Σ, P ) vérifie la condition (4.2) alors φ
est injective. En effet, si φ(f ) = 0, alors f (ai ) = 0, ∀i = 1, . . . , Nℓ , et donc par hypothèse, f = 0. Donc φ est
une application linéaire, φ est injective de P dans IR Nℓ avec dimP = Nℓ . On en déduit que φ est bijective. Donc
toute fonction de P est entièrement déterminée par ses valeurs aux noeuds : (K, Σ, P ) est donc un élément fini de
Lagrange.
On montre facilement que si la condition (4.3) est vérifiée alors φ est surjective. Donc φ est bijective, et (K, Σ, P )
est un élément fini de Lagrange.

Proposition 4.6 Soit (K̄, Σ̄, P̄ ), un élément fini de Lagrange, où Σ̄ est l’ensemble des noeuds de K̄ et P̄ un espace
de fonctions de dimension finie, et soit F une bijection de K̄ dans K, où K est une maille d’un maillage éléments
finis. On pose Σ = F (Σ̄) et P = {f : K → IR; f ◦ F ∈ P̄ } (voir figure 4.4). Alors le triplet (K, Σ, P ) est un
élément fini de Lagrange.

a2 = F (â2 )
â3
F
K
(x, y) a3 = F (â3 )
a1 = F (â1 )
K̂
(x̂, ŷ)

â1 â2

F IGURE 4.4 – Transformation F

Analyse numérique des EDP, M1 138 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.1. ESPACE D’APPROXIMATION CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

Démonstration : Supposons que les hypothèses de la proposition sont réalisées. On veut donc montrer que (Σ, P )
est unisolvant. Soit Σ = (a1 , . . . , aNℓ ), et soit (α1 , . . . , αNℓ ) ∈ IR Nℓ . On veut montrer qu’il existe une unique
fonction f ∈ P telle que
f (ai ) = αi , ∀i = 1, . . . , Nℓ .
Or par hypothèse, (Σ̄, P̄ ) est unisolvant. Donc il existe une unique fonction f¯ ∈ P̄ telle que

f¯(āi ) = αi , ∀i = 1, . . . , Nℓ ,

(où (āi )i=1,...,Nℓ désignent les noeuds de K̄). Soit F la bijection de K̄ sur K, on pose f = f¯ ◦ F −1 . Or par
hypothèse, ai = F (āi ). On a donc : f (ai ) = f¯ ◦ F −1 (ai ) = f¯(āi ) = αi . On a ainsi montré l’existence de f telle
que f (ai ) = αi .
Montrons maintenant que f est unique. Supposons qu’il existe f et g ∈ P telles que :

f (ai ) = g(ai ) = αi , ∀i = 1, . . . , Nℓ .

Soit h = f − g on a donc :
h(ai ) = 0 ∀i = 1 . . . Nℓ .
On a donc h ◦ F (āi ) = h(ai ) = 0. Or h ◦ F ∈ P̄ , et comme (Σ̄, P̄ ) est unisolvant, on en déduit que h ◦ F = 0.
Comme, pour tout x ∈ K, on a h(x) = h ◦ F ◦ F −1 (x) = h ◦ F (F −1 (x)) = 0, on en conclut que h = 0.

Définition 4.7 (Eléments affine-équivalents) . Sous les hypothèses de la proposition 4.6, si la bijection F est
affine, on dit que les éléments finis (K̄, Σ̄, P̄ ) et (K, Σ, P ) sont affine–équivalents.

Remarque 4.8 Soient (K̄, Σ̄, P̄ ) et (K, Σ, P ) deux éléments finis afffine–équivalents. Si les fonctions de base
locales de K, . (resp. de K, , P ) sont affines, alors celles de (resp. ) le sont aussi, et on a :
Remarque 4.8 Soient (K, Σ, P ) et (K, Σ, P ) deux éléments finis afffine–équivalents. Si les fonctions de base
locales de (K̄, Σ̄, P̄ ). (resp. de (K, Σ, P )) sont affines, alors celles de K (resp. K̄) le sont aussi, et on a :

 f¯i = fi ◦ F,
i = 1, . . . , cardΣ

fi = f¯i ◦ F −1 ,

La preuve de cette remarque fait l’objet de l’exercice 55.

Proposition 4.9 (Interpolation) Sous les hypothèses de la proposition 4.10 page 140, soient ΠK̄ et ΠK les
opérateurs d’interpolation respectifs sur K̄ et K, voir définition 4.3 page 137. Soient v ∈ C(K, IR), ΠK̄ v et
ΠK v les interpolées respectives de v sur (K̄, P̄ ) et (K, P ), alors on a :

ΠK v ◦ F = ΠK̄ (v ◦ F )

Démonstration : Remarquons tout d’abord que ΠK v ◦ F et ΠK̄ (v ◦ F ) sont toutes deux des fonctions définies de
K̄ à valeurs dans IR, voir figure 4.5. Remarquons ensuite que, par définition de l’interpolée, ΠK v ∈ P . Comme

F
K K

ΠK v

ΠK v
F ◦ ΠK v

IR

F IGURE 4.5 – Opérateurs d’interpolation ΠK̄ etΠK

Analyse numérique des EDP, M1 139 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.1. ESPACE D’APPROXIMATION CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

(K̄, Σ̄, P̄ ) est l’élément de référence, on a donc :

ΠK v ◦ F ∈ P̄

On a aussi, par définition de l’interpolée : ΠK̄ (v ◦ F ) ∈ P̄ . On en déduit que ΠK v ◦ F et ΠK̄ (v ◦ F ) sont toutes
deux des fonctions de P̄ . Comme l’élément (K̄, P̄ , Σ̄) est unisolvant (car c’est un élément fini de Lagrange), toute
fonction de P̄ est uniquement déterminée par ses valeurs aux noeuds de Σ̄. Pour montrer l’égalité de ΠK v ◦ F et
ΠK̄ (v ◦ F ), il suffit donc de montrer que :

ΠK̄ (v ◦ F )(āi ) = ΠK v ◦ F (āi ), i = 1, . . . , Nℓ ,

où Nℓ = cardΣ̄. Décomposons ΠK̄ (v ◦ F ) sur les fonctions de base locales (f¯j ), j = 1, . . . , Nℓ . On obtient :
Nℓ
%
ΠK̄ (v ◦ F )(āi ) = v ◦ F (āj )f¯j (āi ).
j=1

On a donc :  
Nℓ
%
ΠK̄ (v ◦ F )(āi ) = v ◦  F (āj )f¯j  (āi ) = v ◦ F (āi ) = v(ai ).
j=1

Mais on a aussi :
ΠK v ◦ F (āi ) = ΠK v(F (āi )) = ΠK v(ai ) = v(ai ).
D’où l’égalité.

4.1.2 Construction de H et conformité

4.1.2 Construction de HN et conformité
Nous allons considérer deux cas : le cas où l’espace H est l’espace H 1 tout entier, et le cas où l’espace H est
l’espace H01

Cas H = H 1 (Ω)
Plaçons-nous ici dans le cas où H = H 1 (Ω), où Ω ⊂ IR d est un ouvert borné polygonal (si d = 2, polyèdrique
si d = 3). Soit T un maillage éléments finis, avec T = (Kℓ )ℓ=1,...,L , où les éléments finis Kℓ sont fermés et
tels que ∪Lℓ=1 Kℓ = Ω̄. Soit S = (Si )i=1,...,M l’ensemble des noeuds du maillage éléments finis, avec Si ∈
Ω̄, ∀i = 1, . . . , M.. On cherche à construire une méthode d’éléments finis de Lagrange ; donc à chaque élément
Kℓ , ℓ = 1, . . . , L, est associé un ensemble de noeuds Σℓ = S ∩ Kℓ , et un espace Pℓ de polynômes. On veut que
chaque triplet (Kℓ , Σℓ , Pℓ ) soit un élément fini de Lagrange. On définit les fonctions de base globales (φi )i=1,...,M ,
par :
φi |Kℓ ∈ Pℓ ∀i = 1, . . . , M ; ∀ℓ = 1; . . . , L, (4.4)
et
φi (Sj ) = δij ∀i = 1, . . . , M, ∀j = 1, . . . , M. (4.5)
Chaque fonction φi est définie de manière unique, grâce au caractère unisolvant de (Kℓ , Σℓ , Pℓ ), ℓ = 1, . . . , M .
On pose HN = V ect(φ1 , . . . , φM ). Pour obtenir une méthode d’éléments finis conforme, il reste à s’assurer que
HN ⊂ H 1 .
Une manière de construire l’espace HN est de construire un maillage à partir d’un élément de référence, grâce à la
proposition suivante, qui se déduit facilement de la proposition 4.6 page 138

Proposition 4.10 (Elément fini de référence) Soit T un maillage constitué d’éléments K. On appelle élément
fini de référence un élément fini de Lagrange (K̄, Σ̄, P̄ ), où Σ̄ est l’ensemble des noeuds de K̄ et P̄ un espace de
fonctions, de dimension finie, tel que, pour tout autre élément K ∈ T , il existe une bijection F : K̄ → K telle
que Σ = F (Σ̄) et P = {f : K → IR; f ◦ F ∈ P̄ } (voir figure 4.4). Le triplet (K, Σ, P ) est un élément fini de
Lagrange.

Analyse numérique des EDP, M1 140 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.1. ESPACE D’APPROXIMATION CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

Proposition 4.11 (Critère de conformité, cas H 1 ) Soit Ω un ouvert polygonal (ou polyèdrique) de IR d , d = 2 ou
3. Soit T = (Kℓ )ℓ=1,...,L , un maillage éléments finis de Ω, S = (Si )i=1,...,M l’ensemble des noeuds de maillage.
On se place sous les hypothèses de la proposition 4.10 ; soient (φi )i=1,...,M les fonctions de base globales, vérifiant
(4.4) et (4.5), et on suppose de plus que les hypothèses suivantes sont vérifiées :
Pour toute arête (ou face si d = 3) ǫ = Kℓ1 ∩ Kℓ2 , on a : Σℓ1 ∩ ǫ = Σℓ2 ∩ ǫ et Pℓ1 |ǫ = Pℓ2 |ǫ , (4.6)
où Pℓ1 |ǫ (resp. Pℓ2 |ǫ ) désigne l’ensemble des restrictions des fonctions de Pℓ1 (resp. Pℓ2 ) à ǫ),
Si ǫ est un côté de Kℓ , (Σℓ ∩ ǫ, Pℓ |ǫ ) est unisolvant. (4.7)
1
Alors on a : HN ⊂ C(Ω̄) et HN ⊂ H (Ω). On a donc ainsi construit une méthode d’éléments finis conformes.
(Notons que les côtés de Kℓ sont des arêtes en 2D et des faces en 3D.)
Démonstration : Pour montrer que HN ⊂ C(Ω̄) et HN ⊂ H 1 (Ω), il suffit de montrer que pour chaque fonction
de base globale φi , on a φi ∈ C(Ω̄) et φi ∈ H 1 (Ω). Or par hypothèse, (4.4), chaque fonction φi est polynômiale
par morceaux. De plus, grâce à l’hypothèse (4.6), on a raccord des polynômes sur les interfaces des éléments, ce
qui assure la continuité de φi . Il reste à montrer que φi ∈ H 1 (Ω) pour tout i = 1, . . . , M . Comme φi ∈ C(Ω̄), il
est évident que φi ∈ L2 (Ω) (car Ω est un ouvert borné, donc φi ∈ L∞ (Ω) ⊂ L2 (Ω).
Montrons maintenant que les dérivées faibles Dj φi , j = 1, . . . , d, appartiennent à L2 (Ω). Par définition, la
fonction φi admet une dérivée faible dans L2 (Ω) s’il existe une fonction ψi,j ∈ L2 (Ω) telle que :
( (
φi (x)∂j ϕ(x)dx = − ψij (x)ϕ(x)dx, (4.8)
Ω Ω

pour toute fonction ϕ ∈ Cc1 (Ω) (on rappelle que Cc1 (Ω) désigne l’ensemble des fonctions de classe C 1 à support
L

compact, et que ∂j désigne la dérivée classique par rapport à la j-ème variable). Or, comme Ω̄ = Kℓ , on a :
ℓ=1
 =1
( L (
%
φi (x)Dj ϕ(x)dx = φi (x)Dj ϕ(x)dx.
Ω ℓ=1 Kℓ

Sur chaque élément Kℓ , la fonction φi est polynômiale. On peut donc appliquer la formule de Green, et on a :
( ( (
φi (x)∂j ϕ(x)dx = φi (x)ϕ(x)nj (x)dγ(x) − ∂j φi (x)ϕ(x)dx,
KL ∂Kℓ Kℓ

où nj (x) est la j-ième composante du vecteur unitaire normal à ∂Kℓ en x, extérieur à Kℓ . Mais, si on note Eint
l’ensemble des arêtes intérieures du maillage (i.e. celles qui ne sont pas sur le bord), on a :
%L ( (
X= φi (x)ϕ(x)nj (x)dγ(x) = φi (x)ϕ(x)nj (x)dγ(x)
ℓ=1 ∂Kℓ ∂ Ω̄

% (   
+ (φi (x)ϕ(x)nj (x)) Kℓ1 + (φi (x)ϕ(x)nj (x))Kℓ2 dγ(x).
ǫ∈Eint
où Kℓ1 et Kℓ2 désignent les deux éléments dont ǫ est l’interface.
Comme ϕ est à support compact, (
φi (x)ϕ(x)nj (x)dγ(x) = 0.
∂ Ω̄
 
Comme φi et ϕ sont continues et comme nj (x) Kℓ1 = −nj (x)K pour tout x ∈ ǫ, on en déduit que X = 0. En
ℓ2
reportant dans (4.1.2), on obtient donc que :
( %L (
φi (x)∂j ϕ(x)dx = − ∂j φi (x)ϕ(x)dx.
Ω ℓ=1 Kℓ

Soit ψi,j la fonction de Ω dans IR définie presque partout par



ψij  ◦ = −∂j φi .
Kℓ

Comme ∂j φi est une fonction polynômiale par morceaux, on a ψi,j ∈ L2 (Ω) qui vérifie (4.8), ce qui termine la
démonstration.

Analyse numérique des EDP, M1 141 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.2. EXEMPLES CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

Cas H = H01 (Ω)

Plaçons-nous mainteant dans le cas où H = H01 (Ω). On décompose alors l’ensemble S des noeuds du maillage :

S = Sint ∪ Sext

où
Sint = {Si , i = 1, . . . , N } ⊂ Ω
est l’ensemble des noeuds intérieurs à Ω et

Sext = {Si , i = N + 1, . . . , M } ⊂ ∂Ω

est l’ensemble des noeuds de la frontière. Les fonctions de base globales sont alors les fonctions φi , i = 1, . . . , N
telles que
φi |Kℓ ∈ Pℓ , ∀i = 1, . . . , N, ∀ℓ = 1, . . . , L (4.9)
φi (Sj ) = δij , ∀j = 1, . . . , N, (4.10)
et on pose là encore HN = V ect{φ1 , . . . , φN }. On a alors encore le résultat suivant :

Proposition 4.12 (Critère de conformité, cas H01 ) Soit Ω un ouvert polygonal (ou polyèdrique) de IR d , d = 2 ou
3. Soit T = (Kℓ )ℓ=1,...,L un maillage éléments finis de Ω, S = (Si )i=1,...,M = Sint ∪ Sext l’ensemble des noeuds
du maillage. On se place sous les hypothèses de la proposition 4.6. On suppose que les fonctions de base globale
(φi )i=1,...,M vérifient (4.9) et (4.10), et que les conditions (4.6) et (4.7) sont vérifiées. Alors on a : HN ⊂ C(Ω̄) et
HN ⊂ H01 (Ω)

Démonstration : La preuve de cette proposition est laissée à titre d’exercice.

 emonstration : La preuve de cette proposition est laissée à titre d’exercice.

Remarque 4.13 (Eléments finis conformes dans H 2 (Ω)) On a construit un espace d’approximation HN inclus
dans C(Ω̄). En général, on n’a pas HN ⊂ C 1 (Ω̄), et donc on n’a pas non plus HN ⊂ H 2 (Ω) (en dimension
1 d’espace, H 2 (Ω) ⊂ C 1 (Ω)). Même si on augmente le degré de l’espace des polynômes, on n’obtiendra pas
l’inclusion HN ⊂ C 1 (Ω̄). Si on prend par exemple les polynômes de degré 2 sur les éléments, on n’a pas de
condition pour assurer le raccord, des dérivées aux interfaces. Pour obtenir ce raccord, les éléments finis de
Lagrange ne suffisent pas : il faut prendre des éléments de type Hermite, pour lesquels les degrés de liberté ne
sont plus seulement les valeurs de la fonction aux noeuds, mais aussi les valeurs de ses dérivées aux noeuds. Les
éléments finis de Hermite seront par exemple bien adaptés à l’approximation des problèmes elliptiques d’ordre 4,
dont un exemple est l’équation :
∆2 u = f dans Ω
où Ω est un ouvert borné de IR 2 , ∆2 u = ∆(∆u), et avec des conditions aux limites adéquates, que nous ne
détaillerons pas ici. On peut, en fonction de ces conditions aux limites, trouver un espace de Hilbert H et une
formulation faible de (4.13), qui s’écrit :
 ( (

 ∆u(x)∆ϕ(x)dx = f (x)ϕ(x)dx
Ω Ω

 u ∈ H, ∀ϕ ∈ H.

Pour que cette formulation ait un sens, il faut que ∆u ∈ L2 (Ω) et ∆ϕ ∈ L2 (Ω), et donc que H ⊂ H 2 (Ω). Pour
construire une approximation par éléments finis conforme de ce problème, il faut donc choisir HN ⊂ H 2 (Ω), et le
choix des éléments finis de Hermite semble donc indiqué.

4.2 Exemples
Pour chaque méthode d’élément fini de Lagrange, on définit :
1. un élément de référence K̄
2. des fonctions de base locales sur K̄
3. une bijection Fℓ de K sur Kℓ , pour ℓ = 1, . . . , L, où L est le nombre déléments du maillage.

Analyse numérique des EDP, M1 142 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.2. EXEMPLES CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

4.2.1 Elément fini de Lagrange P 1 sur triangle (d = 2)

Le maillage du domaine est constitué de L triangles (Kℓ )ℓ=1,...,L , et les polynômes d’approximation sont de degré
1.
Elément fini de référence : on choisit le triangle K̄ de sommets (0, 0), (1, 0) et (0, 1), et P̄ = {ψ : K →
IR(x, y) !→ ax + by + c, (a, b, c) ∈ IR 3 }.
Proposition 4.14 (Unisolvance) Soit Σ̄ = (āi )i=1,2,3 avec ā1 = (0, 0), ā2 = (1, 0) et ā3 = (0, 1), et
P̄ = {ψ; K → IR; (x, y) !→ a + bx + cy, (a, b, c) ∈ IR 3 }
Alors le couple (Σ̄, P̄ ) est unisolvant.
Démonstration : Soit (α1 , α2 , α3 ) ∈ IR 3 , et ψ ∈ P̄ . On suppose que ψ(āi ) = αi , i = 1, 2, 3. La fonction ψ est
de la forme ψ(x, y) = a + bx + cy et on a donc :

 a = α1
a + b = α2

a + c = α3
d’où c = α1 , b1 = α2 − α1 et b2 = α3 − α2 . La connaissance de ψ aux noeuds (āi )i=1,2,3 détermine donc
entièrement la fonction ψ.
Fonctions de bases locales.
Les fonctions de base locales sur l’élément fini de référence K̄ sont définies par φ̄i ∈ P̄ φ̄i (āj ) = δij , ce qui
détermine les φ̄ ; de manière unique, comme on vient de le voir. Et on a donc

 φ̄1 (x̄, ȳ) = 1 − x̄ − ȳ
φ̄2 (x̄, ȳ) = x̄

φ̄3 (x̄, ȳ) = ȳ.
 
(x̄, ) = ȳ.
Transformation Fℓ
On construit ici une bijection affine qui transforme K̄ le triangle de référence en un autre triangle K du maillage.
On cherche donc ℓ : K̄ → K, telle que
Fℓ (āi ) = ai i = 1, . . . , 3
où Σ = (ai )i=1,2,3 est l’ensemble des sommets de K. Notons (xi , yi ) les coordonnées de ai , i = 1, 2, 3. Comme
Fℓ est une fonction affine de IR 2 dans IR 2 , elle s’écrit sous la forme.
Fℓ (x̄, ȳ) = (β1 + γ1 x̄ + δ1 ȳ, β2 + γ2 x̄ + δ2 ȳ)
et on cherche βi , γi , δi , i = 1, 2 tels que :

 Fℓ ((0, 0)) = (x1 , y1 )
Fℓ ((1, 0)) = (x2 , y2 )

Fℓ ((0, 1)) = (x3 , y3 ).
Une résolution de système élémentaire amène alors à :
$ '
x1 + (x2 − x1 )x̄ + (x3 − x1 )ȳ
Fℓ (x̄, ȳ) =
y1 + (y2 − y1 )x̄ + (y3 − y1 )ȳ
D’après la remarque 4.8 page 139, si on note φ̄k , k = 1, 2, 3 les fonctions de base locales de l’élément de référence
(ℓ) (ℓ)
(K̄, Σ̄, P̄ ), et φk , k = 1, 2, 3 les fonctions de base locales de l’élément (Kℓ , Σℓ , Pℓ ), on a φk = φ̄k ◦ Fℓ−1
Si on note maintenant (φi )i=1,...,N les fonctions de base globales, on a :

 (ℓ)
φi Kℓ = φk ,

où i = ng(ℓ, k) est l’indice du k-ième noeud de l’élément ℓ dans la numérotation globale. Notons que l’élément
fini de Lagrange ainsi défini vérifie les critères de cohérence 4.6 page 141 et (4.7) page 141. Pour compléter la
définition de l’espace d’approximation HN , il ne reste qu’à déterminer les “noeuds liés”, de la façon dont on a
traité le cas de l’espace H01 (Ω).
Il faut également insister sur le fait que cet élément est très souvent utilisé, en raison de sa facilité d’implantation
et de la structure creuse des systèmes linéaires qu’il génère. Il est particulièrement bien adapté lorsqu’on cherche
des solutions dans l’espace H 1 (Ω). Il se généralise facilement en trois dimensions d’espace, où on utilise alors des
tétraèdres, avec toujours comme espace de polynôme l’espace des fonctions affines.

Analyse numérique des EDP, M1 143 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.2. EXEMPLES CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

4.2.2 Elément fini triangulaire P 2

Comme le titre du paragraphe l’indique, on considère un maillage triangulaire, et un espace de polynômes de degré
2 pour construire l’espace d’approximation.
Elément fini de référence On choisit comme élément fini de référence le triangle de sommets (0,0), (1,0) et (0,1),
voir Figure 4.6 et on prend pour Σ :
1 1 1 1
Σ̄ = {(0, 0), (1, 0), (0, 1), ( , ), (0, ), ( , 0)}
2 2 2 2

3
(0,1)
(0,0,1)

(1/2,1/2)
5 4 (0,1/2,1/2)
(0,1/2)
(1/2,0,1/2)

1 6 2
(0,0) (1/2,0) (1,0)
(0,1,0)
 (0,0) (1/2,0) (1,0)
(1,0,0) (1/2,1/2,0) (0,1,0)

F IGURE 4.6 – Elément de référence pour les éléments finis P2, avec coordonnées cartésiennes et barycentriques
des noeuds

Fonctions de base locales Les fonctions de base locales sont définies à partir des coordonnées barycentriques. On
rappelle que les coordonnées barycentriques d’un point x du triangle K de sommets a1 , a2 et a3 sont les réels
λ1 , λ2 , λ3 tels que :
x = λ1 a1 + λ2 a2 + λ3 a3 .
Dans le cas du triangle de référence K̄ de sommets (0,0), (1,0) et (0,1), les coordonnées barycentriques d’un point
{x} de coordonnées cartésiennes x et y sont donc : λ1 = 1 − x − y, λ2 = x, λ3 = y. Par définition, on a
3
i=1 λi = 1 et λi ≥ 0 (car le triangle K est l’enveloppe convexe de l’ensemble de ses sommets). On peut alors
déterminer les fonctions de base en fonction des coordonnées barycentriques des six noeuds de K̄ exprimés par
leurs coordonnées barycentriques : a1 = (1, 0, 0), a2 = (0, 1, 0), a3 = (0, 0, 1), a4 = (0, 12 , 12 ), a5 = ( 12 , 0, 12 ),
a6 = ( 12 , 12 , 0). Les fonctions de base sont telles que φi ∈ P2 , et φi (aj ) = δij , ∀i = 1, . . . , 6, f orallj = 1, . . . , 6.
Commençons par φ1 ; on veut φ1 (a1 ) = 1, et φi (ai ) = 0, ∀i = 2, . . . , 6. La fonction φ1 définie par
1
φ1 (x, y) = 2λ1 (λ1 − )
2
convient, et comme le couple (Σ̄, P 2) est unisolvant, c’est la seule fonction qui convient. Par symétrie, on définit
1
φ2 (x, y) = 2λ2 (λ2 − ),
2
et
1
φ3 (x, y) = 2λ3 (λ3 − ).
2
Les fonctions de base associées aux noeuds a4 , a5 , a6 sont alors

φ4 (x, y) = 4λ2 λ3 ,

Analyse numérique des EDP, M1 144 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.2. EXEMPLES CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

φ5 (x, y) = 4λ1 λ3 ,
et φ6 (x, y) = 4λ1 λ2 .
Il est facile de voir que ces fonctions forment une famille libre d’éléments de P 2 et comme card Σ̄ = card P 2, le
couple (Σ̄, P 2) est bien unisolvant.
Transformation Fℓ La bijection Fℓ qui permet de passer de l’élément fini de référence K̄ à l’élément Kℓ a déjà été
vue dans le cas de l’élément fini P1 c’est la fonction affine définie par :
- .
x1 + (x2 − x1 )x + (x3 − x1 )y
Fℓ (x, y) =
y1 + (y2 − y1 )x + (y3 − y1 )y
où (xi , yi ), i = 1, 2, 3 sont les coordonnées respectives des trois sommets du triangle Kℓ . Comme cette transfor-
mation est affine, les coordonnées barycentriques restent inchangées par cette transformation.
On peut montrer (ce n’est pas facile) que l’erreur d’interpolation "u − uN "H 1 est contrôlée, en éléments finis P1
et P2 par les inégalités suivantes :
P 1 : si u ∈ H 2 (Ω), on a "u − uN "H 1 (Ω) ≤ Ch"u"H 2 (Ω)
P 2 : si u ∈ H 3 (Ω), on a "u − uN "H 1 (Ω) ≤ Ch2 "u"H 3 (Ω) .
On peut généraliser les éléments finis P1 et P2 aux éléments finis Pk sur triangles, pour k ≥ 1. On prend toujours
le même élément de référence, dont on divise chaque côté en k intervalles. Les extrémités de ces intervalles sont
les noeuds du mailage. On a donc 3k noeuds, qu’on peut repérer par leurs coordonnés barycentriques, qui prennent
les valeurs 0, k1 , k2 , . . . , 1. On peut montrer que si u ∈ H k+1 , alors
"uN − u"H 1 (Ω) ≤ Chk "u"H k+1 (Ω)

4.2.3 Eléments finis sur quadrangles

4.2.3 Eléments finis sur quadrangles
Le cas rectangulaire
On prend comme élément fini de référence le carré K̄ = [−1, 1] × [−1, 1], et comme noeuds les coins de ce carré :
a1 = (1, −1), a2 = (1, 1), a3 = (−1, 1), et a4 = (−1, −1).
On prend comme espace de polynômes
P = {f : K̄ → IR; f ∈ Q1 }
où Q1 = {f : IR 2 → IR; f (x, y) = a + bx + cy + dxy, (a, b, c, d) ∈ IR 4 } Le couple (Σ, P ) est unisolvant. Les
fonctions de base locales sont les fonctions :
1
φ1 (x, y) = − (x + 1)(y − 1)
4
1
φ2 (x, y) = (x + 1)(y + 1)
4
1
φ3 (x, y) = − (x − 1)(y + 1)
4
1
φ4 (x, y) = (x − 1)(y − 1).
4
La transformation Fℓ permet de passer de l’élément de référence carré K̄ à un rectangle quelconque du maillage
Kℓ . Si on considère un rectangle Kℓ parallèle aux axes, dont les noeuds sont notés (x1 , y1 ), (x2 , y1 ), (x2 , y2 ),
(x1 , y2 ), les noeuds du rectangle Kℓ , la bijection Fℓ s’écrit :
- .
1 (x2 − x1 )x + x2 + x1
Fℓ (x, y) = .
2 (y2 − y1 )y + y2 + y1
Considérons maintenant le cas d’un maillage quadrangulaire quelconque. Dans ce cas, on choisit toujours comme
élément de référence le carré unité. La transformation Fℓ qui transforme l’élément de référence en un quadrangle
Kℓ est toujours affine, mais par contre, les composantes de Fℓ ((x, y)) dépendent maintenant de x et de y voir
exercice 54 page 163. En conséquence, le fait que f ∈ Q1 n’entraı̂ne plus que f ◦ Fℓ ∈ Q1 . Les fonctions de base
seront donc des polynômes Q1 sur l’élément de référence K̄, mais pas sur les éléments “courants” Kℓ .

Analyse numérique des EDP, M1 145 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.3. CONSTRUCTION DU SYSTÈME LINÉAIRE CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

Eléments finis d’ordre supérieur Comme dans le cas d’un maillage triangulaire, on peut choisir un espace de
polynômes d’ordre supérieur, Qk , pour les fonctions de base de l’élément de référence K̄ = [−1, 1] × [−1, 1]. On
choisit alors comme ensemble de noeuds : Σ̄ = Σ̄k = {(x, y) ∈ K̄, (x, y) ∈ {−1, −1 + k1 , −1 + k1 , . . . , 1}2 . On
peut montrer facilement que (Σ̄k Qk ) est unisolvant. Là encore, si la solution exacte de problème continu est
suffisamment régulière, on peut démontrer l’estimation d’erreur suivante (voir [3]) :

"u − uN "H ′ (Ω) ≥ C"u"H k+1 (Ω) hk .

Exprimons par exemple l’espace des polynômes Q2 . On a :

Q2 = {f : IR → IR; f (x) = a1 +a2 x+a3 y+a4xy+a5 x2 +a6 y 2 +a7 xy 2 +a8 x2 y+a9x2 y 2 , ai ∈ IR, i = 1, . . . , 9}

L’espace Q2 comporte donc neuf degrés de liberté. On a donc besoin de neuf noeuds dans Σ̄ pour que le couple
(Σ̄, Q2 ) soit unisolvant (voir exercice 59 page 165). On peut alors utiliser comme noeuds sur le carré de référence
[−1, 1] × [−1, 1] :

Σ̄ = {(−1, −1), (−1, 0), (−1, −1), (0, −1), (0, 0), (0, 1), (1, −1), (1, 0), (1, 1)}

En général, on préfère pourtant supprimer le noeud central (0,0)et choisir :

Σ∗ = Σ̄ \ {(0, 0)}.

Il faut donc un degré de liberté en moins pour l’espace des polynômes. On définit alors :

Q∗2 = {f : IR → IR; f (x) = a1 + a2 x + a3 y + a4 xy + a5 x2 + a6 y 2 + a7 xy 2 + a8 x2 y, ai ∈ IR, i = 1, . . . , 8}

Le couple (Σ∗ , Q∗2 ) est unisolvant (voir exercice 60 page 165), et on peut montrer que l’élément Q∗2 est aussi précis
(et plus facile à mettre en oeuvre que l’élément Q ).
Le couple , Q est unisolvant (voir exercice 60 page 165), et on peut montrer que l’élément est aussi précis
(et plus facile à mettre en oeuvre que l’élément Q2 ).

4.3 Construction du système linéaire

On construit ici le système linéaire pour un problème à conditions aux limites mixtes de manière ‘a envisager
plusieurs types de conditions aux limites. Soit Ω un ouvert polygonal 1 , on suppose que ∂Ω : Γ0 ∪ Γ1 avec
mes(Γ0 ) = 0. On va imposer des conditions de Dirichlet sur Γ0 et des conditions de Fourier sur Γ1 ; c’est ce qu’on
appelle des conditions “mixtes”. On se donne donc des fonctions p : Ω → IR, g0 : Γ0 → IR et g1 : Γ1 → IR, et on
cherche à approcher u solution de :


 −div(p(x)∇u(x)) + q(x)u(x) = f (x), x ∈ Ω,

u = g0 sur Γ0 , (4.11)



p(x)∇u(x).n(x) + σu(x) = g1 (x), x ∈ Γ1 ,

où n désigne le vecteur unitaire normal à ∂Ω extérieure à Ω. Pour assurer l’existence et unicité du problème (4.11),
(voir exercice 42), on se place sous les hypothèses suivantes :


 p(x) ≥ α > 0, p.p. x ∈ Ω



 q≥0
(4.12)

 σ≥0



 mes(Γ ) > 0.
0

Pour obtenir une formulation variationnelle, on introduit l’espace

HΓ10 ,g0 = {u ∈ H 1 (Ω); u = g0 sur Γ0 }

et l’espace vectoriel associé :

H = HΓ10 ,0 = {u ∈ H 1 (Ω); u = 0 sur Γ0 }
1. Dans le cas où la frontière ∂Ω de Ω n’est pas polygônale mais courbe, il faut considérer des éléments finis dits “isoparamétriques” que
nous verrons plus loin

Analyse numérique des EDP, M1 146 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.3. CONSTRUCTION DU SYSTÈME CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

Notons que H est un espace de Hilbert. Par contre, attention, l’espace HΓ10 ,g0 n’est pas un espace vectoriel. On va
chercher u solution de (4.11) sous la forme u = u + u0 , avec u0 ∈ HΓ10 ,g0 et u
 ∈ HΓ10 ,0 . Soit v ∈ H, on multiplie
(4.11) par v et on intègre sur Ω. On obtient :
( ( (
−div(p(x)∇u(x))v(x)dx + q(x)u(x)v(x)dx = f (x)v(x)dx, ∀v ∈ H.
Ω Ω Ω

En appliquant la formule de Green, il vient alors :

( ( ( (
p(x)∇u(x)∇v(x)dx − p(x)∇u(x)nv(x)dγ(x) + qu(x)v(x)dx = f (x)v(x)dx, ∀v ∈ H.
Ω ∂Ω Ω Ω

Comme v = 0 sur Γ0 on a :
( (
p(x)∇u(x)nv(x)dγ(x) = p∇u(x)nv(x)dγ(x).
∂Ω Γ1

Mais sur Γ1 , la condition de Fourier s’écrit : ∇u.n = −σu + g1 , et on a donc

( ( ( ( (
p(x)∇u(x)∇v(x)dx+ p(x)σu(x)v(x)dγ(x)+ qu(x)v(x)dx = f (x)v(x)dx+ g1 (x)v(x)dγ(x).
Ω Γ1 Ω Ω Γ1

On peut écrire cette égalité sous la forme : a(u, v) = T̃ (v), avec a(u, v) = aΩ (u, v) + aΓ1 (u, v), où :
 ( (

 a (u, v) = p(x)∇u(x)∇v(x)dx + q(x)u(x)v(x)dx,
 Ω
Ω Ω
(


 aΓ (u, v) = p(x)σ(x)u(x)v(x)dγ(x),
 

 aΓ (u, v) = p(x)σ(x)u(x)v(x)dγ(x),
Γ1

et T̃ (v) = TΩ (v) + TΓ1 (v), avec

( (
TΩ (v) = f (x)v(x)dx et TΓ1 = g1 (x)v(x)dγ(x).
Ω Γ1

On en déduit une formulation faible associée à (4.11) :

1
chercher u ∈H
(4.13)
a(u0 + u
, v) = T̃ (v), ∀v ∈ H,

où u0 ∈ H 1 (Ω) est un révèlement de g0 , c’est à dire une fonction de H 1 (Ω) telle que u0 = g0 sur Γ. Le problème
(4.13) peut aussi s’écrire sous la forme :
1
u
∈H
(4.14)
a( u, v) = T (v), ∀v ∈ H.

où T (v) = T̃ (v) − a(u0 , v). Sous les hypothèses (4.12), on peut alors appliquer le théorème de Lax Milgram (voir
théorème 3.6 page 100) au problème (4.14) pour déduire l’existence et l’unicité de la solution de (4.13) ; notons
que, comme la forme bilinéaire a est symétrique, ce problème admet aussi une formulation variationnelle :

 J(u) = min J(v),
 v∈HΓ1
0 ,g0
(4.15)

 J(v) = 1 a(v, v) + T (v), ∀v ∈ H 1 .
Γ0 ,g0
2
Dans ce cas, les méthodes de Ritz et Galerkin sont équivalentes. Remarquons que l’on peut choisir u0 de manière
abstraite, tant que u0 vérifie u0 = g0 sur Γ0 et u0 ∈ H 1 . Intéressons nous maintenant à la méthode d’approximation
variationnelle. On approche l’espace H par HN = V ect{φ1 , . . . , φN } et on remplace (4.14) par :
1
N ∈ HN
u
(4.16)
a(
uN , φi ) = T (φi ) − a(u0 , φi ), ∀i = 1, . . . , N.

Analyse numérique des EDP, M1 147 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.3. CONSTRUCTION DU SYSTÈME CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

N
%
On pose maintenant u
N = u
j φj . Le problème (4.16) est alors équivalent au système linéaire :
j=1

 = G,
KU
avec 

 Kij = a(φj , φi ), i, j = 1, . . . , N,

 = (
U u1 , . . . , u
N ),



Gi = T (φi ) − a(u0 , φi ), i = 1, . . . , N.
L’implantation numérique de la méthode d’approximation nécessite donc de :
1. construire K et G
 = G.
2. résoudre KU
Commençons par la construction de l’espace HN et des fonctions de base pour une discrétisation par éléments
finis de Lagrange du problème (4.14).

4.3.1 Construction de HN et Φi
On considère une discrétisation à l’aide déléments finis de Lagrange, qu’on note : (Kℓ , Σℓ , Pℓ ) ℓ = 1, . . . , L, où
L est le nombre d’éléments. On note Si , i = 1, . . . , M, les noeuds du maillage, et φ1 , . . . , φN , les fonctions de
base, avec N ≤ M . On peut avoir deux types de noeuds :
– les noeuds libres : Si ∈ Γ0 . On a N noeuds libres
– les noeuds liés : Si ∈ Γ0 . On a M − N noeuds liés.
Notons qu’on a intérêt à mettre des noeuds à l’intersection de Γ0 et Γ1 (ce seront des noeuds liés). Grâce à ceci, et
à la cohérence globale et locale des éléments finis de Lagrange, on a HN ⊂ H. On a donc bien des éléments finis
conformes. Récapitulons alors les notations :
à la cohérence globale et locale des éléments finis de Lagrange, on a HN ⊂ H. On a donc bien des éléments finis
conformes. Récapitulons alors les notations :
• M : nombre de noeuds total
• N : nombre de noeuds libres
• M0 = M − N : nombre de noeuds liés
• J0 = { indices des noeuds liés} ⊂ {1, . . . , M }. On a cardJ0 = M0
• J = { indices des noeuds libres } = {1 . . . M } J0 . On a cardJ = N.
Pour la programmation des éléments finis, on a besoin de connaı̂tre, pour chaque noeud (local) de chaque élément,
son numéro dans la numérotation globale. Pour cela on introduit un tableau ng(L, Nℓ ), où L est le nombre
d’éléments et Nℓ est le nombre de noeuds par élément. (on le suppose constant par souci de simplicité, Nℓ peut en
fait dépendre de L. Exemple : triangle - quadrangle). Pour tout ℓ ∈ {1, . . . , L} et tout r ∈ {1, . . . , Nℓ }, ng(ℓ, r) est
alors le numéro global du r-ième noeud du ℓ-ième élément. On a également besoin de connaı̂tre les coordonnées
de chaque noeud. On a donc deux tableaux x et y de dimension M , où x(i), y(i) représentent les coordonnées
du i-ème noeud. Notons que les tableaux ng, x et y sont des données du mailleur (qui est un module externe par
rapport au calcul éléments finis proprement dit). Pour les conditions aux limites, on se donne deux tableaux :
• CF : conditions de Fourier
• CD : conditions de Dirichlet
(on verra plus tard le format de ces deux tableaux)

4.3.2 Construction de K et G
On cherche à construire la matrice K d’ordre (N × N ), définie par :
Kij = a(φj , φi ) i, j ∈ J
Ainsi que le vecteur G, défini par :
Gi = T (φi ) − a(u0 , φi ) i∈J cardJ = N
La première question à résoudre est le choix de u0 . En effet, contrairement au cas unidimensionnel (voir exercice
38 page 119), il n’est pas toujours évident de trouver u0 ∈ HΓ10 ,g0 . Pour se faciliter la tâche, on commet un “crime
variationnel”, en remplaçant u0 par
%N
u0,N = u0 (Sj )φj .
j∈J0

Analyse numérique des EDP, M1 148 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.3. CONSTRUCTION DU SYSTÈME CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

Notons qu’on a pas forcément : u0,N ∈ HΓ10 ,g0 ; c’est en ce sens que l’on commet un “crime”. Mais par contre,
on a bien u0,N (Sj ) = u0 (Sj ) pour tout j ∈ J. On peut voir la fonction u0,N comme une approximation non
conforme de u0 ∈ HΓ10 ,g0 . On remplace donc Gi par :
%
Gi = T (φi ) − g0 (Sj )a(φj , φi ).
j∈J0

Calculons maintenant a(φj , φi ) pour j = 1, . . . , M, et i = 1, . . . , M. On se sert donc pour l’implatation pratique

de la méthode, des fonctions de forme associées aux noeuds “liés”, même si dans l’écriture du problème discret
théorique, on n’en avait pas besoin.

Calcul de K et G
1. Calcul des contributions intérieures : on initialise les coefficients de la matrice K et les composantes par les
contributions provenant de aΩ et TΩ .
2
Kij = aΩ (φj , φi ) i = 1, . . . , N,
Gi = TΩ (φi ) j = 1, . . . , N.

2. Calcul des termes de bord de Fourier. On ajoute maintenant à la matrice K les contributions de bord :
Kij ←− Kij + aΓi (φj , φi ), i = 1, . . . , N, j = 1, . . . , N.
Gi ←− Gi + TΓi (φi ) i = 1 . . . M.

3. Calcul des termes de bord de Dirichlet. On doit tenir compte ici du relèvement de la condition de bord :
%
Gi ←− Gi − g0 (Ni )Kij ∀i ∈ J
j∈J0
 ←− ij
j∈J0

Après cette affectation, les égalités suivantes sont vérifiées :

Kij = a(φj , φi ) i, j ∈ J(∪J0 )
Gi = T (φi ) − a(u0,N , φi ).
Il ne reste plus qu’à résoudre le système linéaire
%
Kij αj = Gi , ∀i ∈ J. (4.17)
j∈J

4. Prise en compte des noeuds liés. Pour des questions de structure de données, on inclut en général les noeuds
liés dans la résolution du système, et on résout donc le système linéaire d’ordre M ≥ N suivant :
%
K̃ij αj = Gi . ∀i = 1, . . . , N. (4.18)
j=1,...,N

avec K̃ij = Kij pour i, j ∈ J, K̃ij = 0 si (i, j) ∈ J 2 , et i = j, et K̃ii = 1 si i ∈ J. Ces deux systèmes sont
équivalents, puisque les valeurs aux noeuds liées sont fixées.
Si par chance on a numéroté les noeuds de manière à ce que tous les noeuds liés soient en fin de numérotation,
c.à.d. si J = {1, . . . , N } et J0 = {N + 1, . . . , M }, le système (4.18) est de la forme :
   
  α1 G1
|  ..   .. 
 K | 0   .   . 
     
 |     
   αuN   GN 
 − − − | − − −  , U =  −−  , et G =  −− 
     
 |   αN +1   GN +1 
     
 0 | IdM   ..   .. 
 .   . 
|
αM GM

Dans le cas où la numérotation est quelconque, les noeuds liés ne sont pas forcément à la fin, et pour obtenir
le système linéaire d’ordre M (4.18) (donc incluant les inconnues αi , i ∈ J0 , qui n’en sont pas vraiment) on
peut adopter deux méthodes :

Analyse numérique des EDP, M1 149 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.3. CONSTRUCTION DU SYSTÈME CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

(a) Première méthode : on force les valeurs aux noeuds liés de la manière suivante :
Kii ←− 1 pour tout i ∈ J0
Kij ←− 0 pour tout i ∈ J0 j ∈ {1 . . . M } i = j
Gi ←− g0 (Si ) pour tout i ∈ J0
(b) Deuxième méthode : on force les valeurs aux noeuds liés de la manière suivante :
Kii ←− 1020 ∀i ∈ J0
Gi ←− 1020 g0 (Si ) ∀i ∈ J0
La deuxième méthode permet d’éviter l’affectation à 0 de coefficients extra-diagonaux de la matrice. Elle
est donc un peu moins chère en temps de calcul.

Conclusion Après les calculs 1, 2, 3, 4, on a obtenu une matrice K d’ordre M × M et le vecteur G de IR M .

Soit α ∈ IR M la solution du système Kα = G. Rappelons qu’on a alors :
N
%
uN = αi φi ,
i=1
% %
= αi φi + αi φi
i∈J i∈J0

uN =u
N + u0

Remarque 4.15 (Numérotation des noeuds) Si on utilise une méthode itérative sans préconditionnement, la numérotation
des noeuds n’est pas cruciale. Elle l’est par contre dans le cas d’une méthode directe et si on utilise une méthode
itérative avec préconditionnement. Le choix de la numérotation s’effectue pour essayer de minimiser la largeur
des noeuds n’est pas cruciale. Elle l’est par contre dans le cas d’une methode directe et si on utilise une methode
itérative avec préconditionnement. Le choix de la numérotation s’effectue pour essayer de minimiser la largeur
de bande. On pourra à ce sujet étudier l’influence de la numérotation sur deux cas simples sur la structure de la
matrice.

4.3.3 Calcul de aΩ et TΩ , matrices élémentaires.

Détaillons maintenant le calcul des contributions intérieures, c’est à dire aΩ (φi , φj ) i = 1, . . . M, j = 1, . . . , M
et TΩ (φi ) i = 1, . . . , M. Par définition,
( (
aΩ (φi , φj ) = p(x)∇φi (x)∇φj (x)dx + q(x)φi (x)φj (x)dx.
Ω Ω

Décomposons Ω à l’aide du maillage éléments finis.

L

Ω= Kℓ .
ℓ=1

En notant θ(φi , φj )(x) = p(x)∇φi (x)∇φj (x) + q(x)φi (x)φj (x),

On a donc :
L (
%
aΩ (φi , φj ) = θ(φi , φj )dx.
ℓ=1 Kℓ

Pour r et s numéros locaux de l’élément Kℓ , on pose :

(
ℓ
kr,s = θ(φs , φr )dx.
ℓ

ℓ
On va calculer kr,s puis on calcule aΩ (φi , φj ), en effectuant un parcours sur les éléments, ce qui s’exprime par
l’algorithme suivant :
Initialisation : Kij ←− 0, i = 1, . . . M, j ≤ i.
Boucle sur les éléments
Pour ℓ = 1 à L faire
Pour r = 1 à Nℓ faire

Analyse numérique des EDP, M1 150 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.3. CONSTRUCTION DU SYSTÈME CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

i = ng(ℓ, r) numéro global du noeud r de l’élément ℓ

Pour s = 1 à r faire
ℓ
calcul de kr,s
j = ng(ℓ, s)
si i ≥ j
ℓ
Kij ←− Kij + kr,s
sinon
ℓ
Kji ←− Kji + krs
Fin pour
Fin pour
On a ainsi construit complètement la matrice de rigidité K. Il reste à savoir comment calculer
(
ℓ
kr,s = θ(φs , φr )(x)dx.
ℓe
ℓ
Ce calcul s’effectue sur l’élément de référence, et non sur les éléments Kℓ . On calcule ensuite la valeur de kr,s par
des changements de variable à l’aide de la transformation Fℓ (voir Figure 4.4 page 138). Notons :
Fℓ (x̄, ȳ) = (x, y) = (aℓ0 + aℓ1 x̄ + aℓ2 ȳ, bℓ0 + bℓ1 x̄ + bℓ2 y) (4.19)
Notons que les coefficients aℓi et bℓi sont déterminés à partir des la connaissances des coordonnées (x(i), y(i)) où
i = ng(ℓ, r). En effet, on peut déduire les coordonnées locales x(r), y(r), r = 1, Nℓ , des noeuds de l’élément ℓ, à
partir des coordonnées globales des noeuds (x(i), y(i)), et du tableau ng(ℓ, r) = i. Sur l’élément courant Kℓ , le
ℓ
terme élémentaire kr,s s’écrit donc
(
ℓ
kr,s = θ(φs (x, y), φr (x, y))dxdy.
ℓ

Or, x, y) = (x̄, ; donc par changement de variables , on a :

Or, (x, y) = Fℓ (x̄, ȳ) ; donc par changement de variables , on a :
(
ℓ
kr,s = θ(φs ◦ Fℓ (x̄, ȳ), φr ◦ Fℓ (x̄, ȳ)) Jacx̄,ȳ (Fℓ ) dx̄dȳ
ē

ou Jacx̄,ȳ (Fℓ ) désigne le Jacobien de Fℓ en (x̄, ȳ). Or, φs ◦ Fℓ = φ̄s , et, puisque Fℓ est définie par (4.19), on a :
 ℓ 
 a1 bℓ1   ℓ ℓ 
Jac(Fℓ ) = Det(DFℓ ) =  ℓ   = a1 b2 − aℓ2 bℓ1 
a2 bℓ2 
ℓ
donc kr,s = Jac(Fℓ )k̄r,s , où (
k̄r,s = θ(φ̄s (x̄, ȳ), φ̄r (x̄, ȳ))dx̄dȳ
ē
Etudions maintenant ce qu’on obtient pour k̄r,s dans le cas du problème modèle (4.11), on a :
(
 
k̄r,s = p(x̄, ȳ)∇φ̄s (x̄, ȳ)∇φ̄r (x̄, ȳ) + q(x̄, ȳ)φ̄s (x̄, ȳ)φ̄r (x̄, ȳ) dx̄dȳ.
ℓ̄

Les fonctions de base φ̄s et φ̄r sont connues ; on peut donc calculer k̄r,s explicitement si p et q sont faciles à
intégrer. Si les fonctions p et q ou les fonctions de base φ̄, sont plus compliquées, on calcule k̄r,s en effectuant une
intégration numérique. Rappelons que le principe d’une intégration numérique est d’approcher l’intégrale d ’une
fonction continue donnée ψ,
( N
% PI
I = ψ(x̄, ȳ)dx̄dȳ, par I = ωi (Pi )ψ(Pi ),
ℓ̄ i=1

où N PI est le nombre de points d’intégration, notés Pi , qu’on appelle souvent points d’intégration de Gauss, et les
coefficients ωi sont les poids associés. Notons que les points Pi et les poids ωi sont indépendants de ψ. Prenons
par exemple, dans le cas unidimensionnel, K̄ = [0, 1], p1 = 0, p2 = 1, et ω1 = ω2 = 12 . On approche alors
( 1
1
I= ψ(x)dx par I = (ψ(0) + ψ(1)).
0 2
C’est la formule (bien connue) des trapèzes. Notons que dans le cadre d’une méthode, il est nécessaire de s’assurer
que la méthode d’intégration numérique choisie soit suffisamment précise pour que :

Analyse numérique des EDP, M1 151 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.3. CONSTRUCTION DU SYSTÈME CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

1. le système Kα = G(N × N ) reste inversible,

2. l’ordre de convergence de la méthode reste le même.
Examinons maintenant des éléments en deux dimensions d’espace.
1 1

1. Elément fini P1 sur triangle Prenons N PI = 1 (on a donc un seul point de Gauss), choisissons p1 = 3, 3 ,
le centre de gravité du triangle K̄, et ω1 = 1. On approche alors
(
I = ψ(x̄)dx̄ par ψ(p1 ).
ℓ̄

On vérifiera que cette intégration numérique est exacte pour les polynômes d’ordre 1 (exercice 58 page 165).
2. P2 sur triangles. On prend maintenant N PI = 3, et on choisit comme points de Gauss :
$ ' $ ' $ '
1 1 1 1
p1 = , 0 , p2 = , , p3 = 0,
2 2 2 2

et les poids d’intégration ω1 = ω2 = ω3 = 16 . On peut montrer que cette intégration numérique est exacte
pour les polynômes d’ordre 2 (voir exercice 58 page 165).
Remarquons que, lors de l’intégration numérique du terme élémentaire
(
ℓ
 
kr,s = p(x̄, ȳ)(Fℓ (x̄, ȳ))∇φ̄r (x̄, ȳ) · ∇φ̄s (x̄, ȳ) + q(x̄, ȳ)(Fℓ (x̄, ȳ))φ̄r (x̄, ȳ)φ̄s dx̄dȳ,
ℓ̄

ℓ
on approche kr,s par

N
% PI
 
k̄r,s ≃ ωi p(Fℓ (Pi ))∇φ̄r (Pi ) · ∇φ̄s (Pi ) + q(Fℓ (Pi ))φ̄r (Pi )φ̄s (Pi ) .
 N
% P
 
k̄r,s ≃ ωi p(Fℓ (Pi ))∇φ̄r (Pi ) · ∇φ̄s (Pi ) + q(Fℓ (Pi ))φ̄r (Pi )φ̄s (Pi ) .
i=1

Les valeurs ∇φ̄r (Pi ), ∇φ̄s (Pi ), φ̄r (Pi ) et φ̄s (Pi ) sont calculées une fois pour toutes, et dans la boucle sur ℓ, il ne
reste donc plus qu’à évaluer les fonctions p et q aux points Fℓ (Pi ). Donnons maintenant un résumé de la mise en
oeuvre de la procédure d’intégration numérique (indépendante de ℓ). Les données de la procédure sont :
– les coefficients ωi , i = 1, . . . , N PI ,
– les coordonnées (xpg(i), ypg(i)), i = 1, . . . , N PI des points de Gauss,
– les valeurs de φr , ∂φ ∂φr
∂x et ∂y aux points de Gauss, notées φ(r, i), φx (r, i) et φy (r, i), r = 1 . . . Nℓ , i = 1, . . . , N PI .
Pour ℓ donné, on cherche à calculer :
( (
∂φr ∂φs
I= p(Fℓ (x̄, ȳ)) (x̄, ȳ) (x̄, y)dx̄dȳ + q(Fℓ (x̄, ȳ))φr (x̄, y)dx̄dȳφs (x̄, ȳ).
K̄ ∂ x̄ ∂ ȳ e

On propose l’algorithme suivant :

Initialisation : I ←− 0
Pour i = 1 à N PI , faire :
pi = p(Fe (Pi ))
qi = q(Fe (Pi ))
I ←− I + ωi (pi φx (r, i)φy (s, i) + qi φ(r, i)φ(s, i))
Fin pour
On procède de même pour le calcul du second membre
( L
% (
TΩ (φi ) = f (x, y)φ′i x, y)dxdy = gℓ , où gℓ = f (x, y)φi (x, y)dxdy.
Ω ℓ=1 ℓe

L’algorithme sécrit :
Initialisation de G à 0 : Gi ←− 0 i = 1 à M
Pour ℓ = 1 à L
Pour r = 1 à Nℓ "
Calcul de gℓr = ℓe f (x, y)φr (x, y)dxdy
i = ng(ℓ, r)
Gi ←− Gi + gℓr

Analyse numérique des EDP, M1 152 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.3. CONSTRUCTION DU SYSTÈME CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

Fin pour
Fin pour
Il reste le calcul de gℓr qui se ramène au calcul de l’élément de référence par changement de variable.
On a : ( (
gℓr = f (x, y)φr (x, y)dxdy = f ◦ Fℓ (x̄, ȳ)φ̄r (x̄, ȳ)Jacx̄,ȳ (Fℓ )dx̄dȳ.
Kℓ K̄

L’intégration numérique est identique à celle effectuée pour k̄r,s .

4.3.4 Calcul de aΓ1 et TΓ1 (contributions des arêtes de bord “Fourier”.

Détaillons maintenant le calcul des contributions des bords où s’applique la condition de Fourier, c’est à dire
aΓ1 (φi , φj ) i = 1, . . . M, j = 1, . . . , M et TΓ1 (φi ) i = 1, . . . , M. Par définition,
( (
aΓ1 (φi , φj ) = p(x)∇φi (x) · ∇φj (x)dx + q(x)φi (x)φj (x)dx.
Γ1 Γ1

Notons que aΓ1 (φi , φj ) = 0 si φi et φj sont associées à des noeuds Si , Sj de d’un élément sans arête commune avec
les arêtes de la frontière. Soit L1 le nombre d’arêtes ǫk , k = 1, . . . , L1 du maillage incluses dans Γ1 . Rappelons
que les noeuds soumis aux conditions de Fourier sont repertoriés dans un tableau CF , de dimensions (L1, 2), qui
donne les informations suivantes
1. CF (k, 1) contient le numéro ℓ de l’élément Kℓ auquel appartient l’arête ǫk .
2. CF (k, 2) contient le premier numéro des noeuds de l’arête ǫk dans l’élément Kℓ . On suppose que la
numérotation des noeuds locaux a été effectuée de manière “adroite”, par exemple dans le sens trigo-
nométrique. Dans ce cas, CF (k, 2) détermine tous les noeuds de l’arête ǫk dans l’ordre, puisqu’on connait
le nombre de noeuds par arête et le sens de numérotation des noeuds. Donnons des exemples pour trois cas
 nométrique. Dans ce cas, CF k, 2) détermine tous les noeuds de l’arête dans l’ordre, puisqu’on connait
le nombre de noeuds par arête et le sens de numérotation des noeuds. Donnons des exemples pour trois cas
différents, représentés sur la figure 4.7.
(a) Dans le premier cas (à droite sur la figure), qui représente un élément fini P 1, on a CF (k, 2) = 3 et le
noeud suivant sur l’arête est 1.
(b) Dans le second cas (au centre sur la figure), qui représente un élément fini P 2, on a CF (k, 2) = 3 et
les noeuds suivants sur l’arête sont 4 et 5.
(c) Enfin dans l’élément P 1 “de coin” représenté à gauche sur la figure, on a CF (k, 1) = ℓ, CF (k ′ , 1) =
ℓ, CF (k, 2) = 1, CF (k ′ , 2) = 2.

1 5 ǫk′
2
3 Kℓ
ǫk 6 4 ǫk
Kℓ Kℓ ǫk

2 3 1 2 3 1

P1 P1 P 1 en coin

F IGURE 4.7 – Exemples de numérotation d’arête du bord

Pour k = 1, . . . , L1, on note Ŝk l’ensemble des noeuds locaux de ǫk , donnés par CF (k, 2) en appliquant la règle
ad hoc (par exemple le sens trigonomt́rique). On peut alors définir :

Sk = {(r, s) ∈ (Ŝk )2 /r < s}

L’algorithme de prise en compte des conditions de Fourier s’écrit alors :

Pour k = 1 . . . L1

Analyse numérique des EDP, M1 153 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.3. CONSTRUCTION DU SYSTÈME CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

ℓ = CF (k, 1).
Pour chaque (r, (
s) ∈ Sk faire
ℓ
calcul de Irs = p(x)σ(x)φℓr (x)φℓs (x)dx (éventuellement avec intégration numérique)
Ck
i = ng(ℓ, r)
j = ng(ℓ, s)
si j ≤ i
ℓ
Kij ←− Kij + Irs
sinon
ℓ
Kij ←− Kji + Irs
Fin si
Fin pour
Fin pour
ℓ
Le calcul de Irs s’effectue sur l’élément de référence
( (avec éventuellement intégration numérique). De même, on
a une procédure similaire pour le calcul de TΓ1 = p(x)g1 (x)v(x)dγ(x).
Γ1
(
Gi ←− Gi + p(x)g1 (x)φi (x)dγ(x)
Γ2

4.3.5 Prise en compte des noeuds liés dans le second membre

Aprés les calculs précédents, on a maintenant dans Gi :
( (
Gi = f (x)φi (x)dx + p(x)g1 (x)φi (x)dγ(x)
Ω Γ1
 G dx
Ω Γ1

Il faut maintenant retirer du second membre, les combinaisons venant des noeuds liés :
%
Gi ←− Gi − g0 (Sj )a(φj , φi )
j∈J0

où J0 est l’ensemble des indices des noeuds liés. On utilise pour cela le tableau CD qui donne les conditions, de
Dirichlet, de dimension M0 où M0 = cardJ0 . Pour i0 = 1, . . . , M0 , CD(i0 ) = j0 ∈ J0 est le numéro du noeud
lié dans la numérotation globale. La procédure est donc la suivante.
Pour i0 = 1, . . . , M0 , faire
j = CD(i0 )
a = g0 (Sj )
si (i ≤ j) Gi ←− Gi − aKij sinon Gi ←−
Gi − aKji sinon Fin si Fin pour

4.3.6 Stockage de la matrice K

Remarquons que la matrice K est creuse (et même très creuse), en effet a(φj , φi ) = 0 dès que

supp(φi ) ∩ supp(φj ) = φ

Examinons une possibilité de stockage de la matrice K. Soit N K le nombre d’éléments non nuls de la matrice K On
peut stocker la matrice dans un seul tableau KM AT en mettant bout à bout les coefficients non nuls de la première
ligne, puis ceux de la deuxième ligne, etc... jusqu’à ceux de la dernieère ligne. Pour repérer les éléments de K dans
le tableau KM AT , on a alors besoin de pointeurs. Le premier pointeur, nommé, IC est de dimension N K. La
valeur de IC(k) est le numéro de la colonne de K(k). On introduit alors le pointeur IL(ℓ), ℓ = 1, . . . , N L, où
NL est le nombre de lignes, où IL(ℓ) est l’indice dans KM AT du début de la ℓ-ième ligne. L’identification entre
KM AT et K se fait alors par la procédure suivante :
Pour k = 1 . . . N K
si IL(m) ≤ k < IL(m + 1) alors
KM AT (k) = Km,IC(k)
Fin si
Fin pour

Analyse numérique des EDP, M1 154 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.4. ELÉMENTS FINIS ISOPARAMÉTRIQUES CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

Kℓ

Fℓ

K̂

F IGURE 4.8 – Transformation isoparamétrique

La matrice K est symétrique définie positive, on peut donc utiliser une méthode de type gradient conjugué préconditionné
(voir cours de Licence). Notons que la structure de la matrice dépend de la numérotation des noeuds. Il est donc
important d’utiliser des algorithmes performants de maillage et de numérotation.

4.4 Eléments finis isoparamétriques

4.4 Eléments finis isoparamétriques
Dans le cas ou Ω est polygonal, si on utilise des éléments finis de type P2 , les noeuds de la frontière sont effec-
tivement sur la frontière même si on les calcule à partir de l’élément fini de référence. Par contre, si le bord est
courbe, ce n’est plus vrai. L’utilisation d’éléments finis “isoparamétriques” va permettre de faire en sorte que tous
les noeuds frontières soient effectivement sur le bord, comme sur la figure 4.8. Pour obtenir une transformation
isoparamétrique, on définit
Fℓ : K → Kℓ
(x̄, ȳ) !→ (x, y)
à partir des fonctions de base de l’élément fini de référence :
Nℓ
% Nℓ
%
x= xr φ̄r (x̄, ȳ), y = yr φ̄r (x̄, ȳ),
r=1 r=1

où Nℓ est le nombre de noeuds de l’élément et (xr , yr ) sont les coordonnées du r-ième noeud de Kℓ . Remarquons
que la transformation Fℓ isoparamétrique P1 est identique à celle des éléments finis classiques. Par contre, la
transformation isoparamétrique P2 n’est plus affine, alors qu’elle l’est en éléments finis classiques. Notons que les
fonctions de base locales vérifient toujours

φℓr ◦ Fℓ = φr , ∀ℓ = 1, . . . , L, ∀r = 1, . . . , Nℓ .

On peut alors se poser le problème de l’inversibilité de Fℓ . On ne peut pas malheureusement démontrer que Fℓ
est inversible dans tous les cas, toutefois, cela s’avère être le cas dans la plupart des cas pratiques. L’intérêt de
la transformation isoparamétrique est de pouvoir traiter les bords courbes, ainsi que les éléments finis Q1 sur
quadrilatères. Notons que le calcul de φℓr est toujours inutile, car on se ramène encore à l’élément de référence.

4.5 Analyse d’erreur

4.5.1 Erreurs de discrétisation et d’interpolation
On considère toujours le problème modèle (4.11) page 146 sur lequel on a étudié la mise en oeuvre de la méthode
des éléments finis. On rappelle que la formulation faible de ce problème est donnée en (4.14) page 147, et que

Analyse numérique des EDP, M1 155 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.5. ANALYSE D’ERREUR CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

sous les hypothèses (4.12) page 146, le problème (4.14) admet une unique solution u  ∈ H = HΓ10 n = {u ∈
1
H (Ω); u = 0 sur Γ0 }. La méthode d’approximation variationnelle du problème (4.14) consiste à chercher u N ∈
HN = V ect{φ1 , . . . , φN } solution de (4.16) page 147, où les fonctions φ1 , . . . , φN sont les fonctions de base
éléments finis associés aux noeuds x1 , . . . , xN . Comme les hypothèses (3.21) page 107 sont vérifiées, l’estimation
(3.29) page 110 entre u  solution de (4.14) et u (N ) solution de (4.16) est donc vérifiée. On a donc :

M
" u−u N "H 1 ≤ d(
u, HN ),
α

(resp.α) est la constante de continuité (resp. de coercivité) de a. Comme u = u0 + u

où M , on a, en posant
M
c= α ,
"u − uN " ≤ C"u − w"∀w ∈ HN , (4.20)
où uN = u N + u0 . Notons que dans l’implantation pratique de la méthode d’éléments finis, lorsqu’on calcule
T (v) = T (v) − a(u0 , v), on remplace u0 par u0,N ∈ HN , donc on commet une légère erreur sur T . De plus,
on calcule a(φi , φj ) à l’aide d’intégrations numériques : l’inégalité (4.20) n’est donc vérifiée en pratique que
de manière approchée. On supposera cependant, dans la suite de ce paragraphe, que les erreurs commises sont
négligeables et que l’inégalité (4.20) est bien vérifiée. De la même manière qu’on a défini l’interpolée sur un
élément K, (voir définition 4.3 page 137, on va maintenant définir l’interpolée sur H 1 (Ω) tout entier, de manière
à établir une majoration de l’erreur de discrétisation grâce à (4.20).

Définition 4.16 (Interpolée dans HN ) . Soit u ∈ H 1 (Ω) et HN = V ect{φ1 , . . . , φN } où les fonctions φ1 . . . φN
sont des fonctions de base éléments finis associées aux noeuds S1 . . . SN d’un maillage éléments finis de Ω. Alors
on définit l’interpolée de u dans HN , uI ∈ HN par :
N
%
uI = u(Si )φi .
 %
uI = u(Si )φi .
i=1

Comme uI ∈ HN , on peut prendre W = uI dans (4.20), ce qui fournit un majorant de l’erreur de discrétisation :

"u − uN "H 1 ≤ C"u − uI "H 1

On appelle erreur d’interpolation le terme "u − uI "H 1

4.5.2 Erreur d’interpolation en dimension 1

Soit Ω =]0, 1[, on considère un maillage classique, défini par les N + 2 points (xi )i=0...N +1 , avec x0 = 0 et
xN +1 = 1, et on note

hi = xi+1 − xi , i = 0, . . . , N + 1, et h = max{|hi |, i = 0, . . . , N + 1}

On va montrer que si u ∈ H 2 (]0, 1[), alors on peut obtenir une majoration de l’erreur d’interpolation "u − uI "H 1 .
On rappelle (voir exercice 34 page 119) que si Si u ∈ H 1 (]0, 1[) alors u est continue.
En particulier, on a donc H 2 (]0, 1[) ⊂ C 1 ([0, 1]). Remarquons que ce résultat est lié à la dimension 1, voir injection
de Sobolev, cours d’analyse fonctionnelle ou [1]. On va démontrer le résultat suivant sur l’erreur d’interpolation.

Théorème 4.17 (majoration de l’erreur d’interpolation, dimension 1) Soit u ∈ H 2 (]0, 1[), et soit uI son in-
terpolée sur HN = V ect{φi , i = 1, . . . , N }, où φi désigne la i-ème fonction de base élément fini P 1 associée au
noeud xi d’un maillage élément fini de ]0, 1[. Alors il existe C ∈ IR ne dépendant que de u, tel que

"u − uI "H 1 ≤ Ch. (4.21)

Démonstration : On veut estimer

"u − uI "2H 1 = |u − uI |20 + |u − uI |21

où |v|0 = "v"L2 et |v|1 = "Dv"L2 . Calculons |u − uI |21 :

( 1 N (
% xi+1
|u − uI |21 = |u′ − u′I |2 dx = |u′ (x) − u′I (x)|2 dx.
0 i=0 xi

Analyse numérique des EDP, M1 156 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.5. ANALYSE D’ERREUR CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

Or pour x ∈]xi , xi+1 [ on a

u(xi+1 ) − u(xi )
u′I = = u′ (ξi ),
hi
pour un certain ξi ∈]xi , xi+1 [. On a donc :
( xi+1 ( xi+1
|u′ (x) − u′I (x)|2 dx = |u′ (x) − u′ (ξi )|2 dx.
xi xi

On en déduit que : ( ( (
xi+1 xi+1 x
|u′ (x) − u′I (x)|2 dx = | u′′ (t)dt|2 dx,
xi xi ξi
et donc, par l’inégalité de Cauchy-Schwarz,
( xi+1
" xi+1 " x
|u′ (x) − u′I (x)|2 dx ≤ xi ξi |u′′ (t)|2 dt|x − ξi |dx
xi
( xi+1 $( x '
′′ 2
≤ hi |u (t)| dt dx,
xi ξi

car |x − ξi | ≤ hi . En réappliquant l’inégalité de Cauchy-Schwarz, on obtient :

( xi+1 ( xi+1
2 2
′ ′
|u (x) − uI (x)| dx ≤ hi |u′′ (t)|2 dt
xi xi

En sommant sur i, ceci entraine : ( 1

|u − uI |21 ≤ h2 |u′′ (t)|2 dt. (4.22)
0
"1
Il reste maintenant à majorer |u − uI |20 = |u − uI |2 dx. Pour x ∈ [xi , xi+1 ]
 "
Il reste maintenant à majorer |u − uI |20 = |u − uI |2 dx. Pour x ∈ [xi , xi+1 ]
0
$( x '2
2 ′ ′
|u(x) − uI (x)| = (u (t) − uI (t))dt .
xi

Par l’inégalité de Cauchy-Schwarz, on a donc :

( x
2
|u(x) − uI (x)| ≤ (u′ (t) − u′I (t))2 dt |x − xi |
xi & #! )
≤hi

Par des calculs similaires aux précédents, on obtient donc :

( x $( xi+1 '
|u(x) − uI (x)|2 ≤ hi |u′′ (t)|2 dt dxhi
xi xi
( xi+1
≤ h3i |u′′ (t)|2 dt.
xi

En intégrant sur [xi , xi+1 ], il vient :

( xi+1 ( xi+1
2
|u(x) − uI (x)| dx ≤ h4i (u′′ (t))2 dt,
xi xi

et en sommant sur i = 1, . . . , N :
( 1 ( 1
2 4
(u(x) − uI (x)) dx ≤ h (u′′ (t))2 dt.
0 0

On a donc :
|u − uI |0 ≤ h2 |u|2
ce qui entraine, avec (4.22) :
"u − uI "2 ≤ h4 |u|22 + h2 |u|22
≤ (1 + h2 )h2 |u|22
On en déduit le résultat annoncé.
On en déduit le résultat d’estimation d’erreur suivant :

Analyse numérique des EDP, M1 157 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.5. ANALYSE D’ERREUR CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

Corollaire 4.18 (Estimation d’erreur, P1, dimension 1) Soit Ω =]0, 1[ ; soit f ∈ L2 (Ω) et u ∈ H01 (Ω) l’unique
solution du problème 
 u ∈ H01 (Ω)

( ( ,

 a(u, v) = ∇u(x)∇v(x)dx = f (x)v(x)dx,
Ω

et uT L’approximation éléments finis P1 obtenue sur un maillage admissible T de pas hT = max {|hi |}. Alors
i=1,...,N
il existe C ∈ IR ne dépendant que de Ω et f tel que "u − uT " < Ch.

Ce résultat se généralise au cas de plusieurs dimensions d’espace (voir Ciarlet) ; si Ω un ouvert polygonal convexe
de IR d , d ≥ 1, le résultat des conditions géométriques sur le maillage, nécessaires pour obtenir le résultat d’interpo-
lation. Par exemple pour un maillage triangulaire en deux dimensions d’espace, intervient une condition d’angle :
on demande que la famille de maillages considérée soit telle qu’il existe β > 0 tel que β ≤ θ ≤ π − β pour tout
angle θ du maillage.
Remarque 4.19 (Sur les techniques d’estimation d’erreur) Lorsqu’on a voulu montrer des estimations d’er-
reur pour la méthode des différences finies, on a utilisé le principe de positivité, la consistance et la stabilité
en norme L∞ . En volumes finis et éléments finis, on n’utilise pas le principe de positivité. En volumes finis, la
stabilité en norme L2 est obtenue grâce à l’inégalité de Poincaré discrète, et la consistance est en fait la consis-
tance des flux. Notons qu’en volumes finis on se sert aussi de la conservativité des flux numériques pour la preuve
de convergence. Enfin, en éléments finis, la stabilité est obtenue grâce à la coercivité de la forme bilinéaire, et la
consistance provient du contrôle de l’erreur d’interpolation.
Même si le principe de positivité n’est pas explicitement utilisé pour les preuves de convergence des éléments finis
et volumes finis, il est toutefois intéressant de voir à quelles conditions ce principe est respecté, car il est parfois
très important en pratique.
Reprenons d’abord le cas du schéma volumes finis sur un maillage T admissible pour la discrétisation de l’équation
tres important en pratique.
Reprenons d’abord le cas du schéma volumes finis sur un maillage T admissible pour la discrétisation de l’équation
(3.1). 1
−∆u = f dans Ω
u=0 sur ∂Ω.
Rappelons que le schéma volumes finis s’écrit :
 
% % %
 τK,L (uK − uL ) + τK,σ uK  = |K|fK , (4.23)
K∈C σ∈ξK ∩ξint σ∈ξK ∩ξext

avec
|K|L| |σ|
τK,L = et τK,σ = ,
d(xK , xL ) d(xK , ∂Ω)
où |K|, (resp. |σ|) désigne la mesure de Lebesque en dimension d (resp. d − 1) de K (resp. σ).
Notons que les coefficients τK,L et τK,σ sont positifs, grâce au fait que le maillage est admissible (et donc
XKHXL = d(XK , XL )nKL H , où XKHXL désigne le vecteur d’extrémités XK et XL et nKL H la normale unitaire
à K|L sortante de K.
Notons que le schéma (4.23) s’écrit comme une somme de termes d’échange entre les mailles K et L, avec des
coefficients τKL positifs. C’est grâce à cette propriété que l’on montre facilement que le principe de positivité
est vérifié. Considérons maintenant la méthode des éléments finis P1, pour la résolution du problème (3.1) sur
maillage triangulaire. On sait (voir par exemple Ciarlet) que si le maillage satisfait la condition faible de Delau-
nay (qui stipule que la somme de deux angles opposés à une même arête doit être inférieure à π), alors le principe
du maximum est vérifié. Ce résultat peut se retrouver en écrivant le schéma éléments finis sous la forme d’un
schéma volumes finis.

4.5.3 Super convergence

On considère ici un ouvert Ω polygonal convexe de IR d , d ≥ 1, et on suppose que f ∈ L2 (Ω). On s’intéresse à
l’approximation par éléments finis P 1 de la solution u ∈ H01 (Ω) du problème (3.5). On a vu dans le paragraphe
précédent (corollaire 4.18) qu’on peut estimer l’erreur en norme L2 entre la solution exacte u et la solution ap-
prochée par éléments finis P 1 ; en effet, comme l’erreur d’interpolation est d’ordre h, on déduit une estimation

Analyse numérique des EDP, M1 158 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.5. ANALYSE D’ERREUR CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

sur l’erreur de discrétisation, également d’ordre h. En fait, si la solution u de (3.1) est dans H 2 , il se produit un
“petit miracle”, car on peut montrer grâce à une technique astucieuse, dite“truc d’Aubin-Nitsche”, que l’erreur de
discrétisation en norme L2 est en fait d’ordre 2.

Théorème 4.20 (Super convergence des éléments finis P 1) Soit Ω un ouvert polygonal convexe de IR d , d ≥ 1 ;
soit f ∈ L2 (Ω), u solution de (3.5), uT la solution apporchée obtenue par éléments finis P1, sur un maillage
éléments finis T . Soit
hT = max diamK.
K∈T

Alors il existe C ∈ IR ne dépendant que de Ω et f tel que :

"u − uT "H 1 (Ω) ≤ Ch et "u − uT "L2 (Ω) ≤ Ch2 .

Démonstration : Par le théorème de régularité 3.9 page 102, il existe C1 ∈ IR + ne dépendant que de Ω tel que

"u"H 2 (Ω) ≤ C1 "f "L2 (Ω) .

Grâce à ce résultat, on a obtenu (voir le théorème 4.18) qu’il existe C2 ne dépendant que de Ω, β et tel que

"u − uT "H 1 (Ω) ≤ C2 "f "L2 h

Soit maintenant eT = u − uT et ϕ ∈ H01 (Ω) vérifiant

( (
∇ϕ(x).∇ψ(x)dx = eT (x)ψ(x)dx, ∀ψ ∈ H01 (Ω). (4.24)
Ω Ω

On peut aussi dire que ϕ est la solution faible du problème

On peut aussi dire que est la solution faible du problème
,
−∆ϕ = eT dans Ω
ϕ = 0 sur ∂Ω.

Comme e ∈ L2 (Ω), par le théorème 3.9, il existe C3 ∈ IR + ne dépendant que Ω tel que

"ϕ"H 2 (Ω) ≤ C3 "eT "L2 (Ω) .

( (
Or "eT "2L2 (Ω) = eT (x)eT (x)dx = ∇ϕ(x).∇e(x)dx, en prenant ψ = eT dans (4.24).
Ω Ω
Soit ϕT la solution approchée par éléments finis P1 du problème (4.24), c.à.d solution de :

 ϕT ∈ VT ,0 = {v ∈ C(Ω̄); v|K ∈ P1 , ∀K ∈ T , v|∂Ω = 0}

( (
(4.25)

 ∇ϕT (x)∇v(x)dx = e(x)v(x)dx, ∀v ∈ VT ,0
Ω Ω

On sait que uT vérifie : (

∇ϕT (x).∇(u − uT )(x)dx = 0;
Ω
on peut donc écrire que :
(
"eT "2L2 (Ω = ∇(ϕ − ϕT )(x).∇(u − uT )(x)dx ≥ "ϕ − ϕT "H 1 (Ω) "u − uT "H 1 (Ω)
Ω

D’après le théorème 4.18, on a :

"ϕ − ϕT "H 1 (Ω) ≤ C2 "e"L2(Ω) hT et "u − uT "H 1 (Ω) ≤ C2 "f "L2(Ω) hT .

On en déduit que :
"eT "L2 (Ω) ≤ C22 "f "L2 (Ω) h2T .
Ce qui démontre le théorème.

Analyse numérique des EDP, M1 159 Université Aix-Marseille 1, R. Herbin, 2 février 2012
4.6. EXERCICES CHAPITRE 4. ELÉMÉNTS FINIS DE LAGRANGE

4.5.4 Traitement des singularités

Les estimations d’erreur obtenues au paragraphe précédent reposent sur la régularité H 2 de u. Que se passe-t-il si
cette régularité n’est plus vérifiée ? Par exemple, si le domaine Ω possède un coin rentrant, on sait que dans ce cas,
la solution u du problème (3.5) n’est plus dans H 2 (Ω), mais dans un espace H 1+s (Ω), où s dépend de l’angle du
coin rentrant. Considérons donc pour fixer les idées le problème
1
−∆u = f dans Ω, où Ω est un ouvert
u = 0 sur ∂Ω polygônal avec un coin rentrant.

Pour approcher correctement la singularité, on peut raffiner le maillage dans le voisinage du coin. On peut également,
lorsque cela est possible, modifier l’espace d’approximation pour tenir compte de la singularité. Dans le cas d’un
polygône avec un coin rentrant par exemple, on sait trouver ψ ∈ H01 (Ω) (et ψ ∈ H 2 (Ω)) telle que si u est solution
de (3.5) avec f ∈ L2 (Ω), alors il existe un unique α ∈ IR tel que u − αψ ∈ H 2 (Ω).
Examinons le cas d’une approximation par éléments finis de Lagrange. Dans le cas où u est régulière, l’espace
d’approximation est
VT = V ect{φi , i = 1, NT },
où NT est le nombre de noeuds internes du maillage T de Ω considéré et (φi )i=1,NT la famille des fonctions de
forme associées aux noeuds.
Dans le cas d’une singularité portée par la fonction ψ introduite ci-dessus, on modifie l’espace V et on prend
maintenant : VT = V ect{φi , i = 1, NT } ⊕ IRψ Notons que VT ⊂ H01 (Ω), car ψ ∈ H01 (Ω). Reprenons maintenant
l’estimation d’erreur. Grâce au lemme de Céa, on a toujours
M
"u − uT "H 1 ≤ "u − w"H 1 (Ω) , ∀w ∈ VT .
α
On a donc également :
On a donc également :
M
"u − uT "H 1 ≤ "u − αψ − w"H 1 (Ω) , ∀w ∈ VT .
α
puisque αψ + w ∈ VT . Or, u − αψ = u  ∈ H 2 (Ω).
Donc "u − uT "H 1 ≤ M α "
u − w" H (Ω) , ∀w ∈ VT .
1

Et grâce aux résultats d’interpolation qu’on a admis, si on note u

I l’interpolée de u
 dans VT , on a :
M M
"u − uT "H 1 (Ω) ≤ "
u−u
I "H 1 (Ω) ≤ C2 "
u"H 2 (Ω) h.
α α
On obtient donc encore une estimation d’erreur en h.
Examinons maintenant le système linéaire obtenu avec cette nouvelle approximation. On effectue un développement
de Galerkin sur la base de VT . On pose %
uT = ui φi + γψ.
i=1,NT

Le problème discrétisé revient donc à chercher



 (us )i=1,NT ⊂ IR N et γ ∈ IR t.q.
  " " "
j=1,NT uj Ω ∇φj (x) · ∇φi (x)dx + γ Ω ∇ψ(x) · ∇φi (x)dx = Ω f (x)φi (x)dx, ∀i = 1, NT

 " " "
 
j=1,NT uj Ω ∇φi (x) · ∇ψ(x)dx + γ Ω ∇ψ(x) · ∇ψ(x)dx = f (x)ψ(x)dx.

On obtient donc un système linéaire de NT + 1 équations à NT + 1 inconnues.

4.6 Exercices
Exercice 48 (Eléments finis P1 pour le problème de Dirichlet) Corrigé en page 165

Soit f ∈ L2 (]0, 1[). On s’intéresse au problème suivant :

−u′′ (x) = f (x), x ∈]0, 1[,

u(0) = 0, u(1) = 0.

Analyse numérique des EDP, M1 160 Université Aix-Marseille 1, R. Herbin, 2 février 2012