Méthode des éléments finis

En analyse numérique, la méthode des éléments finis est utilisée pour résoudre numériquement des équations aux dérivées partielles.

Sommaire

1 Introduction

1.1 Historique
1.2 Domaines d'application
1.3 Principe

1.3.1 Les difficultés

2 Exemple de problème discret

3 Définition d'un élément fini

4 Formulation d'un élément fini

4.1 Définitions et notations
4.2 Équations fondamentales
4.3 Exemple de formulation : Barre en traction
4.4 Formulation générale (méthode directe)
4.5 Étude des fonctions N(x)
4.6 Les éléments isoparamétriques

4.6.1 Problème
4.6.2 D'où l'idée
4.6.3 Interpolation de la géométrie
4.6.4 Coordonnées locales (cas 2D)
4.6.5 Élément isoparamétrique
4.6.6 Autres classes d'éléments
4.6.7 Évaluation de K
4.6.8 Bénéfice de l'approche

4.7 Cas particulier : les éléments axisymétriques

5 Processus de calcul (cas statique)

6 Références

7 Approche mathématique/théorique

7.1 Formulation faible
7.2 Fonctions tests
7.3 Un ensemble minimal de fonctions tests
7.4 Première discrétisation
7.5 Seconde discrétisation
7.6 Problème sous forme matriciel
7.7 Algorithme

Introduction

Historique

Analyse des structures née vers 1850
- RDM (Résistance des matériaux) ⇒ calculs « manuels »
  Maxwell, Castigliano, Mohr
- Concept d'éléments finis né vers [1940]
  Newmark, Hrenikoff, Mc Henry
- Développement réel depuis 1960
  Calcul numérique sur ordinateur

Domaines d'application

Calcul de structures, électricité, hydraulique, aérodynamique, ...

Principe

Le milieu continu est « idéalisé » par la subdivision en un nombre fini d'éléments dont le comportement est représenté par un nombre finis de paramètres
La résolution du problème global, obtenu par assemblage des éléments, suit les règles qui régissent les structures discrètes.

Les difficultés

D'ordre théorique : formulation des éléments
D'ordre pratique :
- Discrétisation du milieu continu (maillage)
- Qualité des résultats (convergence de la méthode)

Exemple de problème discret

Image manquante
Pb_discret.png

Équation locale du composant e :

$\left. \begin{matrix} I_i^e = {1 \over R^e} (U_i^e - U_j^e) \\ I_j^e = {1 \over R^e} (U_j^e - U_i^e) \end{matrix} \right\} \Longrightarrow \begin{Bmatrix} I_i^e \\ I_j^e \end{Bmatrix} = {1 \over R^e} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}\begin{Bmatrix} U_i^e \\ U_j^e \end{Bmatrix} \,$

soit $\left\{ I^e \right\} = \left[ K^e \right] \left\{ U^e \right\} \,$
On écrit :
- La continuité des potentiels en chaque connexion $i \,$
- L'équilibre des courants à chaque connexion
- L'adjonction des courants externes $P^i \,$
On obtient l'équation globale du système assemblé :

$P^i = \sum_{j=1}^m \sum_{e=1}^m K_{ij}^e U_j^e \,$
soit $\left\{ P \right\} = \left[ K \right] \left\{ U \right\} \,$

Définition d'un élément fini

* Problème vrai :		* Problème idéalisé :
Image manquante Milieu_continu.png image:milieu_continu.png		Image manquante Milieu_discretise.png image:milieu_discretise.png

En calcul de structures, un élément fini est caractérisé par deux matrices :

La matrice de raideur $\left[ K \right] \,$
La matrice de masse $\left[ M \right] \,$

Formulation d'un élément fini

Définitions et notations

Déplacement : $\left\{ U \right\} \,$

C'est un vecteur dont chaque composante est également appelée degré de liberté

3 ddl de translation : $U_x,U_y,U_z \,$
3 ddl de rotation : $\theta_x,\theta_y,\theta_z \,$

Déformation : (voir Tenseur des déformations) $\left[ \epsilon \right] \,$

C'est le rapport de l'allongement à la longueur initiale.

En petites déformations, on a $(A'B')^2 = \left(\ dx + {\partial u \over \partial x} \ dx\right)^2 + \left({\partial v \over \partial x} \ dx\right)^2 \,$

Image manquante
Deformation_EF.png

$\epsilon_x = {{A'B'-AB} \over AB} \,$

Comme $AB=\ dx \,$ , on a $A'B'=(1+\epsilon_x)\ dx \,$

$\Rightarrow 2\epsilon_x+\epsilon_x^2 = 2{\partial u \over \partial x}+{\partial u \over \partial x}^2+{\partial v \over \partial x}^2 \,$

On néglige les termes d'ordre 2 :

$\epsilon_x = {\partial u \over \partial x} \,$

Remarque : $\epsilon \,$ est sans dimension

Contrainte : (voir Tenseur des contraintes) $\left[ \sigma \right] \,$

Elle représente les efforts internes qui s'appliquent dans la structure.

$\sigma_{11},\sigma_{22},\sigma_{33} \,$ : contraintes normales

$\sigma_{12},\sigma_{13},\sigma_{23} \,$ : contraintes de cisaillement

Une contrainte est homogène à une pression (N/m²)
Il existe un système de coordonées dans lequel $\left[ \sigma \right] \,$ est une matrice diagonale : $\left[ \sigma \right] = \begin{bmatrix}\sigma_X & 0 & 0 \\ 0 & \sigma_Y & 0 \\ 0 & 0 & \sigma_Z\end{bmatrix} \,$

$\sigma_X,\sigma_Y,\sigma_Z \,$ sont appelées les contraintes principales

Relations Contraintes-Déformations : (voir Loi de Hooke) ≡ lois de comportement

$\begin{Bmatrix}\sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{12} \\ \sigma_{23} \\ \sigma_{31}\end{Bmatrix} = {{E(1-\nu)} \over {(1+\nu)(1-2\nu)}} \begin{bmatrix}1 & {\nu \over {1-\nu}} & {\nu \over {1-\nu}} & 0 & 0 & 0 \\ & 1 & {\nu \over {1-\nu}} & 0 & 0 & 0 \\ & & 1 & 0 & 0 & 0 \\ & & & {(1-2\nu) \over {2(1-\nu)}} & 0 & 0 \\ & \mbox{(sym)} & & & {(1-2\nu) \over {2(1-\nu)}} & 0 \\ & & & & & {(1-2\nu) \over {2(1-\nu)}}\end{bmatrix} \begin{Bmatrix}\epsilon_{11} \\ \epsilon_{22} \\ \epsilon_{33} \\ \gamma_{12} \\ \gamma_{23} \\ \gamma_{31}\end{Bmatrix}\,$

$E \,$ : module de Young (N/m²)

$\nu \,$ : coefficient de Poisson (sans dimension)

On utilise parfois le module de cisaillement : $G = {E \over {2(1+\nu)}} \,$
Pour un matériau isotrope, il n'y a que 2 paramètres indépendants. Il y en a 6 pour un matériau orthotrope et 21 pour un matériau anisotrope
En notation matricielle, on écrit : $\left\{ \sigma \right\} = \left[D\right]\left\{ \epsilon \right\} \,$

$\left[D\right] \,$ est appelée matrice d'élasticité du matériau.

Énergie de déformation : $W \,$

$W = {1 \over 2}\int_{v}\sigma.\epsilon.\ dv \,$

Travail d'une force :

C'est le produit de la force par le déplacement de son point d'application : $\tau_F = {1 \over 2}F.U_F \,$

Moment :

C'est une force appliquée sur un ddl de type rotation

Équations fondamentales

Équations de l'elasticité linéaire :

${\partial \sigma_x \over \partial x} + {\partial \sigma_{xy} \over \partial y} + {\partial \sigma_{xz} \over \partial z} + F_x = 0\,$

Image manquante
Elasticite_lineaire

${\partial \sigma_y \over \partial y} + {\partial \sigma_{xy} \over \partial x} + {\partial \sigma_{yz} \over \partial z} + F_y = 0\,$

${\partial \sigma_z \over \partial z} + {\partial \sigma_{zx} \over \partial x} + {\partial \sigma_{zy} \over \partial y} + F_z = 0\,$

Relations déformations-déplacements :

${\epsilon_x= {\partial u \over \partial x}} \qquad {\epsilon_y= {\partial v \over \partial y}} \qquad {\epsilon_z= {\partial w \over \partial z}} \,$

$\gamma_{xz}={\partial w \over \partial x} + {\partial u \over \partial z} \,$

$\gamma_{xy}={\partial u \over \partial y} + {\partial v \over \partial x} \,$

$\gamma_{yz}={\partial w \over \partial y} + {\partial v \over \partial z} \,$

Symboliquement, on écrit $\{\epsilon\} = [S]\{U\} \,$

Signification du coefficient de Poisson :

Image manquante
Poisson_coef.png

Si on applique au barreau une contrainte $\sigma_x \,$ , on observe un rétrécissement dans la direction y correspondant à une déformation $-\nu \sigma_x \over E \,$
Quelques valeurs usuelles :

Acier : E = 2,1E11 N/m² = 210 000 MPa $\nu \,$ = 0,29
Aluminium : E = 7,41E10 N/m² = 74 000 MPa $\nu \,$ = 0,32

Remarques : On a toujours $\nu \le \,$ 0,5
Quand $\nu \to \,$ 0,5, les relations précédentes ne sont plus valables. Le matériau est dit viscoélastique.

Exemple de formulation : Barre en traction

Image manquante
Barre_traction.png

On suppose que le déplacement en tout point de la barre est donné par un polynôme du 1^er degré : $u(x) = a_1 + a_2 x \,$

On a $u(0) = u_1 \,$ et $u(L) = u_2 \,$

d'où $u(x) = u_1\left(1-{x \over L} \right) + {x \over L} u_2 \,$

qu'on écrit symboliquement : $u(x) = [N(x)]\{U\} \,$ avec $\left|\begin{matrix}[N(x)] & = & \left[ \left(1-{x \over L}\right) ; \left({x \over L}\right) \right] \\ \{U\} & = & \begin{Bmatrix}u_1 \\ u_2 \end{Bmatrix} \end{matrix}\right. \,$
On en déduit :

$\epsilon = \left[-{1 \over L} ; {1 \over L} \right]\begin{Bmatrix}u_1 \\ u_2\end{Bmatrix} = [B]\{U\} \,$

$\sigma = E\left[-{1 \over L} ; {1 \over L} \right]\begin{Bmatrix}u_1 \\ u_2\end{Bmatrix} = [D]\epsilon \,$

D'autre part, on a par définition :

$\begin{Bmatrix}F_1 \\ F_2 \end{Bmatrix} = \begin{Bmatrix}-\sigma S \\ \sigma S \end{Bmatrix} \,$ où S est l'aire de la section de la barre.

On pose :

$\{F\} = \begin{Bmatrix}F_1 \\ F_2 \end{Bmatrix} = [A]\sigma \qquad \mbox {avec} \qquad \{A\} = S \begin{Bmatrix}-1 \\ 1 \end{Bmatrix} \,$

On obtient finalement :

$\{F\} = [A][D][B]\{U\} \,$

Soit une relation du type :

$\{F\} = [K]\{U\} \qquad \mbox {avec} \qquad [K] = [A][D][B] \,$

En explicitant :

$[K] = {ES \over L} \begin{bmatrix}1 & -1 \\ -1 & 1 \end{bmatrix} \,$

On voit que la matrice de rigidité se calcule comme le produit de 3 matrice :

$[B] \,$ : Transformation des déplacements aux déformations

$[D] \,$ : Matrice d'élasticité du matériau

$[A] \,$ : Transformation des contraintes aux forces

Formulation générale (méthode directe)

La démarche est la suivante :

On exprime le déplacement $\{u(x)\}\,$ en tout point de l'élément en fonction des déplacements aux noeuds $\{u(x)\} = [N(x)]\{U\}\,$
On exprime les déformations en fonction des déplacements $\{\epsilon(x)\}=[S]\{u(x)\}\,$

d'où $\{\epsilon(x)\}=[S]\{N(x)\}\{U\}=[B(x)]\{U\}\,$

On écrit la loi de comportement du matériau qui relie les contraintes aux déformations : $\{\sigma(x)\}=[D]\{\epsilon(x)\}\,$
On écrit que le travail des forces externes appliquées à la structure pour un déplacement virtuel $\delta U \,$ est égal au travail interne des contraintes pour ce même déplacement :

$\{\delta U\}^T\{F\}=\int_V \{\delta \epsilon\}^T\{\sigma\} dv \,$

En explicitant, on a :

$\{\delta U\}^T\{F\}=\{\delta U\}^T \left (\int_V [B]^T[D][B] dv \right ) \{U\} \,$

Comme cette relation est vraie pour tout déplacement virtuel, on en déduit :

$\{F\}=[K]\{U\} \,$

avec $[K] \,$ sous sa forme plus générale : $[K]=\int_V [B]^T[D][B] dv \,$

Remarques :

La relation ci-dessus montre que $[K] \,$ est symétrique.
Le terme courant $K_{i,j} \,$ de la matrice correspond à la force qui s'exerce sur le noeud $j \,$ lorsqu'on impose un déplacement unitaire du noeud $i \,$ .

La symétrie de $[K] \,$ qui s'écrit $K_{i,j} = K_{j,i} \,$ correspond mécaniquement au théorème de réciprocité de Maxwell-Betti.

Qu'est-ce qui va différencier les différents types d'éléments finis ?

Le choix des fonctions N(x)
La nature de l'opérateur S (reliant déformations et déplacements) qui dépend du type de théorie élastique utilisée :

théorie des poutres
théorie des plaques
théorie de la contrainte ou de la déformation plane
théorie des coques
théorie des corps de révolution
théorie de l'elasticité 3D

Remarques :
Nous avons décrit le processus de formulation d'un élément fini dans le cadre de la méthode directe.(dite aussi méthode des déplacements) Il existe d'autres approches :

La méthode des résidus pondérés
L'application du principe des travaux virtuels
La minimisation de l'énergie potentielle

Toutes ces approches sont équivalentes et aboutissent à la construction de la même matrice de rigidité.

Étude des fonctions N(x)

Dans le cas général de l'élasticité tridimensionnelle, ce sont en fait des fonctions de x, y, z.
Les fonctions les plus couramment utilisées sont des polynômes.

Polynôme de degré 1 : élément linéaire (2 noeuds par arête)
Polynôme de degré 2 : élément parabolique (3 noeuds par arête)
Polynôme de degré 3 : élément cubique (4 noeuds par arête)

etc...

Les fonctions N(x) sont appelées fonctions de forme ou fonctions d'interpolation de l'élément.

Les éléments isoparamétriques

$[K]=\int_V [B]^T[D][B] dv \,$

Problème

Soit on construit $[K] \,$ pour un certain nombre d'éléments de forme et de géométrie figée ⇒ nécessité, pour mailler une structure complexe, d'utiliser un grand nombre d'éléments.
Soit on utilise des éléments à géométries variables ⇒ il faut reconstruire $[K] \,$ à chaque fois.

D'où l'idée

Pour construire la matrice de raideur d'un élément à géométrie variable, on va utiliser des fonctions d'interpolation pour décrire non seulement le champ de déplacement de l'élément mais également sa géométrie. De plus, on va travailler en coordonnées locales.

Interpolation de la géométrie

$\begin{matrix}x &=&\ &\tilde N_1(x)x_1 + \tilde N_2(x)x_2 \\ \ &\ &+&\tilde N_3(x)x_3 + \tilde N_4(x)x_4\end{matrix} \,$

Idem pour les autres coordonnées.

Image manquante
EL_isoparam_champ.png
Image:EL_isoparam_champ.png

Coordonnées locales (cas 2D)

Image manquante
EL_isoparam_coord.png

Élément isoparamétrique

Un élément est dit isoparamétrique si on prend les mêmes fonctions d'interpolation pour le déplacement et la géométrie.
$\tilde N(x) = N(x) \,$

Autres classes d'éléments

Image manquante
EL_super_param.png

Évaluation de K

La forme générale s'écrit :
$[K]=\int_x \int_y \int_z G(x,y,z) dxdydz \,$
On passe en variables locales $\xi,\eta,\zeta \,$
On a $dxdydz = \det J\ d\xi d\eta d\zeta \,$

$J \,$ s'appelle la matrice Jacobienne.

On est alors amené à calculer des intégrales du type :
$\int_{-1}^{+1} \int_{-1}^{+1} \int_{-1}^{+1} G(\xi ,\eta ,\zeta) \det J\ d\xi d\eta d\zeta \,$

Bénéfice de l'approche

On s'est ramené à un domaine d'intégration simple et invariant pour lequel on peut appliquer les formules de quadrature de gauss :
$\int_{-1}^{+1} f(x) dx = \sum_{k=1}^n H_k f(a_k) \,$ les $a_k \,$ et $H_k \,$ étant tabulés.
Les $a_k \,$ sont appelés points d'intégration de l'élément ou encore points de gauss de l'élément.

Cas particulier : les éléments axisymétriques

Image manquante
EL_axisymetrique.png

Décomposition en série de Fourier :

$u(r,\theta,z) = \sum_n u_n^s(r,z)cos (n\theta) + \sum_n u_n^s(r,z)sin (n\theta) \,$

L'axisymétrie correspond à la restriction $n=0 \,$ de cette décomposition.

$u(r,\theta,z) = u_n^s (r,z) \,$

Remarque importante : Pour utiliser ce type d'élément, le problème doit être globalement axisymétrique :

la géométrie
les conditions limites
le chargement

Processus de calcul (cas statique)

Maillage
Construction de la matrice de raideur de chaque élément $[K^e] \,$
Assemblage de la matrice globale $[K] \,$
Construction du vecteur chargement $\{F\} \,$
Élimination de certains ddl (si besoin)
Résolution : $\{U\} = [K^{-1}]\{F\} \,$
Calcul des quantités dérivées de $\{U\} \,$

$\{\epsilon^e\} = [B^e] \{U^e\} \,$

$\{\sigma^e\} = [D^e] \{\epsilon^e\} \,$

$W^e = {1 \over 2} \{\epsilon^e\}^T \{\sigma^e\} \,$

etc ...

Références

O. C. Zienkiewicz, R. L. Taylor, J.Z. Zhu : The Finite Element Method: Its Basis and Fundamentals

 Butterworth-Heinemann; 6 edition (March 21, 2005)ISBN: 0750663200

K. J. Bathe : Numerical methods in finite element analysis

 Prentice-Hall (1976) ISBN: 0136271901

Chu-Kia Wang : Matrix methods of structural analysis

 International Textbook Co; 2d ed edition (1970) ISBN: 0700222677

R. H Gallagher : Introduction aux éléments finis

 Pluralis (1977)

N. Willems : Matrix analysis for structural engineers

 Prentice-Hall (1968) ASIN: B0006BUP26

Approche mathématique/théorique

Durant cet exemple, si f est une fonction, la dérivée partielle de f par rapport à x sera notée f_x.

La meilleure façon d'introduire ce sujet est de donner un exemple simple (un « problème modèle »). Nous allons prendre l'équation de Laplace sur l'intervalle [0,1[x[0,1[. Tout d'abord, notons

$\mathrm{\Omega} := \left[ 0,1\right[ \times\left[ 0,1\right[=\{\left( x,y\right); 0 \leq x,y < 1\}.$

Soit g une fonction de $Ω$ vers F, l'ensemble des scalaires (l'ensemble des réels R ou le plan complexe C). Notre problème est de trouver une fonction u du plan R² vers F qui vérifie

u est deux fois différentiable sur R²
u_xx+u_yy=g dans $Ω$
u(x,y)=u(x+1,y)=u(x,y+1), c'est-à-dire, u est périodique en x et en y, de période 1.

Le problème peut être réécrit de la manière suivante :

u appartient à $H 2 (Ω)$ ( $H 2 (Ω)$ est un espace de Hilbert, $H^2(\Omega)=\left\{ v \in L^2(\Omega) , v,v' \in L^2(\Omega)\right\}$ )
Lu=g avec g un élément de $V=H_0^1(\Omega)=\left\{ v \in L^2(\Omega) , v \in L^2(\Omega) ,v(0)=v(1)=0 \right\}$

L étant l'opérateur différentiel défini ainsi

L f = f x x + f y y .

L est appelé opérateur de Laplace. Nous cherchons maintenant un u dans V tel que Lu=g.

Formulation faible

Soit $φ$ une forme linéaire sur V.

$\phi\left(tu+v\right)=t\phi\left(u\right)+\phi\left(v\right)$

Notons par V^* l'ensemble de ces fonctions. Ainsi, les propriétés suivantes sont équivalents

L u = g

$\forall \phi \in \mathrm{V}^*: \phi\left(\mathrm{L}u\right)=\phi\left(g\right).$

Fonctions tests

L'ensemble V^* est plus grand que ce dont on a besoin. Par conséquent, nous choisissons la fonction de la sorte

$\phi_v\left(f\right)=\int_0^1\int_0^1 v\left(x,y\right) f\left(x,y\right) dxdy,$

où $v \in V$

Maintenant, nous écrirons

\int Ω

pour l'intégrale double $\int_0^1\int_0^1$ . En intégrant par parties,

$\int_\Omega v u_{xx}\;dxdy=\int_0^1\left[v u_x\right]_0^1\;dy - \int_\Omega \left(v_x u_x\right)\;dxdy$

$\int_\Omega v u_{yy}\;dxdy=\int_0^1\left[v u_y\right]_0^1\;dx - \int_\Omega \left(v_y u_y\right)\;dxdy$

et du fait des conditions aux limites, les termes $\int_0^1 \left[ v u_x \right]_0^1\;dy$ et $\int_0^1 \left[ v u_y \right]_0^1\;dx$ s'annulent.

$\int_\Omega v\left(u_{xx}+u_{yy}\right)\;dxdy = -\int_\Omega \left(v_x u_x + v_y u_y \right)\;dxdy.$

L'écriture précédente est équivalente à :

$\psi\left(u,v\right) := \phi_v\left(\mathrm{L}u\right) = -\int_\Omega \left(u_x v_x + u_y v_y\right) = \int_\Omega g\cdot v = \phi_v \left(g\right).$

La fonction ψ de u est de v est en fait bilinéaire sur $H^2(\Omega) \times V$ et c'est la forme bilinéaire associée à L. Les fonctions v sont appelées fonctions tests.

Le problème variationnel abstrait s'écrit :

Trouver $u \in V$ tel que, quelque soit $v \in V$ , on ait

$\psi \left( u,v \right) = \phi_v \left( g \right)$

Un ensemble minimal de fonctions tests

En analyse fonctionnelle, nous savons qu'il est inutile d'essayer toutes les fonctions tests v dans V. En fait, si E={e₁,e₂,...} est un sous ensemble de V, l'ensemble des combinaisons linéaires est dense dans V. Par conséquent, on peut utiliser ces fonctions comme fonctions tests. Maintenant, nous allons rechercher une solution au problème

$\textrm{(*)}\ \ \ \forall e_j \in \mathrm{E} : \psi\left(u,e_j\right)=\int_\Omega g\cdot e_j.$

Première discrétisation

De manière à rendre ce processus implémentable, nous devons choisir un sous ensemble fini de E. Prenons F=F_n={e₁,e₂,...,e_n} un ensemble fini, F est un sous ensemble de E, et nous voulons résoudre le problème

$\textrm{(**)}\ \ \forall e_j \in \mathrm{F}_n : \psi\left(u_n,e_j\right)=\int_\Omega g\cdot e_j$

avec l'idée que quand n tend vers l'infini (et F_n croit vers E), les solutions u_n devront converger vers la solution u de (*).

Ce problème est indéterminé, par conséquent, une autre discrétisation est nécessaire.

Seconde discrétisation

Le problème variationnel dans F a une solution et est unique (conséquence du théorème de Lax Milgram sur V, $ψ$ est une forme bilinéaire sur $V \times V$ , V est elliptique et $φ$ est une forme bilinéaire continue.

$u_n = \sum_{k=1}^n a_k e_k$

$g_n = \sum_{k=1}^n b_k e_k.$

Remplaçons dans (**) et en développant, on obtient :

$\textrm{(***) } \sum_{k=1}^n a_k \psi\left(e_k, e_j\right) = \sum_{k=1}^n b_k \int_\Omega e_k e_j,\; j=1,\dots,n.$

Se pose maintenant la question : comment obtenir un g_n acceptable ? Il existe plusieurs approches qui dépendent du choix de l'ensemble E. Si E est une base de Fourier, g_n est la projection de g sur l'ensemble des combinaisons linéaires de F, mais d'autres approches existent.

Une autre approche peut être de calculer directement $\int_\Omega g\cdot e_j$ à partir d'une méthode de quadrature numérique.

Problème sous forme matriciel

Le lecteur astucieux aura noté que le problème précédent peut être formulé sous forme matricielle :

P a = Q b

où a, b sont des vecteurs colonnes

$a = \begin{bmatrix} a_1\\ a_2\\ \vdots\\ a_n \end{bmatrix}, b = \begin{bmatrix} b_1\\ b_2\\ \vdots\\ b_n \end{bmatrix}.$

Les matrices P et Q sont données par (***):

$\mathrm{P}_{jk} = \psi\left(e_k, e_j\right),$

Q j k = \int e k e j . Ω

P est appelé matrice de rigidité et Q, matrice de masse.

Algorithme

La méthode des éléments finis doit être conduite ainsi

On calcule la matrice de rigidité P
On détermine le membre de droite, soit par l'intermédiaire de la matrice de masse, soit directement par une méthode de quadrature
On résoud le problème Pa=Qb. Selon la base qui a été choisie et selon les données du problème, il faut choisir la méthode d'inversion la plus efficace.
On peut écrire u grâce à a.

Note : Dans notre problème,

si l'on choisit une base lagrangienne pour F, la matrice P sera tridiagonale.
étant donné que notre forme bilinéaire $ψ$ est symétrique, P est symétrique
$ψ$ étant V-elliptique, P est symétrique définie positive

on peut utiliser une méthode spécifique pour les systèmes linéaires à matrice tridiagonale ou une factorisation de Cholesky par bande.

Q_jk =	∫	e_ke_j.
	Ω