Éléments finis 1d (séances 6-7)

Laurent Monasse∗

∗ Université Côte d’Azur, Inria, CNRS, LJAD, Nice

27 février 2023

L. Monasse (Inria & LJAD, Nice) Méthodes numériques pour les EDPs 27/02/2023 1 / 12
 Introduction

Rappel : méthode de Galerkin

Soit Vh un sev de dimension finie de (V , k · kV espace de Hilbert, sur lequel

est écrite la formulation variationnelle : trouver u ∈ V tel que

∀v ∈ V , a(u, v ) = b(v ).

Méthode de Galerkin : trouver uh ∈ Vh tel que

∀vh ∈ Vh , a(uh , vh ) = b(vh ).

Choix de la base : prendre une suite d’espaces d’approximation Vh qui soit

dense dans V (convergence par le lemme de Céa)
Méthode des Éléments Finis : choix d’une base polynômiale par morceaux
PN
Si on prend une base (Ψi )1≤i≤N de Vh et qu’on écrit uh = i=1 ui Ψi , alors
le problème à résoudre est de la forme Ah Uh = Bh , avec (Ah )ij = a(Ψj , Ψi ).

L. Monasse (Inria & LJAD, Nice) Méthodes numériques pour les EDPs 27/02/2023 2 / 12
 Introduction

Rappel : méthode de Galerkin

Choix de la base : prendre une suite d’espaces d’approximation Vh qui soit

dense dans V (convergence par le lemme de Céa)
Méthode des Éléments Finis : choix d’une base polynômiale par morceaux
Quels sont ces “morceaux” ? ○1 Maillage
2 Interpolation
Quelle base ? ○
PN
Si on prend une base (Ψi )1≤i≤N de Vh et qu’on écrit uh = i=1 ui Ψi , alors
le problème à résoudre est de la forme Ah Uh = Bh , avec (Ah )ij = a(Ψj , Ψi ).
Comment calculer (Ah )ij en pratique ?
3
○ Réduction au calcul de tableaux élémentaires
4
○ Assemblage de la matrice Ah
5
○ Prise en compte des conditions aux limites ?
6
○ Résolution du système linéaire Ah Uh = Bh ?

L. Monasse (Inria & LJAD, Nice) Méthodes numériques pour les EDPs 27/02/2023 2 / 12
 Introduction

Rappel : méthode de Galerkin

Choix de la base : prendre une suite d’espaces d’approximation Vh qui soit

dense dans V (convergence par le lemme de Céa)
Méthode des Éléments Finis : choix d’une base polynômiale par morceaux
Quels sont ces “morceaux” ? ○1 Maillage
2 Interpolation
Quelle base ? ○
PN
Si on prend une base (Ψi )1≤i≤N de Vh et qu’on écrit uh = i=1 ui Ψi , alors
le problème à résoudre est de la forme Ah Uh = Bh , avec (Ah )ij = a(Ψj , Ψi ).
Comment calculer (Ah )ij en pratique ?
3
○ Réduction au calcul de tableaux élémentaires
4
○ Assemblage de la matrice Ah
5
○ Prise en compte des conditions aux limites ?
6
○ Résolution du système linéaire Ah Uh = Bh ?
Ce sera le plan des deux séances, sur l’exemple le plus simple : éléments finis
P1 en 1d.

L. Monasse (Inria & LJAD, Nice) Méthodes numériques pour les EDPs 27/02/2023 2 / 12
 Maillage

1 Maillage
○
Les éléments finis s’appuient sur un maillage pour pouvoir définir la base
d’approximation : c’est un découpage de l’espace en sous-domaines de forme
prescrite
En 3d, un maillage peut être assez complexe

Courtesy Gmsh, (C) C. Geuzaine et J.F. Remacle

En 1d, c’est juste un découpage en segments ei = [xi−1 , xi ], avec

xi = xi−1 + hi
xi =noeud, ei =élement, hi =taille de l’élément, Nv =nombre de noeuds,
Ne=nombre d’éléments
e1 e2 e3 e4
◦ ◦ ◦ ◦ ◦
x0 x1 x2 x3 x4
L. Monasse (Inria & LJAD, Nice) Méthodes numériques pour les EDPs 27/02/2023 3 / 12
 Interpolation

2 Choix de la base
○
On prend sur chaque élément ei = [xi , xi+1 ] des fonctions affines :
vh (x) = ai x + bi .
1
De plus, Vh ⊂ H ⇒ les fonctions de base doivent être continues
Base : Fonctions de forme ϕi valant 1 sur le noeud xi et 0 sur les autres :
x − x
i−1
 sur ei ,
h



 i
ϕi (x) = xi+1 − x sur e ,
i+1
 hi+1



0 sinon.

PNv
Interpolation : vh (x) = i=0 vi ϕi (x).
ϕ

ϕi

• • • • •
xi−1 xi xi+1
L. Monasse (Inria & LJAD, Nice) Méthodes numériques pour les EDPs 27/02/2023 4 / 12
 Tableaux élémentaires

3 Calcul des tableaux élémentaires

○
Il faut calculer les Aij = a(ϕj , ϕi )
Dans le cas de l’équation de la chaleur,
Z XZ
a(uh , vh ) = uh0 (x)vh0 (x)dx = uh0 (x)vh0 (x)dx.
Ω i ei

Sur chaque ei , vh (x) = vi−1 ϕi−1 (x) + vi ϕi (x)

On se ramène donc à calculer les intégrales des fonctions de forme sur
l’élément ei , donc uniquement pour celles non nulles sur l’élément, ϕi−1 et ϕi .
Calcul de la matrice de rigidité élémentaire K e qui donne la dépendance de
l’intégrale sur ei en fonction des u e = (ui−1 , ui )T et v e = (vi−1 , vi )T :
Z
uh0 (x)vh0 (x)dx = (K e u e |v e ).
ei

Dans notre cas, c’est simple : sur ei

1 1
ϕ0i−1 (x) = − , ϕ0i (x) =
hi hi
L. Monasse (Inria & LJAD, Nice) Méthodes numériques pour les EDPs 27/02/2023 5 / 12
 Tableaux élémentaires

Matrice de rigidité élémentaire pour l’équation de la chaleur

Z Z
uh0 vh0 = (ui−1 ϕ0i−1 + ui ϕ0i )(vi−1 ϕ0i−1 + vi ϕ0i )
ei ei
Z   
ui−1 ui vi−1 vi ui − ui−1 vi − vi−1
= − + − + = hi ·
ei hi hi hi hi hi hi
    
1 +1 −1 ui−1 vi−1
=
hi −1 +1 ui vi

Expression de la matrice de rigidité élémentaire K e sur ei :

 
1 +1 −1
Ke =
hi −1 +1

L. Monasse (Inria & LJAD, Nice) Méthodes numériques pour les EDPs 27/02/2023 6 / 12
 Tableaux élémentaires

Calcul du second membre

R
On a également besoin d’évaluer b(vh ) = Ω fvh .
P R
On découpe sur les éléments : b(vh ) = i ei fvh , puis on calcule le vecteur
élémentaire F e = (fi−1 , fi )T tel que
Z
f (x)vh (x)dx = hF e |v e i
ei
Z Z
= vi−1 f (x)ϕi−1 (x)dx +vi f (x)ϕi (x)dx
ei ei
| {z } | {z }
fi−1 fi

En général, pas de formule exacte pour fi et fi−1 : quadrature numérique

Par exemple, formule des trapèzes :
hi hi
fi−1 ≈ f (xi−1 ), fi ≈ f (xi )
2 2

L. Monasse (Inria & LJAD, Nice) Méthodes numériques pour les EDPs 27/02/2023 7 / 12
 Assemblage

Assemblage de la matrice de rigidité

Comme son nom l’indique : on rassemble les matrices élémentaires...

On a écrit les intégrales par élément. Si on revient à l’intégrale complète :
Z X
uh0 vh0 = hK e u e |v e i
Ω e
Z X
fvh = hF e |v e i
Ω e

Si on revient à Uh = (ui ) et Vh = (vi ) les vecteurs complets, on doit écrire

Z Z
0 0
uh vh = hAh Uh |Vh i, fvh = hF |Vh i
Ω Ω

avec Ah la matrice de rigidité globale, F le second membre global.

On peut alors faire le lien entre les matrices 2 × 2 locales K e et la matrice
N × N globale Ah .

L. Monasse (Inria & LJAD, Nice) Méthodes numériques pour les EDPs 27/02/2023 8 / 12
 Assemblage

Exemple : assemblage pour l’équation de la chaleur

Si on regarde une ligne de la matrice Ah , il s’agit des coefficients devant les
ui de hAh Uh |Ei i où Ei est le vecteur dont tous les termes sont nuls sauf le
ième égal à 1.
Z XZ
0 0
hAh Uh |Ei i = a(uh , ϕi ) = uh ϕi = uh0 ϕ0i
Ω j ej

−1 ui − ui−1 ui − ui+1
Z Z
1
= uh0 + uh0 = +
ei hi ei+1 hi+1 hi hi+1
Donc la matrice Ah est de la forme :
 1
− h11

h1 0 ... ... 0

 −1 1 1
.. 
 h1 h1 + h2 − h12 0 . 

 .. .. 
 0 − h12 1
h2 + h3
1
− h13 . . 
Ah =  .
 
 . .. .. .. .. 
 . . . . . 0


 . 
 . 1 1 1
 . 0 − hN−1 + − h1N

hN−1 hN 
0 ... ... 0 − h1N 1
hN
L. Monasse (Inria & LJAD, Nice) Méthodes numériques pour les EDPs 27/02/2023 9 / 12
 Assemblage

Algorithme d’assemblage

Initialiser Ah(1 :N,1 :N) et F (1 :N) à zéro ;

for ie=1 à Ne do . Boucle sur les éléments :
for i=1 à NDLE do . Boucle sur les DDL par élément :
in=loc(ie,i) . in numéro du ième noeud de l’élément ie
Calculer F e (i)
Faire F (in) = F (in) + F e (i)
for j=1 à NDLE do
jn=loc(ie,j) . jn numéro du jème noeud de l’élément ie
Calculer K e (i, j)
Faire Ah(in, jn) = Ah(in, jn) + K e (i, j)
end for
end for
end for
Termes en bleus : utilisation du tableau de numérotation de degrés de liberté
de l’élément dans le maillage ; récupéré du mailleur

L. Monasse (Inria & LJAD, Nice) Méthodes numériques pour les EDPs 27/02/2023 10 / 12
 Conditions aux limites

Imposition des conditions aux limites

On sait que les conditions de Neumann sont prises en compte de façon

naturelle par la formulation variationnelle ; il faut cependant des conditions de
Dirichlet pour qu’il y ait unicité de la solution
Comment imposer la condition uh (0) = 0 ? Plusieurs possibilités :
Ajouter cette équation supplémentaire sur le DDL u1 : peut compliquer la
résolution du système
Pénaliser la condition de Dirichlet : on remplace par exemple Ah(1, 1) par
1012 , et on a X
u1 = 10−12 (f1 − Ah(1, j)uj ) ≈ 0
j

Supprimer le DDL sur la condition de Dirichlet : naturel puisqu’on en fait une

condition essentielle, mais nécessite de savoir manipuler la taille du système

L. Monasse (Inria & LJAD, Nice) Méthodes numériques pour les EDPs 27/02/2023 11 / 12
 Résolution

Résolution du système linéaire

On n’inverse jamais la matrice Ah : trop coûteux, et Ah est creuse alors que
A−1
h est généralement pleine, donc il faut la stocker
Il y a deux principales façons de résoudre un système linéaire Ah Uh = Bh :
Méthodes directes (en deux étapes)
Factorisation (le plus coûteux) : Ah = LU (Gauss) ou K = B T B (Cholesky),
avec L triangulaire inférieure, U et B triangulaires supérieures
Descente/remontée :
Gauss : Ah x = LUx = f ⇐⇒ Lw = f puis Ux = w
Cholesky : Ah x = B T Bx = f ⇐⇒ B T w = f puis Bx = w
Systèmes assez petits (actuellement, pas plus du million de DDL)
Méthodes itératives (convergent sous conditions sur Ah )
Décomposition : Ah = M − N avec M facile à inverser. Choix possibles :
M = diag(Ah ) (Jacobi) ou M = triang_inf(Ah ) (Gauss-Seidel)
Itérations : initialiser x (0) , puis résoudre
Mx (k+1) = Nx (k) + f
Souvent, le problème converge lentement : il est mal conditionné (rapport
entre plus grande et plus petite valeur propre de Ah ). On peut chercher un
préconditionneur D (matrice inversible “proche” de A−1 h ) et résoudre le
problème équivalent DAh = Df .
L. Monasse (Inria & LJAD, Nice) Méthodes numériques pour les EDPs 27/02/2023 12 / 12
 Résolution

L. Monasse (Inria & LJAD, Nice) Méthodes numériques pour les EDPs 27/02/2023 12 / 12