2019/2020

Master 2 Recherche fondamentale en mathématiques

Méthode des éléments finis

Professeur : Nicolas Seguin

Rédigé par Pierre Le Barbenchon

Table des matières

1 Introduction 1

2 Motivation et prérequis 2

3 Méthode de Lagrange pour des polynômes de degré au plus 1 3

3.1 Théorie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
3.2 Implémentation numérique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3.3 Analyse numérique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.3.1 En théorie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.3.2 En pratique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

4 Méthode de Lagrange pour des polynômes de degré au plus 2 11

4.1 Théorie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
4.2 Implémentation numérique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
4.3 Analyse numérique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4.3.1 En théorie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4.3.2 En pratique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

5 Annexe 16
5.1 Code pour le cas k = 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
5.2 Code pour l’erreur dans le cas k=1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
5.3 Code pour le cas k = 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
5.4 Code pour l’erreur dans le cas k=2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
5.5 Preuve du théorème 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.6 Résultats plus généraux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

1 Introduction
On veut étudier la méthode des éléments finis en dimension 1 pour l’équation
(
−u00 = f sur [0, 10]
u(0) = u(10) = 0

On commencera, en partie 2, par des aspects théoriques sur la méthode des éléments finis qui
motiveront cette étude.

Ensuite, on traitera dans deux grandes parties 3 et 4, les méthodes liées aux polynômes de
Lagrange de degrés 1 et 2 sur des maillages uniforme et non uniforme. L’exemple de maillage non
uniforme qu’on donnera est donné par

[x*b/(n+1)* sinh(x/n)/sinh((n+1)/n) for x in range (n+2)]

Dans chacune des deux parties 3 et 4, on aura une implémentation de la méthode avec, pour
fonction f , la fonction x 7→ sin(πx) sur [0, 10]. Puis pour ce problème, on étudiera l’ordre de
convergences des méthodes pour les normes k · kL2 et k · kH 1 (que l’on abrègera en norme L2 et
norme H 1 ) pour les deux maillages (uniforme et non uniforme).

1
 Enfin, en annexe 5, on donnera le code qui permet d’effectuer les simulations, ainsi que la preuve
du théorème 2.

L’implémentation se fera dans le langage Python, c’est pourquoi les indices du code commencent
à 0.

La plupart des résultats que l’on présente ici sont des restrictions de théorèmes plus généraux
(par exemple, le lemme de Céa, les majorations d’erreur), mon objectif étant de bien comprendre
les résultats dans notre cadre d’application (Lagrange 1 et 2, pour les normes L2 et H 1 avec des
maillages uniforme et non uniforme). On présente en annexe l’énoncé général de certains résultats
(voir 5.6).

2 Motivation et prérequis
Soit I = [0, 10] et l’équation
(
−u00 = f sur I
(1)
u(0) = u(10) = 0

On pose une formule variationnelle associée au problème (1) :

(
Trouver u ∈ H01 (I) tel que
(2)
∀v ∈ H01 (I), a(u, v) = b(v)
Z Z
où a(u, v) = u0 v 0 et b(v) = f v.
I I

On veut utiliser une méthode de Galerkin, en se donnant un sous espace de dimension finie Vh
(pour h > 0) de V = H01 (I) qui tend vers V quand h tend vers 0 dans un certain sens.
On veut donc résoudre le problème suivant :
(
Trouver uh ∈ Vh tel que
(3)
∀vh ∈ Vh , a(uh , vh ) = b(vh )

Soit N la dimension de l’espace Vh , on se donne une base (φi )N

i=1 de l’espace Vh .
Pour résoudre l’équation (3), on procède par analyse-synthèse. Pour l’analyse, on prend une
solution uh de (3), on peut alors écrire
N
X
uh = uj φj
j=1

pour (uj )N N
j=1 ∈ R .
N
X
On a donc a(uh , φi ) = uj a(φj , φi ) = b(φi ) pour tout i ∈ {1, . . . , N }, ce qui peut se réécrire
j=1
en
Ah Uh = Fh
avec (Ah )i,j = a(φj , φi ), (Uh )i = ui et (Fh )i = b(φi ).

2
 Dans notre cas, pour résoudre (1), on a
Z
(Ah )i,j = φ0i (x)φ0j (x)dx
ZI
(Fh )i = f (x)φi (x)dx
I

Pour la synthèse, en calculant Uh = A−1

h Fh , on trouve aussi que uh , associé à Uh , est bien
solution de (3).

Tout le but des parties suivantes sera de résoudre cette équation avec des Vh différents (et donc
des bases (φi ) différentes) et d’étudier les propriétés de convergence entre la solution uh de (3) et
la solution u de (2) quand h tend vers 0.

3 Méthode de Lagrange pour des polynômes de degré au plus 1

3.1 Théorie
Soit I l’intervalle [0, 10].
Définition 1. Soit N ∈ N. On définit un maillage de I relatif aux N points distincts (xi )N
i=1 ∈
N
]0, 10[ par l’ensemble
τh = {[xi , xi+1 ], ∀i ∈ {0, · · · , N }} ,
où x0 = 0 et xN +1 = 10.
Remarque :
On dira qu’un maillage est uniforme si les points (xi )N +1
i=0 sont uniformément répartis dans I.

Définition 2. Soit τh un maillage, on appelle cellule un élément de τh .

On se fixe un maillage τh de I.
Définition 3. L’espace d’approximation Vh,1 que l’on va utiliser ici sera
n o
Vh,1 = v ∈ C 0 (I), v|[xj ,xj+1 ] ∈ P1 [X], j = 0, . . . , N avec v(0) = v(10) = 0

où P1 [X] désigne les polynômes réels de degré au plus 1.

Proposition 4. L’espace Vh,1 est de dimension N .
Exemple 5. La fonction v dessinée ci-dessous est dans Vh,1 . (pour N = 7)

x0 x1 x2 x3 x4 x5 x6 x7 x8

3
 Il faut se donner une base de Vh,1 .

Définition 6. On prend la base B = {φi ∈ Vh,1 , i ∈ {1, . . . , N }, où φi (xj ) = δi,j }.

Les fonctions φi sont de la forme suivante

1 φi

x0 x1 xi−1 xi xi+1 xN xN +1

On pose les fonctions φ0 et φN +1 de la même façon.

1 φ0

x0 x1 xN xN +1

1 φN +1

x0 x1 xN xN +1

Comme vu dans la partie 2, on pose alors les matrices

Z
(Ah )i,j = φ0i (x)φ0j (x)dx pour tout i, j ∈ {1, . . . , N }
I Z

(Fh )i = f (x)φi (x)dx pour tout i ∈ {1, . . . , N }

I

La solution est alors Uh = A−1 h Fh .

Pour calculer les matrices Ah et Fh de manière pratique, on se donner une cellule de référence
K̂ = [0, 1] sur laquelle on se donne deux fonctions de base ϕ0 et ϕ1 .

1 1

ϕ0 ϕ1

0 1 0 1

4
 Ainsi, pour une cellule K = [α, β] de τ , on peut définir la transformation affine
(
K̂ → K
FK : .
x̂ 7→ x̂β + (1 − x̂)α

Cette transformation envoie K̂ sur K.

Ainsi, pour K = [xj , xj+1 ], on note φK,0 = φj et φK,1 = φj+1 , on a alors, pour l, m ∈ {0, 1},
Z Z
−1 0 −1 0
φ0K,l (x)φ0K,m (x)dx = (ϕl ◦ FK ) (x)(ϕm ◦ FK ) (x)dx
K ZK
−1 0 −1 −1 0 −1
= (FK ) (x)ϕ0l (FK (x))(FK ) (x)ϕ0m (FK (x))dx
ZK
1 −1 1 −1
= 0 −1 ϕ0l (FK (x)) 0 −1 ϕ0m (FK (x))dx
K FK (FK (x)) FK (FK (x))
Z
1
= ϕ0 (F −1 (x))ϕ0m (FK −1
(x))dx
(xj+1 − xj )2 K l K
Z
1
= ϕ0 (x̂)ϕ0m (x̂)(xj+1 − xj )dx̂
(xj+1 − xj )2 K̂ l
Z
1
= ϕ0 (x̂)ϕ0m (x̂)dx̂
xj+1 − xj K̂ l

On va donc pouvoir se servir de cette formule sur la cellule de référence pour calculer toutes les
intégrales de la matrices Ah . De même, on peut faire le même type de calcul pour calculer Fh .

3.2 Implémentation numérique

On veut résoudre l’équation
(
−u00 (x) = sin(πx) sur [0, 10]
(4)
u(0) = u(10) = 0

La solution exacte est

sin(πx)
∀x ∈ [0, 10], u(x) =
π2
J’ai choisi d’implémenter le maillage par une liste qui contient les N + 2 points rangés par ordre
croissant (donc les 2 bords sont aux extrémités de la liste).
Pour un maillage uniforme, on aura

Nodes_1_uni = [ b*i/(n+1) + (1-i/(n+1))*a for i in range(n+2)]

x0 x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16

Exemple pour N = 15

On peut se donner un maillage non uniforme comme

Nodes_1 = [x*b/(n+1)* sinh(x/n)/sinh((n+1)/n) for x in range (n+2)]

5
 x0 x5x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16
Exemple pour N = 15

On représente ensuite les cellules par une liste de tableaux de la forme [i, i + 1] où i varie entre
0 et N .

Cells_1 = [Link]([[i,i+1] for i in range (n+1)] )

Pour la matrice Ah , j’ai choisi de créer une matrice A de taille (N + 2) × (N + 2) pour traiter
les cas aux bords, puis d’extraire la sous matrice Ah de taille N × N en enlevant
Z les bords de A.
J’ai choisi d’utiliser la méthode de Simpson pour approcher l’intégrale ϕ0l (x̂)ϕ0m (x̂)dx̂.
K̂
On a alors
Z
1 0
ϕ0l (x̂)ϕ0m (x̂)dx̂ ' ϕl (0)ϕ0m (0) + 4ϕ0l (1/2)ϕ0m (1/2) + ϕ0l (1)ϕ0m (1)


K̂ 6

On peut donc calculer la matrice A en faisant une boucle

Z sur les cellules, puis en faisant varier
1
l et m dans l’ensemble {0, 1} et en ajoutant ϕ (x̂)ϕ0m (x̂)dx̂ à Al,m .
0
xj+1 − xj K̂ l
De même, on calcule le vecteur colonne F de taille N + 2, pour ensuite extraire Fh qui sont les
N coordonnées de F en retirant les bords.
Enfin on résout le système Ah Uh = Fh pour obtenir notre solution à l’aide la bibliothèque numpy
et de la commande Uh = [Link](Ah, Fh).

On présente maintenant les résultats pour les deux maillages présentés ci-dessus, pour N ∈
{20, 60, 200}.

Maillage uniforme Maillage non uniforme

N = 20 N = 20

6
 N = 60 N = 60

N = 200 N = 200

La méthode semble donc converger. Dans la partie suivante, on va quantifier cette convergence
et calculer l’ordre de la convergence pour les normes L2 et H 1 .

3.3 Analyse numérique

On veut donner une estimation de l’erreur de notre approximation.

3.3.1 En théorie

Lemme 1 (de Céa). Pour u la solution de (1), uh solution de (3) et Πh (u) la projection de u sur
l’espace Vh,1 , on a
M
ku − uh kH 1 6 ku − Πh (u)kH 1
α
où M est la constante de continuité de a, α est la constante de coercivité de a avec a la forme
bilinéaire de (2).

Démonstration. On a uh − Πh (u) ∈ Vh,1 par définition de uh et Πh (u).

Donc on a
a(u, uh − Πh (u)) = b(uh − Πh (u)) et a(uh , uh − Πh (u)) = b(uh − Πh (u)),

7
puisque u et uh sont solutions de (1) et (3).
En soustrayant, on a donc
a(u − uh , uh − Πh (u)) = 0.
Par coercivité de a, on a
αku − uh k2H 1 6 a(u − uh , u − uh ) + a(u − uh , uh − Πh (u))
| {z }
=0
6 a(u − uh , u − Πh (u))
6 M ku − uh kH 1 ku − Πh (u)kH 1
Ainsi on a
M
ku − uh kH 1 6 ku − Πh (u)kH 1 .
α

Théorème 1. Pour u solution H 2 (I) de (1) et Πh (u) la projection de u sur l’espace Vh,1 , on a

ku − Πh (u)kH 1 (I) 6 Chku00 kL2 (I)

où h = max |xi+1 − xi |.

[xi ,xi+1 ]∈τh

Démonstration.
Étape 1 : Montrons que ku0 − Πh (u)0 kL2 6 hku00 kL2 .
On va raisonner par densité, on suppose ici que u est C ∞ .
XN Z xi+1
0 0 2
ku − Πh (u) kL2 = (u0 (x) − Πh (u)0 (x))2 dx (5)
i=0 xi

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

u(xi+1 ) − u(xi )
Πh (u)0 (x) = = u0 (yi )
hi
avec hi = xi+1 − xi et yi ∈]xi , xi+1 [ défini par le théorème des accroissements finis.
Z xi+1 Z xi+1
0 0
(u (x) − Πh (u) (x)) dx = 2
(u0 (x) − u0 (yi ))2 dx
xi xi
Z xi+1 Z x 2
00
= u (t)dt dx
xi yi
Z xi+1 Z x  Z x 
00 2
6 u (t) dt 1dt dx par Cauchy–Schwarz
xi yi yi
Z xi+1 Z x 
00 2
6 hi u (t) dt dx
xi yi
Z xi+1 Z xi+1 
00 2
6 hi u (t) dt dx
xi xi
Z xi+1 Z xi+1 
6 hi u00 (t)2 dx dt par Fubini
xi xi
Z xi+1
2
6 hi u00 (t)2 dt (6)
xi

8
Ainsi, en injectant dans l’égalité (5), on a
 2
0
ku − Πh (u)0 k2L2 6 max hi ku00 k2L2 .
i∈J0,N K

En passant à la racine, on a
ku0 − Πh (u)0 kL2 6 hku00 kL2 .

Étape 2 : Montrons que ku − Πh (u)kL2 6 h2 ku00 kL2 .

N Z
X xi+1
ku − Πh (u)k2L2 = |u(x) − Πh (u)(x)|2 dx
i=0 xi

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

Z x 2
2 0 0
|u(x) − Πh (u)(x)| = (u (t) − Πh (u) (t))dt car u(xi ) = Πh (u)(xi )
xi
Z x
6 hi (u0 (t) − Πh (u)0 (t))2 dt par Cauchy–Schwarz
x
Z ixi+1
6 h3i u00 (t)2 dt par l’inégalité (6)
xi

On intègre sur [xi , xi+1 ] pour obtenir

Z xi+1 Z xi+1
2 4
|u(x) − Πh (u)(x)| dx 6 hi u00 (t)2 dt.
xi xi

Ainsi
ku − Πh (u)k2L2 6 h4 ku00 k2L2 .
En passant à la racine carrée, on obtient

ku − Πh (u)kL2 6 h2 ku00 kL2 .

Étape 3 : Conclusion.

ku − Πh (u)kH 1 = ku − Πh (u)kL2 + ku0 − Πh (u)0 kL2

6 h2 ku00 kL2 + hku00 kL2
6 (1 + h)hku00 kL2
6 Chku00 kL2 avec C = |I| + 1

Corollaire 1. Dans notre cadre (k = 1), on a

ku − uh kL2 (I) 6 Ch2

ku − uh kH 1 (I) 6 Ch

où C dépend uniquement du domaine I et de la fonction f .

9
Démonstration. On a, par le lemme de Céa et par le fait que −u00 = f , les inégalités suivantes :
M
ku − uh kL2 (I) 6 ku − Πh (u)kL2 (I)
α
6 C 0 h2 ku00 kL2 en utilisant la preuve du théorème 1
6 C 0 h2 kf kL2
6 C 00 h2 avec C 00 = C 0 kf kL2
et
M
ku − uh kH 1 (I) 6 ku − Πh (u)kH 1 (I)
α
MC
6 C 0 hku00 kL2 avec C 0 =
α
6 C 0 hkf kL2
6 C 00 h avec C 00 = C 0 kf kL2

3.3.2 En pratique
On va maintenant calculer numériquement les erreurs entre la solution exacte et la solution
approchée en normes L2 et H 1 . Pour cela, on calcule la norme L2 de u − uh 1 en utilisant la fonction
quad du module [Link] de Python. On va donner la fonction uh explicitement grâce à la
fonction uh_k1 (voir annexe 5.2). De plus, pour calculer la norme H 1 , il faut calculer u0h qui sera une
fonction constante par morceaux qui vaut sur chaque [xi , xi+1 ] du maillage le coefficient directeur
de la pente entre uh (xi ) et uh (xi+1 ). On donne cette fonction explicitement grâce à uhprim_k1 (voir
annexe 5.2).

Ensuite pour calculer l’erreur, on va tracer la fonction log(ku − uh k) en fonction de log h. On va

trouver une pente qui aura pour coefficient directeur l’exposant de h dans les inégalités que l’on a
démontré dans la partie 3.3.1. On va tracer cette courbe pour N ∈ {100, 200, 500, 800, 1000, 3000}.

On va de nouveau faire une distinction entre les cas où le maillage sera uniforme ou non uniforme.

Les tracés sont de la forme :

Erreur en norme H 1 (I) pour le maillage non uniforme.

1. où u est la solution exacte et uh la solution approchée calculée par la méthode des éléments finis

10
 On résume les résultats en donnant le coefficient directeur moyen des pentes pour les différentes
simulations.

Coefficient directeur Maillage uniforme Maillage non uniforme Théorique

Norme L2 1.9827103 2.0032377 2
Norme H 1 1.0016067 1.0153032 1

Ainsi, la simulation numérique coı̈ncide bien avec les résultats théoriques.

4 Méthode de Lagrange pour des polynômes de degré au plus 2

4.1 Théorie
On va maintenant changer d’espace d’approximation, en élevant le degré des polynômes que l’on
considère.
On se fixe un maillage τh de I.

Définition 7. L’espace d’approximation Vh,2 que l’on va utiliser ici sera

n o
Vh,2 = v ∈ C 0 (I), v|[xj ,xj+1 ] ∈ P2 [X], j = 0, . . . , N avec v(0) = v(10) = 0

où P2 [X] désigne les polynômes réels de degré au plus 2.

Proposition 8. L’espace Vh,2 est de dimension 2N + 1.

Exemple 9. La fonction v dessinée ci-dessous est dans Vh,2 . (pour N = 7)

x0 x1 x2 x3 x4 x5 x6 x7 x8

Il faut se donner une base de Vh,2 . Pour chaque cellule K = [xj , xj+1 ], on va rajouter un point
au milieu xej , afin d’avoir trois points sur chaque cellule. Ainsi, comme il existe un unique polynôme
de degré au plus deux qui passe par trois points distincts, on va pouvoir trouver 3 polynômes pour
chaque cellule tels que chaque polynôme vaut 1 sur un des trois points et 0 sur les deux autres. On
décrit les trois fonctions ϕ0 , ϕ01 et ϕ1 pour la cellule de référence K̂ = [0, 1].

11
 1 1 1
ϕ01

ϕ0 ϕ1

0 1/2 1 0 1/2 1 0 1/2 1

On peut ensuite utiliser les mêmes calculs que pour la méthode de Lagrange avec des polynômes
de degré au plus 1, il faut juste faire varier m, l dans {0, 1, 2} dans les intégrales que l’on considère.

Si on veut représenter la base de Vh,2 que l’on s’est donnée, on tracerait les 2N + 1 fonctions
suivantes (car Vh,2 est de dimension 2N + 1).
On a les N fonctions (φi )Ni=1 :

1 φi

x0 x1 xi−1 xi xi+1 xN xN +1

Puis on a les N + 1 fonctions (φei )N

i=0 :

1 φei

x0 x1 xi−1 xi xi+1 xN xN +1

4.2 Implémentation numérique

On veut toujours résoudre l’équation
(
−u00 (x) = sin(πx) sur [0, 10]
(7)
u(0) = u(10) = 0

J’ai choisi d’implémenter le maillage par une liste qui contient les 2N + 3 points rangés telle que
les N points du maillage sont au début (rangés dans l’ordre croissant), puis on met les N + 1 points
milieux (par ordre croissant) et enfin les deux bords.
Pour un maillage uniforme, on aura

Nodes_2_uni = Nodes_1_uni[1:-1]
+ [(Nodes_1_uni[i] + Nodes_1_uni[i+1])/2 for i in range (n+1)] + [0,10]

x0 x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16

Exemple pour N = 15

12
Nodes = [1.66, 3.33, 5.0, 6.66, 8.33, 0.83, 2.5, 4.16, 5.83, 7.5, 9.16, 0, 10]
Exemple pour N = 5

On peut se donner un maillage non uniforme comme

Nodes_2 = Nodes_1[1:-1] + [(Nodes_1[i] + Nodes_1[i+1])/2 for i in range (n+1)]

+ [0,10]

x0 x5x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16

Exemple pour N = 15

Nodes = [0.22, 0.90, 2.10, 3.92, 6.48, 0.11, 0.56, 1.50, 3.01, 5.20, 8.24, 0, 10]
Exemple pour N = 5

On représente ensuite les cellules par une liste de tableaux de la forme [i, i + 1, i + n + 1] où i
varie entre 0 et N − 2. Ainsi on place le début de la cellule, la fin de la cellule puis le milieu de
la cellule. Les deux cellules du bords sont donc de la forme [31, 0, 15] et [14, 32, 30] dans
l’exemple où N = 15.

Cells_2 = [[2*n+1,0,n]] + [[i,i+1,i+n+1] for i in range(0,n-1)]+[[n-1,2*n+2,2*n]]

Je me suis donné un tableau Bords qui contient les deux bords :

Bords = [2*n+1, 2*n+2]

Pour la matrice Ah , j’ai choisi de créer une matrice Ah de taille N × N cette fois-ci et de gérer
les bords dans la boucle de calcul des coefficients.

On peut donc calculer la matrice Ah en faisant une boucle sur les cellules, puis en faisant varier
1
ϕ0 (x̂)ϕ0m (x̂)dx̂ à Al,m si les indices
R
l et m dans l’ensemble {0, 1, 2} et en ajoutant
xj+1 − xj K̂ l
correspondant à la cellule ne sont pas dans le tableau Bords.
De même, on calcule le vecteur colonne Fh de taille N en faisant attention aux bords.

Enfin on résout le système Ah Uh = Fh pour obtenir notre solution à l’aide la bibliothèque numpy
et de la commande Uh = [Link](Ah, Fh).

On présente maintenant les résultats pour les deux maillages présentés ci-dessus, pour N ∈
{20, 60, 200}.

13
 Maillage uniforme Maillage non uniforme

N = 20 N = 20

N = 60 N = 60

N = 200 N = 200
Comme dans la partie précédente, la méthode semble converger. Dans la section suivante, on va
quantifier cette convergence et calculer l’ordre de la convergence pour les normes L2 et H 1 .

14
4.3 Analyse numérique
On veut donner une estimation de l’erreur de notre approximation.

4.3.1 En théorie

Théorème 2. Pour u solution H 3 (I) de (1) et Πh (u) la projection de u sur l’espace Vh,2 , on a

ku − Πh (u)kH 1 (I) 6 Ch2 ku(3) kL2 (I)

où h = max |xi+1 − xi |.

[xi ,xi+1 ]∈τh

Démonstration. Semblable à la preuve du théorème 1. Pour la preuve détaillée, voir l’annexe 5.5.

Corollaire 2. Dans notre cadre (k = 2), on a

ku − uh kL2 (I) 6 Ch3

ku − uh kH 1 (I) 6 Ch2

où C dépend uniquement du domaine I et de la fonction f .

Démonstration. Semblable à la preuve du corollaire 1.

4.3.2 En pratique
On va maintenant calculer numériquement les erreurs entre la solution exacte et la solution
approchée en normes L2 et H 1 . Pour cela, on calcule la norme L2 de u − uh 2 en utilisant toujours la
fonction quad du module [Link] de Python. On va donner la fonction uh explicitement
grâce à la fonction uh_k2 (voir annexe 5.4) en interpolant entre les 3 points xi , xei et xi+1 . De plus,
pour calculer la norme H 1 , il faut calculer u0h qui sera une fonction affine par morceaux. On donne
cette fonction explicitement grâce à uhprim_k2 (voir annexe 5.4).

Ensuite pour calculer l’erreur, on va tracer la fonction log(ku − uh k) en fonction de log h. On va

trouver une pente qui aura pour coefficient directeur l’exposant de h dans les inégalités que l’on a
démontré dans la partie 3.3.1. On va tracer cette courbe pour N ∈ {100, 200, 500, 800, 1000, 3000}.

On va continuer de faire une distinction entre les cas où le maillage sera uniforme ou non
uniforme.

On résume les résultats en donnant le coefficient directeur moyen des pentes pour les différentes
simulations.

Coefficient directeur Maillage uniforme Maillage non uniforme Théorique

Norme L2 1.9922769 2.0006769 3
Norme H 1 2.0075953 2.0160268 2

2. où u est la solution exacte et uh la solution approchée calculée par la méthode des éléments finis

15
 On a donc une différence entre les résultats théoriques et les simulations pour la norme L2 .
Cependant, pour la norme H 1 , les simulations coı̈ncident avec les résultats théoriques.
La différence de résultats pour la norme L2 vient sans doute d’une erreur d’approximation
d’intégrale que l’on effectue dans la méthode des éléments finis. En effet, dans la partie théorique,
on calcule l’erreur pour la fonction uh qui est donnée par la résolution du système Ah Uh = Fh . Or,
dans la pratique, pour calculer Ah , on a utilisé une méthode de Simpson qui introduit une erreur
de calcul numérique. Cette erreur a été négligée dans la partie théorique.

5 Annexe
5.1 Code pour le cas k = 1
from math import *
import numpy as np
from [Link] import quad
import [Link] as plt

n=20
a=0
b=10

def f(x):
return sin(pi*x)

def solution_exacte(x):
return f(x)/(pi*pi)

Nodes_1 = [b*i/(n+1) + (1-i/(n+1))*a for i in range(n+2)]

#Nodes_1 = [x*b/(n+1)* sinh(x/n)/sinh((n+1)/n) for x in range (n+2)]

Cells_1 = [Link]([[i,i+1] for i in range (n+1)] )

def RefBasicFct(x, k):

if k == 1:
return [Link]([x, 1-x])
if k == 2:
return [Link]([2*(x-1/2)*(x-1), 2*(x-1/2)*x, -4*x*(x-1)])

def DerRefBasisFct(x, k):

if k == 1:
return [Link]([1,-1])
if k == 2:
return [Link]([4*x - 3, 4*x - 1, -8*x + 4])

16
#Calcul de A_h
A = [Link]((n+2,n+2))
for k in range (len(Cells_1)):
for i in range(2):
for j in range(2):
A[Cells_1[k][i]][Cells_1[k][j]] += 1/(abs(Nodes_1[Cells_1[k][1]]-
Nodes_1[Cells_1[k][0]]))*6*(DerRefBasisFct(0, 1)[i]*
DerRefBasisFct(0, 1)[j]+4*DerRefBasisFct(1/2, 1)[i]*
DerRefBasisFct(1/2, 1)[j]+DerRefBasisFct(1, 1)[i]*DerRefBasisFct(1, 1)[j])
A= [Link](A)
A = A[1:-1,1:-1]

#Calcul de F_h
F = [Link](n+2)
for k in range (len(Cells_1)):
for i in range(2):
F[Cells_1[k][i]] += (Nodes_1[Cells_1[k][1]]-Nodes_1[Cells_1[k][0]])*6*
(RefBasicFct(0,1)[i]*f(Nodes_1[Cells_1[k][1]])+4*RefBasicFct(1/2,1)[i]*
f(1/2*Nodes_1[Cells_1[k][0]]+(1/2)*Nodes_1[Cells_1[k][1]])+
RefBasicFct(1,1)[i]*f(Nodes_1[Cells_1[k][0]]))
F = [Link](F)
F = F[1:-1]

#solution approchée
U = [Link](A, F)

sol = [Link](n+2)
for i in range (1,n+1):
sol[i] = U[i-1]

#tracer les courbes

[Link]()
abscisse = [Link](a,b,200)
[Link](abscisse,[solution_exacte(x) for x in abscisse], label="solution_exacte")
[Link](Nodes_1,sol, label="solution_approchée")
[Link](loc=’upper left’)
[Link](-0.11, 0.15)
[Link]()

5.2 Code pour l’erreur dans le cas k = 1

def moy_coef(a,b):
S = 0
for i in range (len(a)-1):

17
 S+= ((b[i+1]-b[i])/(a[i+1]-a[i]))
return S/(len(a)-1)

def uh_k1(L, Nodes, x):

i = 0
while Nodes[i] < x:
i += 1
a = Nodes[i-1]
b = Nodes[i]
fa = L[i-1]
fb = L[i]
return fb/(b-a)*(x-a) + fa/(a-b)*(x-b)

#resol_k1 donne le vecteur U_h

def erreur_k1_L2(n):
Nodes_1 = [b*i/(n+1) + (1-i/(n+1))*a for i in range(n+2)]
# Nodes_1 = [x*b/(n+1)* sinh(x/n)/sinh((n+1)/n) for x in range (n+2)]
sa = resol_k1(n)
def g(x):
return uh_k1(sa,Nodes_1,x)
return sqrt(quad(lambda x : (g(x)-f(x)/(pi*pi))**2,a,b)[0])

def uhprim_k1(L, Nodes, x):

i = 0
while Nodes[i] < x:
i += 1
a = Nodes[i-1]
b = Nodes[i]
fa = L[i-1]
fb = L[i]
return (fb - fa) / (b - a)

def erreur_k1_H1(n):
Nodes_1 = [b*i/(n+1) + (1-i/(n+1))*a for i in range(n + 2)]
#Nodes_1 = [x*b/(n+1)* sinh(x/n)/sinh((n+1)/n) for x in range (n + 2)]
sa = resol_k1(n)
def g(x):
return uh_k1(sa, Nodes_1, x)
def gprim(x):
return uhprim_k1(sa, Nodes_1, x)
return sqrt(quad(lambda x : (g(x)-f(x)/(pi*pi))**2, a, b)[0])
+sqrt(quad(lambda x : (gprim(x)-cos(x*pi)/(pi))**2, a, b)[0])

def trace_erreur_k1(norm):
xlinspace = [Link]([100,200,500,800, 1000, 3000])
def maxliste(L):
max = L[1] - L[0]

18
 for i in range (len(L) - 1):
if (L[i+1] - L[i]) > max :
max = L[i+1] - L[i]
return max
#H = [maxliste([x*b/(n+1)* sinh(x/n)/sinh((n+1)/n) for x in range (n+2)])
for n in xlinspace]
H = [maxliste([b*i/(n+1) + (1-i/(n+1))*a for i in range(n+2)]) for n in xlinspace]
if norm == "L2":
Erreur = [erreur_k1_L2(n) for n in xlinspace]
[Link]()
[Link]([log(h) for h in H], [log(x) for x in Erreur])
[Link]()
if norm == "H1":
Erreur = [erreur_k1_H1(n) for n in xlinspace]
[Link]()
[Link]([log(h) for h in H], [log(x) for x in Erreur])
[Link]()
return moy_coef([log(h) for h in H], [log(x) for x in Erreur])

5.3 Code pour le cas k = 2

from math import *
import numpy as np
from [Link] import quad
import [Link] as plt

n=20
a=0
b=10

def f(x):
return sin(pi*x)

def solution_exacte(x):
return f(x)/(pi*pi)

Nodes_1 = [b*i/(n+1) + (1-i/(n+1))*a for i in range(n+2)]

#Nodes_1 = [x*b/(n+1)* sinh(x/n)/sinh((n+1)/n) for x in range (n+2)]

Nodes_2 = Nodes_1[1:-1] + [(Nodes_1[i] + Nodes_1[i+1])/2 for i in range (n+1)] + [a,b]

Cells_2 = [[2*n+1,0,n]]+[[i,i+1,i+n+1] for i in range(0,n-1)]+[[n-1,2*n+2,2*n]]

Bords = [2*n+1, 2*n+2]

19
def RefBasicFct(x, k):
if k == 1:
return [Link]([x, 1-x])
if k == 2:
return [Link]([2*(x-1/2)*(x-1), 2*(x-1/2)*x, -4*x*(x-1)])

def DerRefBasisFct(x, k):

if k == 1:
return [Link]([1,-1])
if k == 2:
return [Link]([4*x - 3, 4*x - 1, -8*x + 4])

#Calcul de A_h
A_2 = [Link]((2*n+1,2*n+1))
for k in range (len(Cells_2)):
for i in range(3):
for j in range(3):
I = Cells_2[k][i]
J = Cells_2[k][j]
if I not in Bords and J not in Bords:
A_2[I][J] += 1/(abs(Nodes_2[Cells_2[k][1]]-Nodes_2[Cells_2[k][0]]))*6*
(DerRefBasisFct(0, 2)[i]*DerRefBasisFct(0, 2)[j]+4*
DerRefBasisFct(1/2, 2)[i]*DerRefBasisFct(1/2, 2)[j]+
DerRefBasisFct(1, 2)[i]*DerRefBasisFct(1, 2)[j])

A_2 = [Link](A_2)

#Calcul de F_h
F_2 = [Link](2*n+1)
for k in range (len(Cells_2)):
for i in range(3):
I = Cells_2[k][i]
if I not in Bords:
F_2[I] += (Nodes_2[Cells_2[k][1]]-Nodes_2[Cells_2[k][0]])*6*
(RefBasicFct(0,2)[i]*f(Nodes_2[Cells_2[k][1]])+
4*RefBasicFct(1/2,2)[i]*
f(1/2*Nodes_2[Cells_2[k][0]]+1/2*Nodes_2[Cells_2[k][1]])
+RefBasicFct(1,2)[i]* f(Nodes_2[Cells_2[k][0]]))

F_2 = [Link](F_2)

20
#Calcul de U_h
U_2 = [Link](A_2, F_2)

def reordonne(a,n):
V = [Link](len(a))
for i in range(n):
V[2*i+1] = a[i]
for i in range(n+1):
V[2*i] = a[i+n]
return V

#réordonne et ajoute des zeros sur les bords

segment = [Link](([Link](([Link]([a]),
reordonne(Nodes_2[:-2],n))),[Link]([b])))
solution = [Link](([Link](([Link]([0]),
reordonne(U_2,n))),[Link]([0])))

#tracer les courbes

[Link]()
abscisse = [Link](a,b,200)
[Link](abscisse,[solution_exacte(x) for x in abscisse],label="solution_exacte")
[Link](segment, solution, label="solution_approchée")
[Link](loc=’upper left’)
[Link](-0.11, 0.15)
[Link]()

5.4 Code pour l’erreur dans le cas k = 2

def uh_k2(L, Nodes, x):
i = 0
while i < len(Nodes) and Nodes[i] < x:
i += 1
if i\%2==0:
a = Nodes[i-2]
b = Nodes[i-1]
c = Nodes[i]
fa = L[i-2]
fb = L[i-1]
fc = L[i]
else:
a = Nodes[i-1]
b = Nodes[i]
c = Nodes[i+1]
fa = L[i-1]
fb = L[i]
fc = L[i+1]

21
 return fc/((c-a)*(c-b))*(x-a)*(x-b) + fa/((a-b)*(a-c))*(x-b)*(x-c)
+ fb/((b-a)*(b-c))*(x-a)*(x-c)

def erreur_k2_L2(n):
Nodes_1 = [b*i/(n+1) + (1-i/(n+1))*a for i in range(n+2)]
Nodes_2 = Nodes_1[1:-1] +
[(Nodes_1[i] + Nodes_1[i+1])/2 for i in range (n+1)] + [a,b]
sa = resol_k2(n)
solution = [Link](([Link](
([Link]([0]),sa)),[Link]([0])))
segment = [Link](([Link](
([Link]([a]),reordonne(Nodes_2[:-2],n))),[Link]([b])))
def g(x):
return uh_k2(solution,segment,x)
return sqrt(quad(lambda x : (g(x) - f(x)/(pi*pi))**2,a,b)[0])

def uhprim_k2(L, Nodes, x):

i = 0
while i<len(Nodes) and Nodes[i] < x:
i += 1
if i\%2==0:
a = Nodes[i-2]
b = Nodes[i-1]
c = Nodes[i]
fa = L[i-2]
fb = L[i-1]
fc = L[i]
else:
a = Nodes[i-1]
b = Nodes[i]
c = Nodes[i+1]
fa = L[i-1]
fb = L[i]
fc = L[i+1]
return fc/((c-a)*(c-b))*(2*x-a-b) + fa/((a-b)*(a-c))*(2*x-b-c)
+ fb/((b-a)*(b-c))*(2*x-a-c)

def erreur_k2_H1(n):
Nodes_1 = [b*i/(n+1) + (1-i/(n+1))*a for i in range(n+2)]
Nodes_2 = Nodes_1[1:-1] +
[(Nodes_1[i] + Nodes_1[i+1])/2 for i in range (n+1)] + [a,b]
sa = resol_k2(n)
solution = [Link](([Link](
([Link]([0]),sa)),[Link]([0])))
segment = [Link](([Link](
([Link]([a]),reordonne(Nodes_2[:-2],n))),[Link]([b])))
def g(x):
return uh_k2(solution,segment,x)

22
 def gprim(x):
return uhprim_k2(solution,segment,x)
return sqrt(quad(lambda x : (g(x)-f(x)/(pi*pi))**2,a,b)[0])
+ sqrt(quad(lambda x : (gprim(x)-cos(x*pi)/(pi))**2,a,b)[0])

def trace_erreur_k2(norm):
xlinspace = [Link]([100,200,500,800, 1000, 3000])
if norm == "infini":
Erreur = [erreur_k2_infini(n) for n in xlinspace]
[Link]()
[Link]([log(1/n) for n in xlinspace],[log(x) for x in Erreur])
[Link]()

if norm == "L2":
Erreur = [erreur_k2_L2(n) for n in xlinspace]
[Link]()
[Link]([log(1/n) for n in xlinspace],[log(x) for x in Erreur])
[Link]()

if norm == "H1":
Erreur = [erreur_k2_H1(n) for n in xlinspace]
[Link]()
[Link]([log(1/n) for n in xlinspace],[log(x) for x in Erreur])
[Link]()
return moy_coef([log(1/n)for n in xlinspace],[log(x) for x in Erreur])

5.5 Preuve du théorème 2

On veut montrer que
ku − Πh (u)kH 1 6 Ch2 ku(3) kL2 .
Étape 1 : Montrons que ku00 − Πh (u)00 kL2 6 hku(3) kL2 .
On va raisonner par densité, on suppose ici que u est C ∞ .
N Z
X xi+1
00
ku − Πh (u)00 k2L2 = (u00 (x) − Πh (u)00 (x))2 dx (8)
i=0 xi

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

2u(xi ) 4u(xei ) 2u(xi+1 )
Πh (u)(x) = 2 (x − xi+1 )(x − xei ) − 2 (x − xi )(x − xi+1 ) + (x − xi )(x − xei )
hi hi h2i
Ainsi, on a
4
Πh (u)00 (x) = (u(xi ) − 2u(xei ) + u(xi+1 )) = u00 (yi )
h2i
où yi ∈]xi , xi+1 [ 3 .
u(t + h) − 2u(t) + u(t − h)
3. considérer et appliquer le théorème de Rolle à une fonction bien choisie (similaire à
h2
la preuve du théorème des accroissements finis)

23
 Z xi+1 Z xi+1
00 00
(u (x) − Πh (u) (x)) dx = 2
(u00 (x) − u00 (yi ))2 dx
xi xi
Z xi+1 Z x 2
(3)
= u (t)dt dx
x y
Z ixi+1 Z ix  Z x 
(3) 2
6 u (t) dt 1dt dx par Cauchy–Schwarz
x yi yi
Z ixi+1 Z  x
6 hiu(3) (t)2 dt dx
x y
Z ixi+1 Z ixi+1 
(3) 2
6 hi u (t) dt dx
xi xi
Z xi+1 Z xi+1 
(3) 2
6 hi u (t) dx dt par Fubini
xi xi
Z xi+1
6 h2i u(3) (t)2 dt (9)
xi

Ainsi, en injectant dans l’égalité (8), on a

ku00 − Πh (u)00 k2L2 6 h2 ku(3) k2L2 .

En passant à la racine, on a
ku00 − Πh (u)00 kL2 6 hku(3) kL2 .

Étape 2 : Montrons que ku0 − Πh (u)0 kL2 6 h2 ku(3) kL2 .

N Z
X xi+1
0
ku − Πh (u)0 k2L2 = |u0 (x) − Πh (u)0 (x)|2 dx
i=0 xi

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

Z x 2
|u0 (x) − Πh (u)0 (x)|2 = (u00 (t) − Πh (u)00 (t))dt
x
Z ix
6 hi (u00 (t) − Πh (u)00 (t))2 dt par Cauchy–Schwarz
x
Z ixi+1
6 h3i u(3) (t)2 dt par l’inégalité (9) (10)
xi

On intègre sur [xi , xi+1 ] pour obtenir

Z xi+1 Z xi+1
|u0 (x) − Πh (u)0 (x)|2 dx 6 h4i u(3) (t)2 dt.
xi xi

Ainsi
ku0 − Πh (u)0 k2L2 6 h4 ku(3) k2L2 .
En passant à la racine carrée, on obtient

ku0 − Πh (u)0 kL2 6 h2 ku(3) kL2 .

24
Étape 3 : Montrons que ku − Πh (u)kL2 6 h3 ku(3) kL2 .
N Z xi+1
X
ku − Πh (u)k2L2 = |u(x) − Πh (u)(x)|2 dx
i=0 xi

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

Z x 2
2 0 0
|u(x) − Πh (u)(x)| = (u (t) − Πh (u) (t))dt
xi
Z x
6 hi (u0 (t) − Πh (u)0 (t))2 dt par Cauchy–Schwarz
x
Z ixi+1 Z xi+1
6 hi h3i u(3) (s)2 dsdt par l’inégalité (10)
xi xi
Z xi+1
6 h5i u(3) (s)2 ds
xi
On intègre sur [xi , xi+1 ] pour obtenir
Z xi+1 Z xi+1
|u(x) − Πh (u)(x)|2 dx 6 h6i u(3) (t)2 dt.
xi xi
Ainsi
ku − Πh (u)k2L2 6 h6 ku(3) k2L2 .
En passant à la racine carrée, on obtient
ku − Πh (u)kL2 6 h3 ku(3) kL2 .
Étape 3 : Conclusion.
ku − Πh (u)kH 1 = ku − Πh (u)kL2 + ku0 − Πh (u)0 kL2
6 h3 ku(3) kL2 + h2 ku(3) kL2
6 (1 + h)h2 ku(3) kL2
6 Ch2 ku(3) kL2 avec C = |I| + 1

5.6 Résultats plus généraux

Lemme 2 (Céa). Pour u la solution de (1) dans un Hilbert V , uh solution de (3), on a

M
ku − uh kV 6 inf ku − vh kV
α vh ∈Vh
où M est la constante de continuité de a, α est la constante de coercivité de a avec a la forme
bilinéaire de (2).

Théorème 3. Pour u solution H k+1 (I) de (1) et Πh (u) la projection de u sur l’espace Vh des
polynômes de degré au plus k par morceaux sur un maillage τh , on a

ku − Πh (u)kH 1 (I) 6 Chk ku(k+1) kL2 (I)

où h = max |xi+1 − xi |.

[xi ,xi+1 ]∈τh

25