UNIVERSITÉ MOHAMMED PREMIER

Ecole Nationale des Sciences Appliquées

Oujda

Méthodes des Éléments Finis

Filière Génie Civil

Deuxième année

2023-2024

Imad EL MAHI

Adresse : Complexe Universitaire, Hay El Qods, B.P 669, 60000 Oujda, Maroc
Table des matières

1 Quelques notions sur les espaces fonctionnels 4

1.1 Espaces vectoriels normés . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.1.1 Produits scalaires . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.1.2 Normes sur E . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2 Espaces de Banach et de Hilbert . . . . . . . . . . . . . . . . . . . . . . . . 9
1.2.1 Espace complet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.2.2 Espace de Banach . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.2.3 Espace de Hilbert . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3 Espaces fonctionnels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3.1 Espaces de fonctions intégrables . . . . . . . . . . . . . . . . . . . . 11
1.3.2 Le cas particulier de l'espace L2 (Ω) . . . . . . . . . . . . . . . . . . 11
1.4 Fonctions tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.5 Espaces de Sobolev . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.6 Inégalité de Poincaré . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.7 Trace d'une fonction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.8 Formules de Green . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2 Méthode des Éléments nis 18

2.1 Formulation variationnelle . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.1.1 Exemple en 1D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.1.2 Exemple en 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2 Existence et unicité de la solution du problème variationnel . . . . . . . . . 20
2.3 Principe de la méthode des éléments nis . . . . . . . . . . . . . . . . . . . 21
2.4 Mise en ÷uvre de la méthode en dimension 1 . . . . . . . . . . . . . . . . . 23
2.4.1 Maillage et fonctions de base . . . . . . . . . . . . . . . . . . . . . . 24
2.4.2 Méthode des éléments nis P1 . . . . . . . . . . . . . . . . . . . . . 24
2.5 Résolution du système linéaire : . . . . . . . . . . . . . . . . . . . . . . . . 34
2.5.1 Méthodes directes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.5.2 Méthodes itératives . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.6 Problème 1D avec conditions aux limites de Neumann . . . . . . . . . . . . 38
2.6.1 Formulation variationnelle . . . . . . . . . . . . . . . . . . . . . . . 39
2.6.2 Formulation variationnelle discrète . . . . . . . . . . . . . . . . . . 39
2.6.3 Fonctions de base de Vh . . . . . . . . . . . . . . . . . . . . . . . . 40
2.7 Approximation éléments nis P2 . . . . . . . . . . . . . . . . . . . . . . . . 45
2.8 Méthode des Eléments Finis en dimension 2 . . . . . . . . . . . . . . . . . 50
2.8.1 Formulation variationnelle . . . . . . . . . . . . . . . . . . . . . . . 50

2
2.8.2 Formulation variationnelle discrète . . . . . . . . . . . . . . . . . . 51
2.8.3 Espace d'approximation Vh (Eléments Finis P1 ) . . . . . . . . . . . 52

3
Chapitre 1
Quelques notions sur les espaces
fonctionnels

Introduction
Un grand nombre de problèmes issus de la physique et mécanique peuvent être décrits
par des systèmes d'équations aux dérivées partielles (EDP) auquelles sont associées des
conditions initiales et des conditions limites. On s'intéresse généralement à la convergence
de ces systèmes en réponse à ces conditions initiales et aux limites. Quelle est la réponse
du système à ces conditions ? c'est-à-dire comment se comporte la solution ? Il faut donc
s'assurer que la solution du système reste dans l'espace dans lequel on travaille si les
conditions avec lesquelles on a lancé le système appartiennent à cet espace.
La question qui se pose est donc la suivante : si une famille de solutions appartient à un
espace, la solution limite appartient elle aussi à cet espace ?
Pour garantir ceci il faut que l'espace dans lequel on travaille ne possède pas de trou,
c'est-à-dire qu'il n'y a aucun élément manquant. On parle d'espace complet. D'autre
part, la convergence vers la solution physique du problème sera atteinte lorsque la norme
de l'erreur soit assez petite. L'espace du travail doit donc être muni d'une norme.
Un espace vectoriel normé complet est dit Espace de Banach.
Un espace vectoriel muni d'un produit scalaire et complet pour la norme issue de ce
produit scalaire est dit Espace de Hilbert.
Evidemment, un espace de Hilbert est un espace de Banach.

1.1 Espaces vectoriels normés

On considère un espace vectoriel (E, +, .) sur un corps K .
On prendra ici E l'espace des fonctions et K le corps R ou C.

1.1.1 Produits scalaires

Le produit scalaire est une loi qui se rajoute aux autres opérations dans l'espace vectoriel
E . Il permet de dénir l'orthogonalité dans l'espace, la notion d'angle de deux vecteurs,
et de distance.

4
Dénition 1.1.1 Soit E un espace vectoriel sur R.
On dit qu'une application ϕ : E × E −→ R dénit un produit scalaire sur E si elle
vérie
(x, y) 7−→ ϕ(x, y)
les propriétés suivantes :
i/ ϕ est bilinéaire sur E , c'est-à-dire :

3 ϕ(λx + µy, z) = λϕ(x, z) + µϕ(y, z)
∀ (x, y, z) ∈ E :
ϕ(x, λy + µz) = λϕ(x, y) + µϕ(x, z)

ii/ ϕ est symétrique, c'est-à-dire :

∀ (x, y) ∈ E 2 ; ϕ(x, y) = ϕ(y, x)

iii/ ϕ est dénie positive, c'est-à-dire :

∀ x ∈ E ; ϕ(x, x) ≥ 0 et ϕ(x, x) = 0 ⇒ x = 0

Dénition 1.1.2 (Produit scalaire hermitien)

Soit E un espace vectoriel sur C.
Une application ϕ : E × E −→ C est dite produit scalaire hermitien sur E si elle
vérie :
(x, y) 7−→ ϕ(x, y)

i/ ϕ est sesquilinéaire, c'est-à-dire :


3 ϕ(λx + µy, z) = λϕ(x, z) + µϕ(y, z)
∀ (x, y, z) ∈ E :
ϕ(x, λy + µz) = λϕ(x, y) + µϕ(x, z)

ii/ ϕ est symétrique hermitienne, c'est-à-dire :

∀ (x, y) ∈ E 2 ; ϕ(y, x) = ϕ(x, y)

iii/ ϕ est dénie positive, c'est-à-dire :

∀ x ∈ E ; ϕ(x, x) ≥ 0 et ϕ(x, x) = 0 ⇒ x = 0

Notation: On utilise souvent l'écriture ϕ(x, y) = < x, y >.

Exemples :
1. Dans Rn , on dénit le produit scalaire euclidien usuel par :
x = (x1 , . . . , xn ) ; y = (y1 , . . . , yn )
n
X
< x, y > = xi yi = x1 y1 + x2 y2 + . . . + xn yn
i=1

5
 2. Dans Cn , le produit scalaire hermitien usuel est déni par :
n
X
< x, y >= xi yi = x1 y1 + x2 y2 + . . . + xn yn
i=1

3. Soit E = C 0 ([a, b], R) l'ensemble des fonctions continues de [a, b] vers R.

L'application < ·, · > : E × E −→ R
Z b
(f, g) 7−→ f (x) g(x) dx
a
dénit un produit scalaire sur E .
Dénition 1.1.3 (Espace préhilbertien)
Un espace vectoriel E muni d'un produit scalaire euclidien ou hermitien est appelé espace préhilbertien.

1.1.2 Normes sur E

Dénition 1.1.4 Soit E un espace vectoriel sur K .
Une norme sur E est une application k.k : E −→ R vériant :
x 7−→ kxk
i/ ∀ x ∈ E ; kxk ≥ 0 et kxk = 0 ⇔ x = 0
ii/ ∀ λ ∈ K , ∀ x ∈ E ; kλxk = |λ| kxk
iii/ ∀ (x, y) ∈ E 2 ; kx + yk ≤ kxk + kyk (Inégalité triangulaire)
Un espace vectoriel E muni d'une norme est dit espace vectoriel normé (evn).

Exemples :
ˆ Si E = R, la valeur absolue est une norme.
ˆ Si E = C ; le module dénit une norme.
ˆ Si E = Rn ; on peut le faire munir de l'une des normes suivantes, et qui sont souvent
utilisées :
Pour x ∈ Rn ; x = (x1 , . . . , xn )
n
X
kxk1 = |xi |
i=1
n
!1/2
X
kxk2 = x2i dite norme euclidienne
i=1

En général :
n
!1/p
X
kxkp = |xi |p
i=1
kxk∞ = sup |xi |
1≤i≤n

6
Norme euclidienne
À partir d'un produit scalaire < . , . >, on peut dénir une norme dite norme euclidienne
par : √
kxk = < x, x >

Exercice 1
Montrer que k.k : E −→ R est bien une norme sur E.
√
x 7−→ kxk = < x, x >

Normes équivalentes
On dira que deux normes k.k1 et k.k2 sont équivalentes s'il existe

α, β > 0 tels que ∀x ∈ E ; αkxk1 ≤ kxk2 ≤ βkxk1

Exercice 2
Montrer que les normes k.k1 , k.k2 et k.k∞ dénies sur Rn sont équivalentes.

Remarque :
Dans un espace de dimension nie toutes les normes sont équivalentes.
(Attention ceci n'est pas le cas dans un espace de dimension innie)

Propriétés de la norme
P1 Soit k.k une norme sur E . On a :
1. ∀ (x, y) ∈ E 2 ; |kxk − kyk| ≤ kx − yk

P2 Soit k.k une norme euclidienne associée à un produit scalaire < . , . >.
On a ∀ (x, y) ∈ E 2 ,
1. | < x, y > | ≤ kxk kyk Inégalité de Cauchy-Shwarz.

2. kx + yk ≤ kxk + kyk Inégalité de Minkowski.

3. kx + yk2 + kx − yk2 = 2 (kxk2 + kyk2 ) Identité du parallélogramme.

7
Preuve
P1 Il sut d'écrire les deux égalités :

kxk = k(x − y) + yk ≤ kx − yk + kyk ⇒ kxk − kyk ≤ kx − yk 
=⇒ |kxk − kyk| ≤ kx−yk
kyk = k(y − x) + xk ≤ ky − xk + kxk ⇒ kyk − kxk ≤ kx − yk


P2 1. On commence par montrer l'inégalité dans le cas où le nombre < x, y > est
réel. Il sut d'écrire que :
< x + λy, x + λy > ≥ 0 ∀λ ∈ R

⇐⇒ < x, x > +λ < x, y > +λ < y, x > +λ2 < y, y > ≥ 0

⇐⇒ λ2 kyk2 + 2λ < x, y > +kxk2 ≥ 0 (< y, x > = < x, y >)

Ce trinôme en λ est toujours positif si le discriminant :

∆0 =< x, y >2 −kxk2 kyk2 ≤ 0 c-à-d |< x, y >| ≤ kxk kyk

Dans le cas où < x, y > n'est pas réel, on peut se ramener au cas précédent
en multipliant le vecteur x par un complexe de module 1 (par exemple par
< x, y >
).
| < x, y > |
2. Il sut de montrer que kx + yk2 ≤ (kxk + kyk)2
ce qui équivaut à

kxk2 + kyk2 + 2 < x, y >≤ kxk2 + kyk2 + 2kxkkyk

soit
< x, y > ≤ kxkkyk qui est vrai.

3. Facile

Orthogonalité et angles dans un espace vectoriel réel

Dénition 1.1.5
Deux vecteurs x et y de E sont dits orthogonaux si < x, y > = 0.
Dénition 1.1.6
Soient x et y deux vecteurs non nuls de E .
L'angle entre ces deux vecteurs est déni à l'aide du produit scalaire par :
< x, y > = kxk kyk cos(θ)

8
Cette écriture a un sens. En eet, comme

| < x, y > | ≤ kxkkyk

alors :
< x, y > < x, y >
−1 ≤ ≤1 c-à-d ∃ θ / cos(θ) =
kxkkyk kxkkyk

1.2 Espaces de Banach et de Hilbert

1.2.1 Espace complet
Dénition 1.2.1 (Suite de Cauchy)
Soit E un espace vectoriel normé (evn).
Une suite (un ) d'éléments de (E) est dite de Cauchy si :
∀ ε > 0 ; ∃ N ∈ N ; ∀ n , m ≥ N ; kum − un k < ε

qui s'écrit aussi en posant m = n + p

∀ ε > 0 ; ∃ N ∈ N ; ∀ n ≥ N ; ∀ p ∈ N∗ ; kun+p − un k < ε

Remarque :
On voit donc qu'une suite de Cauchy est une suite d'éléments de E qui se rapprochent
l'une de l'autre à partir d'un certain rang. Les suites de Cauchy peuvent donc converger.
En fait, elles convergent si leur limite reste dans l'espace E .

Dénition 1.2.2 (Espace complet)

Un espace vectoriel E normé est dit complet si toute suite de Cauchy d'éléments de E est
convergente dans E .

Exemples :
1. R et C sont des espaces complets.
2. Q n'est pas complet. En eet :
Du fait de la densité√de Q dans R, on peut trouver une suite (un ) d'éléments de
Q qui converge vers 2. (u√ n ) est donc une suite de Cauchy dans Q. Mais (un ) ne
converge pas dans Q car 2 ∈ / Q.

Exercice 3
Montrer que l'ensemble C ([a, b]) des fonctions continues sur [a, b] n'est pas complet.

9
1.2.2 Espace de Banach
Dénition 1.2.3 Un espace vectoriel normé et complet est dit espace de Banach.

Exemple 1.2.1 Soit Ω un ouvert de Rn .

On dénit l'espace des fonctions absolument intégrables sur Ω par :
 Z 
1
L (Ω) = f : Ω → R / |f (x)| dx < +∞
Ω

N.B : L'intégrale sur Ω est ici une intégrale multiple si n ≥ 2.

L1 (Ω) peut être muni de la norme k.k1 suivante :

Z
1
∀ f ∈ L (Ω) ; kf k1 = |f (x)| dx
Ω

Et on a la propriété suivante :

Propriété
L'espace L1 (Ω) muni de la norme k.k1 est un espace de Banach.

1.2.3 Espace de Hilbert

Dénition 1.2.4 Un espace de Hilbert est un espace vectoriel muni d'un produit scalaire
et qui est complet pour la norme associée à ce produit scalaire.

Exemple 1.2.2 L'ensemble Rn muni du produit scalaire usuel est un espace de Hilbert.
Proposition 1.2.1
Tout espace de Hilbert est un espace de Banach.

Preuve
Evident. A partir d'un produit scalaire sur E , on peut dénir une norme.

1.3 Espaces fonctionnels

Dénition 1.3.1 Un espace fonctionnel est un espace vectoriel dont les éléments sont
des fonctions.
Exemple 1.3.1 C n ([a, b]) est l'espace des fonctions dénies sur [a, b] dont toutes les
dérivées d'ordre p (1 ≤ p ≤ n) existent et sont continues sur [a, b].

10
1.3.1 Espaces de fonctions intégrables
Soit Ω un ouvert de Rn .
ˆ On note par Lp (Ω) (p ∈ N∗ ) l'espace des fonctions réelles u dénies sur Ω telles
que : Z
|u(x)|p dx < +∞ ( c-à-d existe )
Ω
Ici l'intégrale est dans Rn . C'est une intégrale simple si n = 1, intégrale double si
n = 2, intégrale triple si n = 3.
ˆ L'espace Lp (Ω) peut être muni de la norme dénie par :
Z 1/p
p p
∀ u ∈ L (Ω) ; kukLp (Ω) = |u(x)| dx
Ω

En particulier :
pour p = 1 : L1 (Ω) est nommé espace des fonctions absolument intégrables (ou som-
mables) :  
Z
1
L (Ω) = u : Ω → R / |u(x)| dx < +∞
Ω

L (Ω) peut être muni de la norme dénie par :

1

Z
1
∀ u ∈ L (Ω) ; kukL1 (Ω) = |u(x)| dx
Ω

pour p = +∞ : On dénit aussi l'espace L∞ (Ω) qui est l'espace des fonctions bornées sur
Ω c'est-à-dire :

L∞ (Ω) = {u : Ω → R / ∃ M ∈ R ; |u(x)| < M ; ∀ x ∈ Ω}

À cet espace, on associe une norme dite norme innie dénie par :

∀ u ∈ L∞ (Ω) : kukL∞ (Ω) = sup |u(x)|

x∈Ω

Théorème 1.3.1 Lp (Ω) muni de la norme k.kLp (Ω) est un espace de Banach (c'est-à-dire
qu'il est complet).

1.3.2 Le cas particulier de l'espace L2(Ω)

ˆ Soit Ω un ouvert de Rn .
On dénit l'espace L2 (Ω) des fonctions de carrés intégrables sur Ω par :
( Z 1/2 )
L2 (Ω) = u : Ω → R / |u(x)|2 dx < +∞
Ω

11
 ˆ L2 (Ω) peut être muni du produit scalaire < . , . >L2 (Ω) déni par :
Z
2
∀ u, v ∈ L (Ω) ; < u, v >L2 (Ω) = u(x)v(x) dx
Ω

et donc d'une norme notée k.k2 dénie par :

Z 1/2
2 2
p
∀ u ∈ L (Ω) ; kukL2 (Ω) = < u, u >L2 (Ω) = |(u(x)| dx
Ω

et On a le théorème suivant :

Théorème 1.3.2 (Fischer - Riesz) Muni du produit scalaire ci dessus, l'espace

L2 (Ω) est un espace de Hilbert (c'est-à-dire complet pour la norme associée à ce
produit scalaire).

Propriétés
1. Si u, v ∈ L2 (Ω) alors u v ∈ L1 (Ω).
2. Si u, v ∈ L2 (Ω) alors u + v ∈ L2 (Ω).
3. Si λ ∈ R ; u ∈ L2 (Ω) alors λu ∈ L2 (Ω).

1.4 Fonctions tests

Dénition 1.4.1 Soit Ω un ouvert de Rn et ϕ : Ω −→ R.
On appelle support de ϕ le plus petit fermé en dehors duquel ϕ est nulle :
supp (ϕ) = {x ∈ Ω / ϕ(x) 6= 0}
Exemple 1.4.1
La fonction nulle : θ : R −→ R
x 7−→ θ(x) = 0
supp (θ) = ∅ = ∅
Exemple 1.4.2
ϕ : R −→ R
x 7−→ 1
supp (ϕ) = R = R
Exemple 1.4.3
La fonction de Schwartz : ϕ : R −→ R
exp( x21−1 ) si |x| < 1
x 7−→
0 sinon
supp (ϕ) = ] − 1, 1[ = [−1, 1]

12
Dénition 1.4.2 On note par D(Ω) l'ensemble des fonctions ϕ de Ω dans R, de classe
C ∞ , et à support borné (donc compact).

D(Ω) = ϕ : Ω → R / ϕ ∈ C ∞ (Ω) et supp (ϕ) est borné



On dit aussi que D(Ω) est l'ensemble des fonctions tests.

Exemple 1.4.4
La fonction nulle : θ : R −→ R est une fonction test c-à-d θ ∈ D(R)
x 7−→ 0
∞

θ ∈ C (R)
En eet : supp (θ) = ∅ qui est fermé borné donc compact.

Exemple 1.4.5
La fonction de Schwartz : supp (ϕ) = [−1, 1] bornée.
On peut montrer aussi que ϕ ∈ C ∞ (R) et donc ϕ ∈ D(R).
La fonction ϕ peut être représentée graphiquement comme suit :

φ( x )

−1 1 x

1.5 Espaces de Sobolev

Soit Ω un ouvert de Rn .
On s'intéresse à l'espace de Sobolev d'ordre 1 déni par :
 
1 2 ∂u 2
H (Ω) = u : Ω → R / u ∈ L (Ω) ; ∈ L (Ω) ; ∀ i = 1, . . . , n
∂xi

L'espace de Sobolev d'ordre 2 est déni par :

∂ 2u
 
2 2 ∂u 2 2
H (Ω) = u ∈ L (Ω) ; ∈ L (Ω) ; ∈ L (Ω) ; ∀ i, j = 1, . . . , n
∂xi ∂xi ∂xj

13
En particulier, en dimension 1, c'est-à-dire si Ω est un ouvert de R alors :

H 1 (Ω) = u ∈ L2 (Ω) ; u0 ∈ L2 (Ω)



H 2 (Ω) = u ∈ L2 (Ω) ; u0 , u00 ∈ L2 (Ω)



Dans cette extension, on voit que :

H 0 (Ω) = L2 (Ω)

Et on a le théorème suivant :

Théorème 1.5.1 L'espace H 1 (Ω) peut être muni du produit scalaire suivant :
∀u, v ∈ H 1 (Ω) ; < u, v >H 1 (Ω) = < u, v >L2 (Ω) + < ∇u, ∇v >L2 (Ω)

et dénit un espace de Hilbert pour ce produit scalaire.

A partir de ce produit scalaire, on dénit une norme dans H 1 (Ω) par :

Z Z
1
∀u ∈ H (Ω) ; kuk2H 1 (Ω) = 2
|u(x)| dx + |∇u(x)|2 dx
Ω Ω
= kuk2L2 (Ω) + 2
k∇ukL2 (Ω)

On dénit aussi l'espace H01 (Ω) par :

H01 (Ω) = u ∈ H 1 (Ω) , u = 0 sur ∂Ω


∂Ω étant la frontière de Ω.
Par exemple si n = 1 et Ω =]a, b[ ouvert de R alors :
H01 (]a, b[) = u ∈ H 1 (]a, b[) , u(a) = u(b) = 0


1.6 Inégalité de Poincaré

Théorème 1.6.1 (Inégalité de Poincaré)
Soit u ∈ H01 (a, b) ; alors il existe une constante C (fonction de (b − a)), telle que :
kukL2 (a,b) ≤ C ku0 kL2 (a,b)

Preuve
Soit x ∈ (a, b), on a :
Z x
u0 (t) dt car u ∈ H01 (a, b)


u(x) =
a

14
En élevant au carré et en utilisant l'inégalité de Cauchy-Schwarz, on a :
Z x 2
2 0
u (x) = u (t) dt
a
Z x  Z x 
02
≤ u (t) dt 1 dt
a a
2
= (x − a) ku0 kL2 (a,b)

On intègre entre a et b :
Z b Z b
2
2
u (x) dx = (x − a) ku0 kL2 (a,b) dx
a a
(b − a)2 0 2
= ku kL2 (a,b)
2
On obtient alors :
b−a
kukL2 (a,b) ≤ √ ku0 kL2 (a,b)
2

En dimension supérieure à 1, on a :

kukL2 (Ω) ≤ CΩ k∇ukL2 (Ω) ∀ u ∈ H01 (Ω)

A partir de l'inégalité de Poicaré, on peut démontrer le lemme suivant :

Lemme 1.6.1
L'expression dénie par :
∀ u ∈ H01 (Ω) , | u |1 = k∇ukL2 (Ω)

est une norme sur l'espace H01 (Ω), équivalente à la norme k · k sur H01 (Ω) induite de celle
de H 1 (Ω), c.a.d qu'il existe deux constantes C1 et C2 telles que

C1 | u |1 ≤ kukH 1 (Ω) ≤ C2 | u |1

Preuve
 1/2
kukH 1 (Ω) = kuk2L2 (Ω) + k∇uk2L2 (Ω)
| u |1 = k∇ukL2 (Ω)

Il est clair tout d'abord que

| u |1 ≤ kukH 1 (Ω) (1)

D'autre part, on a d'après l'inégalité de Poincaré :

kuk2L2 (Ω) ≤ C 2 k∇uk2L2 (Ω)

15
 =⇒ kuk2H 1 (Ω) ≤ (1 + C 2 ) | u |21

√
=⇒ kukH 1 (Ω) ≤ 1 + C 2 | u |1 (2)

(1) et (2) =⇒ l'équivalence des normes.

1.7 Trace d'une fonction

Dénition 1.7.1
Soit u ∈ H 1 (Ω). Notons par ∂Ω la frontière de Ω.
La restriction au bord ∂Ω de la fonction u est appelée la trace au bord de u. Elle est notée
u| ∂Ω .

On peut montrer que :

Lemme 1.7.1
Il existe une constante C telle que :
ku| ∂Ω kL2 (∂Ω) ≤ C k∇ukH 1 (Ω)

L'ensemble des traces au bord des fonctions de H 1 (Ω) est noté H 1/2 (∂Ω).

1.8 Formules de Green

Nous rappelons ici les formules de Green qui sont tout simplement des extensions en
dimension supérieure de l'intégration par parties :
Soit Ω un ouvert de Rn . Notons par ∂Ω sa frontière et n la normale unitaire extérieure
à ∂Ω.

Ω

Formule de Stokes
Z Z Z
div(u) v dΩ = − u · ∇v dΩ + (u · n) v dσ
Ω Ω ∂Ω
ici u est un vecteur et v une fonction.

16
La formule de Stokes est une généralisation à un espace à plusieurs dimensions de la
formule d'intégration par parties :
Z b Z b
0
u v dx = − u v 0 dx + [uv]ba
a a

En particulier :
ˆ Si v=1, on obtient :
Z Z
div(u) dΩ = u · n dσ
Ω ∂Ω

ˆ Si on remplace u par ∇u dans la formule de Stokes, on aura :

Z Z Z
div(∇u) v dΩ = − ∇u · ∇v dΩ + (∇u · n) v dσ
Ω Ω ∂Ω

qui s'écrit aussi :

Z Z Z
∂u
∆u v dΩ = − ∇u · ∇v dΩ + v dσ (Formule de Green)
Ω Ω ∂Ω ∂n
∂u
où = ∇u · n
∂n

17
Chapitre 2
Méthode des Éléments nis
La méthode des éléments nis est une méthode d'approximation des équations aux dérivées
partielles et est l'une des trois méthodes les plus utilisées, avec les méthodes des diérences
nies et de volumes nis. Cette méthode est apparue dans les années 50 et a été appliquées
au départ par les ingénieurs et les mécaniciens pour le calcul de l'aérodynamique autour
des avions. La méthode des éléments nis est actuelement la méthode de base en calcul
de structures et en génie civil.

2.1 Formulation variationnelle

An d'établir la formulation variationnelle d'une EDP avec conditions aux limites, nous
allons considérer deux exemples modèles d'EDP.

2.1.1 Exemple en 1D
On considère l'équation diérentielle (dite de poisson) avec conditions aux limites du type
Dirichlet :
−u00 (x) = f (x) pour a < x < b

(P ) (2.1)
u(a) = u(b) = 0 f ∈ L2 ([a, b])

Dénition 2.1.1 On appelle solution classique (ou solution forte) du problème (P ) une
fonction u de classe C 2 ([a, b]) vériant :
∀ x ∈]a, b[ ; −u00 (x) = f (x)


u(a) = u(b) = 0

Multiplions l'équation (2.1) par une fonction v assez régulière et intégrons entre a et b,
on aura : Z b Z b
00
− u (x) v(x) dx = f (x) v(x) dx
a a
On intègre par partie, on obtient :
Z b Z b
0 0 0 b
u (x) v (x) dx − [u (x) v(x)]a = f (x) v(x) dx
a a

18
Si on prend la fonction v de sorte à vérier la condition initiale, c'est-à-dire telle que
v(a) = v(b) = 0, alors on obtient :
Z b Z b
0 0
u (x) v (x) dx = f (x) v(x) dx
a a

On voit que les intégrales dans cette équation ont un sens dès lors que u et v ∈ H01 ([a, b]).

Le problème revient donc à :

 Trouver Z u ∈ H01 ([a, b]) tel que


b Z b
(Q) 0 0
 u (x) v (x) dx = f (x) v(x) ∀ v ∈ H01 ([a, b])
a a

Le problème (Q) est appelé formulation variationnelle (ou formulation faible) du

problème (P ).

Toute solution du problème (Q) est appelée solution faible.

Et on a bien sûr, si u est solution forte de (P ) alors u est solution faible.

Remarque 2.1.1 Un avantage important de la formulation faible (Q) par rapport au

problème (P ) est que nous faisons apparaître seulement les dérivées premières dans (Q)
alors que l'on a des dérivées secondes dans (P ). La régularité de la solution a été donc
aaiblie.

2.1.2 Exemple en 2D
Soit Ω un domaine régulier de R2 .
On note par ∂Ω la frontière de Ω et n la normale unitaire extérieure à ∂Ω.

Ω

On considère le problème d'équilibre suivant :


 −∆u + u = f pour (x, y) ∈ Ω̊ ; f ∈ L2 (Ω)
(P ) ∂u
 =0 pour (x, y) ∈ ∂Ω
∂n
Soit v une fonction assez régulière. Multiplions par v et intégrons sur Ω :
ZZ ZZ
(−∆u + u) v dxdy = f v dxdy
Ω Ω

19
En utilisant la formule de Green, on aura :
ZZ Z ZZ ZZ
∂u
∇u · ∇v dxdy − v dσ + u v dxdy = f v dxdy
Ω ∂Ω ∂n Ω Ω

∂u
or la CL : = 0 sur ∂Ω =⇒
∂n
ZZ ZZ ZZ
∇u · ∇v dxdy + uv dxdy = f v dxdy
Ω Ω Ω

Ici les termes de cette équation ont un sens dès lors que u, v ∈ H 1 (Ω).
La formulation variationnelle du problème (P ) s'écrit alors :

 Trouver u Z Z∈ H (Ω) vériant Z Z

 1
ZZ
(Q)
 ∇u · ∇v dxdy + uv dxdy = f v dxdy ∀ v ∈ H 1 (Ω)
Ω Ω Ω

2.2 Existence et unicité de la solution du problème va-

riationnel
Les exemples (1D) et (2D) présentés montrent que la formulation variationnelle s'écrit
d'une manière générale sous la forme :
Trouver u ∈ V telle que

(Q)
a(u, v) = L(v) ∀v ∈ V
Dans l'exemple (1D), on a :
 Z b
u0 v 0 dx

 a(u, v) =


Za b
f v dx et V = H01 ([a, b])

 L(v) =


a

Dans l'exemple (2D), on a :

 ZZ ZZ
 a(u, v) = ∇u · ∇v dxdy + u v dxdy


Z ZΩ Ω

 L(v) = f v dxdy et V = H 1 (Ω)



Ω

L'existence et l'unicité de la solution u du problème variationnel (Q) sont assurées par le

théorème de Lax-Milgram qui est le suivant :
Théorème 2.2.1 (Lax-Milgram) Soit V un espace de Hilbert.
a : V ×V −→ R une forme bilinéaire, continue et coercive
(u, v) 7−→ a(u, v)
et
L : V −→ R une forme linéaire et continue sur V
u 7−→ L(u)

20
Alors le problème variationnel :
Trouver u ∈ V telle que

(Q)
a(u, v) = L(v) ∀v ∈ V
admet une solution unique u dans V .

Rappels

a(λu + µv, w) = λa(u, w) + µa(v, w)
ˆ a est bilinéaire ⇐⇒
a(u, λv + µw) = λa(u, v) + µa(u, w)

ˆ a est continue sur V ⇐⇒ ∃ M > 0 ; |a(u, v)| ≤ M kukV · kvkV ∀ u, v ∈ V

ˆ a coercive sur V ⇐⇒ ∃ α > 0 ; a(u, u) ≥ αkuk2V ∀u ∈ V

ˆ L linéaire sur V ⇐⇒ L(λu + µv) = λL(u) + µL(v)

ˆ L continue sur V ⇐⇒ ∃ K > 0 , |L(u)| ≤ K kukV

2.3 Principe de la méthode des éléments nis

Dans le but d'obtenir l'approximation numérique uh de la solution u du problème va-
riationnel (Q), la méthode des éléments nis consiste à remplacer l'espace V (qui est en
général de dimension innie) par un sous espace Vh de V de dimension nie.
Le problème variationnel continue Q sera donc remplacé par le problème variationnel
approché ou discret (Qh ) :
Trouver uh ∈ Vh

(Qh )
a(uh , vh ) = L(vh ) ∀ vh ∈ Vh

Cette méthode est appelée méthode d'éléments nis Galerkin.

En fait, la recherche de la solution approchée uh dans l'espace Vh de dimension nie va
permettre d'avoir un nombre ni d'inconnues (appelées aussi degrés de liberté) et qui
seront en fait les composantes de la solution approchée uh dans une base (ϕi )1≤i≤N de Vh .
N
X
uh (x) = ui ϕi (x)
i=1

Ces composantes ui de uh dans la base (ϕi )1≤i≤N seront calculées en résolvant un système
linéaire.
Expliquons ceci :
Soit (ϕ1 , ϕ2 , . . . , ϕN ) une base de Vh .
Si l'on décompose la solution approchée uh dans cette base :
N
X
uh (x) = ui ϕi (x)
i=1

21
le problème approché (Qh ) s'écrit :

 Trouver u1 , u2 , . . . , uN !tels que



XN

 a ui ϕi , vh = L(vh ) ∀ vh ∈ V h
i=1

La forme a étant bilinéaire, alors on a :

N
X
ui a(ϕi , vh ) = L(vh ) ∀ vh ∈ Vh
i=1

Comme L est linéaire alors, pour que cette relation soit vraie ∀ vh ∈ Vh , il faut et il sut
qu'elle le soit pour chacune des fonctions de base de Vh .
Le problème variationnel approché devient alors :

 Trouver u1 , u2 , . . . , uN tels que



N
(Qh ) X

 ui a(ϕi , ϕj ) = L(ϕj ) ∀ j = 1, 2, . . . , N
i=1

qui s'écrit sous la forme d'un système linéaire :

    
a(ϕ1 , ϕ1 ) a(ϕ2 , ϕ1 ) . . . a(ϕN , ϕ1 ) u1 L(ϕ1 )
    
    
 a(ϕ1 , ϕ2 ) a(ϕ2 , ϕ2 ) . . . a(ϕN , ϕ2 )   u2   L(ϕ2 ) 
    
  = 
    
.   ..   .. 
 

 .
.  .   . 
    
    
a(ϕ1 , ϕN ) a(ϕ2 , ϕN ) . . . a(ϕN , ϕN ) uN L(ϕN )
ou encore :
AU = b

Remarque 2.3.1 On peut montrer que la matrice A est dénie positive c'est-à-dire :
→
−
(X t AX ≥ 0 et X t AX = 0 ⇒ X = 0 )

Elle est donc inversible.

22
Preuve
  
a(ϕ1 , ϕ1 ) a(ϕ2 , ϕ1 ) . . . a(ϕN , ϕ1) x1
  
  
 a(ϕ1 , ϕ2 ) a(ϕ2 , ϕ2 ) . . . a(ϕN , ϕ2)   x2 
  
X t AX = (x1 , x2 , . . . , xN ) 
  
..   .. 
 
.  . 


  
  
a(ϕ1 , ϕN ) a(ϕ2 , ϕN ) . . . a(ϕN , ϕN ) xN
N N
!
X X
= xj a(ϕi , ϕj )xi a bilinéaire ⇒
j=1 i=1
N N
!
X X
= a xi ϕ i , ϕ j xj
i=1 j=1
= a(u, u) ≥ αkuk2 (car a est coercive)
On déduit donc que
X t AX ≥ 0
De plus :
N
X
t 2
X AX = 0 =⇒ αkuk = 0 =⇒ u = 0 =⇒ xi ϕi = 0
i=1
=⇒ xi = 0 ∀ i = 1, . . . , N car (ϕi )1≤i≤N est une base de Vh
→
−
=⇒ X = 0 D'où le résultat.

Quelques commentaires sur le choix des fonctions de base et de l'espace d'ap-

proximation Vh
1. La matrice A = (aij ) est à priori pleine. Dans la pratique, et an de limiter le
stockage dans la matrice A, les fonctions de base (ϕi ) sont choisies à support
assez petit c'est-à-dire que la fonction de base ϕi sera choisie nulle partout sauf sur
quelques mailles. De cette manière, la matrice A du système sera creuse c'est-à-dire
contient beaucoup de 0.
2. De même, l'espace Vh doit être facile à construire : on peut par exemple choisir
l'espace des fonctions anes ou des polynômes de degrés 2 par morceaux.

2.4 Mise en ÷uvre de la méthode en dimension 1

Reprenons l'équation de Poisson en 1D avec conditions aux limites de Dirichlet homo-
gènes :
−u00 (x) = f (x)

∀ x ∈]a, b[
(P )
u(a) = u(b) = 0
Nous avons vu que la formulation variationnelle du problème (P ) s'écrit :
Trouver u ∈ H01 ([a, b]) telle que

(Q)
a(u, v) = L(v) ∀ v ∈ H01 ([a, b])

23
avec  Z b
u0 (x) v 0 (x) dx

 a(u, v) =


Za b

 L(v)

 = f (x) v(x) dx
a
L'espace d'approximation Vh est tel que Vh ⊂ H01 ([a, b]), on peut alors choisir Vh ⊂
C 0 ([a, b]) (fonctions continues sur [a, b]).

2.4.1 Maillage et fonctions de base

On commence par découper le domaine [a, b] en N + 2 points x0 , x1 , . . . , xN +1 tels que
a = x0 < x1 < x2 < . . . < xN +1 = b
b−a
xi = a + ih avec h =
N +1

x 0=a x1 x2 xN x N +1 =b

Cette subdivision est appelée maillage de [a, b].

2.4.2 Méthode des éléments nis P1

Dénissons l'espace Vh par :
vh : [a, b] −→ R ; vh ∈ C 0 ([a, b]) ; vh ∈ P1 sur chaque segment [xi , xi+1 ]
 
Vh =
vh (a) = vh (b) = 0
P1 est l'espace des polynômes de degrés ≤ 1.

Exemple de vh :
vh

v h ( a)=0 v h ( b)=0
a b

24
xN
x2
Le problème variationnel discret s'écrit alors :

Trouver uh ∈ Vh telle que


(Qh )
a(uh , vh ) = L(vh ) ∀ vh ∈ Vh

Et nous avons le théorème :

Théorème 2.4.1
ˆ Les fonctions de l'espace Vh sont entièrement déterminées par leurs valeurs aux
n÷uds x1 , x2 , . . . , xN .
ˆ dim Vh = N est une base de Vh (ϕi )1≤i≤N est dénie par

 ϕi ∈ Vh
1 si i = j

 ϕi (xj ) = δij =
0 si i 6= j
(Symbole de Kronecker)

Les (ϕi ) sont appelées fonctions chapeaux.

Preuve
i/ Soit vh ∈ Vh alors vh ∈ P1 sur chaque [xi , xi+1 ].
c'est-à-dire : vh (x) = αx + β pour x ∈ [xi , xi+1 ].
A partir des valeurs v(xi ) et v(xi+1 ), on peut déterminer d'une manière unique α
et β et donc vh (x).
Il s'en suit que la fonction vh est entièrement déterminée sur [a, b] par ses valeurs aux
points x1 , x2 , . . . , xN . Ces valeurs vh (xi ) aux n÷uds sont appelées degrés de liberté
de l'espace Vh .
ii/ La fonction de base (ϕi ) peut être représentée comme suit :

x 0=a x1 x2 x i−1 xi x i+1 xN x N +1 =b

Montrons que (ϕi )1≤i≤N forme bien une base de Vh .

(ϕi )1≤i≤N est libre. En eet :
si
N
X
λj ϕj (x) = 0 ∀ x ∈ [a, b]
j=1

alors
N
X
pour x = xi on a : λj ϕj (xi ) = 0 ⇒ λi = 0 (i = 1, . . . , N )
j=1

25
Montrons que (ϕi )1≤i≤N est génératrice de Vh .
N
X
Soit vh ∈ Vh , et posons ω(x) = vh (xj ) ϕj (x)
j=1
on a bien sûr ω ∈ Vh .
De plus ∀ i = 1, . . . , N ; ω(xi ) = vh (xi )
or d'après i/ les fonctions de Vh sont entièrement déterminées d'une manière unique par
leurs valeurs aux points xi (i = 1, . . . , N ).
Il s'en suit que vh = ω .
c'est-à-dire :
XN
vh (x) = vh (xj ) ϕj (x)
j=1

La famille (ϕi ) est génératrice de Vh . C'est donc une base de Vh .

→ En décomposant la solution approchée uh sur cette base de fonctions chapeaux :

N
X
uh (x) = ui ϕi (x)
i=1

on se ramène à la résolution d'un système linéaire :

AU = b
dont les inconnues sont les degrés de liberté (u1 , u2 , . . . , uN )t = U .
- La matrice A = (aij )1≤i,j≤N est dénie par
Z b
aij = ϕ0i (x) ϕ0j (x) dx est appelée matrice de rigidité du système
a

- Le second membre b = (bj )1≤j≤N

Z b
bj = f (x) ϕj (x) dx est appelé vecteur de charge.
a

Evaluation de la matrice A et du second membre b

Pour calculer la matrice A, il sut d'évaluer explicitement les ϕi (x) pour 1 ≤ i ≤ N .

a x i−1 xi x i+1 b

26
Il sut d'écrire
sur chaque [xj , xj+1 ]

 ϕi (x) = αx + β
1 si i = j
x 0=a  ϕi (xj ) = δij = x N =b
0 si i 6= j
On obtient facilement
si x ≤ xi−1 ou x ≥ xi+1

 0
 x − xi−1


ϕi (x) = si xi−1 ≤ x ≤ xi
h
 xi+1 − x
si xi ≤ x ≤ xi+1



h
on a : Z b
aij = ϕ0i (x) ϕ0j (x) dx
a

φi−2 φi−1 φi φi+1 φi+2

i−2 i−1 i i+1 i+2

Comme le support de la fonction de base ϕi , qui est l'ensemble des points où ϕi n'est pas
nulle, est égal à [xi−1 , xi+1 ]
alors on a :
aij = 0 si j ≤ i − 2 ou j ≥ i + 2 c'est-à-dire si |i − j| ≥ 2
Pour déterminer la matrice A, il reste donc à évaluer les termes aii , ai i+1 , ai−1 i . Cette
matrice A est donc tridiagonale.
Z b
aii = ϕ0i (x) ϕ0i (x) dx
Za xi+1
= ϕ02
i (x) dx
xi−1
Z xi Z xi+1
1 1 1 1
= · dx + (− ) (− ) dx
xi−1 h h xi h h
1 1
= 2 · h+ 2 h
h h
2
=
h
Z b
ai i+1 = ϕ0i (x) ϕ0i+1 (x) dx
Za xi+1
1 1 1
= (− ) . dx = − 2 . h
xi h h h
1
= −
h

27
 Z b
ai−1 i = ϕ0i−1 (x) ϕ0i (x) dx
Za xi
1 1
= (− ) . dx
xi−1 h h
1
= −
h
La matrice de rigidité A du système est tridiagonale et s'écrit :
 
2 −1 0 ... 0
−1 2 −1 0 ... 0 
 
 0 −1 2 −1 0 . . . 0 
 
1 0 
A= 
h  .. .. . . . . . . . .

. . . . . .
 
0 
.
 

 . . −1 2 −1

0 0 . . . 0 −1 2

Assemblage de la matrice du système

En général, pour déterminer la matrice de rigidité A du système, on est amené à calculer
des coecients aij = a(ϕi , ϕj ) qui sont donnés par des intégrales de la forme :
Z b
aij = H(ϕi , ϕj ) dx
a

(Dans notre cas par exemple H(ϕi , ϕj ) = ϕ0i (x) ϕ0j (x))
Les aij peuvent être calculés en sommant sur chaque élément [xk , xk+1 ] du maillage, c'est-
à-dire :
XN Z xk+1
aij = H(ϕi , ϕj ) dx
k=0 xk

XN
(k)
= aij (2.2)
k=0

Une manière assez simple mais très coûteuse de déterminer la matrice A, est de calculer
chaque coecient aij en utilisant (2.2), par l'algorithme :

pour i = 1 à N faire

←− Boucles sur les n÷uds
pour j = 1 à N faire

aij = 0
pour k = 0 à N faire ←− Boucle sur les éléments
(k)
calculer aij
(k)
aij = aij + aij
Fin pour
Fin pour
Fin pour

28
 (k)
Cet algorithme nécessite un nombre de calcul des aij de l'ordre de N 2 × (N + 1). Si le
nombre de n÷uds N = 100 =⇒ ' 106 calcul.
(k)
Toutefois, l'algorithme tel qu'il est écrit ci-dessus calcule une grande perte des aij qui
(k)
sont nuls. Or les aij sont nuls si et seulement si les n÷uds i et j appartiennent à l'élément
k.

Une manière plus simple et moins coûteuse est de boucler sur les éléments k (ici les
(k)
segments) pour ne calculer que les termes aij utiles, mettre ces termes dans des matrices
de rigidité élémentaire associées à chaque élément k , puis procéder à l'assemblage de ces
matrices élémentaires dans une matrice A globale.

Élément 0 Élément 1 Élément N

x 0=a x1 x2 xN x N +1 =b

En pratique
On se place sur l'élément 0.

r1

x0 x1

Z x1 Z x1  2
(0) 1 1
On calcule a11 = ϕ01 (x)ϕ01 (x) dx = dx = .
x0 x0 h h
Et on a la matrice de rigidité élémentaire associée à l'élément 0 :

(0) 1
(a11 ) =
h
Dans la matrice globale, elle est vue comme suit :
 (0)   
a11 0 . . . 0 1 0 ... 0
   
   
 1 0 0 . . . 0
 0 0 . . . 0 
 
A(0) =  = 
   
 ... .. 
 . ..  h 

 .. .

 .  
   
0 0 ... 0 0 0 ... 0

Sur l'élément 1 :

29
 ρ1 ρ2

x1 x2

(1) (1) (1) (1)

On calcule a11 , a12 , a21 , a22
Z x2
(1)
a11 = ϕ01 (x)ϕ01 (x) dx
x1
Z x2  2
1
= − dx
x1 h
1
=
Zh x2
(1)
a12 = ϕ01 (x)ϕ02 (x) dx
Zx1x2
1 1
= − dx
x1 hh
1
= −
h
(1) 1
a21 = −
h
(1) 1
a22 =
h
=⇒ Matrice de rigidité élémentaire associée à l'élément 1 :
 
(1) (1)  
a11 a12 1 −1
 1
=
 
h

(1)
a21 a22
(1) −1 1

Dans la matrice globale, elle est vue comme suit :

 
1 −1 0 ... 0
−1 1 0 . . . 0
1 
A(1) =  0 0 0 . . . 0

h  .. .. .. .. 

 . . . .
0 0 0 ... 0
Sur l'élément 2 :

r2 r3

x2 x3

30
 (2) (2) (2) (2)
On calcule a22 , a23 , a32 , a33
=⇒ Matrice de rigidité élémentaire associée à l'élément 2 :
 
(2) (2)  
a22 a23 1 −1
 1
=
 
h

(2) (2)
a32 a33 −1 1

Dans la matrice globale, elle est vue comme suit :

 
0 0 0 0 ... 0
0 1 −1 0 ... 0
 
1 0 −1 1 0 ... 0
A(2) = 0 0
 
h 0 0 ... 0
 .. .. .. .. .. 

. . . . .
0 0 0 0 ... 0

On obtient alors la matrice de rigidité globale en assemblant les matrices de rigidité dans
chaque élément, par :
N
X
A = A(k)
k=0
 
2 −1 0 0 ... 0
−1 2 −1 0 . . . 0
.. 
 
.

1 0 −1 2 −1
= . .. 
 .. . . . . . . . . .

h .
 
 
0 . . . 0 −1 2

Calcul du vecteur de charge b

L'évaluation du second membre b nécessite le calcul de N intégrales
Z b
bj = f (x) ϕj (x) dx
a

Si la fonction f est simple alors ces intégrales peuvent être calculées d'une manière expli-
cite. Par contre si f a une forme compliquée alors le calcul de ces intégrales doit se faire
par une méthode numérique.
Z b
On utilise les méthodes dites de quadrature numérique qui consistent à approcher f (x) dx
a
par :
Z b Xn
f (x) dx ' ωi f (xi )
a i=0

où les xi sont des points de [a, b] distincts (ce ne sont pas les noeuds du maillage) et les
ωi sont les poids d'interpolation.

31
On commence par utiliser la formule de Chasles :
Z b n−1 Z
X xi+1
f (x) dx = f (x) dx
a i=0 xi
Rx
puis on évalue xii+1 f (x) dx sur chaque segment [xi , xi+1 ] en utilisant une des formules de
quadrature suivantes :

Formule des rectangles à gauche

Z xi+1 Z b n−1
X
f (x) dx ' (xi+1 − xi ) f (xi ) =⇒ f (x) dx ' (xi+1 − xi ) f (xi )
xi a i=0

Formule des rectangles à droite

32
 Z xi+1 Z b n−1
X
f (x) dx ' (xi+1 − xi ) f (xi+1 ) =⇒ f (x) dx ' (xi+1 − xi ) f (xi+1 )
xi a i=0

Ces deux formules sont sont exactes pour les fonctions f constantes.

Formule du point milieu

x i+x i +1
f( )
2
f (x i+1 )
f (x i )

xi x i+x i+1 x i+1

2

Z xi+1 Z b n−1
xi + xi+1 X xi + xi+1
f (x) dx ' (xi+1 − xi ) f ( ) =⇒ f (x) dx ' (xi+1 − xi ) f ( )
xi 2 a i=0
2
Cette formule est exacte pour les polynômes de degré ≤ 1.

Formule des Trapèzes

xi+1
xi+1 − xi
Z
f (x) dx ' ( ) (f (xi ) + f (xi+1 ))
xi | 2 {z }

rectangle gauche + rectangle droite

'
2
n−1
!
Z b
∆x X
=⇒ f (x) dx ' f (x0 ) + 2 f (xi ) + f (xn ) si xi+1 − xi = ∆x = cte
a 2 i=1
Cette formule est aussi exacte pour les polynômes de degré ≤ 1.

33
Formule de Simpson
Cette formule utilise une interpolation P2 aux points xi , xi +x2 i+1 , xi+1 c'est-à-dire fait passer
un polynôme de degré 2 par les points (xi , f (xi )), xi +x2 i+1 , f ( xi +x2 i+1 ) et (xi+1 , f (xi+1 ))

La formule s'écrit :
Z xi+1  
xi+1 − xi xi + xi+1
f (x) dx = ( ) f (xi ) + 4f ( ) + f (xi+1 )
xi 6 2
Cette formule est exacte pour les polynômes de degré ≤ 2.
Pour retrouver la formule de Simpson, il sut d'écrire :
Z xi+1 Z xi+1 xi+1
ax3 bx2

2
f (x) dx ' (ax + bx + c) dx = ( + + cx)
xi xi 3 2 xi

a, b et c sont calculés en écrivant :



 ax2i + bxi + c = f (xi )
 ax2 + bx + c = f (x )

i+1 i+1
i+1 2    
 xi + xi+1 xi + xi+1 xi + xi+1
 a +b +c=f


2 2 2
Remarquons que la formule de quadrature de Simpson s'écrit aussi :
Z xi+1  
1 xi+1 − xi xi + xi+1
f (x) dx ' (f (xi ) + f (xi+1 )) + 2(xi+1 − xi ) f ( )
xi 3 2 2
1
= {F. Trapèze + 2 × F. point milieu}
3

Exercice 4 Programmer ces méthodes en Matlab et comparer.

2.5 Résolution du système linéaire :

Nous avons vu que la MEF aboutit à la résolution d'un système linéaire A U = b où A
est une matrice creuse.
Il existe 2 grandes familles de résolution des systèmes linéaires :
ˆ Méthodes directes : Gauss, Cholesky, Householder.
ˆ Méthodes itératives : Jacobi, Gausse-Seidel, relaxation, gradient conjugué.

34
2.5.1 Méthodes directes
Méthode de Gauss
Cette méthode consiste à écrire la matrice A sous la forme :
Lower
 
A = LU
Upper

où L est triangulaire inférieure à diagonale unité

U est triangulaire supérieure

Le système s'écrit alors : L |{z}

UX = b
z
On résout tout d'abord le système L z = b par une méthode de remontée.
La matrice U s'obtient par le procédé d'élimination de Gauss.
C'est-à-dire si on pose :
 (1) (1) (1)
  (1) 
a11 a12 . . . a1n b1
   
   
(1)
 (1) (1) (1)   (1) 
A =A= a 21 a 22 . . . a2n
 b =  b2 
 ..  .. 
   
 .  . 


(1) (1) (1) (1)
an1 an2 . . . ann bn

Étape 1 :
Mettre à 0 tous les coecients se trouvant sur la 1ère colonne à partie de la ligne 2 c'est-
à-dire pour i ≥ 2
(1)
(2) (1) ai1 (1)
aij = aij − (1) a1j pour j = 1, . . . , n
a11
(1)
(2) (1) ai1 (1)
bi = bi − (1) b1
a11
On obtient :
(1) (1) (1) (1)
   
a11 a12 . . . a1n b1
   
   
(2) (2)   (2) 
0 a22 ... a2n   b2 


(2)
   
A =  b(2) = 
(2) (2) 
   (2) 
 0 a32 . . . a3n   b3 
.. .. .. ..   .. 
  
. . . .   . 


(2) (2) (2)
0 an2 . . . ann bn

⇐⇒ A(2) X = b(2)

35
Étapes suivantes :
On répète le procédé sur les autres colonnes :
(k)
(k+1) (k) aik (k)
aij = aij − (k)
· akj
akk
(k)
(k+1) (k) aik (k)
bi = bi − (k)
· bk
akk
pour i = k + 1, . . . , n
j = k, . . . , n

À la n du procédé, on obtient la matrice U :

 (1) (1) (1)

a11 a12 . . . a1n
 
 
(2) (2) 
 0 a22 . . . a2n 

 
U = 
(3) (3) 
 
 0 0 a33 . . . a3n 
 .. .. .. .. 

 . . . . 
(n)
0 0 0 . . . ann
La matrice L quant à elle, elle s'écrit sous la forme :
 
1 0 ... 0
 
 
 l21 1 . . . 0
 
L=
 

 l31 l32 1 0
 .. .. .. .. 
 
 . . . .
ln1 ln2 . . . lnn−1 1
où
(k)
aik
pour i ≥ 2 ; lik = (k)
akk

2.5.2 Méthodes itératives

On veut résoudre le système
AX = b (2.3)
Les méthodes itératives consistent à construire une suite de vecteurs X (k) qui converge
vers un vecteur X solution du système (2.3).


 a11 x1 + a12 x2 + . . . + a1n xn = b1
 a21 x1 + a22 x2 + . . . + a2n xn = b2

(2.3) ⇐⇒ .. ..


 . .
an1 x1 + an2 x2 + . . . + ann xn = bn


36
 n

 1 X

 x1 = (b1 − a1j xj )
a11



 j=2

 n
1

 X
 x2 = a (b2 −


 a2j xj )
22
⇐⇒ j=1
j6=2
.

..






 n−1
1

 X
(bn −

 x = anj xj )
 n

 ann j=1

Méthode de Jacobi
1. Initialisation du vecteur  
(0)
x1
 
 
 (0) 
X (0) =  x2 

 .. 

 . 
(0)
xn
2. Si on appelle X (k) le vecteur à l'itération k , alors le vecteur X (k+1) à l'itération
(k + 1) est obtenu par :
 n
(k+1) 1 X (k)
(b1 −


 x1 = a1j xj )



 a11 j=2
..

.
n



 (k+1) 1 X (k)
 xi

 = (bi − aij xj )
aii

j=1,j6=i

Le procédé est arrêté lorsque :

Résidu = kb − A X (k+1) k < Seuil xé

Méthode de Gauss-Seidel
1. Initialisation du vecteur  
(0)
x1
 
 
 (0) 
X (0)  x2 
= 
.
 . 
 . 
(0)
xn
(k+1)
2. Dans le procédé de Gauss-Seidel, une fois une valeur xi à l'itération (k + 1) est
calculée, on l'utilise pour calculer les autres valeurs à la même itération.

37
 L'algorithme s'écrit :

pour i = 1, . . . n



 i−1 n

(k+1) 1 X (k+1)
X (k)
xi = (bi − aij xj − aij xj )

 aii j=1 j=i+1

n pour


Méthode de Gauss-Seidel avec relaxation

La méthode consiste à calculer X (k+1) en faisant une pondération entre le résultat X (k) à
(k+1)
l'itération (k) et celui de Gauss-Seidel XGS à l'itération (k + 1)
(k+1)
X (k+1) = (1 − ω)X (k) + ωXGS

où ω est un paramètre compris entre 0 et 2.

Si ω = 1 : on obtient la méthode de Gauss-Seidel.
Si ω > 1 : la méthode est dite sur-relaxation.
Si ω < 1 : la méthode est dite sous-relaxation.

En général, ω ' 1, 5 marche bien (par expérience).

Remarque 2.5.1 Si la matrice A est à diagonale strictement dominante, alors ces mé-
thodes itératives convergent quelque soit le choix du vecteur initial X (0) .
Une matrice A est à diagonale dominante si :
n
X
|aii | ≥ |aij | ∀ i(1 ≤ i ≤ n)
j=1
j6=i

Lorsque l'inégalité est stricte, A est dite à diagonale strictement dominante.

Si la matrice A est symétrique dénie positive, alors la méthode de Gauss-Seidel converge.

2.6 Problème 1D avec conditions aux limites de Neu-

mann
Nous allons reprendre le problème d'équilibre en 1D avec cette fois-ci en imposant des
conditions aux limites du type Neumann :

−u00 + u = f dans ]a, b[


(P ) 0 0
u (a) = α ; u (b) = β f ∈ L2 ([a, b])

38
2.6.1 Formulation variationnelle
On multiplie par une fonction v assez régulière et on intègre sur [a, b] :
Z b Z b Z b
00
−u v dx + uv dx = f v dx
a a a

En intégrant par parties, on aura :

Z b Z b Z b
0 0 0 b
u v dx − [u v]a + uv dx = f v dx
a a a

soit Z b Z b Z b
0 0
u v dx + uv dx = f v dx + βv(b) − αv(a)
a a a
Le problème variationnel s'écrit alors :
Trouver u ∈ H 1 ([a, b])/

(Q)
a(u, v) = L(v) ∀v ∈ H 1 ([a, b])
avec Z b Z b
0 0
a(u, v) = u v dx + uv dx
a a
Z b
L(v) = f v dx + βv(b) − αv(a)
a

Exercice 5 Montrer à l'aide du théorème de Lax-Milgram que le problème variationnel

(Q) admet une solution unique u dans H 1 ([a, b]).
Remarque 2.6.1 Remarquons que lorsqu'on avait une condition aux limites du type Di-
richlet, nous l'avons imposée sur la fonction v dans l'espace de travail alors que la condi-
tion de Neumann sur la dérivée n'est pas imposée. On dira que la condition aux limites
de Neumann est implicite c'est-à-dire qu'on la trouvera implicitement dans la forme va-
riationnelle.

2.6.2 Formulation variationnelle discrète

Maillage de [a,b]

x 0=a x1 x2 xN x N +1 =b

b−a
∆x = ; xi = xi−1 + ∆x
N +1
= a + i ∆x (i = 0, . . . , N + 1)
On introduit l'espace d'approximation Vh (Eléments nis P1 ) par :
Vh = vh : [a, b] → R ; vh ∈ C 0 ([a, b])/vh |[xi ,xi+1 ] ∈ P1 ∀ i = 0, . . . , N


39
2.6.3 Fonctions de base de Vh
De la même manière que précédemment, on montre que :

dim Vh = N + 2 (ni)

et les fonctions de base de Vh sont les fonctions chapeaux (ϕi ) 0≤i≤N +1 . Il y en a (N + 2)

cette fois-ci.
Les ϕi sont telles que :
 ϕi ∈ P1 sur chaque [xi , xi+1 ]


1 si i = j
 ϕi (xj ) =
0 sinon
ϕ0 :

x 0=a x1 x2 xN x N +1 =b

ϕi (1 ≤ i ≤ N ) :

x 0=a x1 x2 x i−1 xi x i+1 xN x N +1 =b

ϕN +1 :

x 0=a x1 x2 xN x N +1 =b

Le problème variationnel discret s'écrit alors :

Trouver uh ∈ Vh /




 Z b Z b Z b
0 0
(Qh ) uh vh dx + uh vh dx = f vh dx + βvh (b) − αvh (a) ∀ vh ∈ Vh

 a a a

 | {z } | {z }
a(uh ,vh ) L(vh )

En décomposant uh dans la base ϕi :

N
X +1
uh (x) = ui ϕi (x)
i=0

40
et en exprimant le fait que le problème discret (Qh ) serait vrai pour tous les vh ∈ Vh si et
seulement si il l'est pour toutes les fonctions de base ϕj (j = 0, . . . , N + 1), le problème
variationnel discret s'écrit alors :

 Trouver U = (u0 , u1 , . . . , uN +1 ) /
t


NX+1

 ui a(ϕi , ϕj ) = L(ϕj ) ∀ j = 0, . . . , N + 1
i=0
ou encore :
AU = b qui est un système linéaire (N + 2) × (N + 2)
La matrice A = (aij ) 0≤i,j≤N +1 s'écrit :
Z b Z b
0 0
aij = a(ϕi , ϕj ) = ϕi (x)ϕj (x) dx + ϕi (x)ϕj (x) dx
|a } |a
matrice de rigidité Matrice de masse
{z {z }

Le second membre b s'écrit :

Z b
bj = L(ϕj ) = f (x)ϕj (x) dx + βϕj (b) − αϕj (a)
a
Z b
car ϕ0 (a) = ϕ0 (x0 ) = 1
 
b0 = L(ϕ0 ) = f (x)ϕ0 (x) dx − α
a ϕ0 (b) = ϕ0 (xN +1 ) = 0
Z b
car ϕN +1 (a) = ϕN +1 (x0 )
 
=0
bN +1 = L(ϕN +1 ) = f (x)ϕN +1 (x) dx + β
a ϕN +1 (b) = ϕN +1 (xN +1 ) = 1
et pour 1 ≤ j ≤ N :
Z b
car ϕj (a) = ϕj (x0 )
 
=0
bj = L(ϕj ) = f (x)ϕj (x) dx
a ϕj (b) = ϕj (xN +1 ) = 0

Exercice 6 De la même manière que précédemment, calculer explicitement la matrice A

en procédant par :
ˆ Un calcul direct de A.
ˆ Un calcul de A par assemblage des matrices élémentaires sur chaque élément.

Réponse
On va procéder ici par assemblage :
On a : Z b Z b
0 0
aij = ϕi (x)ϕj (x) dx + ϕi (x)ϕj (x) dx
a a

Élément 0 Élément 1 Élément N

x 0=a x1 x2 xN x N +1 =b

41
On se place sur l'élément 0 :
ρ0 ρ1

x0 x1
 
(0) (0)
a00 a01
=⇒ La matrice sur l'élément 0 est : 
 

(0) (0)
a11 a10
Z x1 Z x1
(0) 0 0
a00 = ϕ0 (x)ϕ0 (x) dx + ϕ0 (x)ϕ0 (x) dx
x0 x0
 2
Z x1 Z x1  2
1 x1 − x
= − dx + dx
x0 h x0 h
x
(x1 − x)3 1

1 1
= + −
h h2 3 x0

1 h
= +
h 3

Z x1 Z x1
(0)
a01 = ϕ00 (x)ϕ01 (x) dx + ϕ0 (x)ϕ1 (x) dx
Zx0x1 Z x1 x0
1 1 x 1 − x x − x0
= − dx + . dx
x0 hh x0 h h
Z x1
1 1
= − + 2 (x1 − x)(x − x0 ) dx
h h x0
Z x1
1 1
= − + 2 (x1 − x0 + x0 − x)(x − x0 ) dx
h h x0
( x x )
(x − x0 )2 1 (x − x0 )3 1
 
1 1
= − + 2 (x1 − x0 ) −
h h 2 x0 3 x0
 3 3

1 1 h h
= − + 2 −
h h 2 3
1 h
= − +
h 6

Z x1 Z x1
(0)
a10 = ϕ01 (x)ϕ00 (x) dx
+ ϕ1 (x)ϕ0 (x) dx
x0 x0
x1   Z x1
x − x0 x1 − x
Z
1 1
= − dx + . dx
x0 h h x0 h h
1 h
= − +
h 6

42
 Z x1 Z x1
(0)
a11 = ϕ01 (x)ϕ01 (x) dx + ϕ21 (x) dx
x0 x0
 2
x1 Z x1  2
x − x0
Z
1
= dx + dx
x0 h x0 h
x
1 (x − x0 )3 1

1
= +
h h2 3 x0

1 h
= +
h 3

La matrice pour l'élément 0 est donc :

  
(0) (0) 1 h 1 h

a00 a01 + − +
  h 3 h 6
=
 

(0) (0) 1 h 1 h
a10 a11 −h + 6 h + 3

Dans la matrice globale, elle est vue comme suit :

 1 h
+ 3 − h1 + h6 0 · · ·

h
0
− 1 + h 1 + h 0 · · · 0
 h 6 h 3 
A(0) =  0 0 0 ··· 0


 . . .. .. 
 .. .. . .
0 0 0 ··· 0

On se place sur l'élément 1 :

ρ1 ρ2

x1 x2
 
(1) (1)
a11 a12
=⇒ La matrice sur l'élément 1 est : 
 

(1) (1)
a21 a22
Z x2 Z x1
(1)
a11 = ϕ01 (x)ϕ01 (x) dx + ϕ1 (x)ϕ1 (x) dx
x1 x1
Z x2
 2 Z x2  2
1 x2 − x
= − dx + dx
x1 h x1 h
x
(x2 − x)3 2

1 1
= + −
h h2 3 x1

1 h
= +
h 3

43
 Z x2 Z x1
(1)
a12 = ϕ01 (x)ϕ02 (x) dx + ϕ1 (x)ϕ2 (x) dx
Zx1x2 x2
x1
x2 − x x − x1
Z
1 1
= − dx + . dx
x1 hh x1 h h
1 h
= − +
h 6
(1) 1 h
a21 = − +
h 6
(1) 1 h
a22 = +
h 3
La matrice de l'élément 1 est :
1 h
− h1 + h
 
h
+ 3 6
 
− h1 + h
6
1
h
+ h
3

Dans la matrice globale, elle est vue :

0 ··· 0
 
0 0 0
0 1 + h − 1 + h 0 · · · 0
 h 3 h 6 
(1) 0 − 1 + h 1 + h 0 · · · 0
A = h 6 h 3 
. .. .. .. .. 
 .. . . . .
0 0 0 0 ··· 0
Ainsi de suite et après assemblage, on obtient la matrice globale A :
 
1
h
+ h3 − h1 + h6 0 ··· ··· 0
 1 h 2 2h
− h + 6 h + 3 − h1 + h6 · · · ··· 0 

− h1 + h6 h2 + 2h ··· ···
 
XN  0 3
0 
A= A(k) =  .. .. .. .. 
 
 . . . . 
k=0
 ..
 
 . 2
h
2h 1
+ 3 −h + 6 h

0 ··· ··· 0 − h + h6 h1 + h3
1

Exercice 7 Soit le problème aux limites :

− (ν(x)u0 (x))0 + c(x) u(x) = f (x) dans ]0, 1[

(P )
u0 (0) = α ; u0 (1) = β
où ν , c et f des fonctions continues sur [0, 1] et ν ≥ 0
1. Écrire la formulation faible (Q) du problème (P ).
2. Montrer que le problème (Q) admet une solution unique u dans H 1 ([a, b]).
3. Écrire la formulation faible discrète (Qh ), éléments nis p1 .
4. Écrire le système linéaire obtenu.
5. Proposer une méthode de quadrature pour le calcul du second membre.

44
2.7 Approximation éléments nis P2
Reprenons l'exemple de l'équation de poisson (1D) avec CL de Dirichlet :
−u00 = f dans ]a, b[

(P )
u(a) = u(b) = 0
Nous allons choisir une approximation d'ordre plus élevé que celle P1 .
Considérons alors le sous-espace d'approximation Vh :
Vh = vh ∈ C 0 ([a, b])/vh |[xi ,xi+1 ] ∈ P2 ∀ i = 0, . . . , N et vh (a) = vh (b) = 0


P2 étant l'espace des polynômes de degrés ≤ 2.

Dans la MEF P1 , les fonctions vh étaient supposées anes sur chaque [xi , xi+1 ], et donc
il susait de connaître les valeurs en xi et xi+1 pour déterminer vh sur [xi , xi+1 ].

Dans la MEF P2 , pour déterminer vh sur [xi , xi+1 ] on doit avoir la valeur de vh en 3 points
de l'intervalle [xi , xi+1 ].
On introduit alors un point milieu de [xi , xi+1 ] noté xi+1/2 :
xi + xi+1 h
xi+1/2 = = xi + ∀ i = 0, . . . , N
2 2

Dénition 2.7.1 Le triplet formé par ([xi , xi+1 ], P2 , {xi , xi+1/2 , xi+1 }) est appelé élément
ni.
[xi , xi+1 ] : l'élément.
P2 : le type d'interpolation.
{xi , xi+1/2 , xi+1 } : les n÷uds constituant les degrés de liberté.

Par exemple, dans le cas de la MEF P1 , l'élément ni est :

([xi , xi+1 ], P1 , {xi , xi+1 })

Fonctions de base
Comme pour la MEF P1 , nous allons choisir des fonctions de base (ϕ1/2 , ϕ1 , ϕ3/2 , ϕ2 , . . . , ϕN , ϕN +1/2 )
telles que :

 ϕi (xj ) = δij
ϕi+1/2 (xj+1/2 ) = δij
ϕi , ϕi+1/2 ∈ C 0 ([a, b])


x 0=a x 1/ 2 x1 x 3 /2 x2 xN x N +1 / 2 x N +1 =b

Pour ce choix de fonctions de base, on parle d'éléments nis de Lagrange, et on dira que
l'élément ni ([xi , xi+1 ], P2 , {xi , xi+1/2 , xi+1 }) est un élément ni de Lagrange.

On peut montrer aussi le résultat suivant :

45
Théorème 2.7.1
dim Vh = N + N + 1 = 2N + 1
et la base éléments nis de Vh est donnée par :
(ϕ1/2 , ϕ1 , ϕ3/2 , ϕ2 , . . . , ϕN , ϕN +1/2 )

avec 
2 (x − xi−1 )(x − xi−1/2 )
si


 xi−1 ≤ x ≤ xi
h2


ϕi (x) = 2 (xi+1 − x)(xi+1/2 − x)
 si xi ≤ x ≤ xi+1
 h2
ailleurs


 0

et 
 4 (xi+1 − x)(x − xi )
ϕi+1/2 (x) = si xi ≤ x ≤ xi+1
h2
 0 ailleurs

Représentation des fonctions de base éléments nis P2

1 1

x i−1 x i−1/ 2 xi x i+1 /2 x i+1 x i−1 x i−1/ 2 xi x i+1 /2 x i+1

Le problème variationnel approché s'écrit alors :

Trouver uh ∈ Vh /

(Qh )
a(uh , vh ) = L(vh ) ∀ vh ∈ V h
avec
Z b
a(uh , vh ) = u0h (x) vh0 (x) dx
a
Z b
L(vh ) = f (x) vh (x) dx
a

On décompose uh dans la base de Vh :

N
X N
X
uh (x) = ui ϕi (x) + ui+1/2 ϕi+1/2 (x)
i=1 i=0
2N
X +1
= ui/2 ϕi/2 (x)
i=1

46
En exprimant le fait que le problème variationnel discret est vrai si et seulement si il l'est
pour toutes les fonctions de base (ϕj/2 ) 0≤j≤2N +1 , ce problème discret s'écrit :

 Trouver ui/2 (i = 1, . . . , 2N + 1) tels que :



2N
X +1

 ui/2 a(ϕi/2 , ϕj/2 ) = L(ϕj/2 ) ∀ j = 1, . . . , 2N + 1
i=1

qui représente un système linéaire (2N + 1) × (2N + 1) :

AU = b

où le vecteur inconnu (des degrés de liberté) est :

U = (u1/2 , u1 , u3/2 , . . . , uN +1/2 )

La matrice de rigidité est :

A = (aij ) 1≤i,j≤2N +1 avec aij = a(ϕi/2 , ϕj/2 )

Z b
= ϕ0i/2 (x) ϕ0j/2 (x) dx
a

Le second membre est b = (bj ) 1≤j≤2N +1 :

Z b
bj = L(ϕj/2 ) = ϕj/2 (x) f (x) dx
a

Exercice 8 Vérier que la matrice de rigidité s'écrit :

16
− 83
 
3
0 ··· ··· ··· ··· ··· 0
− 8 14
− 83 1
··· ··· ··· ··· 0 
 3 3 3 
0 − 83 16 8
−3 0 · · · · · · · · · 0
 
 3 
1
− 83 14
− 83 13 0 ···
 
 0 3 3
0 
1 .. ... ... ... ... ...

A=  .
 
h


 
− 83 16
− 83
 
 0 3
0 
..
 
. 1
− 83 14 8
−3

3 3

0 ··· ··· 0 − 38 16
3

Indication pour le calcul de A

Fonctions de base : (ϕ1/2 , ϕ1 , ϕ3/2 , . . . , ϕN +1/2 )

Élément 0 Élément N

x 0=a x 1/ 2 x1 x 3 /2 x2 xN x N +1 / 2 x N +1 =b

47
Calcul de la matrice de rigidité A par assemblage
Z b
A = (aij ) 1≤i,j≤2N +1 avec aij = ϕ0i/2 (x) ϕ0j/2 (x) dx
a

Sur l'élément 0 : [x0 , x1 ] :

ρ1 / 2 ρ1

x0 x 1/ 2 x1
 
(0) (0)
a11 a12
on calcule 
 

(0) (0)
a21 a22
Z x1 Z x1
(0) (0)
a11 = ϕ01/2 (x) ϕ01/2 (x) dx ; a12 = ϕ01/2 (x) ϕ01 (x) dx
Zx0x1 Zx0x1
(0) (0)
a21 = ϕ01 (x) ϕ01/2 (x) dx ; a22 = ϕ01 (x) ϕ01 (x) dx
x0 x0

Dans la matrice globale, on a :

 (0) (0) 
a11 a12 0 · · · 0
 
 
 (0) (0)
a21 a22 0 · · · 0

 
A(0) = 
 .. 
.
 0 

 .
 ..


0 ··· 0

Sur l'élément 1 : [x1 , x2 ] :

ρ1 ρ3/ 2 ρ2

x1 x 3 /2 x2

48
  (1) (1) (1)

a22 a23 a24
 
 
on calcule a32 a33 a34 
 (1) (1) (1) 
 
 
(1) (1) (1)
a42 a43 a44
Z x2 Z x2
(1) 0 0 (1)
a22 = ϕ1 (x) ϕ1 (x) dx ; a23 = ϕ01 (x) ϕ03/2 (x) dx
x1 x1
..
.
Z x2 Z x2
(1) (1)
a42 = ϕ02 (x) ϕ01 (x) dx ; a43 = ϕ02 (x) ϕ03/2 (x) dx
x1 x1

Dans la matrice globale, on a :

 
0 0 0 0 ··· 0
 
 
(1) (1) (1)
0 a
 22 a 23 a 24 0 · · · 0 

 
 
 (1) (1) (1) 
A(1) 0 a32 a33 a34 0 · · · 0
= 
 
 (1) (1) (1)

0 a
42 a 43 a 44 0 · · · 0 
. .
 
 .. .. 


0 ··· 0

Ainsi de suite, on calcule A(2) , A(3) , . . . , A(N ) .

La matrice de rigidité globale sera :
N
X
A= A(k)
k=0

En récapitulatif
La résolution d'une EDP ou EDO par la méthode des éléments nis consiste à :
1. Écrire la formulation variationnelle continue.
2. Choisir un maillage du domaine de calcul, puis dénir un espace d'approximation
Vh selon le type d'éléments nis (P1 ou P2 ).
3. Écrire le problème variationnel approché =⇒ Système linéaire A U = b.
4. Calcul de la matrice A et le second membre b du système linéaire.
5. Résolution du système linéaire A U = b.
Le vecteur inconnu U est constitué par les degrés de liberté.
6. Représenter graphiquement la fonction U c'est-à-dire les ui en fonction des xi .

49
2.8 Méthode des Eléments Finis en dimension 2
Soit Ω un domaine régulier en 2D.
On note par ∂Ω la frontière de Ω et n la normale unitaire extérieure à ∂Ω.

Ω

On veut utiliser la méthode des éléments nis pour approximer la solution du problème
aux limites suivant :

−∆u(x, y) + c(x, y) u(x, y) = f (x, y) pour (x, y) ∈ Ω̊
(P )
u(x, y) = 0 pour (x, y) ∈ ∂Ω

On suppose que f ∈ L2 (Ω) et c est une fonction continue et positive sur Ω.

2.8.1 Formulation variationnelle

On multiplie la première équation de (P ) par une fonction v assez régulière et on intègre
sur Ω :
ZZ ZZ ZZ
−∆u v dxdy + c u v dxdy = f v dxdy
Ω Ω Ω

En appliquant la formule de Green, on obtient :

ZZ Z ZZ ZZ
∂u
∇u · ∇v dxdy − v dσ + c u v dxdy = f v dxdy
Ω ∂n
| ∂Ω {z Ω Ω
car on prend v=0 sur ∂Ω
}
=0

Ces intégrales ont un sens dès lors que l'on prend u, v ∈ H01 (Ω). Le problème variationnel
s'écrit alors :

Trouver u ∈ H01 (Ω) telle que


(Q)
a(u, v) = L(v) ∀ v ∈ H01 (Ω)
avec :  ZZ ZZ
 a(u, v) = ∇u · ∇v dxdy + c u v dxdy


Z ZΩ Ω

 L(v) = f v dxdy


Ω

50
2.8.2 Formulation variationnelle discrète
Maillage de Ω :
An de construire un espace d'approximation de dimension nie, on commence tout
d'abord par réaliser un maillage de Ω à l'aide de triangles Ki . On peut aussi utiliser des
rectangles ou même des éléments mixtes "triangles-rectangles". Ces éléments Ki doivent
constituer un "recouvrement" de Ω, c'est-à-dire
Ne
[
Ω̄ = Ki
i=1

où Ne est le nombre d'éléments du maillage.

Dans le cas où les éléments Ki sont des triangles, on parle de triangulation Th de Ω.
Exemples de Maillages :

i,j+1
i,j+1
i-1,j i+1,j i-1,j i,j i+1,j
i,j

i,j-1 i,j-1

Figure 2.1  Exemple de maillages structurés.

2 8 Matrice de connexion
Eléments Noeuds
4
7 3 1 4 3 6
6
5 2 1 6 5
1
6 3 7 1 5
2
5
. . . .
3 4 . . . .
7 1

Figure 2.2  Exemple de maillages non structurés.

Hypothèses sur le maillage :

(H1 ) : On impose au maillage d'être "conforme" ou "admissible", c.a.d que l'intersection
entre deux éléments quelconques Ki et Kj de Th est :
 soit vide,
 soit réduite à un point,
 soit réduite à une arête commune à Ki et Kj .

51
 Figure 2.3  Maillages conformes (à gauche), maillages non conformes (à droite).

(H2 ) : Les éléments du maillage doivent être susament réguliers. S'il s'agit de triangles
par exemple, cela signie qu'ils doivent être les plus proches de l'équilatéralité ou du moins
ne pas avoir d'angles trop grands. Rappelons que la plupart des schémas numériques dé-
diés pour le traitement de la diusion dans les équations de Navier-Stokes nécessitent des
conditions assez restrictives sur les angles des triangles.

Méthodes de génération de Maillage :

En général, on peut classier les algorithmes de génération de maillages en 2D (tiangu-
laires) et 3D (tétraédriques) en trois grandes familles :

1. Les méthodes dites "frontales" ou par avancement de fronts.

2. Les méthodes de décomposition spatiale.

3. Les algorithmes de type Delaunay-Voronoi.

2.8.3 Espace d'approximation Vh (Eléments Finis P1)

Espace polynomial P1 en dimension 2
En dimension 2, l'espace P1 est déni par :

P1 = p : R 2 → R ; p(x, y) = α + βx + γy .


P1 est en espace vectoriel de dimersion 3.

Proposition
Soit K un triangle non dégénéré de sommets A1 , A2 , A3 (c.a.d A1 , A2 et A3 ne sont pas
alignés).
Il existe un unique polynôme p ∈ P1 prenant des valeurs xées aux sommets de K .

Preuve :
A1(x1 ,y1)

A2(x2 ,y2)
A3(x3 ,y3)

52
Notons par f1 , f2 , f3 les valeurs prises par p anx n÷uds A1 , A2 , A3 .
 
 p(x1 , y1 ) = f1  α + βx1 + γy1 = f1
p(x2 , y2 ) = f2 ⇐⇒ α + βx2 + γy2 = f2
p(x3 , y3 ) = f3 α + βx3 + γy3 = f3
 

Le déterminant du système est :

1 x1 y 1
∆= 1 x2 y 2 = det(K)
1 x3 y 3

Or, la surface du triangle K est dénie par :

| det(K)|
S(K) =
2
et comme le triangle K est supposé non dégénéré alors

S(K) 6= 0
Donc ∆ 6= 0

Par conséquent, le système admet une solution unique α, β, γ . D'ou l'existence et l'unicité
du polynôme p ∈ P1 .

Espace d'approximation Vh :
On considere Vh comme suit (E.F P1 ) :
n o
Vh = vh ∈ C (Ω)/ vh|Ki ∈ P1 ; ∀ i = 1, . . . , Ne et vh = 0 sur ∂Ω
0

Base de Vh :
Pour construire une base de Vh , nous allons procéder comme dans le cas 1D.
Notons par :
 NS le nombre de sommets du maillage.
 NSint le nombre de sommets internes du maillage (c.a.d non situés sur ∂Ω).

∀ i = 1, . . . , NSint , notons par la fonction ϕi :

ϕi (Aj ) = ϕi (xj , yj ) = δij

c.a.d que la fonction ϕi vaut 1 aU sommet Ai et 0 sur les autres sommets.

Proposition :
La famille (ϕ1 , ϕ2 , . . . , ϕNSint ) est une base de Vh .

53
 1

ji

Ai

Figure 2.4  Fonction de base ϕi en 2D.

Ecriture du système linéaire :

Là aussi, en décomposant uh dans la base (ϕi )1≤i≤NSint :
N
X Sint

uh (x, y) = ui ϕi (x, y)
i=1

Le problème variationnel approché s'écrit :

 Trouver u1 , u2 , . . . , uNSint /


NSint
(Qh ) X

 ui a (ϕi , ϕj ) = L (ϕj ) ∀ j = 1, . . . , NSint
i=1

ou encore : AU = b
avec :

• A = (aij )1≤i,j≤NSint
ZZ ZZ
aij = a (ϕi , ϕj ) = ∇ϕi · ∇ϕj dxdy + c ϕi ϕj dxdy
Ω Ω
X Z Z ZZ 
= ∇ϕi · ∇ϕj dxdy + c ϕi ϕj dxdy
K∈Th K K

• b = (bj )1≤j≤NSint
ZZ
bj = f ϕj dxdy
Ω ZZ
X
= f ϕi dxdy
K∈Th K

Les intégrales intervenants dans le second membre b sont généralement évaluées numéri-
quement par des formules de quadratures, par exemple celles de Gauss.

54