Analyse numériques des EDP

MNSA S3

Pr. A. RADID
Université Hassan II
Faculté des Sciences Aïn Chock
Département de Mathématiques et Informatique

Année universitaire 2021-2022

Table des matières

Introduction 1

1 Di¤érences …nies 4
1.1 Préliminaires et dé…nitions . . . . . . . . . . . . . . . . . . . . . . . 4
1.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.1.2 Résolution numérique des problèmes physiques . . . . . . . . 4
1.1.3 Dé…nition d’une EDP . . . . . . . . . . . . . . . . . . . . . . 5
1.1.4 Ordre de l’EDP . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.1.5 EDP linéaire ou non linéaire . . . . . . . . . . . . . . . . . . 7
1.1.6 Les di¤érents types d’équations aux dérivées partielles . . . . 8
1.1.7 Problème aux limites . . . . . . . . . . . . . . . . . . . . . . 8
1.1.8 Problème bien posé . . . . . . . . . . . . . . . . . . . . . . . 8
1.1.9 Les principales méthodes de discrétisation d’EDP . . . . . . 8
1.2 Les méthodes de di¤érences …nies . . . . . . . . . . . . . . . . . . . 9
1.2.1 Présentation de la méthode . . . . . . . . . . . . . . . . . . 9
1.2.2 Maillage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.2.3 Formules d’approximation de dérivées par di¤érences divisées
en dimension 1 . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.2.4 Stencil d’un schéma . . . . . . . . . . . . . . . . . . . . . . . 19
1.2.5 Schéma à multipas . . . . . . . . . . . . . . . . . . . . . . . 19
1.2.6 L’ordre d’un schéma . . . . . . . . . . . . . . . . . . . . . . 19
1.2.7 Discrétisation des conditions aux limites . . . . . . . . . . . 20
1.3 Consistance, Stabilité et Convergence . . . . . . . . . . . . . . . . . 25
1.3.1 Erreur de Troncature . . . . . . . . . . . . . . . . . . . . . . 25
1.3.2 Consistance . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
1.3.3 Stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
1.3.4 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
1.4 Equation de Laplace 1D . . . . . . . . . . . . . . . . . . . . . . . . 28
1.5 Equation de la chaleur 1D . . . . . . . . . . . . . . . . . . . . . . . 33
1.5.1 Schéma explicite de l’équation de la chaleur 1D . . . . . . . 33
1.5.2 Schéma implicite de l’équation de la chaleur . . . . . . . . . 41

2
 1.6 Approximation par di¤érences …nies de l’équation de la chaleur en 2D 43

2 Eléments …nis 46
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
2.1.1 Méthode des éléments …nis . . . . . . . . . . . . . . . . . . . 47
2.1.2 Principe de base . . . . . . . . . . . . . . . . . . . . . . . . 47
2.2 Méthode de Galerkin . . . . . . . . . . . . . . . . . . . . . . . . . . 47
2.2.1 Théorème de Lax Milgram . . . . . . . . . . . . . . . . . . . 47
2.2.2 Choix de l’espace V et formulation variationnelle . . . . . . 48
2.2.3 Résolution d’un système linéaire . . . . . . . . . . . . . . . . 48
2.3 Exemples en dimension 1 et introduction des éléments …nis P1 . . . 51
2.3.1 Exemple en dimension 1 . . . . . . . . . . . . . . . . . . . . 51
2.3.2 Éléments …nis P1 . . . . . . . . . . . . . . . . . . . . . . . . 51
2.3.3 Écriture matricielle . . . . . . . . . . . . . . . . . . . . . . . 53
2.3.4 Calcul du second membre . . . . . . . . . . . . . . . . . . . 55
2.3.5 Majoration de l’erreur dans le cas de l’exemple . . . . . . . . 57
2.3.6 Fonctions de forme locales . . . . . . . . . . . . . . . . . . . 62
2.3.7 Technique d’assemblage . . . . . . . . . . . . . . . . . . . . 69
2.4 Eléments …nis P2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
2.5 Eléments …nis Lagrangiens d’ordre 3 . . . . . . . . . . . . . . . . . 77
2.6 Eléments …nis d’Hermite . . . . . . . . . . . . . . . . . . . . . . . . 78
2.6.1 Résolution pratique de l’équation des plaques par la méthode
des éléments …nis d’Hermite P3 . . . . . . . . . . . . . . . . . 81
2.7 Introduction des éléments …nis en dimension N=2 ou 3 . . . . . . . 84
2.7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
2.7.2 Éléments …nis triangulaires . . . . . . . . . . . . . . . . . . . 84
2.7.3 Coordonnées barycentriques . . . . . . . . . . . . . . . . . . 85
2.7.4 Maillage en dimension N=2 ou 3 . . . . . . . . . . . . . . . 85
2.7.5 Treillis d’ordre k . . . . . . . . . . . . . . . . . . . . . . . . 86
2.7.6 Les Mailles . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
2.7.7 Ensemble de polynômes Pk . . . . . . . . . . . . . . . . . . 88
2.7.8 Eléments …nis Pk . . . . . . . . . . . . . . . . . . . . . . . . 89
2.7.9 Technique de l’élément de référence . . . . . . . . . . . . . . 94
2.7.10 Intégrales doubles et changement de variables . . . . . . . . 96
2.7.11 Intégrales curvilignes . . . . . . . . . . . . . . . . . . . . . . 97
2.7.12 Éléments …nis rectangulaires . . . . . . . . . . . . . . . . . . 98

3
Chapitre 1

Di¤érences …nies

1.1 Préliminaires et dé…nitions

1.1.1 Introduction
Les équations aux dérivées partielles (EDP) interviennent dans la description de
très nombreux problèmes de physique, chimie, sciences de la Terre, biologie : mé-
canique des ‡uides, propagation des ondes, électromagnétisme, phénomènes de
di¤usion ...

1.1.2 Résolution numérique des problèmes physiques

Grandes étapes :

Problème physique continu est décrit par un modèle mathématique continue

(mis en équations)

Modèle mathématique continu est discretisé en s’appuyant sur une(des) méth-

ode(s) numérique(s)

Equations discretisées sont approximées à l’aide des schémas numériques

appropriés, l’algorithme de résolution est établie

Algorithme est codé ( Matlab, Gnuplot, FreeFem++...)

Si tout va bien, la solution approchée du problème initial est obtenue.

4
1.1.3 Dé…nition d’une EDP
Une équation aux dérivées partielles est une équation faisant intervenir une fonc-
tion inconnue de plusieurs variables u (x1 ; : : : ; xd ) ainsi que certaines de ses dérivées
partielles. F (x; u; D uj j p ) = 0:
Dans les applications à la physique, on distingue les problèmes stationnaires
où les variables sont des variables d’espaces : x 2 Rd ! u(x) 2 RN des prob-
lèmes d’évolution où l’une des variables est le temps qui joue un rôle particulier
x 2 Rd ; t 0 ! u(x; t) 2 RN .

Quelques exemples
Équation de la chaleur

– Equation de la chaleur monodimensionnelle :

8 @u(x;t) 2 u(x;t)
>
< @t D @ @x 2 = f (x; t); x 2]0; L[; t > 0
u(x; 0) = u0 (x); x 2 ]0; L[ conditions initiales
>
:
u(0; t) = u(L; t) = 0; t > 0 conditions aux limites
D coe¢ cient de di¤usion thermique
– En dimension supérieure d :
8 @u(x;t)
>
< @t 4u(x; t) = f (x; t); x2 ;t>0
u(x; 0) = u0 (x); x2
>
:
u(x; t) = 0; x2@ ; t>0
P @2u
4u = @x2 i
1 i d
ouvert borné de bord @

Equation des ondes

–8
En dimension 1
@ 2 u(x;t) 2 u(x;t)
>
< @t2 c2 @ @x 2 = f (x; t) x 2]0; L[; t > 0
u(x; 0) = u0 (x); @t u (x; 0) = v0 (x); x 2 ]0; L[
>
:
u(0; t) = u(L; t) = 0 t>0
@t2 au lieu de @t pour l’équation de la chaleur et deux conditions initiales
(position, vitesse) au lieu d’une.
u(x; t) est un champ scalaire dépend de la position x et du temps t,
comme la pression d’un gaz dans le cas du son.
La constante c est la vitesse de propagation des ondes.

5
Pr. A. Radid CHAPITRE 1. DIFFÉRENCES FINIES

– En dimension supérieure
8 @ 2 u(x;t)
>
< @t2 4u(x; t) = f (x; t) x2 ;t>0
u(x; 0) = u0 (x); @t u (x; 0) = v0 (x); x 2
>
:
u(x; t) = 0 x2@ ; t>0
propagation du son (acoustique)
u pression acoustique ou potentiel des vitesses

Equation de poisson

– En une dimension
0
(ku0 ) = f

le réel k est un coe¢ cient de di¤usion

Si k = 1, l’équation de Poisson s’écrit u" = f
En deux dimensions d’espace, elle s’écrit : 4u = f
Dans le cas f = 0, on obtient l’équation de Laplace 4u = 0:
– En trois dimension d’espace

divkru = f

le terme kru est un ‡ux de di¤usion

Equation de transport

– En une dimension
ut + cux = 0; 8x 2 R; 8t > 0
u(x; 0) = u0 (x) 8x 2 R
où c > 0 la vitesse de transport

L’équation de Schrödinger :

@u 1
i = 4u + V u
@t 2

qui décrit en mécanique quantique la fonction d’onde d’une particule massive

dans un potentiel V .
Elle représente aussi la propagation d’une onde harmonique à la limite parax-
iale, dans un milieu de permittivité = n2 = V + 1:

6
1.1.4 Ordre de l’EDP
Dé…nition 1.1 on appelle ordre d’une EDP l’ordre de la plus grande dérivée
présente dans l’équation.

Exemple 1.1 L’équation de la chaleur et l’équation des ondes sont d’ordre

2.

L’équation de la poutre élastique est d’ordre 4:

(
u0000 (x) (a(x)u0 (x))0 + c(x)u(x) = f (x) dans ]0; L[
u(0) = u0 (0) = u(L) = u0 (L) = 0

1.1.5 EDP linéaire ou non linéaire

Dé…nition 1.2 On dit que l’EDP est linéaire si elle se met sous la forme Lu = f
, où u ! Lu est une application linéaire par rapport à u . Sinon, elle est non
linéaire .

Exemples 1.1 1) L’équation de Black-Scholes qui donne l’évolution du prix d’une

option dans un marché …nancier est donnée par :
2
@v 1 2 2@ v(t; x) @v(t; x)
(t; x) + x 2
+ rx rv(t; x) = 0;
@t 2 @x @x

où le prix de l’option, v , est une fonction du prix du sous-jacent, x et du temps

t; r est le taux d’intérêt sans risque, et la volatilité du stock.
L’équation de Black-Scholes est une EDP linéaire d’ordre 2.
2) Equation de propagation d’onde

@2u @2u
=0
@t2 @x2

est linéaire d’ordre 2.

3) l’équation de Bürgers

ut + (u2 )x = 0; t 2 R+ ; x 2 R
u (x; 0) = u0 (0)

est une équation EDP non linéaire.

7
Pr. A. Radid CHAPITRE 1. DIFFÉRENCES FINIES

1.1.6 Les di¤érents types d’équations aux dérivées par-

tielles
Considérons une équation aux dérivées partielles linéaire, de degré 2, de la forme :
2 @2u 2
A @@t2u + B @t@x + C @@xu2 + D @u
@t
+ E @u
@x
+ Fu = 0
On distingue trois grandes catégories d’EDP:

Si B 2 4AC < 0, l’équation est dite elliptique

Exemple 1.2 Équation de Laplace : c4u = f

Si B 2 4AC = 0, l’équation est dite parabolique

@u
Exemple 1.3 Équation de la chaleur: @t
c4u = f ou l’équation de Black
scholes.

Si B 2 4AC > 0, l’équation est dite hyperbolique

@2u
Exemple 1.4 Équation des ondes : @t2
c2 4u = f

1.1.7 Problème aux limites

On appelle problème aux limites une équation aux dérivées partielles munie de
conditions aux limites sur la totalité de la frontière du domaine sur lequel elle est
posée.

1.1.8 Problème bien posé

Le problème est bien posé si pour toute donnée (second membre, domaine, données
au bord, etc.), il admet une solution unique et si cette solution dépend continument
de la donnée.

1.1.9 Les principales méthodes de discrétisation d’EDP

Les principales méthodes de discrétisation d’EDP
Les problèmes de type EDP sont in…ni-dimensionnels (il s’agit de detérminer
toute une fonction). Si l’on veut calculer la solution numériquement, il faut se
ramener à un problème en dimension …nie. Les principales méthodes utilisées sont
:

8
 Les méthodes de di¤érences …nies (ou de volumes …nis): discrétisation du
domaine et remplacement de l’opérateur di¤érentiel par un quotient dif-
férentiel.

Les méthodes de volumes …nis sont adaptées aux équations de conservation

et utilisées en mécanique des ‡uides depuis plusieurs décennies. Le principe
consiste à découper le domaine en des ”volumes de contrôle” ; on intègre
ensuite l’équation de conservation sur les volumes de contrôle ; on approche
alors les ‡ux sur les bords du volume de contrôle par exemple par une tech-
nique de di¤érences …nies.

Les méthodes d’éléments …nis : très liée à la formulation variationnelle des

EDP.
Principe: remplacer un espace de Hilbert H par un sous espace de dimension
…nie H N

Les méthodes spectrales : Chercher une solution approchée sous forme d’un
développement sur une certaine famille de fonctions (Fourier, polynômes,
splines, etc.) :

X
n
u= i (u) pi
i=1

Méthodes souvent coûteuses, mais précises.

Note : On ne va pas étudier cette dernière méthode.

1.2 Les méthodes de di¤érences …nies

1.2.1 Présentation de la méthode
La méthode des di¤érences …nies est basée sur la ré-écriture sous forme discrète de
tous les termes de dérivées présents dans les équations et dans les conditions aux
limites et/ou initiales sous forme discrète.
On considère un domaine physique Rd , où d est la dimension de l’espace.
Le principe des méthodes de di¤érences …nies (DF) consiste à se donner un certain
nombre de points du domaine, qu’on notera (x1 ; : : : ; xd ) (Rd )N . On approche
l’opérateur di¤érentiel en espace en chacun des xji par des quotients di¤érentiels.
Il faut alors discrétiser la dérivée en temps : on pourra par exemple considérer un
schéma d’Euler explicite ou implicite pour la discrétisation en temps.

9
Pr. A. Radid CHAPITRE 1. DIFFÉRENCES FINIES

L’équation discrète di¤ère ainsi de l’équation continue par l’erreur de troncature

associée au schéma. L’interprétation correcte des résultats issus d’un calcul par
di¤érence …nies, requiert la véri…cation des points suivants :

les conditions aux limites sont respectées ;

le schéma est consistant ;

le schéma est stable ;

le maillage est choisi su¢ samment ’…n’pour s’assurer d’une quasi-indépendance

de la solution à son égard ;

Avantages : grande simplicité d’écriture et faible coût de calcul.

Inconvénients : limitation à des géométries simples, di¢ cultés de prise en
compte des conditions aux limites de type Neumann.

1.2.2 Maillage
Un maillage est la discrétisation spatiale d’un milieu continu, ou aussi, une mod-
élisation géométrique d’un domaine par des éléments proportionnés …nis et bien
dé…nis.
Exemples :

1. Pour d = 1, si on discrétise que l’espace.

Découpage irrégulier
On se donne une subdivision de [0; 1], c’est à dire une suite de points
(xk )k=0;:::;N +1 tels que x0 = 0 < x1 < : : : < xN < xN +1 = 1:
Pour i = 0; : : : ; N , on note hi+ 1 = xi+1 xi et on dé…nit le pas du
2
maillage par : h = max hi+ 1 :
i=0;:::;N 2

Découpage régulier
Sur un intervalle [a; b], pour un pas du temps constant on prend xi+1
xi = h 8i 2 [0; N ] ; xi = x0 + ih avec x0 = a et h = bNa

10
 2. Pour d = 2; si on discrétise l’espace = ]0; 1[ ]0; 1[ :

Découpage régulier :
On découpe ]0; 1[ en N intervalles sur axe des x de longueur h, xi = ih avec
h = N1
De même sur axe des y de longueur h, y j = jh.

Maillage d’un regtangle en di¤érences …nies

1.2.3 Formules d’approximation de dérivées par di¤érences

divisées en dimension 1
La résolution d’un problème continu est remplacé par la recherche de (N + 1)
valeurs discrètes ui , approchant la valeur de la solution exacte u(xi ) aux points
du maillage, points de grille, xi :

ui u(xi ) avec i = 0; : : : ; N

11
Pr. A. Radid CHAPITRE 1. DIFFÉRENCES FINIES

On identi…e le pas de discrétisation par xi+1 xi = 4xi = hi : Fréquemment,

le pas est constant, alors il est simplement noté par h.
On donne maintenant quelques formules simples d’approximation des dérivées
par des di¤érences divisées

Dérivées premières
u est de classe C 2 au voisinage de x; en utilisant le développement limité à l’ordre
2, on obtient :
2
u(x + h) = u(x) + hu0 (x) + h2 u00 (x) + O(h3 )
2
u(x h) = u(x) hu0 (x) + h2 u00 (x) + O(h3 )
u0 (x) = u(x+h)h u(x) + O(h) ordre 1, décentré avant
0 u(x) u(x h)
u (x) = h
+ O(h) ordre 1, décentré arrière
0 u(x+h) u(x h) 2
u (x) = 2h
+ O(h ) ordre 2, centré
Les deux premières formules décentrées d’ordre 1, sont obtenues directement
à partir des développemts en avant et en arrière respectivement. La troisième
équation qui provient de la soustraction au deuxième développenet du premier.
Au point du maillage xi ; on obtient :
@u ui+1 ui
u0 (xi ) = (xi ) '
@x h
h2 (3)
On a : @u
@x
(xi ) = u(xi +h)2hu(xi h) 6
u ( i)
u(x + h ) u(x h ) h2 (3)
ou bien @u
@x
(xi ) = i 2 h i 2 24
u ( i)
Ce qui conduit, dans le cas de discrétisations uniformes de pas constant h , à :
ui+ 1 ui 1
u0 (xi ) = @u
@x
(xi ) ' ui+12hui 1 ou 2
h
2

ui+ 1 et ui 1 sont les valeurs approchées de u aux points xi + h2 et xi h

2
respec-
2 2
tivement.

Exercice 1.1 1) Soit u une fonction classe C 2 au voisinage de x, Montrer que :

!
@u u(x + h) u(x) ju"(y)j
(x) sup h:
@x h y2[x;x+h[ 2

2) Si u est de classe C 3 au voisinage de x, alors :

!
@u u(x + h) u(x h) u(3) (y)
(x) sup h2 :
@x 2h y2]x h;x+h[ 6

Réponse :
2) u 2 C 3 au voisinage de x. Pour h su¢ samment petit, u 2 C 3 [x h; x + h]

12
 E¤ectuons un développement de Taylor au point x; il existe 1 et 2 dans ]0; 1[
2 3
u(x + h) = u(x) + hu0 (x) + h2 u00 (x) + h6 u(3) (x + 1 h)
2 3
u(x h) = u(x) hu0 (x) + h2 u00 (x) h6 u(3) (x + 2 h)
En soustrayant, on obtient :
u(x+h) u(x h) 2
2h
= u0 (x) + h12 u(3) (x + 1 h) + u(3) (x + 2 h)
u(3) continue, le théorème des valeurs intermédiaires nous assure qu’il existe un
point entre 1 et 2 que nous notons x + h et qui véri…é
(u(3) (x+ 1 h)+u(3) (x+ 2 h))
u(3) (x + h) = 2
ce qui nous permet de conclure :
0 1
(3)
u(x + h) u(x h) u (y)
u0 (x) h2 @sup A
2h 6
y2]x h;x+h[

Di¤érence divisée décentrée avant (aval) d’ordre deux :

Pour des pas réguliers de longueur h, le développement limité :

u(x + 2h) + 4u(x + h) 3u(x) h2 (3)

u0 (x) = + u ( x)
2h 3
donne la formule d’ordre 2 progressive suivante :

@u ui+2 + 4ui+1 3ui

u0 (xi ) = (xi ) ' :
@x 2h
Di¤érence divisée décentrée arrière (amont) d’ordre deux
Pour des pas réguliers de longueur h, le développement limité :

u(x 4u(x h) + 3u(x) h2 (3)

2h)
u0 (x) = + u ( x)
2h 3
On obtient également une formule régressive d’ordre deux :

@u 3ui 4ui 1 + ui 2
u0 (xi ) = (xi ) '
@x 2h

Dérivées secondes
u 2 C 4 au voisinage de x:Pour h su¢ samment petit, u 2 C 4 [x h; x + h] :
En utilisant la formule de Taylor au point x, il existe 1 et 2 dans ]0; 1[ tels
que :
2 3 4
u(x + h) = u(x) + hu0 (x) + h2 u00 (x) + h6 u(3) (x) + h24 u(4) (x + 1 h)
2 3 4
u(x h) = u(x) hu0 (x) + h2 u00 (x) h6 u(3) (x) + h24 u(4) (x 2 h)
En additionnant ces deux égalités, on obtient :
u(x+h) 2u(x)+u(x h) 4
h2
u00 (x) = h24 u(4) (x + 1 h) + u(4) (x 2 h)

13
Pr. A. Radid CHAPITRE 1. DIFFÉRENCES FINIES

Puisque u(4) est continue, le théorème des valeurs intermédiares nous assure
qu’il existe un point entre x 2 h et x + 1 h que nous notons x + h et qui véri…é
:
(4) u(4) (x + 1 h) + u(4) (x 2 h)
u (x + h) =
2
ce qui entraîne

u(x + h) 2u(x) + u(x h) h2 (4)

u00 (x) = u (x + h) (1.1)
h2 12
avec 2 ] 2 ; 1 [ ] 1; 1[
En négligeant le terme en h2 ; on obtient l’approximation de la dérivée seconde
centrée d’ordre 2
@2u ui+1 2ui + ui 1
2
(xi ) (1.2)
@x h2
En utilisant l’équation (1.1) on obtient aussi la majoration suivante :
0 1
(4)
@2u u(x + h) 2u(x) + u(x h) u (y)
2
(x) 2
h2 @sup A
@x h 12
y2]x h;x+h[

@u u(xi + h ) u(xi h
)
On peut retrouver le résultat (1.2) en utilisant : @x
(xi ) = 2
h
2

h2 (3)
24
u ( i)
@2u u0 (x + h ) u0 (x
h
h
) 2
@x2
(xi )
= i 2
h
i 2
24
u(4) ( i )
Approximation de la dérivée seconde de u :
@2u
@x2
(x) = u(x+h) 2u(x)+u(x
h2
h)
+ O(h2 ) ordre 2, centré
2
@ u u(x+2h) 2u(x+h)+u(x h)
@x2
(x) = h2
+ O(h) ordre 1, décentré avant
@2u u(x) 2u(x h)+u(x 2h)
@x2
(x) = h2
+ O(h) ordre 1, décentré arrière

Dérivées troisièmes
Di¤érence divisée décentrée avant (aval)

@3u ui+3 3ui+2 + 3ui+1 ui

u(3) (xi ) =
(x i ) '
@x3 h3
Di¤érence divisée décentrée arrière (amont)

@3u ui 3 + 3ui 2 3ui 1 ui

u(3) (xi ) = (xi ) '
@x3 h3
Di¤érence divisée centrée
On a deux formules possibles. L’une utilisant des points milieux de segments

14
 @3u ui+3=2 3ui+1=2 + 3ui 1=2 ui 3=2
u(3) (xi ) = 3
(xi ) '
@x h3
L’autre utilisant les noeuds du maillages.

@3u ui+2 2ui+1 + 2ui 1 ui 2

u(3) (xi ) =
3
(xi ) '
@x h3
Ces deux formules centrées sont d’ordre deux.

Dérivées quatrièmes
Di¤érence divisée centrée

d4 u ui+2 4ui+1 + 6ui 4ui 1 + ui 2

u(4) (xi ) = 4
(xi ) '
dx h4
Remarque 1.1 On peut observer que l’on retrouve dans ces formules les coe¢ -
cients du développement du binôme.

La dérivée en dimension 2
On considère maintenant u comme une fonction de deux variables de classe C 4 :
Le principe d’approximation des dérivées est exactement le même.
Si 4x et 4y tendent vers 0, nous avons :

@u u (x + 4x; y) u (x; y)
(x; y) = + O (4x)
@x 4x

@u u (x; y) u (x 4x; y)
(x; y) = + O (4x)
@x 4x

@u u (x + 4x; y) u (x 4x; y)
(x; y) = + O 4x2
@x 24x

@2u u (x + 4x; y) 2u (x; y) + u (x 4x; y)

(x; y) = + O 4x2
@x2 4x2
De même pour les dérivées par rapport y:

Exercice 1.2 Dérivées premières en dimension 2.

Donner l’expression des 3 formules usuelles d’approximation de la dérivée pre-
mière d’une fonction u en un point de grille (xi ; yj ); 0 < i < N; 0 < j < M:
Préciser la régularité requise de u et l’ordre d’approximation obtenu.

15
Pr. A. Radid CHAPITRE 1. DIFFÉRENCES FINIES

Exercice 1.3 Démontrer la formule aux di¤érences …nies à 5 points pour ap-
procher le laplacien d’une fonction u en un point de grille (xi ; yj ); 0 i N; 0
j M.
E¤ectuer une …gure pour illustrer les notations.
Préciser la régularité requise de u et l’ordre d’approximation obtenu.

Correction de l’exercice 1.2

On écrit les Développements limités de Taylor suivants:

h2 2 h3 3
u (xi + h; yj ) = u (xi ; yj ) + h@x u (xi ; yj ) + @xx u (xi ; yj ) + @xxx u (xi ; yj ) + O h4
2 3!

h2 2 h3 3
u (xi h; yj ) = u (xi ; yj ) h@x u (xi ; yj ) + @ u (xi ; yj ) @xxx u (xi ; yj ) + O h4
2 xx 3!

k2 2 k3 3
u (xi ; yj + k) = u (xi ; yj ) + k@y u (xi ; yj ) + @yy u (xi ; yj ) + @yyy u (xi ; yj ) + O k 3
2 3!

k2 2 k3 3
u (xi ; yj k) = u (xi ; yj ) k@y u (xi ; yj ) + @ u (xi ; yj ) @ u (xi ; yj ) + O k 3
2 yy 3! yyy
Et on obtient directement les formules aux D.F. recherchées:

Les formules aux di¤érences …nies décentrées à droite :

@u ui+1;j ui;j @u ui;j+1 ui;j
(xi ; yi ) et (xi ; yi ) ; ordre 1
@x 4x @y 4y

Les formules aux di¤érences …nies décentrées à gauche :

@u ui;j ui 1;j @u ui;j ui;j 1

(xi ; yi ) et (xi ; yi ) ; ordre 1
@x 4x @y 4y

Les formules aux di¤érences …nies centrées :

@u ui+1;j ui 1;j @u ui;j+1 ui;j 1
(xi ; yi ) et (xi ; yi ) ; ordre 2
@x 24x @y 24y

Pour les formules centrées, il su¢ t de soustraire les équations deux à deux.

16
 Opérateurs d’ordre deux en dimension 2

@2u u (x + 4x; y) 2u (x; y) + u (x 4x; y)

2
(x; y) = + O 4x2
@x 4x2
De même pour les dérivées par rapport y:
La formule aux di¤érences …nies à 5 points pour approcher le laplacien s’écrit
:
1 1
4u (xi ; yj ) [ui+1;j 2ui;j + ui 1;j ] + [ui;j+1 2ui;j + ui;j 1 ] ; ordre 2
4x2 4y 2
(1.3)

Exercice 1.4 1) Démontrer la formule aux di¤érences …nies à 5 points pour ap-
procher le laplacien (1.3) d’une fonction u en un point de grille (xi ; yj ); 0 i
N; 0 j M .
E¤ectuer une …gure pour illustrer les notations.
Préciser la régularité requise de u et l’ordre d’approximation obtenu.
2) Discrétiser l’opérateur di¤érentiel div ( ru) ; où la fonction est donnée,
a-priori dépendante de x.
h i
1
div ( ru) 1 (ui+1;j ui;j ) + i 1 ;j (ui+1;j ui;j )
4x2
h i+ 2 ;j 2 i
1
+ 4y2 i;j+ 1 (ui;j+1 ui;j ) + i;j 1 (ui;j 1 ui;j )
2 2

Corrigé de l’exercice 1.4

1) On écrit les D.L. de Taylor jusqu’à l’ordre 4 (u doit être de classe C 4 ) :

@ h2 @ 2 h3 @ 3 h4 @ 4
u (xi + h; yj ) = u (xi ; yj ) +h u (xi ; yj ) + u (x ;
i jy ) + u (x ;
i jy ) + u (xi ; yj ) +O h5
@x 2 @x2 3! @x3 4! @x4
@ h2 @ 2 h3 @ 3 h4 @ 4
u (xi h; yj ) = u (xi ; yj ) h
u (xi ; yj ) + u (xi ; yj ) u (x ;
i jy ) + u (xi ; yj ) +O h5
@x 2 @x2 3! @x3 4! @x4
Puis de même selon l’axe des y
@ k2 @ 2 k3 @ 3 k4 @ 4
u (xi ; yj + k) = u (xi ; yj ) +k u (xi ; yj ) + u (x i ; y j ) + u (x i ; y j ) + u (xi ; yj ) +O k 5
@y 2 @y 2 3! @y 3 4! @y 4
@ k2 @ 2 k3 @ 3 k4 @ 4
u (xi ; yj k) = u (xi ; yj ) k u (xi ; yj ) + 2
u (xi ; yj ) 3
u (xi ; yj ) + 4
u (xi ; yj ) +O k 5
@y 2 @y 3! @y 4! @y
Par combinaison linéaire (on ajoute les deux équations), on obtient la formule
D.F. à 5 points qui approche la valeur de 4u (xi ; yj ) ; formule d’ordre 2 car en
O(h2 + k 2 ):
2) Voir Exercice 1 des TD1.

17
Pr. A. Radid CHAPITRE 1. DIFFÉRENCES FINIES

Dérivées croisées
@ u 2
Déterminons une approximation de la dérivée croisée @x@y de la fonction de 2
variables u(x; y).
La discrétisation du domaine de calcul est bidimensionnelle et fait intervenir
deux pas d’espace supposés constants 4x et 4y dans les directions x et y .
La principe est toujours basé sur les développements de Taylor :
@u @u (l4x)2 @2u
u(xi+l ; yj+m ) = u (xi ; yj ) + l4x @x i
+ m4y @y
+ 2 @x2
j i
2
@2u @2u
+ (m4y)
2 @y 2
+ 2ml4x4y
2 @x@y
+
j i;j

Exercice 1.5 Montrer que

@2u ui+1;j+1 ui+1;j 1 ui 1;j+1 + ui 1;j 1

=
@x@y i;j 44x4y
Corrigé:
@u @u
u (xi+1 ; yj+1 ) = u(xi ; yj ) + 4x (xi ; yj ) + 4y (xi ; yj ) + (1.4)
@x @y
4x2 @ 2 u 4y 2 @ 2 u 4x4y @ 2 u
(x ;
i jy ) + (x ;
i jy ) + 2 (xi ; yj ) + :::
2! @x2 2! @y 2 2 @x@y

@u @u
u (xi+1 ; yj 1 ) = u(xi ; yj ) + 4x (xi ; yj ) 4y (xi ; yj ) + (1.5)
@x @y
4x2 @ 2 u 4y 2 @ 2 u 4x4y @ 2 u
(x ;
i jy ) + (x ;
i jy ) 2 (xi ; yj ) + :::
2! @x2 2! @y 2 2 @x@y

@u @u
u (xi 1 ; yj+1 ) = u(xi ; yj ) 4x (xi ; yj ) + 4y (xi ; yj ) + (1.6)
@x @y
4x2 @ 2 u 4y 2 @ 2 u 4x4y @ 2 u
(x ;
i jy ) + (x ;
i jy ) 2 (xi ; yj ) + :::
2! @x2 2! @y 2 2 @x@y

@u @u
u (xi 1 ; yj 1 ) = u(xi ; yj ) 4x (xi ; yj ) 4y (xi ; yj ) + (1.7)
@x @y
4x2 @ 2 u 4y 2 @ 2 u 4x4y @ 2 u
(x ;
i jy ) + (x ;
i jy ) + 2 (xi ; yj ) + :::
2! @x2 2! @y 2 2 @x@y
On fait le calcul de (1.4)-(1.5)-(1.6)+(1.7) et on obtient : ui+1;j+1 ui+1;j 1
@2u
ui 1;j+1 + ui 1;j 1 = 44x4y @x@y
i;j

18
1.2.4 Stencil d’un schéma
Dé…nition 1.3 On appelle stencil d’un schéma le nombre de points nécessaire à
l’écriture du schéma.

1.2.5 Schéma à multipas

Dé…nition 1.4 Un schéma aux di¤érences …nies est dit schéma à p pas en temps
si les valeurs un+1
j des solutions approchées au temps tn+1 sont fonctions des
valeurs aux p instants précédents, soit aux temps tn , tn 1 , ...tn p+1 . En partic-
ulier, un schéma est dit à un pas si les un+1
j ne dépendent que des unj .

1.2.6 L’ordre d’un schéma

Dé…nition 1.5 L’ordre du schéma mesure la précision ou erreur de troncature
mathématique commise en remplaçant les dérivées partielles exactes par leurs ap-
proximations sous formes de di¤érences divisées. L’ordre est déterminé par des
développements de Taylor obtenus en injectant dans l’écriture du schéma numérique
la fonction solution continue exacte du problème di¤érentiel.

Dé…nition 1.6 Un schéma aux di¤érences …nies est d’ordre p en temps et d’ordre
q en espace si la di¤érence entre l’équation et le schéma appliqué à la fonction
solution du problème continu est un in…niment petit d’ordre p en temps et d’ordre
q en espace.

Exercice 1.6 On considère le schéma saute-mouton pour approcher l’équation de

transport
@u @u
+a = 0 (a > 0)
@t @x
sur une grille spatiale régulière de pas 4x et avec un pas de temps 4t :

un+1
j = unj 1
unj+1 unj 1

4t
où n 1 et u0j et u1j sont donnés et où on a posé = a 4x :
Dessiner le stencil du schéma
De quel type de schéma s’agit-il ?
Quel est l’ordre de ce schéma ?

Réponse : Le stencil su schéma :

2) c’est un schéma explicite à deux pas.

3) On a :
un+1
j = unj 1 unj+1 unj 1

19
Pr. A. Radid CHAPITRE 1. DIFFÉRENCES FINIES

4t
or = a 4x donc

un+1
j unj 1
unj+1 unj 1
+a =0
24t 24x
Le schéma étant écrit au temps tn et au point xj . Faisons des développe-
ments limités autour de ces valeurs. On obtient par développement de Taylor les
approximations suivantes :

un+1
j unj 1
@u
= (xj ; tn ) + O 4t2
24t @t
unj+1 unj 1 @u
= (xj ; tn ) + O 4x2
24x @x
Donc le schéma est d’ordre deux en temps et en espace, comme peut le montrer
sa réécriture sous forme de di¤érences centrées

1.2.7 Discrétisation des conditions aux limites

Les di¤érences …nies consistent à approcher les opérateurs de dérivation par des
opérateurs discrets de dérivation et il faut de plus ajouter des conditions aux limites
qui peuvent être de plusieurs types.. Par exemple, pour d = 1 :

20
 Si l’équation est dé…nie sur un domaine borné, par exemple x 2 [ ; ], le
maillage sera restreint à cet intervalle c’est à dire j 2 f0; 1; : : : ; N g avec x0 = et
xN =
Par exemple, si on a des conditions aux limites de Dirichlet

u( ; t) = L; u( ; t) = R; 8t 2 R+

elles se traduisent au niveau discret en

un0 = L; unN = R

Si on a des conditions de Neumann

@x u( ; t) = L; @x u( ; t) = R; 8t 2 R+

elles se traduisent au niveau discret en

un1 un0 un unN
= L; N 1
= R; 8n;
4x 4x

Si on a des conditions périodiques

u (x + ; t) = u (x + ; t) ; 8t

elles se traduisent au niveau discret en

un0 = unN ; 8n

et plus généralement unj = unN +j :

Exercice 1.7 Soit f une fonction continue sur [0; 1] et c 2 C ([0; 1] ; R+ )

1. Déterminer la solution explicite de l’équation

u" + c(x)u = f dans ]0; 1[ (1.8)

associée aux conditions de Dirichlet homogènes au bord :

u (0) = u (1) = 0

2. Déterminer la solution explicite de l’équation (1.8) associée aux conditions

de Dirichlet:
u (0) = ; u (1) = ; pour ; 2 R

21
Pr. A. Radid CHAPITRE 1. DIFFÉRENCES FINIES

3. Même question en considérant les conditions de Newmann:

u0 (0) = ; u0 (1) = :

4. Ecrire le problème 2) sous forme matricielle : Ah uh = bh et montrer que la

matrice Ah est inversible

Corrigé

1. Soit N un entier, on subdivise l’intervalle [0; 1] en seegments de longueur

h = N 1+1 : Soit xi = ih pour i variant de 0 à N + 1.
On va chercher une solution approchée (u0 ; : : : ; uN +1 ) dé…nie en les xi : ui =
u(xi )
les conditions initiales u (0) = 0; u (1) = 0 implique u0 = 0; uN +1 = 0:
On cherche uh = (u1 ; : : : ; uN )T
Sachant que u00 (xi ) + c (xi ) u (xi ) = f (xi ) pour i = 1; : : : ; N
On obtient en utilisant l’approximation de la dérivée seconde de u

u (xi+1 ) 2u(xi ) + u (xi 1 ) h2 (4)

+ c (x i ) u (x i ) + u (xi + i h) = f (xi )
h2 12
ui+1 2ui + ui 1
+ ci ui = fi pour i = 1; : : : ; N
h2
Pour i = 1
1
( u2 + 2u1 ) + c1 u1 = f1
h2
Pour i = N
1
( uN 1 + 2uN ) + cN uN = fN
h2
2. les conditions initiales u (0) = ; u (1) = implique u0 = ; uN +1 = :
Le schéma explicite est le même comme dans 1) :
ui+1 2ui + ui 1
+ ci ui = fi pour i = 1; : : : ; N
h2

Pour i = 1
u2 2u1 + u0
+ c1 u1 = f1
h2
1
( u2 + 2u1 ) + c1 u1 = f1 + 2
h2 h

22
 Pour i = N
uN +1 2uN + uN 1
+ cN uN = fN
h2
1
( uN 1 + 2uN ) + cN uN = fN +
h2 h2
3.
ui+1 2ui + ui 1
+ ci ui = fi pour i = 1; : : : ; N
h2
u0 (0) = ; u0 (1) = :
Approximation d’ordre 1
u(h) u(0)
= u0 (0) = h
+ O (h) ) u1 = u0 + h
u(1) u(1 h)
= u0 (1) = h
+ O (h) ) uN +1 = uN + h
Approximation d’ordre 2
h2 00
u(h) = u(0) + hu0 (0) + 2
u (0) + O (h3 )
u00 (x0 ) + c (x0 ) u (x0 ) = f (x0 ) c’est à dire u00 (0) = c0 u0 f0
2 h2 h2
u1 = u0 + h + h2 (c0 u0 f0 ) + O (h3 ) c’est à dire 1 + 2 0
c u0 u1 = 2 1
f

De même on trouve uN = uN +1 + 12 h2 (cN +1 un+1 fn+1 ) :

4. Matriciellement, le problème 2) s’écrit :

Ah u h = b h (1.9)

où
0 1 0 1
2 + c1 h2 1 f1 + h2
B 1 2 + c2 h2 1 C B C
1 B
B . . ...
C
C
B
B ..
C
C
Ah = 2 B .. .. C ; bh = B . C
h B C B C
@ 1 2 + cN 1 h2 1 A @ A
1 2 + cN h2 fN + h2

et

t
Uh = u1 ; : : : ; uN
Supposons que c 0. Montrer que la matrice Ah est symétrique dé…nie positive.
La matrice Ah est clairement symétrique. Montrons que Ah est dé…nie positive.

23
Pr. A. Radid CHAPITRE 1. DIFFÉRENCES FINIES

Soit v un vecteur de RN ; de composantes v1 ; : : : ; vN : On a :

(0)
X
N
(0)
t t
v Ah v = v Ah v + ci vi2 v t Ah v
i=1

or

(0) 1
P
N
v t Ah v = h2
vi ( vi 1 + 2vi vi+1 )
i=1
Par changement d’indice
" N #
1 X XN X
N +1
t (0)
v Ah v = 2 ( vi 1 vi ) + vi2 vi 1 vi
h i=1 i=1 i=2

et comme on a posé v0 = vN +1 = 0; on peut écrire :

1X 2 1X
N N
t (0)
v Ah v = 2 2v + ( 2vi 1 vi )
h i=1 i h2 i=1

(0) 1
P
N P
N P
N
v t Ah v = h2
vi2 + vi2 2vi 1 vi
i=1 i=1 i=1
1
PN NP+1 PN
= h2
vi2 + vi2 1 2vi 1 vi
i=1 i=1 i=1
1
PN
= h2
vi2 2vi 1 vi + vi2 1 + vN 2
i=1
PN
= 1
h 2 (vi vi 1 )2 + vN 2
i=1
NP +1
= (vi vi 1 )2 0
i=1

Donc 8v = (v1 ; : : : ; vN )t 2 RN ; Ah v:v 0 donc Ah est positive.

Si on suppose Ah v:v = 0; on a alors
X
N
ci h2 vi2 = 0 et vi = vi 1 8i = 1; : : : ; N + 1
i=1

On a donc v1 = v2 = = vN = vN +1 = 0:
Remarquons que ces égalités sont véri…ées même si les ci sont nuls.
Donc la matrice Ah est bien dé…nie.
La matrice Ah est symétrique dé…nie positive donc elle est inversible, ce qui
entraine l’existence et l’unicité de la solution.

24
Remarque 1.2 On peut aussi monter que A est inversible en calculant ses valeurs
propres ou en véri…ant que A est une matrice à diagonale fortement dominante et
irréductible.

1.3 Consistance, Stabilité et Convergence

Un certain nombre de notion est nécessaire lors de la résolution d’équations aux
dérivées partielles (EDP) au moyen de leurs équivalents discrétisés. Les trois prin-
cipales sont la convergence, la stabilité et la consistance. Ces trois propriétés
permettent de relier la solution exacte des équations continues à la solution exacte
des équations discrétisées et à la solution numérique obtenue.
Ces di¤érents liens sont :

la stabilité, c’est la propriété qui assure que la di¤érence entre la solu-

tion numérique obtenue et la solution exacte des équations discrétisées, est
bornée.

la consistance, c’est la propriété qui assure que la solution numérique des

équations discrétisées tende vers la solution exacte des équations continues
lorsque le pas de discrétisation (4t; 4x) tendent vers zéro.

la convergence, c’est la propriété qui assure que la solution numérique tende

vers la (ou une) solution exacte des équations continues. C’est évidemment
la propriété la plus recherchée !

1.3.1 Erreur de Troncature

Notons S x; t u(xj ; tn ) l’application d’un schéma aux di¤érences …nies à la solution
continue u .

Dé…nition 1.7 on appelle erreur de troncature E.T, l’ensemble des termes nég-
ligés dans les développements limités lors de l’obtention d’une équation (ou schéma)
discrète.

Il est en e¤et possible d’écrire :

Équation continue = Équation discrète + E.T

Remarque 1.3 Le calcul des erreurs de troncature est usuellement basé sur les
développements de Taylor.

25
Pr. A. Radid CHAPITRE 1. DIFFÉRENCES FINIES

1.3.2 Consistance
Un schéma aux di¤érences …nies est consistant avec l’équation aux dérivées par-
tielles, si l’erreur de troncature du schéma tend vers zéro, lorsque tous les pas de
discrétisation tendent vers zéro.

lim E:T: = 0
(4x;4t)!0
4t
Des problèmes peuvent se poser si l’erreur de troncature varie comme 4x . Dans
4t
ce sas, on est obligé de ra¢ ner de sorte que 4x ! 0: On dit dans ce cas le schéma
a une consistance conditionnelle.

Remarque 1.4 Si on n’approche pas correctement le problème continu, on n’aura

aucune chance d’obtenir la convergence du schéma. La qualité de la consistance
n’est autre que la précision du schéma.

1.3.3 Stabilité
Autre les outils qui permettent de comparer les performances des di¤érents sché-
mas, on doit également choisir les pas 4x et 4t de sorte que le schéma correspon-
dant donnera une solution approchée correcte, au sens où une petite perturbation
de la donnée initiale n’induira par une perturbation trop grande sur la solution
calculée au temps …nal.
Soit U n = unj 1 j N la solution numérique d’un schéma au temps tn .
Nous pouvons exprimer le vecteur U n+1 au temps tn+1 en fonction de U n . par
:

U n+1 = C (4t; 4x) U n

où C est une matrice caractérisant le schéma et dépendant des pas de temps
et d’espace.
On en déduit :
U n = C nU 0
où U 0 est le vecteur des conditions initiales.
Critère de stabilité au sens des normes
Introduisons les normes L1 et L2 discrètes suivantes. Pour tout v 2 RM on
note
kvk1 = max jvi j
1 i M

!1=2
X p
kvk2;h = h jvi j2 = h kvk2
1 i M

26
 Dé…nition Un schéma aux di¤érences …nies est dit stable pour la norme k:k
s’il existe une constante K > 0 indépendante de 4x et 4t ( lorsque ces valeurs
tendent vers zéro) telle que quelque soit la donnée initiale U 0

kU n k K U0 pour tout n 0

Si cette inégalité n’a lieu que pour des pas 4x et 4t restreints à certaines inégalités,
on dit que le schéma est conditionnellement stable.

Remarque 1.5 La stabilité est un concept qui s’applique sur les schémas numériques
et pas sur les équations continues.

Di¤érents types de stabilité

Un schéma peut être :

Inconditionnellement stable : Quels que soient 4t et 4x les erreurs

causées par le schéma numérique n’explose pas au …l des itérations.

Conditionnellement stable : On doit poser une condtion sur 4t et 4x

pour que la solution n’explose pas.

Inconditionnellement instable : Quels que soient 4t et 4x les erreurs s’ampli…ent

au …l des itérations. Ceci cause des résultats complétement faux.

Stabilité d’un schéma linéaire:

La stabilité d’un schéma linéaire à deux niveaux est facile à interpréter. En ef-
fet, par linéarité tout schéma linéaire à deux niveaux peut s’écrire sous la forme
condensée,
AU n = U n+1
où A est une matrice (dite d’itération) et on obtient An U 0 = U n+1 , (la notation
An désigne ici la puissance n-ème de A) et par conséquent la stabilité du schéma
est équivalente à

An U 0 K U 0 ; 8n 0; 8U 0 2 RN
kM uk
Introduisant la norme matricielle subordonnée kM k = sup kuk
, la sta-
u2RN 1 ;u6=0

bilité du schéma est équivalente à

kAn k K ; 8n 0

qui veut dire que la suite des puissances de A est bornée.

27
Pr. A. Radid CHAPITRE 1. DIFFÉRENCES FINIES

La stabilité de majorations de la norme de la matrice A souvent obtenues par

le calcul de ses valeurs propres.
La condition nécessaire et su¢ sante de stabilité pour les équations aux dif-
férences s’énonce kAk 1;
Remarque : Si kAk 1 alors (A) 1 car (A) kAk
Si la matrice A symètrique réelle alors kAk2 = (A):

1.3.4 Convergence
La convergence d’un schéma aux di¤érences …nies est une propriété naturelle qui
assure que, pour des valeurs su¢ samment petites des pas d’espace et de temps,
la solution numérique calculée sera proche de la solution exacte du problème de
départ.

Dé…nition 1.8 Le schéma aux di¤érence …nies utilisé pour la résolution numérique
de l’équation aux dérivées partielles est convergent si, pour toute solution u , la
suite unj converge vers u(xj ; tn ) avec (4x; 4t) ! (0; 0) :

Théorème 1.1 (Théorème de Lax) Pour un problème linéaire, la consistance et

la stabilité sont nécessaires et su¢ santes pour assurer la convergence.

Remarque 1.6 La consistance et la stabilité d’un schéma sont en général beau-

coup plus facile à étudier que sa convergence.

1.4 Equation de Laplace 1D

On considère le problème aux limites 1D suivant :

u00 = f; x 2 ]0; 1[
(1.10)
u (0) = u (1) = 0:

Si u 2 C 4 alors
2u (x) u (x + h) u (x h) h2 (4)
u00 (x) = + u (x + h) avec 2 ] 1; 1[
h2 12
(1.11)
On note u(x) la fonction inconnue, solution de l’EDP. Le domaine spatial des
EDP qu’on va résoudre est [0; 1] (une dimension spatiale).
Schéma de di¤érences …nies :
Discrétisation du domaine : On le discrétise en utilisant un pas d’espace h ; les
points de discrétisation sont x0 ; :::; xN +1 , avec xi = i h; on a donc 1 = (N +1)h: On
note ui l’approximation numérique de la solution u en xi ; et on note fi = f (xi ) :

28
 On obtient les équations discrètes en approchant u00 (xi ) par le quotient dif-
férentiel par développement de Taylor (1.11),

2u (xi ) u (xi + h) u (xi h) h2 (4)

u00 (xi ) = + u (xi + h) avec 2 ] 1; 1[
h2 12
(1.12)
le problème (1.10) devient :

2u (xi ) u (xi + h) u (xi h) h2 (4)

+ u (xi + h) = f (xi ); i = 1; : : : ; N; (1.13)
h2 12
avec u(x0 ) = 0; u(xN +1 ) = 0:
En négligeant le terme en h2 et en approchant alors u (xi ) par ui ; on obtient
le schéma discrétisé suivant :
1
(ui+1 2ui + ui 1 ) = fi ; i = 1; : : : ; N; (1.14)
h2
avec u0 = uN +1 = 0:
le schéma (1.14) est un schéma à 3 points (ui 1 ; ui ; ui+1 ):
Écriture matricielle :
On peut écrire les équations (1.14) sous forme matricielle :

A(N ) U N = F (N ) (1.15)

avec
0 1 0 1 0 1
2 1 0 u1 f1
B
1 B 1 2 1 C B .. C B .. C
C B C B C
A(N ) = 2B ... C 2 MN (R) ; U N = B . C et F N = B . C
h @ 1 A @ A @ A
0 1 2 uN fN

Exercice 1.8 Montrer que la matrice A(N ) symétrique dé…nie positive.

La matrice A(N ) est évidemment symétrique. Montrons qu’elle est dé…nie pos-
itive.
A est symétique,
Notons < :; : > le produit scalaire euclidien.
Soit v = (v1 ; v2 ; : : : ; vN )t 2 Rn ; v 6= 0;
on pose v0 = vN +1 = 0: Calculons le produit scalaire A(N ) v; v = v t A(N ) v: On
a :

1
P
N
A(N ) v; v = h2
vi ( vi 1 + 2vi vi+1 )
i=1

29
Pr. A. Radid CHAPITRE 1. DIFFÉRENCES FINIES

Par changement d’indice

" N #
1 X XN X
N +1
A(N ) v; v = 2 ( vi 1 vi ) + 2vi2 vi 1 vi
h i=1 i=1 i=2

et comme on a posé v0 = vN +1 = 0; on peut écrire :

1X 2 1X
N N
(N )
A v; v = 2 2v + ( 2vi 1 vi )
h i=1 i h2 i=1

1
P
N P
N P
N
A(N ) v; v = h2
vi2 + vi2 2vi 1 vi
i=1 i=1 i=1
1
PN NP
+1 P
N
= h2
vi2 + vi2 1 2vi 1 vi
i=1 i=1 i=1
PN
= + h12 vi2 2vi 1 vi + vi2 1
2
+ vN
i=1
P
N
= + h12 (vi vi 1 )2 + vN
2
i=1
NP+1
= + h12 (vi vi 1 )2 0
i=1

Donc 8v = (v1 ; : : : ; vN )t 2 RN ; A(N ) v:v 0 donc A(N ) est positive.

Si on suppose A(N ) v:v = 0; on a alors

vi = vi 1 8i = 1; : : : ; N + 1

On a donc v1 = v2 = = vN = vN +1 = 0:
(N )
Donc la matrice A est bien dé…nie.
(N )
La matrice A est symétrique dé…nie positive donc elle est inversible, ce qui
entraine l’existence et l’unicité de la solution.
Etude de la consistance du schéma:
L’erreur de consistante ici est donc l’erreur qu’on commet en remplaçant l’opérateur
u"par le quotient 2u(xi ) u(xih+h)
2
u(xi h)
:
h2 (4)
on note : i = 12 u (xi + h) ; d’après (1.13), on obtient :

h2 (4) h2
j ij = u (xi + h) max u(4) (x) ; 8i = 1; : : : ; N
12 12 0 x 1
lorsque h tend vers 0 l’erreur de troncature i tend aussi vers 0 donc le schéma est
consistant à l’ordre 2.

30
Exercice 1.9 Calculer les valeurs propres et les vecteurs propres de la matrice
0 1
2 1 0
B ..
. C
B 1 2 C
A=B ... .. C:
@ . 1 A
0 1 2
.

Rappel :

Dé…nition 1.9 Soit M = (aij ) 2 Mn (K) : Pour tout i 2 [1; n] on nomme :

! ( )
X
n Xn
Di = Bf aii ; jaij j = x 2 K= jaii xj jaij j
j=1;j6=i j=1;j6=i

la ième ligne de Gerschgorin.

Théorème 1.2 (Théorème de Gerschgorin ) Soit M 2 Mn (K) ; l’ensemble des

valeurs propres est inclus dans la réunion des disques de Gerschgorin. C’est à dire
: n
est une valeur proprede M =) 2 [ Di :
i=1

Corrigé de l’exercice 1.9:

Cherchons les valeurs propres de la matrice A 0 1
0 1 0
B ..
. C
B 1 0 C
En écrivant la matrice A sous la forme : A = 2I B avec B = B ... ... C
@ 1 A
0 1 0
la matrice B est symétrique réelle, donc toutes ses valeurs propres réelles. Avec
le théorème de Gerchgörin, on déduit que pour toute valeur valeur propre 2 R
on a j j 2:Une telle valeur propre peut s’écrire = 2 cos avec 2 [0; ] : Si x
est un vecteur propre non nul associé à la valeur propre ses composantes sont
alors solutions de la récurrence :

xk 1 xk + xk+1 = 0 (1 k N) (1.16)

avec les conditions aux limites x0 = 0 et xN +1 = 0:

Le polynôme caractéristique de (1.16) est :

r2 2 cos r + 1 = r ei r e i

31
Pr. A. Radid CHAPITRE 1. DIFFÉRENCES FINIES

les racines sont donc r1 = ei ; r2 = e i

et on a :

xk = c1 eik + c2 e ik

De x0 = 0; on déduit que c2 = c1 et de xN +1 = 0 avec (c1 ; c2 ) 6= (0; 0) ; on

déduit que sin (N + 1) = 0 et = Nj+1 avec 1 j N:
Les valeurs propres de A sont donc : k = 2 2 cos Nk+1 = 4 sin2 2(Nk +1) pour
1 k N:
Les valeurs propres de A sont : k = 2 1 cos Nk+1 ; avec k = 1; : : : ; N ; et
les vecteurs propres sont : xk = sin Njk+1 j=1;:::;N
Etude de la stabilité
Rappel :
On rappelle que par dé…nition, si M = (mij ) 2 MN (K), v 2 KN
kM vk
Norme matricielle induite ou subordonnée : kM k = sup kvk 1 = sup kM vk1 ;
1
v2RN ;v6=0 v2RN ;kvk1 =1
avec kvk1 = sup jvi j
1 i N
En particulier
X
N
kM k1 = max jmij j
1 i N
j=1
r P
kvkp = jvi jp ; 8p 1:
1 i N
P
kM k1 = max jaij j ;
1 j N1 i N
P
kM k1 = max jaij j ;
1 i N1 j N
p
kM k2 = (M M ):

Pour toute norme matricielle sur MN (C) : (M ) kM k :

Pour tout " > 0, il existe au moine une norme matricielle subordonnée telle
que : kM k (M ) + ":
1
Si une matrice A est inversible et symétrique, alors la matrice A et symétrique.
1
Le spectre de la matrice A est l’ensemble des inverses du spectre de A:

La matrice A(N ) = h12 A donc les valeurs propres de A(N ) sont : h22 1 cos Nk+1 :
la valeur propre la plus petite en valeur absolue est 1 = h22 (1 cos ( h)) et
développement de Taylor, on a
2 2
2 h
1 = 2 + O h2 = 2
+ O h2 donc min = 1 ' 2
:
h 2

32
Par conséquent,
A(N ) v 1
A(N ) 1
= sup 1
= ;
v2RN ;v6=0 kvk1 1

1 1
A(N ) 2
1

1 1
UN 1
= A(N ) F (N ) 2
kF k1
1

D’où la stabilité de la méthode.

Majoration de l’erreur :

1 1
EN A(N ) (N )
2 2
(N )
2
2

(N ) 1 2 (4) (N )
où j la jème composante de (N ) s’écrit : j = 12 hu j
N 2
Donc E Ch ; d’où la convergence du schéma numérique.

1.5 Equation de la chaleur 1D

1.5.1 Schéma explicite de l’équation de la chaleur 1D
Considérons l’équation de la chaleur suivante :

8
@u @2u
>
< @t = D @x2 (x; t) 2]0; 1[ ]0; T [
u(x; 0) = (x) donnée : condition initiale (1.17)
>
:
u(0; t) = u(1; t) = 0 8t conditions aux limites de Dirichlet homogène

D coe¢ cient de di¤usion thermique

On peut considérer qu’il s’agit de l’évolution de la température d’une barre
soumise aux extremités à une température nulle.
Pour résoudre cette équation numériquement, on doit discrétiser les variables
x et t:

Discrétisation de l’équation de la chaleur

Construisons un maillage régulier. on note xj = j4x; j = 0; : : : ; N +1 et tn = n4t:

On note unj l’approximation de u(xj ; tn ).

33
Pr. A. Radid CHAPITRE 1. DIFFÉRENCES FINIES

En utilisant pour la dérivée temporelle une approximation à droite, on obtient

:

@u un+1
j unj
(xj ; tn ) '
@t 4t
et une approximation centré pour la dérivée seconde en espace :

@2u unj+1 2unj + unj 1

(x ;
j nt ) '
@x2 4x2
Le schéma se pose alors :
8 n+1 n
> uj uj un 2un n
j +uj
>
< 4t = D
j+1
4x2
1
j = 1; : : : ; N
u0j = u0 (xj ) j = 1; : : : ; N (1.18)
>
>
:un = un = 0 8n 2 N
0 N +1

Ce schéma est un schéma à un pas, car le vecteur des solutions approchées au

temps tn+1 ne dépend que des solutions approchées au temps tn . C’est un schéma
explicite car il donne une formule explicite de calcul de la solution au temps tn+1
en fonction des valeurs de la solution au temps précédent. Il n’y a pas d’équation
à résoudre pour obtenir la valeur au nouvel instant tn+1

Consistance du schéma :
Exercice : Montrer que le schéma d’Euler explicite centré en espace pour l’équation
de la chaleur est d’ordre un en temps et d’ordre deux en espace.
On condère l’équation de la chaleur suivante :

@u @2u
D =0
@t @x2
Corrigé
En développant un+1
j et unj par développement de Taylor, soustrayant les résul-
tats et divisant le tout par 4t,

34
 @u un+1
j unj 4t @ 2 u (xj ; tn )
= O (4t)2
@t 4t 2 @t2

En développant unj+1 et unj 1 par développement de Taylor,additionnant les

résultats, on obtient

@2u (unj+1 2unj + unj 1 ) 4x2 @ 4 u (xj ; tn )

= O 4x4
@x2 4x2 12 @x4

2
h i
@u u(xj ;tn+1 ) u(xj ;tn ) u(x ;t ) 2u(xj ;tn )+u(xj 1 ;tn )
@t
D @@xu2 4t
D j+1 n 4x2
=
2 2 4
4t @ u(xj ;tn ) @ u(xj ;tn )
2 @t2
+ D 4x12 @x4
+ O (4t)2 + O (4x4 )
u(xj ;tn+1 ) u(xj ;tn ) u(xj+1 ;tn ) 2u(xj ;tn )+u(xj 1 ;tn )
On pose : S x; t u(xj ; tn ) = 4t
D 4x2
donc
@ @2
u (xj ; tn ) u (xj ; tn ) S x; t u(xj ; tn ) = O (4t) + O (4x2 )
D
|@t @x2 {z }
= Erreur de Trancature
lorsque 4t et 4x tendent vers zéro, l’erreur de trancature tend aussi vers 0;
le schéma est donc consistant.
on retrouve l’erreur de troncature sur le temps, O (4t). En suivant le même
principe avec le terme sur l’espace, on trouve une erreur en O (4x2 ) :
L’ordre de l’erreur de troncature est O(4t; 4x2 )
D4t n
un+1
j = u +
4x2 j 1
1 2 D4t
4x2
unj + D4t n
u pour
4x2 j+1
j = 1; : : : ; N
en posant = D4t
4x2
, on obtient :

un+1
j = unj 1 + (1 2 ) unj + unj+1 pour j = 1; : : : ; N (1.19)

On connaît les u0j , on en déduit les u1j et ainsi de suite jusqu’aux unj . Ce
schéma est dit explicite car les un+1
j se déduisent directement des unj . Dans ce
genre de schéma, aucune système d’équation n’est à résoudre (pas de matrice à
inverser numériquement!!!), le nouveau pas se déduit des pas précédement calculés.

Formulation matricielle du schéma

Soit sous forme matricielle :

U n+1 = CU n

35
Pr. A. Radid CHAPITRE 1. DIFFÉRENCES FINIES

0 1 0 10 1
un+1
1
1 2 0 0 un1
B un+1 C B 1 2 0 CB un2 C
B 2 C B .. CB C
B .. C B ..
.
..
.
..
. CB .. C
B . C = B 0 . CB . C
B n+1 C B ... CB n C
@ uN 1 A @ 0 1 2 A @ uN 1 A
un+1
N 0 ::: 0 1 2 unN
n
= (I A) U
avec 0 1
2 1 0
B ... C
B 1 2 C
A=B .. .. C
@ . . 1 A
0 1 2

Analyse de Stabilité du schéma explicite de l’équation de la chaleur

Etude matricielle de la stabilité La matrice C = I D4t 4x2
A est une matrice
symétrique réelle. Ces vecteurs propres sont ceux de A (Voir Exercice 1.9), ses
valeurs propres sont égales à :

D4t k
k =1 k avec k =2 2 cos
4x2 N +1

Les vecteurs propres de C associés aux valeurs propres k sont les vecteurs vk =
(sin Nk+1 ); sin N2k+1 ; : : : ; sin N
Nk
+1
avec 1 k N:
Majorons la norme euclidienne de U n+1

4t
U n+1 2
I D A kU n k2
4x2 2

Comme la matrice C est une matrice symétrique, sa norme euclidienne est égale
à son rayon spectral donc

D4t D4t D4t 2 k

kCk2 = I A = I A = max 1 4 sin ;
4x2 2 4x2 1 k N 4x2 2 (N + 1)
D4t 1
la condition de stabilité kCk2 1 se traduit donc par : 0 4x2 2
donc le
schéma explicite est donc conditionnellement stable.
Cette condition D4t
4x2
1
2
est la condition classique de stabilité du schéma
d’Euler explicite pour l’équation de la chaleur. Elle impose des pas de temps
très petits, ce qui condamne pratiquement l’usage de ce schéma explicite pour les
problèmes paraboliques.

36
 Si D4t
4x2
> 1
2
: les valeurs propres de C sont les k = 1 4 D4t
4x2
sin2 k
2(N +1)
; 1
k N:

D4t 2 N
N =1 4 sin = 1 4 +O(4x2 ) = (4 1)+O(4t); 4 1>1
4x2 2 (N + 1)
T
0 n 4t
T
Si n0 = E 4t
! 1 alors :

C n0 v N 1 n0
=j Nj = (4 1)n0 exp (n0 O (4t)) ! 1
kv N k1

Etude de stabilité par l’analyse de Fourier

La transformé de Fourier est un outil très utile pour l’étude de la stabilité de
schémas numériques dans L2 (R) : L’utilisation de la T.F. dans ce contexte est
encore appelée méthode de von Neumann. Elle se généralise notamment à des EDP
avec d variables d’espace (transformée de Fourier dans L2 Rd ), à des systèmes
d’EDP, ainsi qu’au cas de conditions aux limites périodiques (la transformée de
Fourier est alors remplacée par la décomposition en séries de Fourier).
La solution u de l’équation de la chaleur
8
@u @2u
>
< @t = @x2 (x; t) 2 [0; 1] [0; T ]
u(x; 0) = u0 (x)
>
:
u(0; t) = u(L; t) = 0 8t
sous la forme du développement est donnée par :
X
u(x; t) = u~k (t) exp (ikx)
k2Z

où k est un coe¢ cient réel intégrant le nombre d’onde, le type de conditions

aux limites et la dimension du domaine. Les coe¢ cients de Fourier u~k véri…ent
chacun une équation di¤érentielle en temps dont la solution s’écrit : u~k (t) =
u~k (0) exp ( k 2 t)
Faisons la même analyse dans le cas discret. Injectons dans le schéma numérique
une suite de solutions de la forme : unj = u~nk eikj4x
Ces solutions ont pour conditions initiales u0j = u~0k eikj4x et correspondent cha-
cune a une composante harmonique. L’étude de la stabilité se ramène à l’étude
de l’évolution au cours des pas de temps n des suites u~nk quand n augmente. La

37
Pr. A. Radid CHAPITRE 1. DIFFÉRENCES FINIES

condition minimale de stabilité numérique nécessite que les nombres u~nk restent
bornés 8k et 8n = 0; : : : ; N .
Si l’on veut de plus décroissance de la solution approchée, on devra avoir
décroissance des u~nk quand n augmente.
Dans les schémas à p pas, on obtient les u~n+1
k par multiplication par une matrice
d’ampli…cation G(4t; k) selon :
0 n+1 1 0 1
u~k u~nk
B u~n C B u~n 1 C
B k C B k C
B .. C = G (4t; k) B .. C
@ . A @ . A
u~nk p+2
u~nk p+1

Dans les schémas à un pas, la matrice d’ampli…cation se réduit à un facteur

d’ampli…cation G (4t; k) tel que u~n+1
k = Gk (4t) u~nk .

Condition nécessaire de stabilité de Von Neumann Pour que le schéma

soit stable, il faut qu’il existe > 0 tel que les valeurs propres i de la matrice
d’ampli…cation G(4t; k) soient toutes majorées en module selon :

j ij 1 + c4t 8i = 1; : : : ; p

avec c > 0; quel que soit k et pour tout 0 < 4t <

Dans le cas de schéma à un pas la matrice d’ampli…cation se réduit à un facteur
scalaire et la condition de Von Neuman est su¢ sante.

Conditions su¢ santes de stabilité 1) Si la matrice d’ampli…cation est nor-

male, c’est à dire qu’elle commute avec son adjointe (ou transposée dans le cas
réel)

GG = G G
ou bien, ce qui est équivalent, si elle admet une base de vecteurs propres or-
thonormés, alors la condition de Von Neumann est su¢ sante.
2) La condition précédente étant parfois di¤cile à véri…er, on peut utiliser la
condition su¤sante suivante : le schéma est stable si les coe¤cients de la matrice
G(4t; k) sont bornés et si ses valeurs propres sont toutes de module strictement
inférieur à 1 sauf éventuellement une de module égal à 1.
Analyse du schéma explicite :
On peut écrire le schéma sous la forme :

un+1
j = unj 1 + (1 2 ) unj + unj+1 pour j = 1; : : : ; N

38
 On considère une solution discrète particulière unj = Gn eikxj avec xj = j4x
avec le coe¢ cient d’ampli…cation.
Gn+1 eikj4x = Gn eik(j 1)4x + (1 2 ) Gn eikj4x + Gn eik(j+1)4x
On simpli…e par Gn eikj4x

G = e ik4x + (1 2 ) + e+ik4x
= 1 + e ik4x 2 + e+ik4x
ik4x ik4x 2
= 1+ e+ 2 e 2

k4x
= 1 4 sin2 2
1
donc la condition jGj 1 impose j1 4 j 1 d’où 0 2
:

Stabilité L1
Dé…nition 1.10 On dit qu’un schéma est L1 -stable si la solution approchée est
bornée dans L1 indépendamment du pas du maillage.
Proposition 1.1 Si la condition de stabilité = D 4x
4t
< 1
2
est véri…ée, alors le
schéma (1.18) est L1 stable au sens où :
sup unj = ku0 k1 :
j=1;:::;N
n=1;:::;M

Preuve 1.1 On peut écrire le schéma (1.18) sous la forme :

un+1
j = unj 1 + (1 2 ) unj + unj+1 pour j = 1; : : : ; N
Si 0 1
2
, on a 0 et 1 2 0 et la quantité un+1
j est combinaison convexe
n n n
de uj ; uj 1 et uj+1 :
Soit M n = max unj , on a alors
j=1;:::;N

un+1
j M n + (1 2 ) M n + M n pour tout j = 1; : : : ; N
et donc un+1
j M n , On en déduit en passant au maximum que : M n+1 M n:
On montre de la même manière que
max un+1
j max unj ;
j=1;:::;N j=1;:::;N

On en déduit
max un+1
j max u0j et min un+1
j min u0j
j=1;:::;N j=1;:::;N j=1;:::;N j=1;:::;N

d’où
sup un+1
j ku0 k1 :
j=1;:::;N
1 n M

39
Pr. A. Radid CHAPITRE 1. DIFFÉRENCES FINIES

Convergence
Dé…nition 1.11 Soit u la solution du problème (1.17) et unj j=1;:::;N la suite de
(1.19). On appelle erreur de discrétisation au point (xj ; tn ) la quantité enj = u unj

Théorème 1.3 (Lax) Soit u(x,t) la solution su¢ samment régulière de l’équation
de la chaleur (1.17) (avec des conditions aux limites appripriées). Soit unj la solu-
tion numérique discrète obtenue par un schéma de di¤érences …nies avec la donnée
initiale u0j = u0 (xj ) :On suppose que le schéma est linéaire, à deux niveaux, con-
sistant, et stable pour une norme kk : Alors le schéma est convergent au sens où

8T > 0; lim sup ken k =0

4t;4x !0 tn T

avec en le vecteur "erreut" dé…ni par ses composantes enj = unj u (tn ; xj )
De plus, si le schéma est précis à l’ordre p en espace et à l’ordre q en temps,
alors pour tout T > 0 il existe une constante CT > 0 telle que

sup ken k CT ((4x)p + (4t)q )

tn T

Preuve 1.2 On note unj = u (xj ; tn ). On a donc par dé…nition de l’erreur de

consistance :
un+1
j unj D
2unj unj 1 unj+1 = nj (1.20)
4t 4x2
D’autre part, le schéma numérique s’écrit :
un+1
j unj D
2
2unj unj 1 unj+1 = 0 (1.21)
4t 4x
Retranchons (1.20) et (1.21), on obtient :

en+1
j enj D
2
2enj enj 1 enj+1 = n
j (1.22)
4t 4x
Soit encore :

en+1
j = enj 1 + (1 2 ) enj + enj+1 + 4t n
i pour j = 1; : : : ; N

Or
1
enj 1 + (1 2 ) enj + enj+1 ken k1 ; car ;
2
et donc comme le schéma est consistant, nj C (4t + 4x2 ) ; ce qui entraine
que :
en+1
j ken k1 + 4tC 4t + 4x2

40
 On a donc par récurrence :

ken k1 e0 1
+ n4tC 4t + 4x2

ken k1 e0 1
+ T C 4t + 4x2
Ce qui est démontre le théorème.

1.5.2 Schéma implicite de l’équation de la chaleur

Discrétisation par Euler implicite
Discrétisons la même équation avec le même maillage mais utilisons cette fois une
approximation à gauche pour la dérivée temporelle :

@u un+1
j unj
(xj ; tn ) =
@t 4t
Observons le nouveau schéma :
8 n+1 n
> uj uj un+1 2un+1 +un+1
>
< 4t = D j+1 j
4x2
j 1
j = 1; : : : ; N
> u0j = (xj ) j = 1; : : : ; N
>
:un = un = 0
0 N 8n 2 N

la solution à l’itération n est donné par :

un+1 n+1
j 1 + (1 + 2 ) uj
n+1
uj+1 = unj pour j variant de 1à N (1.23)

On constate que les inconnues à l’itération n sont reliées entre elles par une
relation implicite (d’où le nom de la méthode)

Analyse de stabilité du schéma implicite

Stabilité matricielle : On écrit le schéma (1.23) sous forme matricielle :
0 1 0 n+1 1 0 1
1+2 0 0 u1 un1
B 1+2 0 C B un+1 C B un C
B .. CB 2 C B 2 C
B ..
.
..
.
..
. C B .
. C B .
. C
B 0 . CB . C = B . C
B ... C B n+1 C B n C
@ 0 1+2 A @ uN 1 A @ uN 1 A
0 ::: 0 1+2 un+1
N unN
A chaque itération, le vecteur des inconnues discrètes se détermine par résolu-
tion d’un système linéaire. La matrice du système étant tridiagonale, un algorithme
de Thomas (basé sur la méthode du pivot de Gauss) est très souvent utilisé.

41
Pr. A. Radid CHAPITRE 1. DIFFÉRENCES FINIES

Ce schéma nécessite la résolution d’un système d’équations. Au temps tn+1 le

système d’équations s’écrit :
0 1 0 n+1 1 0 1
u1 un1
B C B un+1 C B un C
B CB 2 C B 2 C
B . . .
.. .. .. C B . C B . C
B C B .. C = B .. C
B C B n+1 C B n C
@ A @ uN 1 A @ uN 1 A
un+1
N unN

avec = 1 + 2 D4t
4x2
= 1 + 2 et = D4t4x2
=
Soit
0 n+1 1 0 1 1 0 1
u1 un1
B un+1 C B C B un2 C
B 2 C B C B C
B .. C B ... ... ... C B .. C 1
B . C=B C B . C=C Un
B n+1 C B C B n C
@ uN 1 A @ A @ uN 1 A
un+1
N unN

avec C = I + 2 D4t
4x2
A; les valeurs propres de la matrice C sont k = 1+
4 D4t
4x2
sin2 k
2(N +1)
1; 81 k N et par suite les valeurs propres de C 1
sont
1 1
(C ) = max K
1; 81 k N; donc le schéma est inconditionnement
1 k N
stable.
Trouver les un+1
j à partir de unj nécessite ici la résolution d’un système i.e.
l’inversion d’une matrice. Le système est ici tridiagonal. On peut le résoudre
facilement par l’algorithme de Thomas, spécialement adapté aux matrices bandes.
Un tel schéma est implicite.

Analyse du schéma implicite par Von Neumann:

un+1 n+1
j 1 + (1 + 2 ) uj un+1 n
j+1 = uj pour j variant de 1à N

Gn+1 eik(j 1)4x

+ (1 + 2 ) Gn+1 eikj4x Gn+1 eik(j+1)4x = Gn eikj4x

G est le facteur d’ampli…cation

On simpli…e par Gn eikj4x on obtient :
G e ik4x + (1 + 2 ) eik4x = 1
G e ik4x + eik4x + (1 + 2 ) = 1
G( (2 cos k4x) + (1 + 2 )) = 1

42
 1
G=
1 + 2 (1 cos k4x)
ou encore
1
G= k4x
1 + 4 sin2 2
On constate que pour ce schéama la condition jGj 1 est toujours respectée, donc
qu’il est inconditionnement stable.

1.6 Approximation par di¤érences …nies de l’équation

de la chaleur en 2D
La méthode des di¤érences …nies s’étend sans di¢ culté aux problèmes en plusieurs
dimensions d’espace. Considérons par exemple l’équation de la chaleur en deux
dimensions d’espace dans le domaine rectangulaire = (0; 1) (0; L) avec des
conditions aux limites de Dirichlet
@u @2u @2u
@t @x2
= 0 pour (x; y; t) 2
@y 2
R+
u (t = 0; x; y) = u0 (x; y) pour (x; y) 2 (1.24)
u (t; x; y) = 0 pour t 2 R+ ; (x; y) 2 @

Pour discrétiser le domaine ; on introduit deux pas d’espace 4x = Nx1+1 > 0

et 4y = Ny1+1 > 0 (avec Nx et Ny deux entiers positifs). Avec le pas de temps
4t > 0, on dé…nit ainsi les noeuds d’un maillage régulier.
(tn ; xj ; yk ) = (n4t; j4x; k4y) pour n 0; 0 j Nx + 1; 0 k Ny + 1:
On note unj;k la valeur d’une solution discrète approchée au point (tn ; xj ; yk ) ; et
(t; x; y) la solution exacte de (1.24).
Les conditions aux limites de Dirichlet se traduisent, pour n > 0; en
un0;k = unNx +1;k = 0; 8k et unj;0 = unj;Ny +1 = 0; 8j
La donnée initiale est discétisée par
u0j;k = u0 (xj ; yk ) 8j; k
La généralisation au cas bidimentionnelle du schéma explicite est évidente :
un+1
j;k unj;k unj 1;k + 2unj;k unj+1;k unj;k 1 + 2unj;k unj;k+1
+ + =0
4t (4x)2 (4y)2
pour n 0; j 2 f1; : : : ; Nx g et k 2 f1; : : : ; Ny g : La seule di¤érence notable
avec le cas unidimensionnel est le caractère deux fois plus sévère de la condition
CFL.

43
Pr. A. Radid CHAPITRE 1. DIFFÉRENCES FINIES

Nous avons le résultat suivant à démontrer :

Exercice : Montrer que le schéma explicite est stable en norme L1 (et même
qu’il véri…e le principe du maximum) sous la condition CFL.
4t 4t 1
2 +
(4x) (4y)2 2
Correction

4t 4t 4t 4t
un+1
j;k = 1 2 2 + unj;k + n
2 uj 1;k + unj+1;k + n
2 uj;k 1 + unj;k+1
(4x) (4y)2 (4x) (4y)
4t 4t
Si (4x)2
+ (4y) 2
1
2
alors un+1
j;k est une combinaison convexe de coordonnées de
n
u et
un+1
j;k kun k1 :
Schéma implicite :
Le schéma implicite est :
un+1
j;k unj;k un+1 n+1
j 1;k + 2uj;k un+1
j+1;k un+1 n+1
j;k 1 + 2uj;k un+1
j;k+1
+ + = 0 (1.25)
4t (4x)2 (4y)2
Remarquons que le schéma implicite nécessaite, pour calculer un+1 en fonction
de un ; la résolution d’un système linéaire sensiblement plus compliqué qu’en une
dimension d’espace (la situation serait encore pire en 3 dimensions).
Rappelons qu’en dimension un, il su¢ t d’inverser une matrice tridiagonale.
Nous allons voir qu’en dimension deux la matrice a une structure moins simple.
L’inconnue discrète unj;k est indicée par deux entiers j et k, mais en pratique on
utilise un seule indice pour stocker un sous la forme d’un vecteur dans l’ordinateur.
Une manière (simple et e¢ cace) de ranger dans un seul vecteur les inconnues
unj;k est d’écrire

un = un1;1 ; : : : ; un1;Ny ; un2;1 ; : : : ; un2;Ny ; : : : ; unNx ;1 ; : : : ; unNx ;Ny

Remarquons qu’on a rangé les inconnues"colonne par colonne", mais qu’on aurait
aussi bien pu le faire " en "ligne par ligne"
Avec cette convention, le schéma implicite (1.25) requiert l’inverse de la matrice
symétrique tridiagonale par "blocs"
0 1
D1 E1
B E1 D2 E2 C
B C
M =B B C
C
@ ENx 2 DNx 1 ENx 1 A
ENx 1 DNx

44
 où les blocs diagonaux Dj sont des matrices carrées de taille Nx :
0 1
1 + 2 (cx + cy ) cy
B ..
. C
B cy 1 + 2 (cx + cy ) C
Dj = B . . C
@ .. .. cy A
cy 1 + 2 (cx + cy )
4t
avec cx = (4x) 2 et cx =
4t
(4x)2
; et les blocs extra-diagonaux Ej = (Ej )t sont des
matrices carrées de taille Ny
0 1
cx 0
B ... C
Ej = @ A
0 cx

Au total la matrice M est pentadiagonale et symétrique.

Cependant les cinq diagonales ne sont pas contigués, ce qui entraîne une aug-
mentation considérable du coût de la résolution d’un système linéaire associé à M .
La situation serait encore plus pire en 3 dimensions.

45
Chapitre 2

Eléments …nis

46
2.1 Introduction
2.1.1 Méthode des éléments …nis
La MEF se propose de discrétiser. La discrétisation intervient à plusieurs niveaux
:

Discrétisation : il est nécessaire de disposer d’une description du do-

maine sur lequel on souhaite travailler. Cette description va se faire en
l’approximant par un maillage, qui sera constitué d’éléments.

Interpolation : il est ensuite nécessaire de disposer d’une manière de

représenter le ou les champs inconnus. Ce que se proposer de faire la MEF,
c’est d’approximer ces champs par des fonctions plus simples (disons polynômi-
ales de degré un ou deux) dé…nies sur chacun des éléments du maillage (le
champ est approximé par des bouts de fonctions qui, elles, ne sont dé…nies
chacune que sur un seul élément).

Approximation : selon le type d’approximation, on remplace non seule-

ment l’espace V de dimension in…nie par des approximations Vh de dimension
…nie, mais parfois également les formes bilinéaires et linéaires dé…nissant le
problème.

2.1.2 Principe de base

L’idée de base de la méthode des éléments …nis est de remplacer l’espace de Hilbert
V sur lequel est posée la formulation variationnelle par un sous-espace Vh de di-
mension …nie. Le problème “approché”posé sur Vh se ramène à la simple résolution
d’un système linéaire, dont la matrice est appelée matrice de rigidité. Par ailleurs,
on peut choisir le mode de construction de Vh de manière à ce que le sous-espace Vh
soit une bonne approximation de V et que la solution uh dans Vh de la formulation
variationnelle soit “proche”de la solution exacte u dans V .

2.2 Méthode de Galerkin

2.2.1 Théorème de Lax Milgram
Théorème 2.1 Soit V un espace de Hilbert et soit a(:; :) une forme bilinéaire de
V V dans R telle que :

a(u; v) M kuk kvk ; pour tout u; v 2 V (Continuité) (2.1)

47
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

a(u; v) kvk2 ; > 0 pour tout v 2 V (V-ellipticité) (2.2)

Soit d’autre part L(:) une forme linéaire et continue de V dans R. Alors le
problème variationnel abstrait : trouver u 2 V tel que :

a(u; v) = L(v); pour tout v 2 V (2.3)

a une solution unique.

2.2.2 Choix de l’espace V et formulation variationnelle

Soit un domaine ouvert de Rn (n = 1; 2 ou 3 en pratique), de frontière = @
et sur lequel on cherche à résoudre une EDP munie de conditions aux limites.
Trouver u 2 V tel que

a (u; v) = L(v); pour tout v 2 V (2.4)

On est amené à chercher une solution approchée uh appartenant à un sous espace
Vh V de dimension …nie, soit
Trouver uh 2 Vh tel que

a (uh ; vh ) = L(vh ); pour tout vh 2 Vh (2.5)

Théorème 2.2 On suppose que l’espace V , la forme bilinéaire a(:; :) et la forme

linéaire L satisfont les hypothèses du Théorème de Lax Milgram. Alors le problème
(2.5) admet une solution unique uh 2 Vh :

Preuve 2.1 L’espace Vh étant de dimension …ni, c’est un sous espace vectoriel
fermé de V et c’est donc un Hilbert pour la norme induite par celle de V (k:k).
On peut donc appliquer le Théorème de Lax Milgram, avec V remplacé par Vh :

2.2.3 Résolution d’un système linéaire

Le problème (2.5) est équivalent à la résolution d’un système linéaire.
En e¤et, soit fwj g ; 1 j N , une base de Vh (avec dim Vh = N )
On décompose uh sur cette base.

X
N
uh = i wi (2.6)
i=1

48
et le problème (2.5) consiste à trouver le vecteur ( 1 ; : : : ; n) tel que :
X
N

i a (wi ; wj ) = L (wj ) ; 1 j N (2.7)

i=1

Ah = Ahij 1 i;j N
= (a (wi ; wj ))1 i;j N

bh = (L (w1 ) ; : : : ; L (wj ) ; : : : ; L (wN ))

La coercivité de la forme bilinéaire a(:; :) entraîne le caractère dé…ni positif de
la matrice Ah , et donc son inversibilité. En e¤et,
2
X
N
Ah uh :uh ui wi C juh j2 avec C > 0;
i=1

car toutes les normes sont équivalentes en dimension …nie (j j désigne la norme
euclidienne dans RN ). De même, la symétrie de a(u; v) implique celle de Ah .
Dans les applications mécaniques la matrice Ah est appelée matrice de rigid-
ité.
Remarque 2.1 Le résultat précédent permet de retrouver l’existence de la solution
uh du problème (2.5). En e¤et, il su¢ t alors de montrer l’unicité de la solution.
(1) (2)
Donc, soient uh et uh deux solutions du problème (2.5), on a :
(1)
a uh ; vh = (f; vh );
(2)
a uh ; vh = (f; vh ); donc
(1) (2) (1) (2) (1) (2)
a uh uh ; vh = 0; pour tout vh 2 Vh ; soit a uh uh ; uh uh = 0:
Puisque la forme bilinéaire a(:; :) est V -elliptique. On a donc
(1) (2) (1) (2)
uh uh = 0; donc uh = uh :

Remarque 2.2 La résolution du système linéaire (2.7) est d’autant plus simple et
rapide que la matrice A = fa (wi ; wj )g a une structure "agréable" (matrice bande
par exemple). Nous verrons au paragraphe prochain sur un exemples simple que
la structure de cette matrice dépend du choix du sous espace Vh et du choix de la
base fwj g N
j=1 pour ce sous espace de dimension …nie.

On a donc maintenant un théorème abstrait de majoration de l’erreur entre la

solution exacte et la solution approchée (2.7):
Théorème 2.3 (Céa) Soit u la solution du problème (2.7) et soit uh la solution
du problème (2.5). On a la majoration :
M
ku uh k inf fku vh k ; vh 2 Vh g (2.8)

49
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

Preuve 2.2 On a par coercivité de a : pour tout vh 2 Vh

ku uh k2 a (u uh ; u uh ) = a (u uh ; u vh ) + a (u uh ; vh uh )
On peut écrire
a (u; vh uh ) = L (vh uh ) ; car vh uh 2 Vh V
et on a a (uh ; vh uh ) = L (vh uh ) ;donc a (u uh ; uh vh ) = 0
Par continuité de a, on obtient:
ku uh k2 a (u uh ; u vh ) M ku uh k ku vh k
d’où on déduit
ku uh k M ku vh k ; pour tout vh 2 Vh ; ce qui entraîne l’inégalité (2.8).

C’est un bon départ pour montrer la convergence de uh vers u lorsque h tend

vers 0, mais il reste à trouver une estimation de

inf ku vh kV ;
vh 2Vh

ce qui est en général assez di¢ cile. Nous verrons ( dans l’exercice 2 de la série
2 pour le problème u" + cu = f dans ]0; 1[) comment le faire.
Lemme. On se place sous les hypothèses précédentes. On suppose qu’il existe
un sous-espace V V dense dans V et une application rh de V dans Vh (appelée
opérateur d’interpolation) tels que :

lim kv rh (v)k = 0 8v 2 V
h!0

Alors la méthode d’approximation variationnelle interne converge, c’est à dire que

lim ku uh k = 0:
h!0

Démonstration 2.1 Soit > 0: Par densité de V, il existe v 2 V tel que ku vk

:
Par ailleurs, il existe un h0 > 0 (dépendant de ) tel que, pour cet élément
v 2 V; on a
kv rh (v)k 8h h0
D’après le lemme du Céa, on a :

ku uh k C ku rh (v)k C (ku vk + kv rh (v)k)

2C ;

d’où l’on déduit le résultat.

50
2.3 Exemples en dimension 1 et introduction des
éléments …nis P1
2.3.1 Exemple en dimension 1
On considère l’équation suivante :

d2 u
= f; x 2 ]0; 1[ (2.9)
dx2
u(0) = u(1) = 0 (2.10)
Sous forme variationnelle, ce problème se formule
Trouver u 2 V = H01 (0; 1) tel que

a(u; v) = (f; v); pour tout v 2 V;

avec
Z 1 Z 1
du dv
a(u; v) = dx; (f; v) = f vdx:
0 dx dx 0

2.3.2 Éléments …nis P1

On note par Pk l’ensemble des polynômes à coe¢ cients réels d’une variable réelle
de degré inférieur ou égal à k.
La méthode des éléments …nis P1 repose sur l’espace discret des fonctions glob-
alement continues et a¢ nes sur chaque maille

Vh = fv 2 C ([0; 1]) tel que vj [xj ; xj+1 ] 2 P1 pour tout 0 j N + 1g (2.11)

et sur son sous-espace

V0h = fv 2 Vh tel que v(0) = v(1) = 0g (2.12)

Il est clair que Vh est un sous espace vectoriel : si on prend v1 et v2 2 Vh ; et

1; 2 2 R; alors 1 v1 + 2 v2 2 Vh : en e¤et, la fonction reste a¢ ne par morceaux
et les conditions aux limites sont véri…ées. Il existe donc une base pour représenter
toutes les fonctions de Vh :
On peut représenter les fonctions de Vh et V0h ; a¢ nes par morceaux à l’aide de
fonctions de base très simples. Introduisons la fonction châpeau dé…nie par :

1 jxj si jxj 1
(x) =
0 si jxj > 1

51
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

Si le maillage est uniforme : pour 0 j N + 1; on dé…nit les fonctions de base

x xj
j (x) = : (2.13)
h

Théorème 2.4 L’espace Vh dé…nie par (2.11) est un sous espace de H 1 (0; 1) de
dimension n + 2; et toute fonction vh 2 Vh est dé…nie de manière unique par ses
valeurs aux sommets (xj )0 j n+1

X
n+1
vh (x) = vh (xj ) j (x) ; 8x 2 [0; 1]
j=0

De même,l’espace V0h ; dé…ni par (2.12), est un sous-espace de H01 (0; 1) de dimen-
sion n; et toute fonction vh 2 V0h est dé…nie de manière unique par ses valeurs
aux sommets (xj )1 j n

X
n
vh (x) = vh (xj ) j (x) ; 8x 2 [0; 1]
j=1

Démonstration 2.2 Rappelons que les fonctions continues et de classe C 1 par

morceaux appartiennent à H 1 (0; 1) : Donc Vh et V0h sont bien des sous espaces

52
de H 1 (0; 1) ; et V0h est dans H01 (0; 1) : Par ailleurs on véri…e que les ( i )0 i N +1
forment une base de Vh : Remarquant que : j (xi ) = ij: ; où ij: est le symbole de
Kronecker qui vaut 1 si i = j et 0 sinon..
Soient n+2 scalaires 0 ; 1 ; : : : ; n+1 tels que la fonction

X
n+1
vh (x) = j j (x) ; 8x 2 [0; 1] ;
j=0

en particulier

X
n+1
8i 2 f0; 1; : : : ; N + 1g ; 0 = vh (xi ) = j j (xi ) = i;
j=0

donc la famille est libre.

C’est une famille génératrice car soit vh 2 Vh et

X
n+1
w (x) = vh (xj ) j (x) ; 8x 2 [0; 1] :
j=0

Evidemment w 2 Vh et w (xi ) = vh (xi ) ; 8i donc w = vh et par suite c’est une

base de Vh :
Soit vh 2 V0h alors vh 2 Vh et

X
n+1 X
n
vh (x) = vh (xj ) j (x) = vh (xj ) j (x)
j=0 j=1

car vh (0) = vh (1) = 0 donc la famille ( i )1 i N est génératrice et elle est libre
donc c’est base de V0h ;

Remarque 2.3 La base j dé…nie par (2.13), permet de caractériser une fonc-
tion de Vh par ses valeurs aux noeuds du maillage. Dans ce cas on parle d’éléments
…nis de Lagrange. Par ailleurs, comme les fonctions sont localement P1 , on dit que
l’espace Vh dé…nie par (2.11) est l’espace des éléments …nis de Lagrange d’ordre 1.

2.3.3 Écriture matricielle

La fomulation variationnelle de la solution approrchée devient :
Trouver uh 2 V0h tel que :
Z 1 Z 1
0 0
uh (x) vh (x) dx = f (x)vh (x) dx 8vh 2 V0h
0 0

53
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

On décompose uh sur la base j 1 j n et on prend vh = i ce qui donne :

X
n Z 1 Z 1
0 0
uh (xj ) j (x) i (x) dx = f (x) i (x) dx:
j=1 0 0

R1
En notant uh = (uh (xj ))1 j n ; bh = 0
f (x) i (x) dx et introduisant la
1 i n
matrice de rigidité
Z 1
0 0
Ah = j (x) i (x) dx
0 1 i;j n

La formulation variationnelle dans V0h revient à résoudre dans Rn le système

linéaire
Ah u h = b h (2.14)
Comme les fonctions de base j ont un petit support, l’intersection des supports
de i et j est souvent vide et la plupart des coe¢ cients de Ah sont nuls. Un calcul
simple montre que :
8 1
Z 1 >
> si j=i 1
< 2h
0 0 h
si j=i
j (x) i (x) dx = 1
0 >
> si j = i+1
: h
0 si j 6= fi 1; i; i + 1g

On a donc à résoudre le système linéaire

Z
1 1 xj+1
(uj 1 2uj + uj+1 ) = f 'j dx; 1 j I 1 (2.15)
h2 h xj 1

u0 = uI = 0; (2.16)
IP1
car uh (x) = uj 'j (x) ; avec uj = uh (xj ) :
j=1
0 1
2 1 0
1B C
...
B 1 2 C
Ah = B .. .. C (2.17)
h@ . . 1 A
0 1 2
et la matrice Ah est tridiagonale et on obtient toujours cette structure tridi-
agonale si on remplace le problème (2.9), (2.10), par le problème (2.18) . Ceci
provient du fait que les fonctions de base ont un support local.

54
Remarque 2.4 Si on considère le problème de Dirichlet homogène suivant :
d2 u
dx2
+ u(x) = f (x) 8x 2 ]0; 1[ =
(P ) (2.18)
u (0) = u (1) = 0

on a au plus de la matrice de rigidité (2.17), la matrice de masse qui s’écrit de la

forme suivante : 2 3
4 1 0 ::: 0
6 1 4 1 0 0 7
6 7
6 .. 7
h6 0 1 4 1 . 7
M= 6 6 . . . . 7
66 .. .. .. .. 0 7
7
6 . 7
4 .. 0 1 4 1 5
0 0 ::: 0 1 4

Remarque 2.5 Si on modi…e la numérotation des fonctions de base 'j (x) ; la

matrice du système linéaire (2.15) perd sa structure diagonale.

2.3.4 Calcul du second membre

Pour obtenir le second membre bh il faut calculer les intégrales
Z xj+1
(bh )j = f (x)'j (x) dx; pour tout 1 j N
xj 1

L’évaluation exacte du second membre bh peut être di¢ cile ou impossible si

la fonction f est compliquée. En pratique, on utilise des formules d’intégrations
numériques qui donnent une approximation des intégrales dé…nissant bh : Par ex-
emple, on peut utiliser :
La formule des rectangles à gauche :
Z xi+1
1
(x) dx (xi )
xi+1 xi xi

La formule des rectangles à droite :

Z xi+1
1
(x) dx (xi+1 )
xi+1 xi xi

La formule du point milieu :

Z xi+1
1 xi+1 + xi
(x) dx
xi+1 xi xi 2

55
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

La formule des trapèzes :

Z xi+1
1 1
(x) dx ( (xi+1 ) + (xi ))
xi+1 xi xi 2

La formule de Simpson :
Z xi+1
1 1 1 2 xi+1 + xi
(x) dx (xi+1 ) + (xi ) +
xi+1 xi xi 6 6 3 2

Les deux premières formules sont exactes pour les fonctions a¢ nes et la
troisième est exacte pour les polynômes de second degré.
Pour les fonctions règulières, ces formules sont approchées avec un reste d’ordre
O(h); O(h2 ); O(h3 ) respectivement.
Par exempe, La formule des trapèzes :
Rx
(bh )j = xjj+11 f (x) j (x) dx
Rx Rx
= xjj 1 f (x) j (x) dx + xjj+1 f (x) j (x) dx
x x
= j 2 j 1 f (xj ) j (xj ) + f (xj 1 ) j (xj 1 )
x x
+ j+12 j (f (xj+1 ) j (xj+1 ) + f (xj ) j (xj ))
= h2 f (xj ) + h2 f (xj )
= hf (xj ) :

Ainsi,
0 le système à résoudre
1 0d’écrit
1: 0 1
2 1 0 u1 f (x1 )
B ... C B u2 C B f (x2 ) C
1 B 1 2 CB C B C
h2 B . . C B .. C = B .. C
@ .. .. 1 A@ . A @ . A
0 1 2 un f (xn )
Avec cette approximation, on retrouve exactement le même système linéaire
que celui obtenu par la méthode des di¤érences …nies.
La matrice A est tridiagonale, symétrique et dé…nie positive, alors le système
ci-dessus admet une unique solution.

Remarque 2.6 Considérons l’expression Rx

Ej (u) = h12 (u (xj 1 ) 2u (xj ) + u (xj+1 )) h1 xjj+11 f 'j dx; où u est la solu-
tion exacte du problème (2.9), (2.10).
Si on considère les développements de Taylor de u (xj 1 ) et u (xj+1 ) autour de
xj ; on obtient si la solution appartient à C 4 (0; 1) :
0
1 d2 u h2 d4 u d4 u
h2
(u (xj 1 ) 2u (xj ) + u (xj+1 )) = dx2
(xj ) 24 dx4
(xi + h) + dx4
xi h ;
0
0< < 1; 0 < < 1;

56
 D’autre part, si on développe la fonction f (x) autour du point xj :
0 2 d2 f
f (x) = f (xj ) + (x xj ) f (xj ) + (x xj ) ( )
dx2
où est un point appartenant à ]xj ; x[
8 x x
> j 1
< xj xj 1 si x 2 [xj 1 ; xj ]
x xj+1
'j (x) = si x 2 [xj ; xj+1 ]
> x x
: j 0 j+1 si x > x
j+1 et x < xj 1
R xj+1
xj 1R
(x xj ) 'j (x) f 0 (xj )dx = 0; En e¤et
xj+1 0
R xj R xj+1
xj 1
(x x j ) f (x j )' j (x) dx = xj 1
(x xj ) f 0 (xj )'j (x) dx+ xj
(x xj ) f 0 (xj )'j (x) dx
R xj R xj+1
= f 0 (xj ) xj 1
(x
x j ) ' j (x) dx + xj
(x xj ) 'j (x) dx
h 2
i xj R
(x xj ) xj (x xj )2
= f 0 (xj ) 2
' j (x) xj 1 2h
dx
xj 1
h ixj+1 R
(x xj )2 x (x x )2
+ 2
' j (x) + xjj+1 2hj dx
x
h ij
3 xj
h ixj+1
(x xj ) (x xj )3
= f 0 (xj ) 6h
+ 6h
xj 1 xj
0 h2 h2
= f (xj ) 6
+ 6
=0
onR obtient alors: R xj+1 2f ( )
1 xj+1
h xj 1
1
f (x) 'j (x) dx = f (xj )+ 2h xj 1
(x xj )2 'j (x) d dx2 dx; Si f
00
2 C 0 (]0; 1[) ; on
a: R n 2 o
1 xj+1 h2 d f
h xj 1
f 'j dx f (xj ) 12
max dx2 ; x 2 ]0; 1[ :
d2 u 4 2f
On a donc, puisque
n dx2
= f: donco ddxu4 = dxd
2
2 d4 u
jEj (u)j h6 max dx4
;x 2 ]0; 1[ , or on peut écrire
1
Rx d2 u
h2
(u (xj 1 ) 2u (xj ) + u (xj+1 )) h1 xjj+11 f 'j dx dx2
(xj ) f (xj ) =
Ej (u) :
La relation (2.15) est donc une approximation de l’équation à un terme d’ordre
h (O (h2 )) près. La quantité Ej (u) ( en terminologie de di¤érences …nies) s’appelle
2

l’erreur de troncature au point xj :

2.3.5 Majoration de l’erreur dans le cas de l’exemple

Dé…nition 2.1 On appelle opérateur d’interpolation P1 l’application linéaire rh de
H01 (0; 1) dans Vh dé…nie pour tout v 2 H 1 (0; 1) par
X
n+1
(rh v) (x) = v (xj ) j (x)
j=0

57
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

Interpolation P1 d’une fonction de H 1 (0; 1)

Remarque 2.7 Cette dé…nition a bien un sens car les fonctions de H 1 (0; 1) sont
continues et leurs valeurs ponctuelles sont donc bien dé…nies.
L’interpolée rh v d’une fonction v est simplement la fonction a¢ ne par morceaux
qui coïncide avec v sur les sommets du maillage xj :

Soit u 2 H01 (0; 1) = V la solution du problème :

Trouver u 2 V tel que a(u; v) = (f; v), avec
Z 1 Z 1
du dv
a(u; v) = dx; (f; v) = f vdx :
0 dx dx 0

Soit uh 2 Vh solution de

a(uh ; v) = (f; vh ); pour tout vh 2 Vh

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

M
ku uh k inf fku vh k ; vh 2 Vh g ; (2.19)

R1 dv 2
1
2
où kvk = 0 dx
dx :

Remarque 2.8 Soit v une fonction de C 0 [0; 1]. On rappelle que rh la fonction Vh -
interpolée de v; c’est à dire l’unique fonction de Vh telle que rh v (xi ) = v(xi ); 0
i N + 1:

58
 L’inégalité (2.19) permet d’écrire, puisque u 2 H01 (0; 1) C 0 [0; 1]

M
ku uh k ku rh uk (2.20)

Il su¢ t alors de majorer l’erreur entre fonction et son interpollation.

Lemme 2.1 (d’interpolation) Soit rh l’opérateur d’interpolation P1 . Pour tout

v 2 H 2 (0; 1) ; alors
00
jv rh vj1;(0;1) Ch v (2.21)
L2 (0;1)
00
kv rh vkL2 (0;1) Ch2 v (2.22)
L2 (0;1)

R1 dv 2
1
2
avec jvj1;(0;1) = 0 dx
dx

Preuve 2.3 rh v(x)j (xi ; xi+1 ) = xxi xxi+1

i+1
x xi
v(xi ) + xi+1 xi
v(xi+1 ) = x xi+1
h
v(xi ) +
x xi
h
v(xi+1 )
N R
P xi+1 2
jv rh vj21;(0;1) = xi
d
dx
(v rh v) dx; et
i=0

d dv v(xi+1 ) v(xi )
(v rh v) j (xi ; xi+1 ) = (2.23)
dx dx h
Utilisant une formule de Taylor avec reste d’intégrale
Z b
dv d2 v (t)
v (b) = v (a) + (b a) (a) + (b t) dt;
dx a dt2

il vient, en développant v (xi+1 ) et v (xi ) autour du point x :

Z xi+1
dv d2 v (t)
v (xi+1 ) = v (x) + (xi+1 x) (x) + (xi+1 t) dt;
dx x dt2
Z xi
dv d2 v (t)
v (xi ) = v (x) + (xi x) (x) + (xi t) dt; (2.24)
dx x dt2
d’où d’après(2.23)
Z xi+1 Z x
d 1 d2 v 1 d2 v
(v rh v) j(xi ;xi+1 ) = (xi+1 t) 2 dt + (xi t) dt (2.25)
dx h x dt h xi dt2

On a en appliquant (a + b)2 2a2 + 2b2

59
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

d 2 2
R xi+1 2
2 Rx 2
2
dx
(v rh v) j(xi ;xi+1 ) h2 x
(xi+1 t) ddt2v dt + h22 xi (xi t) ddt2v dt
R xi+1 R xi+1 d2 v 2
2
h2 x
(xi+1 t)2 dt x dt2
dt +
Rx Rx d2 v
2
2
h2 xi
(xi t)2 dt xi dt2
dt
2 (xi+1 x)3 (x xi )3 R xi+1 d2 v
2
h2 3
+ 3 xi dt2
dt

car x 2 ]xi ; xi+1 [ : On a donc

Z Z ! Z 2
xi+1
d 2 2 xi+1
(xi+1 x)3 (x xi )3 xi+1
d2 v
(v rh v) dx + dx dt
xi dx h2 xi 3 3 xi dt2
" #xi+1 Z
2
2 (xi+1 x)4 (x xi )4 xi+1
d2 v
+ dt
h2 12 12 xi dt2
xi
4 4 Z xi+1 2 2
2 h h dv
+ dt
h2 12 12 xi dt2
Z xi+1 2
h2 d2 v
dt
3 xi dt2

donc
N Z
X xi+1
2 d
kv 0
rh v 0 kL2 (0;1) = (v rh v)2 dx
i=0 xi dx
Z
h2 X
N xi+1 2
d2 v h2 00 2
dt = kv kL2 (0;1)
3 i=0 xi dt2 3

D’où (2.21).
On a
N Z
X xi+1
kv rh vk2L2 (0;1) = (v rh v)2 dx;
i=0 xi

On a :
x xi
rh vj(xi ;xi+1 ) = v(xi ) + (v (xi+1 ) v(xi )) (2.26)
h
et d’après (2.24) on a :
Z xi
dv d2 v (t)
v(x) = v (xi ) (xi x) (x) (xi t) dt (2.27)
dx x dt2

60
 alors (2.27)-(2.26) donne
Z xi
dv d2 v v (xi+1 ) v(xi )
v rh vj(xi ;xi+1 ) = + (x xi ) (x) (xi t) 2 dt (x xi )
dx x dt h
Z xi 2
dv v (xi+1 ) v(xi ) dv
= (xi t) 2 dt (x xi ) (x)
x dt h dx
Z xi
d2 v dv drh v
= (xi t) 2 dt + (x xi ) (x) (x)
x dt dx dx
Z xi Z
d2 v x xi xi+1 d2 v
= (xi t) 2 dt (xi+1 t) 2 dt
x dt h x dt
Z xi 2
x xi dv
+ (xi t) 2 dt (d’après(2.25))
h x dt
Z xi Z
x xi+1 d2 v x xi xi+1 d2 v
= (xi t) 2 dt (xi+1 t) 2 dt
h x dt h x dt

En élevant au carré, il vient

2 x xi+1 2 R xi 2
2 2 R xi+1 2
2
v rh vj(xi ;xi+1 ) 2 h x
(xi t) ddt2v dt +2 x hxi x
(xi+1 t) ddt2v dt
x xi+1 2 (xi x)3 2 3 R xi+1 d2 v 2
2 h 3
+ x hxi (xi+13 x) xi dt2
dt:
Rx 2
2
2
3h
(x xi+1 )2 (x xi )2 xii+1 ddt2v dt:
On intégre sur l’intervalle ]xi ; xi+1 [ :
Z xi+1 Z xi+1 Z xi+1 2
2 2 2 2 d2 v
(v rh v) dx (x xi+1 ) (x xi ) dx dt:
xi xi 3h xi dt2

On applique deux fois l’intégration par parties :

Z " #xi+1 Z
xi+1
2 2 (x xi+1 )3 2
xi+1
(x xi+1 )3
(x xi+1 ) (x xi ) dx = (x xi ) 2 (x xi ) dx
xi 3 xi 3
" #xxi i+1 Z
4
(x
xi+1 ) xi+1
(x xi+1 )4
= (x xi ) + dx
6 xi 6
xi
" #xi+1
(x xi+1 )5 h5
= =
30 30
xi

Z xi+1 Z xi+1 2
2 h4 d2 v
(v rh v) dx dt:
xi 45 xi dt2

61
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

On fait la somme de ces inégalités, on obtient :

N Z
X Z
h4 X
xi+1 N xi+1 2
d2 v h4 00 2
kv rh vk2L2 (0;1) = (v 2
rh v) dx dt = v
i=0 xi 45 i=0 xi dt2 45 L2 (0;1)

On en déduit l’inégalité (2.22).

Ce résultat montre qu’en ra¢ nant le maillage (c’est-à-dire en prenant un pas

de discrétisation h de plus en plus petit), nous pouvons approcher toute fonction
de H 2 ( ) \ H01 ( ) par une fonction de Vh avec une plus grande précision et que
l’erreur d’interpolation v rh v tend vers 0 lorsque h ! 0: En…n, cette erreur tend
vers zéro à l’ordre 1 en h pour la norme H 1 et à l’ordre 2 en h pour la norme L2 .
De l’inégalité (2.20) et du Lemme, on déduit le

Théorème 2.5 Soit u 2 H01 (0; 1) et uh 2 V0h les solutions de (2.9) et (2.14),
respectivement. Alors la méthode des éléments …nis P1 converge, c’est à dire que :
lim ku uh kH 1 (0;1) = 0:
h!0
De plus, si u 2 H 2 (0; 1) (ce qui est vrai si f 2 L2 (0; 1)); alors il existe une
constante C indépendante de h telle que :
00
ku uh kH 1 (0;1) Ch u = Ch kf kL2 (0;1)
L2 (0;1)

2.3.6 Fonctions de forme locales

A…n de généraliser ce qui a été fait à des éléments …nis de Lagrange utilisant des
polynômes de degré élevé, il est utile d’introduire la notion de fonction de forme
locale. Considérons la maille de référence K = [0; 1] sur laquelle nous dé…nissons
les fonctions de forme locales fP0 ; P1 g par

P0 = 1 t; P1 = t:

En introduisant les noeuds de K dé…nis par a0 = 0 et a1 = 1; il vient

Pi (aj ) = ij ; 80 i; j 1:

Par ailleurs, toute maille Ki = [xi ; xi+1 ] ; 0 i n; est l’image de référence

K pour la transformation a¢ ne

Ti : K 3 t 7 ! x = xi + thi 2 Ki (2.28)

avec hi = xi+1 xi ; de transformation inverse

1 x xi
Ti : Ki 3 x 7 ! t = 2K (2.29)
hi

62
On véri…e facilement que :
8
< P0 Ti 1 (x) si x 2 Ki
81 i n; 'i (x) = P1 Ti 11 (x) si x 2 Ki 1
:
0 sinon

si bien qu’en utilisant la régle de la dérivation composée, il vient

8 1 0 1
< hi P0 Ti (x) si x 2 Ki
0 1 0 1
81 i n; 'i (x) = P T (x) si x 2 Ki 1
: hi 1 1 i 1
0 sinon

Introduisons la matrice de rigidité élémentaire d’ordre 2 donnée par :

Z R1 R1 0 !
(P 0 (t))2 dt P (t) P 0
(t) dt
Ae = Pm0 (t) Pn0 (t) dt = R 10 0 0 0 0R 0
1 0
1
2
K 0 m;n 1 0
P1 (t) P0 (t) dt 0
(P 1 (t)) dt
R1 R1 !
( 1)2 dt 1 1dt 1 1
= R 10 0R
1 = :
0
1 ( 1) dt 0
(1)2 dt 1 1

Les termes non-nuls dans la matrice de rigidité Ah peuvent s’évaluer en utilisant

cette matrice. En e¤et, il vient pour tout 1 i n;
R1
(Ah )ii = 0 ('0i (x))2 dx
Pn R
= Kj
('0i (x))2 dx
Rj=0 R
= Ki 1 ('0i (x))2 dx + Ki ('0i (x))2 dx
R 2 R 2
= Ki 1 h21 P10 Ti 11 (x) dx + Ki h12 P00 Ti 1 (x) dx
R i 1 R i
= hi1 1 K P10 (t)2 dt + h1i K P00 (t)2 dt
= hi1 1 Ae22 + h1i Ae11

De même, pour tout 81 i n;

R1
(Ah )ii+1 = R0 '0i (x) '0i+1 (x) dx
= Ki '0i (x) '0i+1 (x) dx
R1
= h1i 0 P0 (t) P1 (t) dt
= h1i Aeii+1
= hi1

Si le maillage est uniforme c’est à dire : 81 i n; hi = h; On obtient :

2 1
(Ah )ii = h ; (Ah )ii+1 = h ; On retrouve bien la matrice (Ah ) donnée par (2.17).

63
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

Exercice 2.1 On considère le problème d’advection-di¤usion suivant :

d2 u
dx2
+ u0 (x) = f (x) 8x 2 ]0; 1[ =
(P )
u (0) = u (1) = 0

où réel positif et f une fonction dans L2 ( ) :

Résoudre cette équation par la méthode des éléments …nis P1

Corrigé de l’exercice 2.1 La formulation variationnelle de ce problème s’écrit

: Trouver la fonction u 2 H01 ([0; 1]) telle que pour tout v 2: H01 ([0; 1]) on ait :
Z 1 Z 1 Z 1
0 0 0
u (x) v (x) dx + u (x) v (x) dx = f (x)v (x) dx
0 0 0

La formulation variationnelle, issue de l’utilisation des éléments …nis P1 ; con-

siste à trouver uh 2 V0h = fv 2 Vh tel que v(0) = v(1) = 0g avec Vh = fv 2 C ([0; 1]) tel que vj [x
1
xj = jh où h = n+1 tel que :
Z 1 Z 1 Z 1
u0h (x) vh0 (x) dx + u0h (x) vh (x) dx = f (x)vh (x) dx 8vh 2 V0h
0 0 0

et
uh (0) = 0; uh (1) = 0:
On note ( i )i=1;:::;n la base de Vh0 dé…nie par i (xj ) = ij : On décompose uh sur

la base j 1 j n et on prend vh = i comme fonction test, ce qui donne :

X
n Z 1 Z 1 Z 1
0 0 0
uh (xj ) j (x) i (x) dx+ j (x) i (x) dx = f (x) i (x) dx: pour tout 1 i n
j=1 0 0 0

Le système correspondant à la discrétisation s’écrit :

KU + CU = F

où la matrice C a pour coe¢ cients

Z 1
0
Cij = i (x) j (x) dx
0
R xi x xi 1 1 R xi+1 xi+1 x 1
Cii = dx + dx
h xi 1 h i h
xi h xi hi
xi+1
h
(x xi 1 )2 (xi+1 x)2
= 2h 2 + 2h 2
xi 1 xi
h2 h2
= 2h2 2h2
= 0

64
 R1 0
Cii+1 = R0 x' i (x) 'i+1 (x) dx
i+1 xi+1 x 1
= dx
h xi h i h
2 xi+1
(xi+1 x)
= 2h2
= 12 :
xi
R1
Ci+1i = R0 'i+1 (x) '0i (x) dx
xi+1 x xi 1
= dx
h xi hi
2 xi+1
h
(x xi )
= 2h2
= 12 :
xi
0 1 0 1 1
2 1 0 2
B .. ..
. . C B 1 ... ... C
1 B 1 C B C
K= h B ... C et C = B 2 . C
@ 1 A @ .. 1 A
2
1
1 2 2
0
ci11 ci12
Ce =
ci21 ci22
R xi+1 R1 R1 h i1
1 x xi x xi (1 t)2
avec ci12 = h xi
P0 h
P10 h
dx = 0
P0 (t) 1dt = 0
(1 t)dt = 2
=
0
1
2
:

1
R xi+1 x xi x xi
R1 R1 t2 1
ci21 = h xi
P1 h
P00 h
dx = 0
tdt = 0 2
dt = 2
:
R xi+1
ci11 = 1
hi xi
P0 x hixi P00 x xi
hi
dx
R1
= P (t)P00 (t)dt
R01 0
= (1 t)dt
h0 2 i1
(t 1)
= 2
0
1
= 2

R xi+1
ci22 = 1
P1 x hixi P10 x xi
dx
Rhi1 xi hi
= P (t)P10 (t)dt
R01 1
= tdt
h0 2 i 1
t
= 2
0
1
= 2

Exercice 2.2 On considère le problème de Neumann non-homogène suivant :

d2 u
u(x) dx2
= f (x) 8x 2 ]a; b[ =
(P )
u0 (a) = et u0 (b) =

65
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

où et sont deux réels positifs et f une fonction dans L2 ( ) :

Trouver la matrice du système linéaire issu de la méthode des éléments …nis
P1 :

Corrigé de l’exercice 2.2 Le problème est mis sous la forme variationnelle faible
:

Z b Z b Z b
0 0 b
u (x) v (x) dx + u (x) v (x) dx = f (x)v (x) dx + [v(x)u0 (x)]a
a a a

Trouver u 2 H 1 (0; 1) telle que :

Z b Z b Z b
0 0
u (x) v (x) dx + u (x) v (x) dx = f (x)v (x) dx + (v(b)u0 (b) v(a)u0 (a)) ;
a a a
8v 2 H 1 (0; 1)

le domaine = [a; b] est découpé en n sous domaines élémentaires Kj ; cor-

respondant à la subdivision a = x0 ; x1 ; : : : ; xn+1 = b avec xi = a + ih pour
i = 0; : : : ; n + 1

Le problème approché est alors

Trouver uh 2 Vh telle que :

Rb Rb Rb
8vh 2 Vh ; a h
u (x) vh (x) dx + a u0h (x) vh0 (x) dx = a f (x)vh (x) dx + (vh (b) vh (a) )
(2.30)
Trouver uh (x0 ) ; ; uh (xn+1 ) 2 R tels que :

P
n+1
uh (x) = ui 'i (x) avec ui = uh (xi ) 8i; 0 i n + 1:
i=0

On prend vh = 'j ; alors (2.30) devient :

trouver le vecteur U = (u1 ; : : : ; un+1 ) tel que
P Rb
n+1 P Rb 0
n+1
ui a 'i (x) 'j (x) dx + ui a 'i (x) '0j (x) dx =
R bi=0 0
i=0

a
f (x)' j (x) dx + ' j (b)u (b) 'j (a)u0 (a) ; 8j; 0 j n+1

On pose S = ( u0 (a) ; : : : ; u0 (b)) = ( ; : : : ; ) et en dé…nissant la matrice de

masse par :
Z b
M = (Mij )1 i;j n+1 = 'i (x) 'j (x) dx
a 1 i;j n+1

66
et la matrice de rigidité
Z b
K = (Kij )1 i;j n+1 = '0i (x) '0j (x) dx
a 1 i;j n+1
Rb
et le vecteur de charge F par : (F1 ; : : : ; Fn+1 )t avec Fj = a
f (x)'j (x) dx
l’équation s’écrit sous la forme matricielle :

( M + K) U = (F + S)

Si on a CL de Neumann homogène c’est à dire :u0 (a) = u0 (b) = 0 alors S =

(0; : : : ; 0) et le système devient :

( M + K) U = F
8 1
Z >
> h
si j=i 1
1 < 2
0 0 h
si j=i
j (x) i (x) dx = 1
0 >
> h
si j =i+1
:
0 si j 6= fi 1; i; i + 1g
R1 0 2 1
( 0 (x)) dx = h
R01 0 2
n+1 (x) dx = h1 h ix
R0b R xi x xi 1 xi x (x xi 1 )2 xi x i Rx 2
'
a i
(x) ' i 1 (x) dx = xi 1 h h
dx = 2h h
+ xii 1 (x x2hi 1 ) h1 dx
h i xi xi 1
(x xi 1 )3 h
= 6h2
= 6 pour tout 1 i n
xi 1 h i xi h ixi+1
Rb 2 R xi (x xi 1 )2 R xi +1 (xi+1 x)2 (x xi 1 )3 (xi+1 x)3
'
a i
(x) dx = xi 1 h2
+ xi h2
dx = 3h2
+ 3h2
=
xi 1 xi
h
3
+ h3 = 2 h3 ;
pour tout 1 i n h i x1
Rb Rx 2
(x1 x)3
a
(x) dx = a 1 (x1h2x) dx = h12
'20 3
= h3 :
Rb 2 R 1 (x xn )2 h i1
a
1 (x xn )3
'
a n+1
(x) dx = xn h2 dx = h2 3
= h3 :
xn
Si on prend par exemple f (x) = 1
Z 1 Z 1
F =( '1 dx; :::; 'n dx)
0 0
R1 R xi R xi+1 h i xi h ixi+1
x xi 1 xi+1 x (x xi 1 )2 (xi+1 x)2
0
'i dx = xi 1 h
dx + xi h
dx = 2h
+ 2h
=
xi 1 xi
h
2
+ h2 = h: h i x1
R1 R x1 x1 x (x1 x)2
0
'0 dx = a h
dx = 2h
= h2 :
R1 R 1 x xn h i1
a
(x xn )2
0
' n+1 dx = xn h
dx = 2h
= h2 :
xn

67
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

Sur la maille Ki = [xi ; xi+1 ] sur cet élément, il n y a que 2 fonctions de base
non nulles : 'i et 'i+1
8
< P0 Ti 1 (x) si x 2 Ki
81 i n; 'i (x) = P1 Ti 11 (x) si x 2 Ki 1
:
0 sinon
avec Ti 1 (x) = x hixi et hi = xi+1 xi P0 (t) = 1 t et P1 (t) = t
l’élément Ki = [xi ; xi+1 ] produira donc e¤ectivement une contribution non nulle
aux 4 coe¢ cients Mii ; Mii+1 ; Mi+1i et Mi+1i+1 de la matrice globale M .
Calculons les contributions élémentaires de Ki et disposons les sous la forme
d’une matrice élémentaire 2 2

ei11 ei12
Me =
ei21 ei22
R xi+1 R1 R1 h i1
x xi (1 t)3
avec ei11 = xi
P02 hi
dx = hi 0
P02 (t)dt = hi 0
(1 t)2 dt = hi 3
=
0
hi
3
:
R xi+1 x xi
R1 R1 hi
ei22 = xi
P12 hi
dx = hi 0
P12 (t)dt = hi 0
t2 dt = 3
:

R xi+1 x xi x xi
ei21 = ei12 = xi
P0 hi
P1 hi
dx
R1
= hi 0 P0 (t)P1 (t)dt
R1
= hi 0 t(1 t)dt
h2 i1
t3
= hi t2 3
0
hi
= 6

Donc

hi 2 1
Me = :
6 1 2
Dé…nissons une expansion de la matrice élémentaire, en remplaçant dans une
matrice de zéros, les coe¢ cients des ieme et i + 1eme lignes et colonnes par les
composantes 0 de la matrice élémentaire
1 M e:
0 ::: 0
B C
B C
B .. .. C
B . C
M (i) = h6i B . 2 1 C
B 1 2 C
B C
@ A
0 ::: 0

68
 La matrice de masse globale s’obtient en assemblant les matrices élémentaires.
Si hi = 0
h 8i = 1; : : : ; n;
1 on obtient
0 alors 1 0 1
2 1
B 1 2 C B 2 1 C B C
B C B C B C
M = h6 BB
C+ hB 1 2
C 6B
C+
C + h B
6 B
C
C
@ A @ A @ 2 1 A
0 1 1 2
2 1
B 1 4 1 C
B C
h B . . C
= 6B . C
B C
@ 1 4 1 A
1 2
De la même manière, le calcul de la matrice de rigidité s’e¤ectue après expan-
sion de la matrice de rigidité élémentaire,
i i 1
k11 k12 1 1
Ke = i i = :
k21 k22 hi 1 1
en posant 0 1
0 ::: 0
B C
B C
1 B .
B .
.. C
. C
K (i) = B . 1 1 C
hi B 1 1 C
B C
@ A
0 ::: 0
La matrice de masse rigidité s’obtient en assemblant les matrices élémentaires.
Si hi = h 8i = 1; : : : ; n; on obtient alors
0 1
1 1 ::: 0
B 1 2 1 C
B C
B
1 B ... .. .. .. C
K= B . . . CC
hB C
B C
@ 1 2 1 A
0 ::: 1 2

2.3.7 Technique d’assemblage

Considérons un maillage à N éléments et notons B la matrice globale à assembler
(matrice de raideur ou de masse globale), bk les vecteurs élémentaires correspon-
dants relatifs à chaque élément Kk et F le vecteur global du second membre, lk les
vecteurs élémentaires correspondants relatifs à chaque élément Kk :

69
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

L’algorithme d’assemblage est très simple dès que l’on dispose d’un tableau
associant les points d’un élément Kk et les noeuds du maillage global.
Son schéma est le suivant :
Pour k=1:N faire %boucle sur les éléments
pour i=1:2 faire %boucle sur les numéros locaux
pour j=1:2 faire %boucle sur les numéros locaux
I=k+i-1 %numéros globaux
J=k+j-1 %numéros globaux
B (I; J) = B (I; J) + b (i; j) %matrcice globale, b:matrice
élémentaire
FIN des 3 boucles

Pour second membre :

Pour k=1:N faire %boucle sur les éléments
pour i=1:2 faire %boucle sur les numéros locaux
I=k+i-1 %numéros globaux
F(I)=F(I)+l(i) %F: vecteur global, l: vecteur élémentaire
Fin des 2 boucles

2.4 Eléments …nis P2

Dans certaines applications, on peut considérer que l’approximation par des droites
sur chaque élément du maillage [xj ; xj+1 ] est trop grossière, c’est à dire qu’il fournit
une fonction approchée trop éloignée de la fonction exacte. On peut alors chercher
à approcher u sur chaque maille par des polynômes de plus haut degré.
La méthode des éléments …nis P2 repose sur l’espace discret

Vh = fv 2 C ([0; 1]) tel que vj [xj ; xj+1 ] 2 P2 pour tout 0 j n + 1g (2.31)

et sur son sous-espace

V0h = fv 2 Vh tel que v(0) = v(1) = 0g (2.32)

Ces deux espaces sont composés de fonctions continues, paraboliques par morceaux
qu’on peut représenter à l’aide de fonctions de base très simples.
Introduisons tout d’abord les points milieux des segments [xj ; xj+1 ] dé…nis par
x +x
xj+ 1 = xj + h2 = j 2 j+1 pour 0 j n: On dé…nit aussi deux fonctions "mères"
2
: 8
< (1 + x) (1 + 2x) si 1 x 0
(x) = (1 x) (1 2x) si 0 x 1 (2.33)
:
0 si jxj > 1

70
 Figure 2.1: Les fonctions de base des éléments …nis P2

et
1
1 4x2 si jxj 2
(x) = 1 (2.34)
0 si jxj > 2

8j; 1 j n; j(xi ) = ij et
8j; 0 j n; j+ 1 (xi ) = ij
2

Les fonctions de base correspondant à un point xj extrémité d’un élément sont

dé…nies par :
8
>
> x xj 1 (xj 1 x)
>
> 2
si x 2 [xj 1 ; xj ]
>
>
> j j 12 (xj 1 xj )
< x x

j (x) = x xj+ 1 (xj+1 x) (2.35)

>
> 2
si x 2 [x ; x ]
>
> j j+1
>
> xj xj+ 1 (xj+1 xj )
>
: 2
0 sinon
Les fonctions de base correspondant à un point milieu d’un élément sont dé…nies
par :
8 (x xj )(xj+1 x)
< si x 2 [xj ; xj+1 ]
xj+ 1 xj xj+1 xj+ 1 (2.36)
j+ 21 (x) =
: 2 2
0 sinon

71
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

L’espace Vh ; dé…ni par (2.31), est un sous-espace de H 1 (0; 1) de dimension

2n + 3; et toute fonction vh 2 Vh est dé…nie de manière unique par ses valeurs aux
sommets (xj )0 j n+1 et aux milieux xj+ 1
2 0 j n

X
n+1 X
n
vh (x) = vh (xj ) j (x) + vh xj+ 1 j+ 12 (x) ; 8x 2 [0; 1]
2
j=0 j=0

De même,l’espace V0h ; dé…ni par (2.32), est un sous-espace de H01 (0; 1) de dimen-
sion 2n + 1; et toute fonction vh 2 V0h est dé…nie de manière unique par ses valeurs
aux sommets (xj )1 j n et aux milieux xj+ 1
2 0 j n

X
n X
n
vh (x) = vh (xj ) j (x) + vh xj+ 1 j+ 12 (x) ; 8x 2 [0; 1]
2
j=1 j=0

Remarque 2.9 On pourra changer l’indexation des points (x0 ; x1=2 ; x1 ; : : : ; xn+ 1 ; xn+1 )
2
en les notant xk=2 0 k 2n+1 et celle des fonctions de base ( 1=2 ; 1 ; : : : ; n+ 1 ; n+1 )
2

Dé…nition 2.2 L’opérateur d’interpolation P2 est l’application rh dé…nie de H 1 (0; 1)

dans Vh par :
X
2N +2
1
8v 2 H (0; 1) rh v = vh xj=2 j=2
j=0

où les fonctions j sont données par (2.33) et (2.34).

Le problème approché par la méthode des éléments …nis P2 du problème de

Dirichlet (2.9), (2.10) est :

Trouver uh 2 V0h telle que :

Z 1 Z 1
0 0
8vh 2 V0h ; uh vh dx = f vh dx
0 0

qui est équivalent à

8
< Trouver uh xj=2 2 R, j = 1; : : : ; 2N + 1 tels que
P+1
2N R1 R1
: 8i = 1; : : : ; 2N + 1; uh xj=2 0 0i=2 0j=2 dx = 0 f i=2 dx
j=1

t
On pose uh = uh x1=2 ; uh (x1 ) ; : : : ; uh (xN ) ; uh xN +1=2
Alors uh est solution de
Kh u h = b h (2.37)

72
où Kh 2 R(2N +1) (2N +1) et bh 2 R2N +1 sont donnés par
Z 1 Z 1
0 0
Kh = i=2 j=2 dx et bh = f i=2 dx
0 1 i;j 2N +1 0 1 i 2N +1

Les fonctions de base k=2 ont un petit support, et la plupart des coe¢ cients de
Kh sont donc nuls.
Fonctions de forme locales
Introduisons les fonctions de forme locales fP0 ; P1 ; P2 g par
P0 = (2t 1) (t 1) ; P1 = 4t (1 t) ; P2 = t (2t 1)
1
En introduisant les noeuds de K dé…nis par a0 = 0; a1 = 2
et a2 = 1, il vient
Pi (aj ) = ij ; 80 i; j 2:
On véri…e facilement que :
8
< P0 Ti 1 (x) si x 2 Ki
81 i n; 'i (x) = P2 Ti 11 (x) si x 2 Ki 1
:
0 sinon
et que les fonctions i :

P1 Ti 1 (x) si x 2 Ki
80 i n; i (x) =
0 sinon
où les Ti sont dé…nies précédemment dans (2.28) et 2.29
8 1 0 1
< hi P0 Ti (x) si x 2 Ki
0 1 0 1
81 i n; 'i (x) = P T (x) si x 2 Ki 1
: hi 1 2 i 1
0 sinon
1 1
0 P0
hi 1
Ti (x) si x 2 Ki
80 i n; i (x) =
0 sinon
Exercice 2.3 Véri…er que la matrice de rigidité Kh est:
0 16 8 1
3 3
0
B 38 14 8 1 C
B 3 3 3 C
B 0 8 16 8
0 C
B 3 3 3 C
1B C
1 8 14 8 1
B 3 3 3 3 3 C
B . .. . .
.. .. .. . C
hB C
B . C
B 0 .. 0 C
B C
@ 1 8 14 8 A
3 3 3 3
8 16
0 3 3

73
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

et la matrice de masse :
2 3
16 2 0 0
6 2 8 2 1 7
6 7
6 0 2 16 2 0 7
6 7
h 66 1
... ... 7
7
M= 6 7
30 6 7
6 7
6 0 16 2 0 7
6 7
4 1 2 8 2 5
0 0 2 16
Corrigé de l’exercice 2.3 Sur chaque sous-domaine élémentaire, on approche
la fonction v(x) par un polynôme de degré 2. Sur chaque intervalle [xk ; xk+1 ] on
prend pour v(x) le polynôme d’interpolation de Lagrange sur les trois valeurs
vk ; vk+ 1 ; vk+1 : Dans ces conditiond on démontre que l’approximation s’écrit :
2
v(x) = vk 0 (x) + vk+ 1 1=2 (x) + vk+1 1 (x) avec i (x) = P2i x hxk et
2

P0 (t) = (2t 1) (t 1) ; P1 (t) = 4t (1 t) ; P2 (t) = t (2t 1)

Z xk+1 Z xk+1
x xk x xk
i j dx = P2i P2j dx
xk xk h h
Z 1
1
= h P2i (t) P2j (t) dt avec i,j=0, ; 1
0 2
La matrice
2 de masse globale pour la condition 3 de Dirichlet :
16 2 0 0
6 2 8 2 1 7
6 7
6 0 2 16 2 0 7
6 7
6 ... ... 7
h 6 1 7
M = 30 6 7
6 7
6 7
6 0 16 2 0 7
6 7
4 1 2 8 2 5
0 0 2 16
La matrice
0 de masse élémentaire
1 vaut :
4 2 1
ck = 1 @ 2 16 2 A
M 30
1 2 4
Z xk+1 Z xk+1
0 0 1 x xk x xk
i j dx = 2
P2i0 0
P2j dx
xk h xk h h
Z
1 1 1
= P2i (t) P2j (t) dt avec i; j = 0; ; 1
h 0 2

74
La matrice de rigidité élémentaire :
0 1
K11 K12 K13
Ke = @ K21 K22 K23 A
K31 K32 K33
0 R1 R1 R 1 0 dP21
dP0 2 dP0 dP1
dt dt 0 dP dt
B R0 dt R 01 dt dt
dP1 2
R 1 dP
dt dt
C
Ke = @ 01 dP1 dP0
dt 1 dP2
dt 0 dt dt dt A
R1 dt dt
dP2 dP0
R01 dt
dP2 dP1
R1 2
0 dt dt
dt 0 dt dt
dt 0 dP dt
2
dt
dP0 dP1 dP2
dt
= 4t 3; dt = 8t + 4; dt = 4t 1
R 1 dP0 dP0 R1 h i1
2 (4t 3)3 1
K11 = 0 dt dt dt = 0 (4t 3) dt = 12
= 12 + 27
12
= 28
12
= 37
h 0 i1
R 1 dP1 dP1 R1 2 ( 8t+4)3
K22 = 0 dt dt dt = 0 ( 8t + 4) dt = 24
= 64
24
64
+ 24 = 16
3
h i 0
R 1 dP2 dP2 R1 3 1
K33 = 0 dt dt dt = 0 (4t 1)2 dt = (4t121) = 27
12
1
+ 12 = 28
12
= 37
R 1 0 dP1 R1 0R1
K12 = 0 dP dt dt
dt = 0
(4t 3) ( 8t + 4) dt = 0
( 32t2 + 40t 12) dt =
h 3 2
i 1
32 t3 + 40 t2 12t = 38 = K21 :
R 1 0 dP2 R0 1 R1
K13 = 0 dP dt dt
dt = 0
(4t 3) (4t 1) dt = 0 (16t2 16t + 3) dt =
h 3 2
i1
16 t3 16 t2 + 3t = 13 = K31 :
R 1 2 dP1 0 R 1 R1
K32 = 0 dP dt dt
dt = 0
(4t 1) ( 8t + 4) dt = 0
( 32t2 + 24t 4) dt =
h 3 2
i 1
32 t3 + 24 t2 4t = 38 = K23 :
0
Donc 0 1
7 8 1
Ke = 13 @ 8 16 8 A
1 8 7
La matrice de rigidité globale est :
Pn
K = h K (k) avec
k=1
K (k) = 2
[] 3
7 8 1
6 8 16 8 7
6 7
6 1 8 7 7
6 7
6 7
6
1 6
7
K1 = 3 6 7
7
6 7
6 7
6 7
6 7
4 5

75
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

2 3
6 7
6 7
6 7 8 1 7
6 7
6 8 16 8 7
6
1 6
7
K2 = 3 6 1 8 7 7
7
6 7
6 7
6 7
6 7
4 5
2 3
6 7
6 7
6 7
6 7
6 7
6
1 6
7
K3 = 3 6 7 8 1 7
7
6 8 16 8 7
6 7
6 1 8 7 7
6 7
4 5
2 3
6 7
6 7
6 7
6 7
6 7
6
1 6
7
K4 = 3 6 7
7
6 7
6 7
6 7 8 1 7
6 7
4 8 16 8 5
1 8 7
la matrice
P de rigidité globale
K= Ki

En adaptant les arguments de la méthode P1 , on peut montrer le résultat de

convergence suivant :

Théorème 2.6 (Convergence de la méthode P2 ). Soit u 2 H01 (0; 1) la solution de

(2.9), (2.10) et uh 2 Vh la solution de (2.37). Alors, on a :

lim ku uh kH 1 (0;1) = 0
h!0

De plus, si u 2 H 3 (0; 1) (ce qui est vrai si f 2 H 1 (0; 1)) alors il existe une
constante C > 0 indépendante de h telle que :

ku uh kH 1 (0;1) Ch2 u(3) L2 (0;1)

76
On dit que la convergence est quadratique.

Daprès le théorème de convergence, on voit que l’avantage de la méthode P2

est une accélération de la convergence (quadratique au lieu de linéaire) lorsque la
solution est régulière (u 2 H 3 (0; 1)): Le désavantage est que le système Kh uh = bh
est plus couteux à résoudre car la matrice Kh n’est pas tridiagonale mais pentadi-
agonale et il y a deux fois plus d’inconnues (exactement 2n + 1 au lieu de n pour
les éléments …nis P1 ) donc la matrice est deux fois plus grande. En particulier, on
remarque que si u 2 = H 3 (0; 1) ; il n’y a aucun intérêt à employer la méthode P2 .

Exercice 2.4 On considère le problème de convection-di¤usion suivant :

d u 2
u(x) dx2
= f (x) 8x 2 ]0; 1[ =
(P )
u (0) = u (1) = 0

où et sont deux réels positifs et f une fonction dans L2 ( ) :

Résoudre cette équation par la méthode des éléments …nis P1 et P2.

Corrigé de l’exercice 2.4 Voir TD.

2.5 Eléments …nis Lagrangiens d’ordre 3

La méthode des éléments …nis P3 repose sur l’espace discret

Vh = fv 2 C ([0; 1]) tel que vj [xj ; xj+1 ] 2 P3 pour tout 0 j n + 1g (2.38)

et sur son sous-espace

V0h = fv 2 Vh tel que v(0) = v(1) = 0g (2.39)

Pour obtenir les éléments lagrangiens d’ordre 3, il su¢ t de chercher sur chaque
sous-domaine éléùentaire, une fonction test polynomiale de degré 3. Sur chaque
maille [xk ; xk+1 ], on peut écrire :

v(x) = vk 0 (x) + vk+ 1 1=3 (x) + +vk+ 2 2=3 (x) + vk+1 1 (x)
3 3

avec
x xk
i (x) = P3i
h

77
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

Les fonctions de forme locales fP0 ; P1 ; P2 ; P3 g sont données par

P0 (t) = 29 t 13 t 2
3
(t 1)
P1 (t) = 27 2t t 32 (t 1)
P2 (t) = P1 (1 t)
P3 (t) = P2 (1 t)

Et les noeuds de K dé…nis par a0 = 0; a1 = 13 ; a2 = 2

3
et a3 = 1, il vient

Pi (aj ) = ij ; 80 i; j 3:

Dans ce qui suit on notera Pk l’ensemble des polynômes à coe¢ cients réels
d’une variable réelle de degré inférieur ou égal à k :
( )
X k
Pk = p(x) = aj xj ; aj 2 R :
j=0

Remarque Généralement, le sous-espace Vh de V est noté :

Vhk = vh 2 C ( ) tel que vhjK 2 Pk (K)

Sa dimension est égale à :

dim Vhk = (k + 1) (N + 1) N

2.6 Eléments …nis d’Hermite

Au lieu d’utiliser l’approximation de Lagrange, nous utilisons ici l’interpolation
d’Hermite en imposant à la fonction v et sa dérivée d’être continues sur tout le
domaine.
On introduit une méthode des éléments …nis d’Hermite P3 qui repose sur
l’espace discret

Vh = v 2 C 1 ([0; 1]) tel que vj [xj ; xj+1 ] 2 P3 pour tout 0 j n (2.40)

On dé…nit deux fonctions "mères"

8 2
< (1 + x) (1 2x) si 1 x 0;
(x) = (1 x)2 (1 + 2x) si 0 x 1;
:
0 si jxj > 1;

78
 Figure 2.2: Les fonctions de base des éléments …nis d’Hermite

et 8 2
< x (1 + x) si 1 x 0;
(x) = x (1 x)2 si 0 x 1;
:
0 si jxj > 1
Si le maillage est uniforme, pour 0 j n + 1; on dé…nit les fonctions de base
(voir la …gure suivante)

x xj
j (x) = pour 0 j n + 1;
h
x xj
j (x) = pour 0 j n + 1;
h

L’espace Vh ; dé…ni par (2.40), est un sous-espace de H 1 (0; 1) de dimension

2 (n + 2) : Toute fonction vh de Vh est dé…nie de manière unique par ses valeurs et
celles de sa dérivée aux sommets (xj )0<j<n+1 ; et on a pour tout x 2 [0; 1]

X
n+1 X
n+1
vh (x) = vh (xj ) j (x) + (vh )0 (xj ) j (x)
j=0 j=0

Démonstration 2.3 Les fonctions de Vh étant de classe C 1 ; il s’agit bien d’un

sous-espace de H 1 (0; 1) : On véri…e facilement que les i ; j forment une base

79
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

de Vh en remarquant que
0 0
j (xi ) = ij ; j (xi ) = 0; j (xi ) = 0; j (xi ) = ij :

Éléments …nis hermitiens

Au lieu d’utiliser l’approximation de Lagrange, nous utilisons ici l’interpolation
d’Hermite en imposant à la fonction v(x) et à sa dérivée d’être continues sur tout
le domaine. Pour une approximation cubique, nous chercherons donc une fonction
polynomiale d’ordre 3 de la forme a3 x3 + a2 x2 + a1 x + a0 véri…ant
8
>
> v (xk ) = vk
<
v 0 (xk ) = vk0
>
> v (xk+1 ) = vk+1
: 0 0
v (xk+1 ) = vk+1
0
où les nombres vk ; vk0 ; vk+1 ; vk+1 sont des nombres quelconques. Une telle fonction
est unique et s’écrit

v(x) = vk N0 (x) + vk0 N1 (x) + vk+1 N2 (x) + vk+1

0
N3 (x)

avec
x xk
Ni (x) = Pi h
si i est pair
x xk
Ni (x) = hPi h
si i est impair
Les polynômes de base étant donnés par

P0 (x) = (x 1)2 (2x + 1)

P1 (x) = x (x 1)2
P2 (x) = P0 (1 x)
P3 (x) = P1 (1 x)
La matrice de masse élémentaire :
0 1
156 32 54 13
B 3 C
Mc = 1 B 32 4 13 C
420 @ 54 13 156 32 A
13 3 32 4
Puis la matrice de rigidité élémentaire :
0 1
36 3 36 3
B 1 C
Kb= 1 B 3 4 3 C
30 @ 36 3 36 3 A
3 1 3 4
L’expansion de la matrice de masse s’écrit :

80
Exemple 2.1 On veut résoude numériquement l’équation des plaques en dimen-
sion 1: 0000
u =f dans ]0; 1[
0 0 (2.41)
u (0) = u (1) = u (0) = u (1) = 0;
qui admet une unique solution u 2 H02 (0; 1) si f 2 L2 (0; 1) :Vh n’est pas seulement
un sous-espace de H 1 (0; 1) ; mais est aussi un sous-espace de H2 (0; 1) (ce qui n’est
pas le cas pour les éléments …nis de Lagrange).
Pour résoudre (2.41) nous considérons le sous espace :

V0h = fv 2 Vh tel que v (0) = v (1) = v 0 (0) = v 0 (1) = 0g : (2.42)

L’espace Vh , et son sous espace V0h dé…ni par (2.42), sont des sous-espaces de
H (0; 1) ; et de H02 (0; 1) respectivement, de dimensions 2(n + 2) ; et 2n respective-
2

ment. Toute fonction vh de V0h est dé…nie de manière unique par ses valeurs et
celles de sa dérivée aux sommets (xj )1 j n ; et on a pour tout x 2 [0; 1]

X
n X
n
vh (x) = vh (xj ) j (x) + (vh )0 (xj ) j (x)
j=1 j=1

Démonstration 2.4 Soit vh 2 Vh : elle est de classe C 1 sur [0; 1] et C 2 par

morceaux. Donc sa dérivée vh0 ; étant continue et C 1 par morceaux, appartient bien
à H 1 (0; 1).
Par conséquent, vh est un élément de H 2 (0; 1) :

2.6.1 Résolution pratique de l’équation des plaques par la

méthode des éléments …nis d’Hermite P3 .
La formulation variationnelle de l’approximation interne est :
trouver uh 2 V0h tel que
Z 1 Z 1
u00h (x) vh00 (x) dx = f (x) vh (x) dx 8vh 2 V0h (2.43)
0 0

0
On décompose uh sur la base des j ; j 1 j n et on note uh = uh (xj ) ; uh (xj ) 1 j n
le vecteur de ses coordonnées dans cette base. La formulation variationnelle (2.43)
revient à résoudre dans R2n un système linéaire.

Kh u h = b h

Exercice 2.5 Calculer explicitement la matrice de rigidité Kh pour (2.43)

81
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

X
N X
N
uh (x) = uh (xj ) j (x) + (uh )0 (xj ) j (x)
j=1 j=1

Trouver uh 2 Vh tel que

Z 1 Z 1
00 00
uh (x) vh (x) dx = f (x) vh (x) dx 8vh 2 V0h
0 0

X
N X
N Z 1
00 00 00 00
uh (xj ) j (x) i (x) + u0h (xj ) j (x) i (x) = f (x) i (x) dx; 1 i N:
j=1 j=1 0

X
N X
N Z 1
00 00 00 00
uh (xj ) j (x) i (x)+ u0h (xj ) j (x) i (x) = f (x) i (x) dx; 1 i N:
j=1 j=1 0

C’est un système linéaire à 2N équations et 2N inconnus uj et u0j 1 j N:

uh (xj ) = ui ; u0h (xj ) = u0i ;
On peut l’écrire sous la forme suivante :

A B u F
=
Bt C u0 F
R 00
avec A = (aij )1 i;j N où aij = (x) 00i (x) dx
R 00 j
B = (bij )1 i;j N où bij = (x) 00i (x) dx
R 00j 00
C = (cij )1 i;j N où bij = j (x) i (x) dx
B t est la transposée de B.
t
u = (u1 ; u2 ; : : : ; uN )t u0 = (u01 ; u02 ; : : : ; u0N )

8 3
> x x x xj 2
2 hj < 3 h
+1 si x 2 [xj 1 xj ]
x x 3 x xj 2
j (x) = 2 hj 3 +1 si x 2 [xj xj+1 ]
>
: h
0 ailleurs
8
> x xj 3 x xj 2 x xj
< h
+2 h
+ h
si x 2 [xj 1 xj ]
(x) = x xj 3 x xj 2 x xj
j
> h
2 h
+ h
si x 2 [xj xj+1 ]
: 0 ailleurs
Z
00 00
aij = j (x) i (x) dx

82
 De même nous introduisons les fonctions de forme locales f 0 ; 1; 0; 1g par

0 (t) = (1 t)2 (1 + 2t)

1 (t) = 2t3 + 3t2
0 (t) = hi t (1 t)2
1 (t) = hi t2 (t 1)

k est dans [0 1] la fonction correspondant à la restriction de i dans [xi xi+1 ]

k est dans [0 1] la fonction correspondant à la restriction de i dans [xi xi+1 ]
hi = xi+1 xi

00 1 00
81 i n; i (x) = k (t)
h2i
00 1 00
81 i n; i (x) = 2 k (t)
hi

R1
(Ah )ii = 0
('00i (x))2 dx
P
n R
= Kj
('00i (x))2 dx
Rj=0 R
= Ki
('00i (x))2 dx +
('00i (x))2 dx Ki
R 1
00 1 2 R 00 1 2
= Ki 1 h21 1 Ti 1 (x) dx + Ki h12 0 Ti (x) dx
R i 1 R i
= hi1 1 K 001 (t)2 dt + h1i K 000 (t)2 dt
= hi1 1 Ae22 + h1i Ae11

R1 00 2
(Ch )ii = 0
( i (x)) dx
P
n R 00 2
= Kj
( i (x)) dx
Rj=0 00 2 R 00 2
= Ki
( i (x)) dx +
(x)) dx Ki
( i
R 1
1
1 2 R 1 2
= Ki 1 h 2 00
1 Ti 1 (x) dx + Ki h12 00
0 Ti (x) dx
R i 1 R i
= hi1 1 K 001 (t)2 dt + h1i K 000 (t)2 dt
= hi1 1 C22
e
+ h1i C11e

On obtient en dé…nitive l’expression suivante de la matrice élémentaire relative

à un élément [xi xi+1 ]
0 1
12 6hi 12 6hi
1 B 6hi 4h2i 6hi 2h2i C
Ki = 2 B C
hi @ 12 6hi 12 6hi A
6hi 2h2i 6hi 4h2i

83
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

2.7 Introduction des éléments …nis en dimension

N=2 ou 3
2.7.1 Introduction
On considère le problème de Dirichlet homogène.
ouvert borné de RN et f 2 L2 ( )

4u = f sur
(2.44)
u=0 sur =@

Il existe une solution unique dans H01 ( )

Dans tout ce qui suit nous supposerons que le domaine est polyèdrique
(polygonal si N = 2), a…n que nous puissions le mailler exactement.

2.7.2 Éléments …nis triangulaires

N-simplexe
Tout commence par la dé…nition d’un maillage du domaine par des triangles en
dimensin N = 2 et des tétraèdres en dimension N = 3. On regroupe les triangles
et les tétraèdres dans la famille plus générale des N -simplexes.

Dé…nition 2.3 On appelle N -simplexe K de RN l’enveloppe convexe de (N + 1)

points (aj )1 j N +1 de RN ; appelés sommets de K:
On dit que le N -simplexe K est non dégénéré si les points (aj )1 j N +1 n’appartiennent
pas à un même hyperplan de RN :

Remarque 2.10 Un 2-simplexe est un triangle et un 3-simplexe est un tétraèdre.

Si on note (ai;j )1 i N les coordonnées du vecteur aj ; le N-simplexe est non

dégénéré si la matrice
0 1
a1;1 a1;2 : : : a1;N +1
B a2;1 a2;2 : : : a2;N +1 C
B C
B .. .. .. C
A=B . . . C (2.45)
B C
@ aN;1 aN;2 : : : aN;N +1 A
1 1 ::: 1

est inversible. Un N-simplexe a autant de faces que de sommets, qui sont

elles-mêmes des (N-1)-simplexes.

84
Maillage : descrption et structure
Dé…nition 2.4 Soit un ouvert connexe polyédrique de RN ; Un maillage triangu-
laire ou une triangulation de est un ensemble Th de N-simplexes (non-dégénérés)
(Ki )1 i n qui véri…ent
n
1. Ki et = [ Ki
i=1

2. L’intersection Ki \Kj de deux N -simplexes distincts est un m-simplexe, avec

0 m N 1; dont tous les sommets sont aussi des sommets de Ki et Kj :

Les sommets ou noeuds du maillage Th sont les sommets des N-simplexes Ki qui
le composent.

2.7.3 Coordonnées barycentriques

Dans un N -simplexe K il est commande d’utiliser des coordonnées barycentriques
au lieu de coordonnées cartésiennes usuelles. Les coordonnées barycentriques
( j )1 j N +1 de x 2 RN sont dé…nies par

X
N +1 X
N +1

j = 1; ai;j j = xi pour 1 i N;
j=1 j=1

qui admet bien une unique solution car la matrice A, dé…nie par (2.45), est in-
versible.
Remarquons que les j sont des fonctions a¢ nes de x. On véri…e alors que

K = x 2 RN tel que j (x) 0 pour 1 j N +1 ;

et que les (N + 1) faces de K sont les intersections de K et des hyperplans j (x) =

0; 1 j N + 1:

2.7.4 Maillage en dimension N=2 ou 3

Soit un ouvert connexe polyèdrique de RN : Un maillage ou une triangulation de
. est un ensemble Th de N-simplexes (non dégénérés) (Ki )1 i n qui véri…ent

1. Ki
2. en dimension N = 2, l’intersection Ki \Kj de deux triangles distincts est soit
vide, soit réduite à un sommet commun, soit à une arête commune entière
(en dimension N = 3, l’intersection est soit vide, soit un sommet commun,
soit une face commune entière, soit une arête commune entière.

85
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

Les sommets (ou noeuds) du maillage Th sont les sommets des N-simplexes
Ki qui le composent. Par convention, le paramètre h désigne le maximum des
diamètres des N-simplexes Ki
Exemple de maillage en dimension N=2 :

Exemple de maillage triangulaire en dimension N=2

Exemples de situations interdites pour un maillage triangulaire

Exemples de maillage non conforme

2.7.5 Treillis d’ordre k

Pour tout entier k 1 on appelle treillis d’ordre k l’ensemble (…ni) :
1 k 1
k = x 2 K tel que j (x) 2 0; ; : : : ; ;1 pour 1 j N (2.46)
k k

86
 dont les points sont notés (b
aj )1 j nk :
Pour k = 1 il s’agit de l’ensemble des sommets de K, et pour k = 2 des sommets
et des points milieux des arêtes reliant deux sommets. Dans le cas général, k est
un ensemble …ni de points (b aj )1 j nk :

2.7.6 Les Mailles

Les mailles sont des N -simplexes (triangles en 2-D, tétraèdres en 3-D) .

Dimension N=2

Treillis d’ordre 1, 2 et 3 pour un tétraèdre

Les ronds représentent les points du treillis

Dimension N=3

Treillis d’ordre 1, 2 et 3 pour un triangle

87
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

2.7.7 Ensemble de polynômes Pk

Soit Pk l’ensemble des polynômes à coe¢ cients réels de RN dans R de degré in-
férieur ou égal à k, c’est à dire que tout p 2 Pk s’écrit sous la forme :
X iN
i1
p (x) = i1 ;:::;iN x1 : : : xN avec x = (x1 ; : : : ; xN )
i1 ;:::;iN 0;i1 +:::+iN k

L’intérêt de la notion de treillis k d’un N-simplexe K est qu’il permet de

caractériser tous les polynômes de Pk (on dit que k est unisolvant pour Pk ).
Degré k=1 : polynômes a¢ nes
X
p(x) = 0 + i xi avec x = (x1 ; : : : ; xN ) :
1 i N

Unisolvance de k pour Pk
Lemme 2.2 Tout polynôme de Pk est déterminé de manière unique par ses valeurs
aux points (baj )1 j nk du treillis k : Autrement dit, il existe une base j 1 j nk
de Pk telle que
j (b
ai ) = ij ; 1 i; j nk :
Démonstration 2.5 (k = 1): On véri…e que tout polynôme p 2 P1 se met sous
la forme
X
N +1
p(x) = p (b
aj ) j (x)
j=1

car les coordonnées barycentriques j (x) sont a¢ nes en x: Comme j (b

ak ) = jk ;
la base recherchée est j = j :
Continuité à l’interface entre 2 mailles
Lemme 2.3 Soit K et K 0 deux N-simplexes ayant une face commune = @K \
@K 0 : Soit un entier k 1. Alors, leurs treillis d’ordre k, k et 0k coïncident sur
cette face . De plus, étant donné pK et pK 0 deux polynômes de Pk ; la fonction v
dé…nie par :
pK (x) si x 2 K
v(x) =
pK 0 (x) si x 2 K 0
est continue sur K [ K 0 ; si et seulement si pK et pK 0 ont des valeurs qui coïncident
aux points du treillis sur la face commune :
Démonstration 2.6 Il est clair que la restriction à une face de K de son treillis
d’ordre k est aussi un treillis d’ordre k dans l’hyperplan contenant cette face,
qui ne dépend que des sommets de cette face, les treillis k et 0k coïncident sur
leur face commune : Si les polynômes pK et pK 0 coïncident aux points de k \ ;
alors par application du lemme précédente 2.2 ils sont égaux sur ; ce qui prouve
la continuité.

88
2.7.8 Eléments …nis Pk

Dé…nition 2.5 Etant donné un maillage Th d’un ouvert ; la méthode des élé-
ments …nis Pk ; ou éléments …nis triangulaires de Lagrange d’ordre k, associée à ce
maillage, est dé…nie par l’espace discret.

Vh = v 2 C tel que vjKi 2 Pk pour tout Ki 2 Th

On appelle noeuds des degrés de liberté l”ensemble des points (distincts) (b

ai )1 i ndl
des treillis d’ordre k de chacun des N-simplexes Ki 2 Th :
On appelle degrés de liberté d’une fonction v 2 Vh l’ensemble des valeurs de v
en ces noeuds (b ai )1 i ndl :
On dé…nit aussi le sous-espace V0h par

V0h = fv 2 Vh tel que v = 0 sur @ g

Proposition 2.1 L’espace Vh est un sous-espace de H 1 ( ) dont la dimension est

le nombre de degrés de liberté,et il existe une base ( i )1 i ndl de Vh dé…nie par

i (b
aj ) = ij 1 i; j ndl ;

telle que
ndl
X
v (x) = v (b
aj ) i (x)
i=1

Remarque 2.11 L’appellation "éléments …nis de Lagrange" veut dire que toute
fonction de l’espace Vh est caractérisé par ses valeurs ponctuelles (ses degrés de
liberté) aux noeuds (b
aj ) :
On parle d’éléments …nis de Hermite si les degrés de liberté sont les valeurs de
la fonction et de ses dérivées partielles d’ordre 1.

Fonction P1 en dimension 2

89
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

Exemples

N=2, k=1,

P1=vectf1; x; yg ; dimension =3

N=2, k=2,

P2=vectf1; x; y; x2 ; xy; y 2 g ; dimension =6

N=2, k=3,

P3=vectf1; x; y; x2 ; xy; y 2 ; x3 ; x2 y; xy 2 ; y 3 g ; dimension =10

N=3, k=1,

P1=vectf1; x; y; zg ; dimension =4

N=3, k=2,

P2=vectf1; x; y; z; x2 ; xy; y 2 ; yz; z 2 g ; dimension =10

90
 Approximation par un domaine polyédrique
h d’un ouvert régulier

Graphe d’une fonction 'A relative à un sommet A

P x
8v 2 Vh ; v(X) = v(A)'A (X); X = 2 R2
A2Th y

On suppose que le domaine est polyédrique (polygonal si N = 2) et on le

découpe en triangles Tl :On considère une trianglulation du domaine : Nous allons
construire un espace Vh de fonctions vh , linéaires en x et y sur chaque triangle (de
la forme a + bx + cy) et continues sur le domaine :

91
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

Lemme 2.4 Soit un triangle T dont les sommets sont points Si de coordonnées
(xi ; yi ) ; 1 i 3: Il existe un polynôme p(x; y) 2 P1 (de la forme a + bx + cy)
unique prenant des valeurs données arbitraires pi ; 1 i 3; aux sommets Si :
Démonstration 2.7 Soit p(x; y) = a + bx + cy. On a :
a + bxi + cyi = pi ; 1 i 3 (2.47)
La matrice du système s’écrit :
0 1
1 x1 y1
@ 1 x2 y2 A
1 x3 y3
et son déterminant est di¤érent de zéro dès que les trois points S1 ; S2 ; S3 sont non
alignés. On peut donc déterminer de façon unique les quantités a, b, et c.
On appelle coordonnées barycentriques d’un point M (x; y) par rapport aux
points Ai (xi ; yi ); 1 i 3;non alignés, les scalaires i (x; y) (ou i (M )) tels que
X
3

i (x; y) = 1; (2.48)
i=1

X
3

i (x; y) xi = x; (2.49)
i=1

X
3

i (x; y) yi = y; (2.50)
i=1
Géométriquement, i (M ) représente l’aire du triangle M Aj Ak (Aj et Ak sont
les deux autres sommets du triangle).
D’après le lemme précédent, on peut toujours trouver les scalaires i (x; y) ; 1
i 3; satisfaisant les relations (2.48) à (2.50). Ces quantités sont linéaires en x et
y. Le polynôme p(x; y) du Lemme 2.2 s’exprime de façon simple en fonction des
coordonnées barycentriques.
X
3
p(x; y) = i (x; y) pi
i=1

On a les relations suivantes

0 si i 6= j
i (Sj ) =
1 si i = j
Soit Vh l’espace des fonctions vh ; dont la restriction à chaque triangle Tl ap-
partient à P1 ; et déterminées par leurs valeurs vl aux sommets Sl des trinagles
Tl :

92
Dé…nition 2.6 Vh = vh 2 C 0 ; vh jTl 2 P1 ; 8l 2 f1; : : : ; Ntri g
où P1 est l’espace vectoriel des polynômes de degré au plus 1.
V0h = fvh 2 Vh ; vh j = 0g

Proposition 2.2 Vh H 1 ( ) :
Vh est un espace vectoriel de dimension Nsom dont une base est donnée par les
fonctions 'i 2 C 0 , 8l 2 [1; Ntri ] ; 'i jTl 2 P1 'i (Sj ) = i;j ; 1 i; j Nsom :
NP
som
8vh 2 Vh ; vh = vh (Si ) 'i
i=0
Les scalaires vh (Si ) ; i 2 f0; 1; : : : ; Nsom g sont les degrés de libertés de la fonc-
tion vh 2 Vh :

Démonstration 2.8 D’après la dé…nition de 'i ; on a : pour toute fonction v 2 Vh

X
N som

v(x) = vh (Si ) 'i (x)

i=0

puisque cette relation est vrai sur chaque triangle Tl 2 Th ; la famille ('i )1 i NS
forme donc un système générateur de Vh :

Preuve 2.4 Pour montrer que cette famille est libre, supposons une combinaison
NP
som
linéaire nulle: i 'i (x) = 0
i=0
En évaluant cette somme au sommet Sj pour j 2 f0; 1; : : : ; Nsom g ; implique
j = 0; 8j 2 f0; 1; : : : ; NS g

Corollaire 2.1 1) Les fonctions de V0h sont entièrement déterminées par les valeurs
qu’elles prennent en chacun des Ni sommets internes du maillage.
PNi
2) dimV0h = Ni et 8vh 2 V0h ; vh = vh (Si ) 'i :
i=0
3) V0h H01 ( ) :

Preuve 2.5 Pour obtenir une base de V0h , il su¢ t de retirer les fonctions de base
correspondant aux sommets situés sur @ , puisque une fonction de Vh s’annule sur
si et seulement si elle s’annule aux sommets situés sur @ ;(les fonctions sont
a¢ nes).

Proposition 2.3 Le support de la fonction 'i est formé de la réunion des triangles
admettant Si comme sommet.
Support 'i = [ Tl :
l;Si 2Tl

93
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

Calcul de uh
Soit uh 2 Vh la solution du problème (2.44) (cas de l’équation de Poisson ).
Le problème discret est :

Trouver u 2 Vh tel que

a (uh ; v) = L (v) 8v 2 V0h

On cherche alors uh sous la forme

Ni
X
uh = uh (Si ) 'i (x)
i=1

On pose
ui = uh (Si )
Les ui sont donnés par la résolution du système linéaire
Ni
X
a 'i ; 'j ui = L('j ); 1 j Ni :
i;j=1

Les coe¢ cients a 'i ; 'j ne sont di¤érents de zéro que si les supports des fonctions
'i et 'j ont une intersection non vide, c’est à dire si les points Si et Sj sont les
sommets d’un même triangle.
Les ui pour i 2 f1; : : : ; Ni g sont déterminés en résolvant le système Au = b
Avec u = (u1 ; : : : ; uNi ) ;
Z Ntri Z
X
A = a ' i ; 'j = r'i r'j = r'i r'j
k=1 Tk

et
Ntri Z
X
bi = f 'i
k=1 Tk

2.7.9 Technique de l’élément de référence

Les calculs des fonctions de forme et des matrices et second-membres élémentaires
peuvent s’e¤ectuer en se ramenant à un élément de référence simple.
Dans le cas d’éléments triangulaires on choisit le triangle rectangle isocèle unité
de sommets : a1 = (0; 0) ; a2 = (1; 0) ; a3 = (0; 1)
Transformation a¢ ne entre l’élément de référence et un triangle

94
 Dans ce triangle de référence les coordonnées barycentriques admettent les
b et yb :
expressions suivantes en fonction des coordonnées cartésiennes notées x
8
< b1 = 1 x
> b yb
b2 = x
b
>
: b = yb
3

Par transformation a¢ ne inversible on fait correspondre à ce triangle de référence

un élément triangulaire droit quelconque T de sommets

A1 (x1 ; y1 ) ; A2 (x2 ; y2 ) ; A3 (x3 ; y3 )

Pour cela il su¢ t de construire une transformation a¢ ne F qui associe respective-

ment les sommets a1 ; a2 ; a3 du triangle quelconque.
Il est très commode pour cela d’utiliser les fonctions coordonnées barycen-
triques.
On obtient les coordonnées x; y du point M de T correspondant au point
b x; yb) du triangle de référence selon :
m(b

x = x1 b1 (b
x; yb) + x2 b2 (b
x; yb) + x3 b3 (b
x; yb)

y = y1 b1 (b
x; yb) + y2 b2 (b
x; yb) + y3 b3 (b
x; yb)
Le jacobien de la transformation est égàal deux fois l’aire algébrique du triangle T
:
@x @x
x2 x1 x3 x1
det ( J) = @y @@yyb =
@b
x
= 2:Aire(T )
@b
x @ yb
y2 y1 y3 y1
car c’est le rapport de l’aire du triangle T et celle du triangle de référence (qui
vaut 12 ):

95
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

Cette transformation a¢ ne qui associe les sommets a1 ; a2 ; a3 du triangle de

référence aux sommets A1 ; A2 ; A3 , est donc inversible si les points A1 ; A2 ; A3 ne
sont pas alignés. Pour tout point M correspondant au point m b on a :

bi (b
x; yb) = i (x; y) 8i = 1; 2; 3

autrement dit, les coordonnées barycentriques sont conservées dans la transforma-

tion a¢ ne ainsi dé…nie.
On ramène alors simplement le calcul des intégrales sur un T quelconque au
calcul homologue sur le triangle de référence Tb selon :
Z Z Z Z
(x; y) dxdy = 2Aire(T ) x; yb) ; y (b
(x (b x; yb)) db
xdb
y
T Tb

@ i @ bi @b
x @ bi @b
y
= +
@x @bx @x @b
y @x
et de même
@ i @ bi @b
x @ bi @b
y
= +
@y @bx @y @b
y @y
Ce qui s’exprime sous fomre matricielle selon :

1 y3 y1 y1 y2
grad i = grad bi
2Aire (T ) x1 x3 x2 x1

2.7.10 Intégrales doubles et changement de variables

On rappelle que pour un changement de variables correspondant à une transfor-
mation (x; y) 2 ! (s; t) 2 D dé…nie par : x = f1 (s; t) y = f2 (s; t)
on appelle déterminant jacobien de la transformation le déterminant :
@f1 @f2
J (f1 ; f2 ) = @s @t
@f2 @f2
@s @t

D(x;y)
que l’on note également : D(s;t)
On a alors l’égalité :
Z Z Z Z
D (x; y)
f (x; y) dxdy = f (x (s; t) ; y (s; t)) dsdt
D D (s; t)

96
2.7.11 Intégrales curvilignes
Soit une courbe C d’extrémités A et B et de longueur L, orientée de A vers B.
On note s(M ) l’abscisse curviligne d’un point M de C. Soit f une fonction qui à
tout point M de C associe le réel f (M ) = f (s). On dé…nit l’intégrale curviligne de
f sur C par l’intégrale simple suivante :
Z Z L
f (M )ds = f (s)ds
C 0

Si C admet une représentation paramétrique de paramètre t, nous aurons, en

posant f (s(t)) = g(t) :
Z Z tB
ds
f (M )ds = g (t) dt
C tA dt
où tA et tB sont respectivement les valeurs du paramètre t associées aux points A
et B de la courbe C.
Exercice
On considère le carré = ] 1; 1[2 maillé suivant la …gure :
Calculer la matrice de rigidité Kh des éléments …nis P1 appliqués au Laplacien
avec condition aux limités de Newmann.
(On utilisera les symétries du maillage).

97
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

2.7.12 Éléments …nis rectangulaires

Éléments …nis rectangulaires
Si le domaine est de type rectangulaire (c’est-à-dire que est un ouvert
polyédrique dont les faces sont perpendiculaires aux axes), on peut le mailler par
des rectangles (voir la Figure ) et utiliser une méthode d’éléments …nis adaptée.
Nous allons dé…nir des éléments …nis de type Lagrange (c’est-à-dire dont les
degrés de liberté sont des valeurs ponctuelles de fonctions), dits éléments …nis
Qk : Commençons par dé…nir un N-rectangle K de Rn comme le pavé (non-dégénéré)
N
[li ; Li ] avec 1 < li < Li < +1. On note (aj )1 j 2N les sommets de K.
i=1

Dé…nition 2.7 Soit un ouvert connexe polyédrique de RN . Un maillage rec-

tangulaire de est un ensemble Th de N -rectangles (non dégénérés) (Ki )1 i n qui
véri…ent
n
1. Ki et = [ Ki
i=1

2. l’intersection Ki \ Kj de deux N-rectangles distincts est un m-rectangle, avec

0 m N 1; dont tous les sommets sont aussi des sommets de Ki et Kj :

(En dimension N = 2, l’intersection de deux rectangles est soit vide, soit un

sommet commun, soit une face commune entière).
Figure –Treillis d’ordre 1, 2, et 3 pour un rectangle (les ronds représentent les
points du treillis).
Les sommets ou noeuds du maillage Th sont les sommets des N-rectangles Ki qui
le composent. Par convention, le paramètre h désigne le maximum des diamètres
des N-rectangles Ki :
Nous dé…nissons l’ensemble Qk des polynômes à coe¢ cients réels de RN dans
R de degré inférieur ou égal à k par rapport à chaque variable, c’est-à-dire que

98
tout p 2 Qk s’écrit sous la forme
X
p (x) = i1
i1 ;:::;iN x1 : : : xiNN avec x = (x1 ; : : : ; xN )
0 i1 k;:::;0 iN k

Remarquons que le degré total de p peut être supérieur à k, ce qui di¤érencie

l’espace Qk de Pk .
Pour tout entier k 1 on dé…nit le treillis d’ordre k du N-rectangle K comme
l’ensemble

xj lj 1 k 1
k = x 2 K tel que 2 0; ; : : : ; ;1 pour 1 j N (2.51)
Lj lj k k
Pour k = 1 il s’agit de l’ensemble des sommets de K, et pour k = 2 et N = 2
des sommets, des points milieux des arêtes reliant deux sommets, et du barycentre
(voir la Figure ).
Le treillis k d’un N-rectangle K est unisolvant pour Qk , c’est-à-dire qu’il
permet de caractériser tous les polynômes de Qk :
Lemme 2.5 Soit K un N-rectangle. Soit un entier k 1. Alors, tout polynôme
de Qk est déterminé de manière unique par ses valeurs aux points du treillis d’ordre
k, k dé…ni par (2.51).
Dé…nition 2.8 Étant donné un maillage rectangulaire Th d’un ouvert , la méth-
ode des éléments …nis Qk est dé…nie par l’espace discret
Vh = v 2 C tel que vjKi 2 Qk pour tout Ki 2 Th (2.52)
On appelle noeuds des degrés de liberté l’ensemble des points (b
ai )1 i ndl des
treillis d’ordre k de chacun des N-rectangles Ki 2 Th .
Comme dans le cas triangulaire, la Dé…nition a un sens grâce à la proposition
suivante (dont nous laissons la démonstration au lecteur en guise d’exercice).
Proposition 2.4 L’espace Vh , dé…ni par (2.52), est un sous-espace de H 1 ( )
dont la dimension est le nombre de degrés de liberté ndl . De plus, il existe une base
de Vh ('i )1 i ndl dé…nie par:
'i (b
aj ) = ij 1 i; j ndl
telle que
ndl
X
v(x) = v(b
ai )'i (x)
i=1
Comme les éléments …nis Qk sont des éléments …nis de type Lagrange, on peut
démontrer les mêmes résultats de convergence que pour la méthode des éléments
…nis Pk :

99
Pr. A. Radid CHAPITRE 2. ELÉMENTS FINIS

Fonction de base Q1 en dimension N=2

Bibliographie

G. Allaire, Analyse numérique et optimisation,Éditions de l’École Polytech-

nique 2ème édition 2012
P. Ciarlet et É. Lunéville. La méthode des éléments …nis - Tome 1: De la
théorie à la pratique. Concepts généraux,
Les Presses de l’ENSTA. 2009.
A. Quarteroni, R. Sacco, F. Saleri, Méthodes Numériques : Algorithmes,
analyse et applications, Springer Milan, 2008.
A. Quarteroni, A. Valli, Numerical Approximation of Partial Di¤erential Equa-
tions, Springer Berlin Heidelberg, 2009.
M. Renardy, R.C. Rogers, An introduction to partial di¤erential equations,
1993.

100