Dossier délivré pour

DOCUMENTATION
19/09/2008

Introduction à la méthode
des éléments finis
par Pierre SPITERI
Docteur ès sciences mathématiques
Professeur à l’École nationale supérieure d’électronique, d’électrotechnique,
d’informatique, d’hydraulique et de télécommunication de Toulouse (ENSEEIHT)

1. Méthode de Galerkin............................................................................... AF 504 – 2

1.1 Exposé général de la méthode ................................................................... — 2
1.2 Exemples en dimension 1 et introduction de la méthode
des éléments finis........................................................................................ — 3
2. Introduction des éléments finis en dimension 2............................. — 7
2.1 Interpolation par des polynômes de degré 1 ............................................ — 7
2.2 Exemples de matrices de discrétisation .................................................... — 10
Références bibliographiques ......................................................................... — 12

n a vu dans l’article [AF 503] qu’une équation aux dérivées partielles ellip-
O tique pouvait être exprimée sous diverses formulations équivalentes, en ce
sens que toute solution d’une formulation est solution d’une autre formulation
et réciproquement. La formulation forte du problème présente un intérêt dans la
mesure où l’utilisation de la méthode des différences finies est envisagée pour
effectuer une approximation du problème. La formulation équivalente du pro-
blème basée sur la formulation d’un problème d’optimisation associé à la fonc-
tionnelle v → J ( v ) , avec J ( v ) définie par :
1
J ( v ) = --- a ( v, v ) – L ( v )
2
nécessite que la forme bilinéaire a ( u, v ) soit symétrique, ce qui en soit est res-
trictif dans la mesure où certains phénomènes sont modélisés à partir d’opéra-
teurs non autoadjoints. Cependant, lorsque a(.,.) est symétrique, cette
formulation du problème conduit à la méthode de Ritz ; numériquement, l’idée
est de chercher à minimiser J(.) non plus sur l’ensemble E tout entier, mais sur
un sous-espace de E construit à partir de fonctions facilement calculables ; la
fonction inconnue qui réalise le minimum est représentée comme combinaison
linéaire de fonction de forme (ou de tout autre famille physiquement admissible)
et les coefficients de cette combinaison linéaire sont les inconnues du problème.
J(.) est alors transformée en une fonctionnelle quadratique et déterminer le
minimum de cette nouvelle fonctionnelle revient alors à annuler les dérivées
partielles de celle-ci par rapport à ces inconnues, ce qui conduit classiquement à
la résolution d’un système linéaire. Nous ne développerons pas cette méthode.
L’autre formulation équivalente mise en évidence est la formulation faible ou
formulation variationnelle basée sur le théorème des travaux virtuels ; cette
expression équivalente du problème contient l’ensemble des informations rela-
tives à ce dernier, c’est-à-dire l’équation aux dérivées partielles et les conditions
aux limites ; elle offre, de plus, une grande possibilité de calculs effectifs des
solutions approchées en choisissant des sous-ensembles de l’espace fonction-
nel E bien adaptés au calcul numérique. La méthode de base de l’approximation
est la méthode de Galerkin, qui consiste à choisir la fonction inconnue sous
forme de combinaison linéaire de fonctions de forme de manière à obtenir un
système discret en choisissant successivement comme fonction test dans la
formulation faible ces mêmes fonctions de forme. Bien évidemment, dans le cas

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
© Techniques de l’Ingénieur, traité Sciences fondamentales AF 504 − 1

Dossier délivré pour

DOCUMENTATION
19/09/2008
 Dossier délivré pour
DOCUMENTATION
19/09/2008
MÉTHODE DES ÉLÉMENTS FINIS __________________________________________________________________________________________________________

où les sous-espaces de E sont identiques, les systèmes linéaires dérivés de la

méthode de Ritz et de la méthode de Galerkin sont les mêmes. Un cas particulier
de la méthode de Galerkin est la méthode des éléments finis que nous allons
développer dans cet article et dont le grand intérêt est constitué par sa grande
souplesse.
Cet article est le second volet d’un ensemble de trois articles traitant de la méthode des élé-
ments finis :
— [AF 503] Approche variationnelle pour la méthode des éléments finis ;
— [AF 504] Introduction à la méthode des éléments finis ;
— [AF 505] Présentation générale de la méthode des éléments finis.

1. Méthode de Galerkin Théorème 1. La forme bilinéaire étant symétrique et coercive,

alors la résolution du problème (3) aboutit à la résolution d’un
système algébrique linéaire dont la matrice A = ( a ( ϕ j , ϕ i ) ) est
symétrique, définie positive.
1.1 Exposé général de la méthode
Preuve ♦
On a vu dans l’article [AF 503] que divers problèmes d’équation
aux dérivées partielles pouvaient se mettre sous la forme variation- Écrivons le problème (3), en remplaçant u h par son expression
nelle suivante : donnée par la relation (2); on a donc, en choisissant
v h = ϕ i ( x ) , ∀i ∈ { 1, …, m } et à cause de la bilinéarité de a (.,.) :
 déterminer u ∈ E tel que : m
 (1)
 a ( u, v ) = L ( v ), ∀ v ∈ E ∑ a ( ϕ j ( x ) , ϕ i ( x ) ) µ j = L ( ϕ i ( x ) ) , ∀ i ∈ { 1, … , m }
j=1
avec a ( u, v ) forme bilinéaire,
L(v) forme linéaire, Cette dernière relation définit un système linéaire dont les incon-
j=m
E espace vectoriel normé de dimension infinie. nues sont les nombres { µ j } j = 1 , qui représentent en fait les compo-
Bien qu’il existe des méthodes d’approximation de la solution j=m
santes de u h dans la base { ϕ j } j = 1 . Or, la forme bilinéaire étant
d’équation aux dérivées partielles dans des espaces de dimension
infinie (en particulier décomposition de la solution en série de Fou- symétrique on a :
rier), la principale difficulté du traitement numérique est liée à cette
dernière caractéristique. Une des techniques pour contourner cette a ( ϕ j ( x ) , ϕ i ( x ) ) = a ( ϕ i ( x ) , ϕ j ( x ) ) , ∀ i, j ∈ { 1 , … , m }
difficulté est de projeter sur un espace noté E h de dimension finie,
inclus dans l’espace fonctionnel E. La matrice A = ( a ( ϕ j , ϕ i ) ) est donc symétrique. De plus soit Π le
Soit donc E h ⊂ E un sous-espace vectoriel de E de dimension j=m
vecteur de composantes { µ j } j = 1 ; formons :
finie m ; on considérera dans la suite que E h est engendré par les
fonctions de base ϕ j ( x ), j = 1, …, m . On cherche la solution du pro-
m m
blème, qu’on notera à présent u h , pour spécifier qu’elle appartient     
à E h , sous la forme d’une combinaison linéaire des fonctions tests AΠ = 

∑ a ( ϕj ( x ), ϕi ( x ) ) µj  = a  ∑ µj ϕj ( x ), ϕi ( x )  = (a ( u h , ϕ i ( x ) ) )
ϕ j ( x ) ∈ E h , soit : j=1 j=1
m
et
uh ( x ) = ∑ µj ϕj ( x ) (2)
j=1 m m
 
Π A Π = ∑ a ( u h , ϕ j ( x ) ) µ j = a u h , ∑ µ j ϕ j ( x ) = a ( u h , u h )  α u h
t 2
et l’on considère le problème en dimension finie associé au pro- Eh
>0
blème (1), que l’on écrit comme suit : j=1  j=1 
t
 déterminer u h ∈ E h tel que : et Π A Π > 0 dès que Π ≠ 0 ; donc la matrice A est bien symétrique,
 (3) définie positive. ♦
 a ( u h , v h ) = L ( v h ), ∀ v h ∈ E h

Remarque
Corollaire 1. La matrice A = ( a ( ϕ i, ϕ j ) ) est inversible.
Si l’on suppose que les hypothèses du théorème de Lax-Milgram
sont vérifiées, c’est-à-dire, E étant un espace de Hilbert, que la forme
bilinéaire a(.,.) est continue et coercive et que la forme linéaire L(.)
Remarque
est continue, alors, par application de ce même résultat, le problème
(3) admet une solution unique ; en effet, l’espace E h étant un sous- La résolution du système linéaire précédent sera d’autant plus
espace vectoriel de dimension finie, c’est un sous-espace fermé de aisée que la matrice A a la structure la plus creuse possible et même
E et, par conséquent, c’est un espace de Hilbert pour la norme possède une structure bande. L’obtention d’une telle structure
induite par celle de E ; on peut donc appliquer le théorème de Lax- dépend essentiellement du choix :
Milgram en remplaçant E par E h . De plus, on peut facilement véri- — du sous-espace E h et, en particulier, de celui de la base
fier que le problème (3) est équivalent à un problème de minimisa- { ϕj } j = m ;
tion de la fonctionnelle J ( u h ) dans E h . On va préciser ce résultat en j=1
donnant des démonstrations plus simples. — de la numérotation des points du maillage.

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
AF 504 − 2 © Techniques de l’Ingénieur, traité Sciences fondamentales

Dossier délivré pour

DOCUMENTATION
19/09/2008
Dossier délivré pour
DOCUMENTATION
19/09/2008
__________________________________________________________________________________________________________ MÉTHODE DES ÉLÉMENTS FINIS

Le fait de projeter la solution de l’espace de dimension infinie E calculer, à titre d’exercice, les coefficients de la matrice A ; tous
dans l’espace de dimension finie E h a pour effet de faire perdre de calculs faits, on trouve :
l’information sur le problème initial ; en fait, dans cette opération,
1 dϕ
i ( x ) dϕ j ( x )
on commet une erreur dont on peut donner, sans démonstration
(cf. [1]), une majoration :
a ( ϕ i ( x ), ϕ j ( x ) ) =
∫ 0
-----------------
dx
- ------------------ dx
dx

∫
1
Théorème 2. Soit u la solution du problème (1) et u h celle du =  ( i + 1 ) ( 1 – x )x i – x i + 1  ( j + 1 ) ( 1 – x )x j – x j + 1 dx
0   
problème (3). On a alors la majoration :
soit :
M  
u – uh  ----- Inf  u – v h , ∀v h ∈ E h 
E α 
E
 2(i + 1)(j + 1)
a ( ϕ i ( x ), ϕ j ( x ) ) = -------------------------------------------------------------------------
(i + j + 1)(i + j + 2)(i + j + 3)

Remarque On remarque que a ( ϕ i ( x ), ϕ j ( x ) ) > 0, et en particulier

Le résultat précédent montre que la qualité de l’approximation a ( ϕ i ( x ), ϕ j ( x ) ) ≠ 0. Par conséquent, bien que la base soit simple, on
des fonctions de l’espace E par celles de l’espace E h dépend de la
richesse de ce dernier espace. obtient une matrice de discrétisation pleine. Ainsi, si l’on résout le
système linéaire obtenu par la méthode de Gauss, on devra effec-
Remarque
m3
Lorsque la forme bilinéaire a (.,.) est symétrique, on peut obtenir tuer, compte tenu de la symétrie de la matrice, de l’ordre de --------
3
une majoration plus fine du type : opérations arithmétiques.

M   ■ Second choix de base

u – uh  ----- Inf  u – v h , ∀v h ∈ E h 
E α 
E
 Soit E h l’espace engendré par les fonctions ϕ j ( x ) = sin ( jπx ) ,
1  j  m . De la même façon que lors de la situation précédente, on
remarque que ϕ j ( 0 ) = ϕ j ( 1 ) , pour tout indice j ∈ [ 0, m – 1 ] , et les
fonctions ϕ j ( x ) vérifient les conditions aux bords. De plus, les fonc-
1.2 Exemples en dimension 1 tions ϕ j ( x ) sont indéfiniment différentiables, de carré intégrable,
et introduction de la méthode ainsi que leur dérivées. Par conséquent ϕ j ( x ) ∈ H 01 ( [ 0, 1 ] ), pour tout
indice j ∈ [ 1, m ] et l’on a bien l’inclusion E h ⊂ E = H 01 ( [ 0, 1 ] ) . De
des éléments finis plus :
1 dϕ
i ( x ) dϕ j ( x )
∫ ∫
2 1
Dans ce paragraphe, nous allons examiner un certain nombre a ( ϕ i ( x ), ϕ j ( x ) ) = - ------------------ dx = ij π
----------------- cos ( iπx ) cos ( jπx ) dx
d’exemples mettant en œuvre la méthode de Galerkin et qui permet- 0 dx dx 0
tront d’introduire la méthode des éléments finis.
Tous calculs faits, on trouve :
Exemple 1. Considérons la situation la plus simple correspondant à
l’équation de Poisson, définie dans le segment [0,1] et munie des 0 si i ≠ j

conditions de Dirichlet homogènes : a ( ϕ i ( x ), ϕ j ( x ) ) =  i 2 π 2
 ------------ si i = j
 d 2u  2
 – ---------2- = f ( x ), x ∈ ]0, 1[
 dx Donc, la matrice A est diagonale. Cela provient du choix parti-
 culier de fonctions de base de l’espace E h . En effet, on vérifie sim-
 u(0) = u(1) = 0
plement que les fonctions ϕ j ( x ) considérées sont les fonctions pro-
On a vu au paragraphe 2 de l’article [AF 503] que la formulation varia-
tionnelle associée à ce problème est donnée par : d2
pres de l’opérateur – ---------- associé aux conditions de Dirichlet
dx 2

 déterminer u ∈ H 0 ( [ 0, 1 ] ) tel que :
1
homogènes. On a donc obtenu la matrice A la plus creuse possible
 et l’on peut penser avoir réglé le problème du choix optimal des
 a ( u, v ) = L ( v ), ∀v ∈ H 01 ( [ 0, 1 ] )
 fonctions de base de l’espace E h , en choisissant les fonctions pro-
avec : pres de l’opérateur aux dérivées partielles considéré. En fait, il n’en
est rien, car si l’on change légèrement le problème aux limites en

∫ ∫
1 du dv 1
a ( u, v ) = ------- ------dx et L ( v ) = f v dx considérant le problème à coefficient variable :
0 dx dx 0
 d  du ( x )
 – -------  c ( x ) ---------------- = f ( x ), x ∈ ]0, 1[
1
pour construire le sous-espace de dimension finie de H 0 ( [ 0, 1 ] ) , on va
 dx dx
considérer plusieurs choix de fonctions de base { ϕ j } j = 1 .
j=m 
 u ( 0 ) = 0, u ( 1 ) = 0
■ Premier choix de base auquel on associe, la forme bilinéaire :
j
Soit E h l’espace engendré par les fonctions ϕ j ( x ) = x ( 1 – x )x ,

∫
0  j  m – 1 . Remarquons que, pour tout indice j ∈ [ 0, m – 1 ] , les
1 du ( x ) d v ( x )
a ( u, v ) = c ( x ) ---------------- --------------- dx
fonctions ϕ j ( x ) sont nulles au bord ; de plus, ces fonctions étant 0 dx dx
polynomiales de degré inférieur ou égal à m + 1, elles sont de carré
intégrable ; en outre, ces fonctions sont indéfiniment différentiables les coefficients de la matrice sont donnés par :
et leurs dérivées étant polynomiales de degré inférieur ou égal à m

∫
1 1
sont également de carré intégrable. Donc ϕ j ( x ) ∈ H 0 ( [ 0, 1 ] ) , pour a ( ϕ i ( x ), ϕ j ( x ) ) = ij π
2
c ( x ) cos ( iπx ) cos ( jπx ) dx
tout indice j ∈ [ 0, m – 1 ] et l’on a bien l’inclusion E h ⊂ E . On peut 0

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
© Techniques de l’Ingénieur, traité Sciences fondamentales AF 504 − 3

Dossier délivré pour

DOCUMENTATION
19/09/2008
 Dossier délivré pour
DOCUMENTATION
19/09/2008
MÉTHODE DES ÉLÉMENTS FINIS __________________________________________________________________________________________________________

En général, pour une fonction c ( x ) quelconque, il n’y a aucune

raison pour que les coefficients a ( ϕ i ( x ), ϕ j ( x ) ) soient nuls pour ϕ j (x )
i ≠ j . Donc le choix particulier de fonctions de base considéré est
mal adapté pour un problème plus général, dans la mesure où dans
ce dernier cas, on obtient encore une matrice pleine. De plus, les
fonctions propres de l’opérateur différentiel à coefficients variables
sont en général peu commodes à déterminer.
■ Troisième choix de base
La définition de l’espace Eh va nous permettre d’introduire une
base élément fini pour la résolution numérique du problème aux
limites considéré, ce choix pouvant être étendu à des situations plus x j –1 xj xj +1
générales. On considère, pour simplifier, un découpage uniforme de
1
l’intervalle [ 0, 1 ] en ( m + 1 ) parties de longueur h, avec h = --------------- Figure 1 – Fonction de base linéaire par morceaux
m+1
et on pose x j = jh , j = 1, …, m, m + 1 , avec x 0 = 0 et x m + 1 = 1 .
On considère l’espace F h ⊂ E h défini comme suit :
ϕj – 1( x ) ϕ j (x ) ϕj + 1 (x )
0
F h = { v h ∈ C ( [ 0,1 ] ), v h linéaire sur chaque segment [ x j , x j + 1 ] }
0
avec C ( [ 0,1 ] ) l’ensemble des fonctions continues sur [ 0, 1 ] .
On pose v h ( x j ) = v j de telle sorte que, sur chaque intervalle
[ x j , x j +1 ] , on puisse interpoler la fonction v h ( x ) par un polynôme
de degré inférieur ou égal à un, ayant pour expression analytique :

x – xj
v h ( x ) = v j + -------------- ( v j + 1 – v j ) , ∀x ∈ [ x j , x j + 1 ]
h
xj – 1 xj xj + 1
De manière équivalente, on peut définir l’espace F h comme suit :
F h = { v h, v h [ x j, x j + 1 ]
∈ P 1 , v h entièrement déterminée par ses valeurs Figure 2 – Fonctions de base associées à 3 points successifs

v j aux points v j, 0  j  m + 1 }

avec P 1 l’ensemble des polynômes de degré inférieur ou égal à 1. Remarque

On pose : Sur ce dernier point, on retrouve une similitude avec la méthode
Eh = { vh ∈ Fh , vh ( 0 ) = vh ( 1 ) = 0 } des différences finies.

pour tenir compte des conditions aux limites de Dirichlet homogè- Il convient de noter une propriété très importante dans la
nes. Une base de l’espace E h est composée des fonctions construction des fonctions de base ; en effet dans la méthode des
ϕ j ( x ) ∈ E h , pour i = 1, …, m continues sur [ 0, 1 ] , linéaires sur cha- éléments finis, les fonctions de forme sont à supports disjoints,
que intervalle [ x j , x j + 1 ] , satisfaisant de plus ϕ j ( x i ) = δ i, j ( δ i, j sym- excepté pour les points strictement contigus, c’est-à-dire d’une part
bole de Kronecker). Ces fonctions de base sont appelées dans la ϕ j ( x ) et ϕ j + 1 ( x ) , ainsi que d’autre part ϕ j ( x ) et ϕ j – 1 ( x ) . Cette pro-
littérature les fonctions chapeaux et sont représentées sur les figu- priété permet d’obtenir une matrice A creuse car, en général, on a :
res 1 et 2. On a, de plus, dim ( E h ) = m . On a bien l’inclusion
E h ⊂ E = H 01 ( [ 0, 1 ] ) car E h ⊂ H 01 ( [ 0, 1 ] )C 0 ( [ 0,1 ] ) ; en effet les a ( ϕ j ( x ), ϕ j – 1 ( x ) ) ≠ 0 , a ( ϕ j ( x ), ϕ j ( x ) ) ≠ 0 , a ( ϕ j ( x ), ϕ j + 1 ( x ) ) ≠ 0
fonctions de E h , étant polynomiales par morceaux, sont des fonc-
tions de carré intégrable ; de plus, les dérivées premières des fonc-
tions ϕj ( x ) sont des fonctions en escalier, donc de carré intégrable. et
Par ailleurs, comme E h ⊂ H 01 ( [ 0, 1 ] ) , on a ϕ j ( 0 ) = ϕ j ( 1 ) = 0 ,
j = 1, … , m . a ( ϕ j ( x ), ϕ i ( x ) ) = 0 dès que i – j  2
Soit u h ∈ E h ; considérons le développement de u h dans la base
Donc, en général, la matrice A = ( a ( ϕ j ( x ), ϕ i ( x ) ) ) est une matrice
j=m
{ ϕ j } j = 1 ; on a : tridiagonale, dans le cas où les fonctions de base sont linéaires par
morceaux. Plus précisément, on peut calculer exactement les coeffi-
m cients de la matrice A ; auparavant il est nécessaire de déterminer
uh ( x ) = ∑ µj ϕj ( x ) ⇒ uh ( xi ) ≡ ui = µ i car ϕ j ( x i ) = δ i, j l’expression des fonctions de base et de leurs dérivées ; on procède
j=1 par identification en cherchant une fonction ϕ j ( x ) de la forme
donc : ϕ j ( x ) = ax + b (a et b inconnues du problème), vérifiant sur chaque
m intervalle la condition ϕ j ( x i ) = δ i, j ; tous calculs faits, on trouve :
uh ( x ) = ∑ ui ϕj ( x )
j=1 x – xj – 1
 --------------------
- sur [ x j – 1 , x j]  --1- sur [ x
 h j – 1 , x j]
h
et les coefficients du développement de u h dans la base
j=m
{ ϕj }j = 1  dϕ j ( x ) 
ϕ j ( x ) =  0 sinon , ------------------ =  0 sinon
 dx 
ne sont rien d’autre que les valeurs de la fonction u h aux points du x x
–
 --------------------
j+1
- sur [ x j , x j + 1 ]  – --1- sur [ x j , x j ]
maillage x j , 0  j  m + 1 .  h  h +1

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
AF 504 − 4 © Techniques de l’Ingénieur, traité Sciences fondamentales

Dossier délivré pour

DOCUMENTATION
19/09/2008
 Dossier délivré pour
DOCUMENTATION
19/09/2008
__________________________________________________________________________________________________________ MÉTHODE DES ÉLÉMENTS FINIS

On déduit aisément les coefficients de la matrice A de dimension – u ( x j + 1 ) + 2u ( x j ) – u ( x j – 1 ) 1

∫
xj + 1
m ; en effet : E j = ---------------------------------------------------------------------------
2
- – --- f ( x ) ϕ j ( x ) dx
h h xj – 1
2 2 2
1  dϕ  dϕ j x j + 1  dϕ
j j
∫ ∫ ∫
xj
2
a ( ϕj , ϕj ) =  -------- dx =  -------- dx +  -------- dx = --- avec u ( x ) solution exacte du problème de Poisson considéré.
0  dx  x j – 1  dx  xj  dx  h
En supposant que la solution exacte du problème est quatre fois
et continûment différentiable sur l’intervalle [ 0,1 ] et en considérant
les développements de u ( x j + 1 ) et u ( x j – 1 ) autour du point x j , on
1 dϕ dϕ j + 1 x j + 1 dϕ dϕ j + 1
∫ ∫
1 obtient :
a ( ϕj , ϕj + 1 ) = --------j ---------------
- dx = --------j ---------------
- dx = – ---
0 d x dx xj dx dx h
2
– u ( x j + 1 ) + 2u ( x j ) – u ( x j – 1 ) d u ( xj ) h2  ( 4 ) ^
- – ------ u ( x j ) + u ( x~j )
1 dϕ (4)
dϕ j – 1 d ϕ j dϕ j – 1
∫ ∫
xj ---------------------------------------------------------------------------
- = – -------------------
1 24  
a ( ϕj , ϕj – 1 ) = --------j --------------
- dx = -------- --------------- dx = – --- h
2
dx 2
0 dx dx x j – 1 dx dx h
(4)
Remarque avec u ( x ) dérivée quatrième au point x,
^ < x et x < x~ < x
xj – 1 < x
On vérifie bien que la matrice est symétrique ; elle est aussi défi- j j j j j + 1.
nie positive par application du résultat établi au théorème 1, ce que D’autre part, en supposant que la fonction f ( x ) est deux fois
l’on peut encore vérifier par un calcul direct simple. continûment différentiable, on peut la développer autour du point
Finalement on a à résoudre le système linéaire suivant : x j , ce qui conduit à :

∫ ∫
2
1 xj + 1 1 xj + 1 d f(ξ)
∫ ∫
 – u j – 1 + 2u j – u j + 1 xj xj + 1
--- f ( x ) ϕ j ( x ) dx = f ( x j ) + ------- ( x – x j ) ϕ j ( x ) ---------------- dx, x j < ξ < x
 -------------------------------------------------- = f ( x ) ϕ j ( x ) dx + f ( x ) ϕ j ( x ) dx, 1  j  m h 2h
 h xj – 1 xj
xj – 1 xj – 1 dx 2
 u0 = um + 1 = 0
 car on vérifie facilement que :

∫
1 xj + 1
Remarque --- ( x – x j ) ϕ j ( x ) dx = 0
Il est important de noter que, aussi bien dans le calcul des coeffi- h xj – 1
cients de la matrice que dans celui du second membre du système
Si la dérivée seconde est continue, on a donc :
linéaire, le calcul des intégrales s’effectue sur chaque intervalle
[ x j , x j + 1 ] , j = 0, 1, …, n , c’est-à-dire élément par élément.
 2f ( x ) 
∫
xj + 1 2
1 h
Remarque --- f ( x ) ϕ j ( x ) dx – f ( x j )  ------ Sup  d ----------------- 
h xj – 1 12 0 < x < 1  dx 2 
Comme prévu, la matrice A est tridiagonale ; de plus quel que soit
l’opérateur aux dérivées partielles considéré, on obtient une telle
structure de matrice, en particulier si l’on considère un problème or :
monodimensionnel à coefficient variable. Cette propriété provient 2 4 2
du fait que les fonctions de base choisies sont à supports disjoints. d u(x) d u(x) d f(x)
------------------ = f ( x ) ⇒ ------------------
4
= -----------------
On verra, par la suite, qu’il est possible de définir l’espace Eh en uti- dx 2 dx dx 2
lisant des fonctions de base polynomiales par morceaux de degré
plus élevé à supports disjoints sauf pour les fonctions de base atta- et
chées à des points contigus ; dans ce cas, cependant, la demi-lar-
geur de bande de la matrice augmente. Un autre choix de la base de h
2  d 4u ( x ) 
E j  ------ Sup  ------------------ 
6 0 < x < 1
dx 4 
l’espace Eh composé de fonctions à supports non disjoints conduit
à une matrice de discrétisation pleine.
Remarque Donc, le schéma numérique obtenu par discrétisation élément fini
Si l’on modifie la numérotation des fonctions de base (c’est-à-dire est donc une approximation de l’équation de Poisson à un terme h 2
des points du maillage), la matrice du système linéaire perd sa près. Cette estimation permet, par des calculs simples, compte tenu
structure tridiagonale. du fait que la matrice A est définie positive, de donner une majora-
tion de la norme discrète de l’erreur. Nous donnerons dans l’article
En général, il n’est pas toujours possible de calculer une primitive [AF 505], § 3.1 des résultats plus généraux de majoration d’erreur.

de l’intégrale
∫ f ( x ) ϕ j ( x ) dx ; dans la pratique, on utilise des formules Exemple 2. Soit le problème aux limites suivant :

d’intégration numérique. Si, par exemple, on utilise la formule des

 d 2u ( x )
trapèzes, on obtient : - + qu ( x ) = f ( x ), x ∈ ]0, 1[, ( q > 0 )
 – -----------------
 dx
2

∫ ∫
1 xj + 1
f ( x ) ϕ j ( x ) dx = f ( x ) ϕ j ( x ) dx 
0 xj – 1 u(0) = u(1) = 0

∫ ∫
xj xj + 1
= f ( x ) ϕ j ( x ) dx + f ( x ) ϕ j ( x ) dx = hf ( x j ) et la forme bilinéaire associée :
xj – 1 xj

∫
1
et, dans ce cas particulier, on retrouve le schéma classique de dis- a ( u, v ) =  du
------- ------ + qu v dx
dv
0  dx dx 
crétisation obtenu par différences finies. C’est un hasard que l’on ne
notera que dans peu de cas. Pour les fonctions chapeaux définies précédemment, on obtient, en
Remarque effectuant le même type de calcul, une matrice de discrétisation par
Dans le schéma défini précédemment, on obtient, en résolvant le élément fini symétrique, tridiagonale de coefficients :
système linéaire, une solution approchée ; il conviendra d’estimer
2 2qh 1 qh
l’erreur commise. La façon la plus classique est de calculer l’erreur a ( ϕ j , ϕ j ) = --- + ---------- , a ( ϕ j , ϕ j + 1 ) = a ( ϕ j , ϕ j – 1 ) = – --- + -------
de troncature définie par : h 3 h 6

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
© Techniques de l’Ingénieur, traité Sciences fondamentales AF 504 − 5

Dossier délivré pour

DOCUMENTATION
19/09/2008
 Dossier délivré pour
DOCUMENTATION
19/09/2008
MÉTHODE DES ÉLÉMENTS FINIS __________________________________________________________________________________________________________

Cet exemple est intéressant dans la mesure où l’on constate que Le calcul des coefficients de la matrice se déduit des calculs pré-
l’on ne retrouve pas la matrice de discrétisation obtenue par la cédents, et l’on trouve finalement :
méthode des différences finies. On verra dans l’article [AF 505], 1 1
§ 2.1 un moyen systématique de construction de la matrice de dis- a ( ϕ m ( x ), ϕ m + 1 ( x ) ) = – --- , a ( ϕ m + 1 ( x ), ϕ m + 1 ( x ) ) = ---
h h
crétisation par élément fini.
la dernière équation s’écrivant :
Exemple 3. Soit le problème aux limites suivant :

∫
um + 1 – um xm + 1
- =
---------------------------- f ( x ) ϕ m + 1 ( x ) dx
 d 2u ( x ) h xm
 – ------------------ = f ( x ) , x ∈ ]0, 1[
 dx 2 Remarque
 L’intégrale intervenant dans la dernière équation peut être évaluée
 du ( 1 )
 u ( 0 ) = 0 , --------------- = 0 en utilisant des formules de quadrature approchée ; par exemple,
 dx l’utilisation de la formule des trapèzes conduit à :

∫
dont la formulation variationnelle est : xm + 1 h h
f ( x ) ϕ m + 1 ( x ) dx ≈ ---f ( x m + 1 ) = ---f ( 1 )
xm 2 2
a ( u, v ) = L ( v ), ∀ v ∈ E
Remarque
avec : La matrice A du système linéaire est symétrique définie positive,
conformément au résultat du théorème 1. Un calcul direct montre
E = { v ∈ H 1 ( [ 0, 1 ] ) , v ( 0 ) = 0 } que :
m
et
∑ ( uj + 1 – uj )
t 2 2 m+1
U AU = u 1 + , ∀U ∈ 

∫ ∫
1 du dv 1 j=1
a ( u, v ) = ------- ------dx , L ( v ) = f v dx
0 dx dx 0 Exemple 4. Dans les exemples précédents, on a considéré une
approximation par des fonctions linéaires par morceaux sur chaque
Pour discrétiser le problème précédent par la méthode des éléments intervalle [ x j , x j + 1 ] ; on peut considérer également des fonctions qua-
finis, on considère le sous-espace E h de E défini par : dratiques par morceaux, et même plus généralement polynomiales de
degré inférieur ou égal à k par morceaux.
0
E h = { v h ∈ C ( [ 0, 1 ] ), v h linéaire sur chaque segment Considérons le cas de l’approximation quadratique ; on réalise une
[ x j , x j + 1 ] et v h ( 0 ) = 0 } approximation dans l’espace E h des fonctions continues dont la res-
triction à chaque intervalle [ x j , x j + 1 ] est un polynôme de degré infé-
L’espace E h est engendré par ( m + 1 ) fonctions de base
rieur ou égal à 2 ; pour construire une telle base de l’espace E h , on
j = m+1
{ ϕj }j = 1 où : introduit le point milieu noté x 1 de chaque segment [ x j , x j + 1 ].
j + ----
2
— pour 1  j  m , ϕ j ( x ) est une fonction chapeau considérée pré-
cédemment ; Dans le cas de problème avec conditions aux limites de Dirichlet, on
— pour j = m + 1 , ϕ j ( x ) est une fonction demi-chapeau. aura donc ( 2m + 1 ) degrés de liberté au lieu de m ; comme précé-
demment, à chaque point de discrétisation, on associe une fonction de
Remarque base ϕ j ( x ) polynomiale de degré 2 sur chaque intervalle qui vérifie
Dans les exemples précédents, la dimension de E h était m, alors ϕ j ( x i ) = δ i, j et on a dim ( E h ) = 2m + 1. La détermination des fonc-
qu’ici elle est m + 1 ; cela provient du fait que, dans le cas du pro-
blème avec conditions aux limites mêlées Dirichlet-Neumann, seule tions de base s’effectue par identification, ϕ j ( x ) ayant une expression
la dérivée au point x = 1 est connue, la valeur u ( 1 ) ne l’étant pas, analytique du type :
ce qui se traduit par un degré de liberté supplémentaire à détermi- 2
ϕ j ( x ) = ax + bx + c
ner.
Dans ces conditions, on peut développer la solution comme suit : On détermine les inconnues a, b et c par les équations données par
ϕ j ( x i ) = δ i, j . On a donc deux types de fonctions de base :
m+1 — celles attachées aux points x j et déterminées par les conditions :
uh ( x ) = ∑ vj ϕj ( x )
ϕ j ( x j ) = 1 , ϕ j  x 1  = ϕ j ( x j + 1 ) = 0
j=1 j + ----
2
et, comme précédemment, le système linéaire se déduit de : qui conduisent à l’expression suivante :
2 2

∫
1 2x 4x j + 3h 2x j + 3hx j
a ( u h ( x ), ϕ j ( x ) ) = f ( x ) ϕ j ( x ) dx , j = 1, … , m + 1 ϕ j ( x ) = ---------
2
- – ---------------------x
2
+ 1 + -----------------------------
2
-
0 h h h
— celles attachées aux points milieux x 1 et déterminées par les
Or, pour 1  j  m , les calculs sont identiques à ceux effectués conditions : j + ----
2
dans les exemples 1 et 2, et l’on retrouve les mêmes équations. Il ne
reste donc plus qu’à déterminer la dernière équation du système ϕ   = 1, ϕ = ϕ
1 x 1 1 ( xj ) 1 ( xj + 1 ) = 0
linéaire, qui se déduit de l’écriture suivante : j + ----  j + ---- j + ---- j + ----
2 2 2 2

a ( u h ( x ), ϕ m +1 ( x ) ) = a ( ϕ m ( x ), ϕ m +1 ( x ) )u m + a ( ϕ m +1 ( x ), ϕ m +1 ( x ) )u m +1 qui conduisent à l’expression suivante :

2 4 ( 2x + h ) 4x j ( x j + h )
∫
xm + 1 4x j
= f ( x ) ϕ m +1 ( x ) dx ϕ 1 (x) 2
- + --------------------------
= – -------- 2
- x – ---------------------------
2
j + ---- h h h
xm 2

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
AF 504 − 6 © Techniques de l’Ingénieur, traité Sciences fondamentales

Dossier délivré pour

DOCUMENTATION
19/09/2008
Dossier délivré pour
DOCUMENTATION
19/09/2008
__________________________________________________________________________________________________________ MÉTHODE DES ÉLÉMENTS FINIS

1,2
18
1 15 16 17
* *
* * *
* * 14
0,8 * * 13
12
* * 11
10
0,6 * * 9
* * 8
6 7
0,4 5
* * 4
2 3
0,2 1
* *

0* *

-0,2 Figure 4 – Exemple de maillage par des triangles

0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1
fonction de base en x = 0
* fonction de base en x = 0,5

Figure 3 – Fonctions de base attachées aux points x = 0 et x = 0,5 sur

l’élément K = [0,1] K1

Tous calculs effectués, la matrice de discrétisation associée est pen-

tadiagonale. Les fonctions ϕ j ( x ) , ϕ j + 1 ( x ) et ϕ 1 ( x ) sont représen- K2 K3
j + ----
tées sur la figure 3. 2

2. Introduction des éléments

finis en dimension 2 Figure 5 – Situation interdite

Dans la suite de ce paragraphe, nous allons construire l’espace

Nous allons généraliser la méthode décrite au paragraphe 1 ; sauf E h de fonctions vh , polynomiales en x et y sur chaque « triangle »
indication contraire spécifiée dans les remarques qui suivront, on de degré inférieur ou égal à k et continues sur le domaine
considérera dans ce paraghraphe le cas d’un problème aux limites Ω = Ω ∂ Ω
avec conditions aux limites de Dirichlet homogènes. On supposera
que le domaine Ω ⊂  2 est polyédrique ; on s’affranchira de cette
hypothèse restrictive au paragraphe 2.4 de l’article [AF 505]. On
effectue une « triangulation » du domaine Ω , c’est-à-dire que l’on 2.1 Interpolation par des polynômes
suppose que l’on peut découper Ω en M triangles K i , tels que les de degré 1
conditions suivantes soient vérifiées :

∅ L’espace E h est ici l’espace des fonctions v h linéaires en x et en y

M  sur chaque triangle K (c’est-à-dire de la forme ϕ ( x, y ) = ax + by + c )
Ω =  K i , avec K i K j =  un sommet, i ≠ j
et continues sur Ω . On a le résultat préliminaire suivant :
i = 1 
 un côté

Cette « triangulation » est représentée sur la figure 4. Lemme 1. Soit un triangle K dont les sommets A i ont les coor-
données ( x i , y i ) , 1  i  3 . Alors il existe un unique polynôme
Remarque p ( x, y ) ∈ P 1 de degré inférieur ou égal à un, prenant des valeurs
On suppose dans la suite que d’un sommet de chaque triangle données arbitraires p i , 1  i  3 , aux sommets A i .
partent toujours plusieurs côtés, la situation représentée sur la
figure 5 étant interdite.
Preuve ♦
Remarque
En effet, il suffit d’écrire p i = ax i + by i + c , 1  i  3 . C’est un
Si Ω n’est pas un polygone, on peut toujours le rendre polygonal système linéaire de trois équations à trois inconnues, dont la
en sous-tendant les arcs de la frontière ∂ Ω par des cordes. On verra matrice C a un déterminant donné par :
ultérieurement que l’on peut « trianguler » le domaine Ω en consi-
dérant les quadrilatères quelconques (carré, rectangle, parallélo- det ( C ) = x 1 ( y 2 – y 3 ) + x 2 ( y 3 – y 1 ) + x 3 ( y 1 – y 2 )
gramme, losange…), ainsi que des triangles et quadrilatères à côté
courbe, ce qui permet de s’affranchir de l’hypothèse concernant sa det ( C ) est différent de zéro dès que les points A i ne sont pas ali-
forme polyédrique et de prendre en compte la courbure de la fron- gnés, ce qui a été implicitement supposé. On peut donc déterminer
tière. de façon unique les réels a, b et c sur chaque triangle.

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
© Techniques de l’Ingénieur, traité Sciences fondamentales AF 504 − 7

Dossier délivré pour

DOCUMENTATION
19/09/2008
 Dossier délivré pour
DOCUMENTATION
19/09/2008
MÉTHODE DES ÉLÉMENTS FINIS __________________________________________________________________________________________________________

Le polynôme p ( x , y ) peut s’exprimer de façon simple en fonction Preuve ♦

des coordonnées barycentriques λ i ( M ) = λ i ( x , y ) , 1  i  3 , d’un t
En effet, on a établi que C Λ = Z ; or C étant inversible, on a
point quelconque M intérieur au triangle K. ♦ –t
Λ = C Z et, par conséquent, les composantes de Λ sont bien des
polynômes de degré un. De plus, si on particularise le point M, en le
choisissant confondu avec le sommet A 1 , on a :
Définition 1. Soit M = ( x , y ) un point quelconque intérieur
au triangle K de sommets non alignés A i , 1  i  3 ; soit 3 3 3
( x i , y i ) , 1  i  3 , les coordonnées des sommets du triangle. On
appelle coordonnée barycentrique du point M par rapport aux
∑ λi ( A1 ) = 1, ∑ xi λi ( A1 ) = x1 , ∑ yi λi ( A1 ) = y1
i=1 i=1 i=1
trois points A i , les scalaires λ i ( M ) = λ i ( x , y ) , 1  i  3 , véri-
fiant : On constate que λ 1 ( A 1 ) = 1 et λ 2 ( A 1 ) = λ 3 ( A 1 ) = 0 est solu-
3 3 3 tion de ce système ; comme, d’après le résultat du lemme 1, ce der-
∑ λi ( x , y ) = 1, ∑ xi λi ( x , y ) = x, ∑ yi λi ( x , y ) = y nier admet une solution unique, on a par permutation circulaire
λ i ( A j ) = δ i, j , 1  i , j  3 . Cette dernière propriété nous permet de
i=1 i=1 i=1
vérifier trivialement que les nombres λ i ( M ) , 1  i  3 , sont linéai-
rement indépendants ; comme ce sont également des polynômes
Remarque de degré inférieur ou égal à un, l’ensemble λ i ( M ) , 1  i  3,
t constitue une base de P 1 , ensemble des polynômes de degré infé-
Posons Λ =  λ 1 ( M ), λ 2 ( M ), λ 3 ( M ) et Z = ( 1, x , y ) . Les coor-
t
rieur ou égal à 1. On peut donc écrire :
 
données barycentriques sont alors solutions du système linéaire 3

C Λ = Z et d’après le résultat du lemme 1, la matrice C est inver-

t t p(M) = ∑ θ i λ i(M)
i=1
sible, ce qui prouve l’existence et l’unicité des coordonnées bary-
centriques. Particularisons le point M au sommet A j ; on a donc :
Remarque 3

Notons S 1 l’aire du triangle A 2 MA 3 , S 2 celle du triangle A 1 MA 3 p ( Aj ) = pj = ∑ θ i λ i ( Aj ) = θ j , car λ i ( A j ) = δ i, j

i=1
et S 3 la surface du triangle A 1 MA 2 . Soit S l’aire du triangle K. On a
3
3
évidemment S = ∑ S i ; de plus, on rappelle que : Donc on peut écrire p ( x , y ) = ∑ pi λ i ( x , y ) . Enfin, soit M un
i=1
i=1
point appartenant au côté A j A k opposé au sommet A i ; on peut
 1 1 1 
  toujours écrire M = θ A k + ( 1 – θ )A j , de sorte que les coordonnées x
1 1 1
S = --- A 1 A 2 ∧ A 1 A 3 = --- det  x 1 x 2 x 3  = --- det ( C ) et y du point M sont également combinaison convexe de celles des
2 2   2
 y1 y2 y3  points A j et A k ; λ i ( M ) étant un polynôme de degré un en x et
en y, on a :
S =  ---  x 1 ( y 2 – y 3 ) + x 2 ( y 3 – y 1 ) + x 3 ( y 1 – y 2 )
1 λi ( M ) = α x + β y + χ
 2  
soit, puisque χ = θχ + ( 1 – θ ) χ :
On peut, de la même façon, déterminer S i pour 1  i  3 . On
λ i ( M ) = α ( θ x k + ( 1 – θ )x j ) + β ( θ y k + ( 1 – θ )y j ) + χ
S
vérifie alors que λ i ( M ) = -----i , 1  i  3 . Dans le cas particulier où K = θ ( α xk + β yk + χ ) + ( 1 – θ ) ( α xj + β yj + χ )
S
est le triangle rectangle isocèle de côté l’unité, noté dans la suite Kˆ , ou encore :
λi ( M ) M ∈ Aj Ak = θλ i ( A k ) + ( 1 – θ ) λ i ( A j ) = 0
les points A i , 1  i  3 , ayant pour coordonnées A 1 = ( 1, 0 ),
A 2 = ( 0, 1 ) et A 3 = ( 0, 0 ) , l’expression des coordonnées barycen- Donc :

triques associée au point M = ( x, y ) s’écrit : p(M) M ∈ Aj Ak = λi ( M ) pi + λj ( M ) pj + λk ( M ) pk = λj ( M ) pj + λk ( M ) pk

♦
λ 1 ( M ) = x , λ 2 ( M ) = y et λ 3 ( M ) = 1 – x – y
À partir de ces résultats techniques préliminaires, on peut généra-
À partir de cette notion on établit le résultat suivant : liser à la dimension 2 les résultats obtenus en dimension 1. Soit
donc :
Fh = { vh , vh K ∈ P 1, v h entièrement déterminée par ses valeurs
Lemme 2. Les coordonnées barycentriques du point M,
λ i ( x , y ) , 1  i  3 , sont des polynômes de degré 1 en x et en y, aux nœuds de la triangulation },
vérifiant de plus λ i ( A j ) = δ i, j , 1  i, j  3 . Le polynôme p ( x, y )
de degré un prenant les valeurs données arbitraires p i , l’espace des fonctions dont la restriction à chaque triangle K est un
1  i  3 , aux sommets A i s’exprime au moyen des coordon- polynôme de degré inférieur ou égal à un, les fonctions de l’espace
nées barycentriques de la façon suivante : F h étant entièrement déterminées par leurs valeurs aux sommets
A j , 1  j  3 des triangles K. On alors le résultat important :
3
p ( x, y ) = ∑ pi λi ( x , y ) 2
i=1 Théorème 3. On a les inclusions suivantes F h ⊂ L ( Ω ) ,
et la restriction de p ( x, y ) pour un point M appartenant au côté Fh ⊂ H 1 ( Ω ) , Fh ⊂ C ( Ω )
A j A k est un polynôme de degré 1 qui ne dépend que des 1
valeurs p j et p k . ⇒ Fh ⊂ C ( Ω ) ∩ H ( Ω )

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
AF 504 − 8 © Techniques de l’Ingénieur, traité Sciences fondamentales

Dossier délivré pour

DOCUMENTATION
19/09/2008
Dossier délivré pour
DOCUMENTATION
19/09/2008
__________________________________________________________________________________________________________ MÉTHODE DES ÉLÉMENTS FINIS

Preuve ♦
2
1) F h est inclus dans L ( Ω ) , car les fonctions de cet espace sont
des polynômes de degré inférieur ou égal à un et définies par mor-
ceaux sur chaque triangle K ; donc, sur chaque triangle, ces fonc-
tions sont de carré sommable.
ϕ 0( x , y )
2) De même, les dérivées des fonctions de l’espace F h sont des
polynômes de degré inférieur ou égal à un et définies par morceaux
sur chaque triangle K ; donc les dérivées des fonctions de F h sont
de carré sommable.
3) Soit K et K’ deux triangles adjacents le long d’un côté A 2 A 3 ,
(K) ( K′ )
tels que K, K′ ⊂ Ω . Soit v h ∈ F h et notons v h , v h les restrictions
(K)
de v h aux triangles K et K’. La restriction vh sur l’arête commune
ne dépend que des valeurs prises aux sommets de cette arête, de
( K′ )
même que la restriction v h sur cette même arête, conformément A0
au résultat du lemme 2 ; donc pour tout point M appartenant à cette
(K) ( K′ )
arête commune, on a l’égalité v h ( M ) = v h ( M ) ; donc les fonc-
tions de l’espace F h sont continues sur K ∪ K′ . ♦
Remarque
Lorsque l’on passe d’un triangle K à un triangle adjacent K’, les
fonctions v h de F h sont continues sur K ∪ K′ ⊂ Ω . Figure 6 – Fonction de base linéaire par morceaux associée
au point A0
Soit à présent E h le sous-espace de F h des fonctions v h nulles
aux sommets situés sur la frontière ∂ Ω ; ces fonctions sont en fait
1
nulles sur toute la frontière et on a l’inclusion E h ⊂ C ( Ω ) ∩ H ( Ω ) .
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
La dimension de E h est égale au nombre N de sommets non situés
sur la frontière ∂ Ω . Soit A 1, …, A N , ces sommets, c’est-à-dire les X X X X
nœuds de la « triangulation ». On considère les fonctions ϕ j ( x, y ) X X X X X
de E h , représentées sur la figure 6, telles que :
X X X X X
ϕ j ( A i ) = δ i, j , 1  i , j  N X X X X
Le support de la fonction ϕ j est formé de la réunion des triangles X X X X
admettant Aj comme sommet. L’ensemble des fonctions { ϕ j } est
une base de E h . À chaque nœud Aj on associe donc une fonction X X X X X X X
de base ϕ j , définie comme en dimension 1 par : X X X X X X X
1) ϕ j ( A i ) = δ i, j , 1  i , j  N ; X X X X X X X
2) la restriction de ϕ j à chacun des triangles possédant Aj
comme sommet est un polynôme de degré inférieur ou égal à un ; X X X X X
3) pour chaque triangle K i , n’appartenant pas à support ( ϕ j ) , on X X X X X
a ϕ j ( x, y ) = 0 .
X X X X X X X
Soit u h la solution approchée de l’équation de Poisson en dimen-
sion 2 avec conditions aux limites de Dirichlet homogènes ; posons X X X X X X X
u i = u h ( A i ) . Si l’on développe u h dans la base ϕ j , on aura comme X X X X X X X
précédemment :
X X X X
N
0
uh ( M ) = ∑ uj ϕj ( M ) , ∀ M ∈ Ω X X X X
j=1 X X X X X
où les composantes u i , valeur de u h aux sommets des triangles, X X X X X
sont données par la résolution du système linéaire :
X X X X
N
Les croix représentent des éléments non nuls
∑ a ( ϕi, ϕj )ui = <f, ϕ j > , 1  j  N
j=1 Figure 7 – Matrice de discrétisation
Lorsque l’on calculera les éléments de matrice a ( ϕ i, ϕ j ) , ces ter-
mes seront nuls lorsque support ( ϕ j ) ∩ support ( ϕ i ) = ∅ . Donc cette gonaux sont également tridiagonaux, conformément à la représen-
tation de la figure 7.
matrice A sera creuse. La matrice A s’appelle matrice de rigidité du
système. Les points A i auxquels sont attachées les valeurs u i de la Remarque
solution u h s’appellent les nœuds de la « triangulation » ; les côtés On cherchera à minimiser la largeur de la bande de la matrice en
adoptant une numérotation des points du maillage adéquate. De
des triangles sont les arêtes de la « triangulation ». manière générale, la largeur de bande est caractérisée par la diffé-
Si l’on considère le maillage représenté sur la figure 4, par des rence maximale entre l’indice de la dernière et de la première fonc-
considérations sur les supports des fonctions de base, on peut don- tion de base dont l’intersection des supports est différente du vide ;
ner la forme de la matrice de discrétisation (figure 7), où, dans cette il existe des algorithmes automatiques de numérotation qui permet-
représentation les croix représentent des éléments non nuls. On tent une telle minimisation, tel l’algorithme de Cuthill-Mac Kee
remarque que cette matrice est tridiagonale par blocs, les blocs dia- (cf. [2]).

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
© Techniques de l’Ingénieur, traité Sciences fondamentales AF 504 − 9

Dossier délivré pour

DOCUMENTATION
19/09/2008
 Dossier délivré pour
DOCUMENTATION
19/09/2008
MÉTHODE DES ÉLÉMENTS FINIS __________________________________________________________________________________________________________

La détermination de ϕ 0 ( x , y ) sur les autres triangles ne présente

pas de difficulté majeure et est laissée à titre d’exercice. Par ailleurs,
comme en dimension un, on devra déterminer, par le même principe,
ϕ0 (x,y) l’expression des fonctions de base ϕ i ( x , y ) , i = 1, …, 6 , attachées
aux points A i , i = 1, …, 6 , voisins du point A 0 . Par exemple, sur le
triangle K 1 , l’expression de la fonction de base ϕ 1 ( x , y ) et celle de
son gradient sont données par :
x–x
=  ---, 0 
1
ϕ1 ( x , y ) = --------------0- et ∇ ϕ 1 ( x , y )
K1
h K1 h 

Il convient de noter que seules les expressions de ces fonctions

ϕ i ( x , y ) , sur les triangles utiles, est nécessaire ; en effet, dans la suite
du calcul, on ne doit prendre en compte que les triangles pour lesquels
A0 support ( ϕ j ) ∩ support ( ϕ 0 ) = ∅ pour j ≠ 0 . Le reste du calcul con-
∂Ω siste à calculer les coefficients de la matrice A, c’est-à-dire un certain
nombre d’intégrables doubles ; par exemple, on aura :
6

Figure 8 – Fonction de base demi-chapeau linéaire par morceaux

a ( ϕ 0, ϕ 0 ) =
∫∫  ∇ ϕ 2 + q ϕ 2 dxdy =
Ω
0 0 ∑
i=1
∫∫ Ki 
 ∇ ϕ 2 + q ϕ 2 dxdy
0 0

associée au point frontière A0

soit, dans la situation considérée ici et puisque les gradients des fonc-
2
Remarque h
tions de base sont constants sur chaque triangle et que aire ( K i ) = ------ ,
Considérons le cas d’un problème aux limites avec condition de 2
type Neumann. La dimension de l’espace d’approximation E h est pour tous les indices i, il vient :
identique à celle de l’espace F h ; elle est égale au nombre de nœuds 6 2
∑  ∇ ϕ0 K ⋅ aire ( Ki ) + q ∫∫
2 2 h
0 a ( ϕ 0, ϕ 0 ) = ϕ dxdy = 4 + q ------
appartenant à Ω = Ω ∪ ∂ Ω (et non Ω intérieur de Ω ), car la valeur Ki 0 2
i=1
de la fonction sur le bord n’est pas à priori connue. On aura donc
deux sortes de fonctions de base : 0 Le calcul des termes hors diagonaux s’effectue de la même façon ;
— les fonctions chapeaux pour les points A i ∈ Ω , représentées par exemple :
sur la figure 6 ;

∫∫
— les fonctions demi-chapeaux pour les points A i ∈ ∂ Ω , ces der-  < ∇ ϕ , ∇ ϕ > + q ϕ ϕ  dxdy
a ( ϕ 0, ϕ 1 ) =
Ω
nières étant représentées sur la figure 8. 0 1 0 1

6
2.2 Exemples de matrices
de discrétisation
=
i=1
∑ ∫∫  < ∇ ϕ , ∇ ϕ > + q ϕ ϕ  dxdy
Ki  0 1 0 1

On examine, dans ce paragraphe, un certain nombre de situations Or la fonction ϕ 1 ( x , y ) n’est différente de zéro que sur les triangles
classiques. K 1 et K 6 et ce terme se réduit à :
2
Exemple 5. On considère le domaine Ω = [ 0, 1 ] × [ 0, 1 ] ⊂  , et
le problème suivant :
0
a ( ϕ 0, ϕ 1 ) =
∫∫  < ∇ ϕ , ∇ ϕ > + q ϕ ϕ  dxdy
K1  0 1 0 1

 – ∆u ( x ) + q u ( x ) = f ( x ), ( q  0 ) sur Ω 2

∫∫  < ∇ ϕ , ∇ ϕ > + q ϕ ϕ  dxdy = – 1 + q -----

 h
+ -
 u ( x ) = 0, sur ∂ Ω K6 0 1 0 1 12
On considère également une « triangulation » du domaine Ω par Finalement, tous calculs effectués on trouve :
des triangles rectangles isocèles de coté h, comme représenté sur la
figure 9 ; on constate que chaque nœud est le centre d’un polygone à 2
h
six côtés rassemblant six triangles K i , i = 1, …, 6 , possédant donc ce a ( ϕ 0, ϕ 1 ) = a ( ϕ 0, ϕ 2 ) = a ( ϕ 0, ϕ 4 ) = a ( ϕ 0, ϕ 5 ) = – 1 + q ------
12
nœud comme sommet. Soit A 0 un point de coordonnées ( x 0, y 0 ) et
A i , i = 1, …, 6 , les voisins de ce point A 0 , les points A 1 , A 2 , A 4 et et
A 5 étant de plus en position de points cardinaux autour du point A 0 .
Soit ϕ 0 ( x, y ) la fonction de base attachée au point A 0 ; le support de la h
2

fonction de base ϕ 0 ( x, y ), représenté sur la figure 10, est constitué a ( ϕ 0, ϕ 3 ) = a ( ϕ 0, ϕ 6 ) = q ------

12
par l’union des triangles K i , i = 1, …, 6 . De plus, dans le cas de
l’approximation par des fonctions de base linéaires par morceaux, on L’évaluation du second membre, s’effectue suivant le même prin-
peut exprimer les restrictions de la fonction de base K i , i = 1, …, 6 , cipe dans la mesure où il suffit de calculer l’intégrale suivante :
sur chaque triangle K i , grâce à la condition ϕ j ( A i ) = δ i, j ; en effet, on
cherche ϕ 0 ( x , y ) de la forme : 6

ϕ0 ( x , y ) K = ai x + bi y + ci
i ∫∫ Ω
f ( x, y ) ϕ 0 ( x , y )dxdy = ∑
i=1
∫∫ Ki
f ( x, y ) ϕ 0 ( x , y )dxdy

Par identification sur chaque triangle K i , on peut déterminer les

coefficients a i , b i , c i ; par exemple, si nous considérons, le triangle En général, il n’est pas toujours possible d’avoir une primitive de la
K 1 de sommets A 0 , A 1 = ( x 0 + h, y 0 ) et A 2 = ( x 0 , y 0 – h ) , on
obtient :
x0 – x y – y0
fonction sous le signe
∫ , ici f ( x, y ) ϕ 0 ( x, y ) ; on utilise donc, pour éva-

- + -------------- + 1 et ∇ ϕ 0 ( x , y ) K =  – ---, ----

1 1 luer cette intégrale, des formules d’intégration approchée. Nous consi-
ϕ 0 ( x , y ) K = ---------------  h h
1 h h 1 dérons ci-dessous deux types d’approximation.

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
AF 504 − 10 © Techniques de l’Ingénieur, traité Sciences fondamentales

Dossier délivré pour

DOCUMENTATION
19/09/2008
Dossier délivré pour
DOCUMENTATION
19/09/2008
__________________________________________________________________________________________________________ MÉTHODE DES ÉLÉMENTS FINIS

Remarque. On constate que, dans cet exemple, pour q ≠ 0 , on

ne trouve pas un schéma de discrétisation identique à celui
obtenu par la méthode des différences finies. Par contre, lorsque
le nombre q est nul, on retrouve le schéma de discrétisation à
cinq points, identique à celui obtenu par la méthode des différen-
ces finies classique ; ce résultat n’est que le fruit du hasard et per-
met de démarquer la méthode des éléments finis par rapport à
celle des différences finies.

Approximation n° 2 : pour calculer

∫∫ Ω
f ( x, y ) ϕ 0 ( x, y ) dx dy , on

i=6
interpole la fonction f ( x, y ) dans la base { ϕ i } i = 0 ,
6
pour ( x, y ) ∈  Ki ; on peut donc écrire :
i=1
6
f ( x, y ) = ∑ f ( A i ) ϕ i ( x, y )
i=0
Figure 9 – Maillage élément fini uniforme
donc
6

A5 A6
∫∫Ω
f ( x, y ) ϕ 0 ( x, y ) dx dy = ∑ f ( Ai ) ⋅
i=0
∫∫ ϕ
Ω
i ( x, y ) ϕ 0 ( x , y ) d x dy

Tous calculs faits, on trouve :

∫∫ ∫∫
2 2
2 h h
K5 ϕ 0 ( x, y ) dx dy = ------ et ϕ i ( x, y ) ϕ 0 ( x, y ) dx dy = ------ , i = 1, 2, …, 6
Ω 2 Ω 12

K4 ce qui permet d’avoir l’estimation :

6
 f ( A0 ) f ( A i ) 2
K6
∫∫ Ω  2
- + ∑ ------------
f ( x, y ) ϕ 0 ( x, y ) dx dy =  -------------
i=1
12 
- h

A0 2
Exemple 6. On considère le domaine Ω = [ 0, 1 ] × [ 0, 1 ] ⊂  , et
A4 A1 le problème suivant :

K3  – ∆u ( x ) + q u ( x ) = f ( x ), ( q > 0 ) sur Ω
0

K1  ∂u ( x )
 -------------- = 0, sur ∂ Ω
 ∂n
K2
Discrétisons ce problème en considérant de nouveau une triangula-
tion du domaine par des triangles rectangles isocèles de côté h. Pour
les points internes à Ω , il n’y a aucun changement par rapport à la
A3 A2 situation envisagée dans l’exemple précédent. Les seuls changements
qui interviennent sont relatifs aux points appartenant à la frontière ;
Figure 10 – Support de la fonction ϕ0 associée au point A0 dans ce cas, il faut attacher à ces points frontières des fonctions demi-
chapeaux, de telle sorte que le support de la fonction de base reste
inclus dans Ω ; par conséquent, certains coefficients de la matrice de
rigidité vont diminuer. On suppose, par exemple, que les points A 0 ,
Approximation n° 1 : pour calculer chacune des six intégrales, on A 1 et A 4 sont situés sur ∂ Ω et que, pour un tel point, le support de
utilise une formule de quadrature approchée qui donne un résultat la fonction demi-chapeau est constitué par le polygone A 4 A 5 A 6
exact pour les polynômes de degré inférieur ou égal à un : A 1 A 0 A 4 . Dans ce cas, les contributions de a ( ϕ 0, ϕ 0 ) , a ( ϕ 0, ϕ 1 ) et
a ( ϕ 0, ϕ 4 ) sont diminuées de moitié, les autres a ( ϕ 0, ϕ 5 ) et a ( ϕ 0, ϕ 6 )
3 restant inchangées. Plus exactement, on obtient en effectuant les
∫∫
aire ( K )
K
g ( x, y ) dx dy ≈ ---------------------
3 ∑ g ( Ai ) calculs :
i=1 2 2
h 1 h
a ( ϕ 0, ϕ 0 ) = 2 + q ------ , a ( ϕ 0, ϕ 1 ) = a ( ϕ 0, ϕ 4 ) = – --- + q ------ ,
4 2 24
Par un calcul très simple on obtient :
2 2
h h
a ( ϕ 0, ϕ 5 ) = – 1 + q ------ et a ( ϕ 0, ϕ 6 ) = q ------
∫∫
2
h 12 12
f ( x, y ) ϕ 0 ( x, y ) dx dy ≈ ------ f ( A 0 )
K 6
de plus, pour ce point, la part contributive du second membre diminue
de moitié, et on obtient dans le cadre du premier type d’approximation
donc : de l’intégrale de surface :
6

∫∫
2 2

∫∫
h
∑ -----6- f ( A0 )
2
f ( x, y ) ϕ 0 ( x, y ) dx dy = = h f ( A0 ) h
f ( x, y ) ϕ 0 ( x , y )dxdy = ------ f ( A 0 )
Ω Ω 2
i=1

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
© Techniques de l’Ingénieur, traité Sciences fondamentales AF 504 − 11

Dossier délivré pour

DOCUMENTATION
19/09/2008
 Dossier délivré pour
DOCUMENTATION
19/09/2008
MÉTHODE DES ÉLÉMENTS FINIS __________________________________________________________________________________________________________

Remarque ∂u
nant de la formule de Green. Donc si ------ = g , on aura à tenir compte
On verra dans l’article [AF 505] un procédé de calcul systématique ∂n
qui permet d’intégrer automatiquement la part contributive appor- de l’expression suivante de la forme bilinéaire :
tée par chaque triangle au calcul des coefficients de la matrice et du
second membre. En fait, il suffit de balayer élément par élément et
non nœud par nœud comme on l’a fait jusqu’à présent.
Remarque
a ( u, v ) =
∫∫
Ω
( ∇u ⋅ ∇v + quv ) dx dy =
∫∫Ω
fv dx dy +
°∫
∂Ω
gv ds

Dans le cas de conditions aux limites de type Neumann non

homogènes, il faut rajouter la contribution du terme de bord prove- et rajouter la part de l’intégrale curviligne.

Références bibliographiques

[1] RAVIART (P.A.) et THOMAS (J.M.). – Introduc- [2] AXELSON (O.) et BARKER (V.A.). – Finite ele-
tion à l’analyse numérique des équations aux ment solution of boundary value problems.
dérivées partielles. Collection Mathématiques Theory and computation. Academic Press
Appliquées, Masson (1983). (1984).

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
AF 504 − 12 © Techniques de l’Ingénieur, traité Sciences fondamentales

Dossier délivré pour

DOCUMENTATION
19/09/2008