Table of Contents
This document is under active development and discussion!

If you find errors or omissions in this document, please don’t hesitate to submit an issue or open a pull request with a fix. We also encourage you to ask questions and discuss any aspects of the project on the Feel++ Gitter forum. New contributors are always welcome!

Functional analysis

1. Outils d’analyse fonctionnelle

include_latex_macros::docs/math/fem/macros.tex[]

1.1. Quelques rappels

1.1.1. Normes et produits scalaires

Soit \(E\) un espace vectoriel.

stem:[\|.\|] : stem:[E \rightarrow \RR] est une *norme* sur stem:[E] ssi elle vérifie :
(N1)

\[ \left( \| x \| = 0 \right) \Longrightarrow (x=0) \]

(N2)

\[\forall\, \lambda\in\RR,\; \forall x\in E, \quad \| \lambda x \| = |\lambda| \; \| x \| \]

(N3)

\[\forall\, x,y \in E, \quad \| x+ y \| \le \|x \| + \|y\|\] (inégalité triangulaire)

Pour \(E=\RR^n\) et \(x=(x_1,\ldots,x_n) \in\RR^n\), on définit les normes

\[ \| x \|1 = \sum{i=1}^n |x_i| \qquad \| x \|2 = \left( \sum{i=1}^n x_i^2 \right)^{1/2} \qquad \| x \|\infty = \sup{i} |x_i| \]

Produit scalaire

On appelle produit scalaire sur \(E\) toute forme bilinéaire symétrique définie positive.

\(\quad<.,.>\) : \(E\times E \rightarrow \RR\) est donc un produit scalaire sur \(E\) ssi il vérifie :

S1
\[\forall\; x,y \in E, \quad <x,y> = <y,x>\]
S2
\[\forall\; x_1,x_2,y \in E, \quad <x_1+x_2,y> = <x_1,y>+ <x_2,y>\]
S3
\[\forall\; x,y \in E, \, \forall\, \lambda\in\RR,\quad <\lambda x,y> = \lambda <x,y>\]
S4

\(\forall\; x \in E, x\ne 0, \quad <x,x>\; > 0\)

A partir d’un produit scalaire, on peut définir une norme induite : \( \| x \| = \sqrt{<x,x>} \)
On a alors, d’après (N3), l’inégalité de Cauchy-Schwarz : \({ | <x,y> | \le \| x \| \; \| y \| }\)

Pour \(E=\RR^n\), on définit le produit scalaire \(<x,y> = \sum_{i=1}^n x_i \, y_i.\) Sa norme induite est \(\| . \|_2\) définie précédemment.

Espace normé

Un espace vectoriel muni d’une norme est appelé espace normé.

Espace préhilbertien

Un espace vectoriel muni d’un produit scalaire est appelé espace préhilbertien. En particulier, c’est donc un espace normé pour la norme induite.

1.1.2. Suites de Cauchy - espaces complets

Soit \(E\) un espace vectoriel et \((x_n)_n\) une suite de \(E\). \((x_n)_n\) est une suite de Cauchy si et seulement si

\[\forall \varepsilon > 0,\;\; \exists N / \forall p>N, \forall q>N, \quad \|x_p - x_q \| < \varepsilon\]

Toute suite convergente est de Cauchy.

La réciproque de [th-conv-cauchy] est fausse.
Espace complet

Un espace vectoriel est complet ssi toute suite de Cauchy y est convergente.

Espace de Banach

Un espace normé complet est un espace de Banach.

Espace de Hilbert

Un espace préhilbertien complet est un espace de Hilbert.

Espace euclidien

Un espace de Hilbert de dimension finie est appelé espace euclidien.

1.1.3. Espaces fonctionnels

Un espace fonctionnel est un espace vectoriel dont les éléments sont des fonctions.

\({\cal C}^p([a;b\))]

\({\cal C}^p([a;b\))] désigne l’espace des fonctions définies sur l’intervalle \([a,b\)], dont toutes les dérivées jusqu’à l’ordre \(p\) existent et sont continues sur \([a,b\)].

Dans la suite, les fonctions seront définies sur un sous-ensemble de \(\RR^n\) (le plus souvent un ouvert noté \(\Omega\)), à valeurs dans \(\RR\) ou \(\RR^p\).

La température \(T(x,y,z,t)\) en tout point d’un objet \(\bar{\Omega}\subset \RR^3\) est une fonction de \( \bar{\Omega} \times \RR \longrightarrow \RR\).

Les normes usuelles les plus simples sur les espaces fonctionnels sont les normes \(\bf L^p\) définies par :

\[\| u \|_{L^p} = \left ( \int_{\Omega } |u|^p \right) ^{1/p} \quad ,\; p\in [1,+ \infty[ , \qquad \hbox{et}\qquad \| u \|_{L^\infty} = {\hbox{Sup}}_{\Omega } |u|\]

Comme on va le voir, ces formes \(L^p\) ne sont pas nécessairement des normes. Et lorsqu’elles le sont, les espaces fonctionnels munis de ces normes ne sont pas nécessairement des espaces de Banach. Par exemple, les formes \(L^\infty\) et \(L^1\) sont bien des normes sur l’espace \({\cal C}^0([a;b\))], et cet espace est complet si on le munit de la norme \(L^\infty\), mais ne l’est pas si on le munit de la norme \(L^1\).

Pour cette raison, on va définir les espaces \({\cal L}^p(\Omega)(p\in [1,+ \infty[\)) par

\[[{\cal L}^p(\Omega) = \left\{ u : \Omega \rightarrow \RR, \hbox{ mesurable, et telle que }\int_\Omega |u|^p<\infty \right\}\]
on rappelle qu’une fonction \(u\) est mesurable ssi \(\{ x / |u(x)|<r \}\) est mesurable \(\forall r>0\).

Sur ces espaces \({\cal L}^p(\Omega)\), les formes \(L^p\) ne sont pas des normes. En effet, \(\| u \|_{L^p} = 0\) implique que \(u\) est nulle presque partout dans \({\cal L}^p(\Omega)\), et non pas \(u=0\). C’est pourquoi on va définir les espaces \(\bf L^p(\Omega)\) :

Égalité presque partout

\(L^p(\Omega)\) est la classe d’équivalence des fonctions de \({\cal L}^p(\Omega)\) pour la relation d’équivalence égalité presque partout. Autrement dit, on confondra deux fonctions dès lors qu’elles sont égales presque partout, c’est à dire qu’elles ne diffèrent que sur un ensemble de mesure nulle.

Les espaces \(L^p(\Omega)\) sont complets

La forme \(L^p\) est une norme sur \(L^p(\Omega)\), et \(L^p(\Omega)\) muni de la norme \(L^p\) est un espace de Banach (c.a.d. est complet).

Un cas particulier très important est \(p=2\). On obtient alors l’espace fonctionnel \(L^2(\Omega)\), c’est à dire l’espace des fonctions de carré sommable sur \(\Omega\) (à la relation d’équivalence égalité presque partout près). A la norme \(L^2\) : \(\| u \|_{L^2} = \left( \int_\Omega u^2 \right)^{1/2} \), on peut associer la forme bilinéaire \((u,v)_{L^2} = \int_\Omega u\, v\). Il s’agit d’un produit scalaire, dont dérive la norme \(L^2\).

D’où le théorème suivant

\(L^2(\Omega)\) est un espace de Hilbert.

1.2. Notion de dérivée généralisée

Nous venons de définir des espaces fonctionnels complets, ce qui sera un bon cadre pour démontrer l’existence et l’unicité de solutions d’équations aux dérivées partielles, comme on le verra plus loin notamment avec le théorème de Lax-Milgram.

Toutefois, on a vu que les éléments de ces espaces \(L^p\) ne sont pas nécessairement des fonctions très régulières.

Dès lors, les dérivées partielles de telles fonctions ne sont pas forcément définies partout.

Pour s’affranchir de ce problème, on va étendre la notion de dérivation.

Le véritable outil à introduire pour cela est la notion de distribution, due à L. Schwartz (1950).

Par manque de temps dans ce cours, on se contentera ici d’en donner une idée très simplifiée, avec la notion de dérivée généralisée.

Cette dernière a des propriétés beaucoup plus limitées que les distributions, mais permet de “sentir" les aspects nécessaires pour mener à la formulation variationnelle.

Dans la suite, \(\Omega\) sera un ouvert (pas nécessairement borné) de \(\RR^n\).

1.2.1. Fonctions tests

Soit \(\varphi : \Omega \rightarrow \RR\).

On appelle support de \(\bf \varphi\) l’adhérence de \(\{ x \in \Omega / \varphi(x) \ne 0 \}\).

Pour \(\Omega = ]-1,1\[\), et \(\varphi\) la fonction constante égale à 1, \(\hbox{Supp}\, \varphi = [-1,1\)].

Espace des fonctions tests

On note \({\cal D}(\Omega)\) l’espace des fonctions de \(\Omega\) vers \(\RR\), de classe \({\cal C}^\infty\), et à support compact inclus dans \(\Omega\).

\({\cal D}(\Omega)\) est parfois appelé espace des fonctions-tests.

L’exemple le plus classique dans le cas 1-D est la fonction

\[ \varphi(x) = \left\{ \begin{array}{ll} { e^{- \frac{1}{1-x^2}} } & \hbox{si } |x|<1\\ 0 & \hbox{si } |x|\ge 1\\ \end{array} \right. \] \(\varphi\) est une fonction de \({\cal D}(]a,b\[)\) pour tous \(a < -1 < 1 < b\).

Cet exemple s’étend aisément au cas multi-dimensionnel (\(n>1\)).

Soit \(a\in\Omega\) et \(r>0\) tel que la boule fermée de centre \(a\) et de rayon \(r\) soit incluse dans \(\Omega\).

On pose alors :

\[ \varphi(x) = \left\{ \begin{array}{ll} { e^{- \frac{1}{r^2-|x-a|^2}} } & \hbox{si } |x-a|<r\\ 0 & \hbox{sinon }\\ \end{array} \right.\]

\(\varphi\) ainsi définie est un élément de \({\cal D}(\Omega)\).

Adhérence de \(\overline{{\cal D}(\Omega)\)

\(\overline{{\cal D}(\Omega) } = L^2(\Omega)\)

1.2.2. Dérivée généralisée

Soit \(u\in {\cal C}^1(\Omega)\) et \(\varphi \in {\cal D}(\Omega)\).

Par intégration par parties (annexe [sec:green]), on a :

\[\int_\Omega \partial_i u\; \varphi = - \int_\Omega u \; \partial_i\varphi + \int_{\partial \Omega} u \; \varphi \; {\bf e}_i.{\bf n}\]

Ce dernier terme (intégrale sur le bord de \(\Omega\)) est nul car \(\varphi\) est à support compact (donc nul sur \(\partial \Omega\)).

Or \(\int_\Omega u \; \partial_i\varphi\) a un sens par exemple dès que \(u\in L^2(\Omega)\).

Donc le terme \(\int_\Omega \partial_i u\; \varphi\) a aussi du sens, sans que \(u\) ne soit nécessairement de classe \({\cal C}^1\).

Ceci permet de définir \(\partial_i u\) même dans ce cas.

cas 1-D \(\quad\) Soit \(I\) un intervalle de \(\RR\), pas forcément borné.

On dit que \(u\in L^2(I)\) admet une dérivée généralisée dans \(L^2(I)\) ssi \(\exists u_1\in L^2(I)\) telle que

\[ \forall \varphi\in {\cal D}(I), \quad \int_I u_1\;\varphi = - \int_I u \varphi' \]

Soit \(I=]a,b[\) un intervalle borné, et \(c\) un point de \(I\). On considère une fonction \(u\) formée de deux branches de classe \({\cal C}^1\), l’une sur \(]a,c[\), l’autre sur \(]c,b[\), et se raccordant de façon continue mais non dérivable en \(c\). Alors \(u\) admet une dérivée généralisée définie par \(u_1(x)=u'(x)\quad \forall x\ne c\). En effet :

\[ \forall \varphi\in {\cal D}(]a,b[)\qquad \int_a^b u \varphi' = \int_a^c + \int_c^b = - \int_a^c u' \varphi - \int_c^b u'\varphi + \underbrace{(u(c-)-u(c+))}_{=0} \, \varphi(c) \]

par intégration par parties. La valeur \(u_1(c)\) n’a pas d’importance: on a de toute façon au final la même fonction de \(L^2(I)\), puisqu’elle est définie comme classe d’équivalence de la relation d’équivalence égalité presque partout.

En itérant, on dit que \(u\) admet une dérivée généralisée d’ordre \(\bf k\) dans \(L^2(I)\), notée \(u_k\), ssi \({\forall \varphi\in {\cal D}(I), \quad \int_I u_k\;\varphi = (- 1)^k \; \int_I u \varphi^{(k)} }\)

Ces définitions s’étendent naturellement pour la définition de dérivées partielles généralisées, dans le cas \(n>1\).

Unicité de la dérivée généralisée

Quand elle existe, la dérivée généralisée est unique.

Quand \(u\) est de classe \({\cal C}^1(\bar{\Omega})\), la dérivée généralisée est égale à la dérivée classique.

1.3. Espaces de Sobolev

1.3.1. Les espaces \(H^m\)

\({ H^1(\Omega) = \left\{ u \in L^2(\Omega)\; / \; \partial_i u \; \in L^2(\Omega), \quad 1 \le i \le n \right\} }\) où \(\partial_i u\) est définie au sens de la dérivée généralisée.

\(H^1(\Omega)\) est appelé espace de Sobolev d’ordre 1.

Pour tout entier \(m\ge 1\), \[ H^m(\Omega) = \left\{ u \in L^2(\Omega) \; / \; \partial^\alpha u \; \in L^2(\Omega) \quad \forall \alpha =(\alpha_1,\ldots,\alpha_n) \in \NN^n\hbox{ tel que}\; |\alpha|= \alpha_1+\cdots+\alpha_n \le m \right\}\]

\(H^m(\Omega)\) est appelé *espace de Sobolev d’ordre \(\bf m\).

Par extension, on voit aussi que \(H^0(\Omega)=L^2(\Omega)\).

Dans le cas de la dimension 1, on écrit plus simplement pour \(I\) ouvert de \(\RR\) :

\[ H^m(I) = \left\{ u \in L^2(I) \; / \; u', \ldots, u^{(m)} \in L^2(I) \right\} \]

\(H^1(\Omega)\) est un espace de Hilbert pour le produit scalaire \[(u,v)1 = \int\Omega u \, v\, + \sum_{i=1}^n \; \int_\Omega \partial_i u \; \partial_i v = (u,v)0 + \sum{i=1}^n (\partial_i u, \partial_i v )_0\]

en notant \((.,.)_0\) le produit scalaire \(L^2\). On notera \(\|.\|_1\) la norme associée à \((.,.)_1\).

On définit de même un produit scalaire et une norme sur \(H^m(\Omega)\) par \(\[(u,v)_m = \sum_{|\alpha| \le m} ( \partial^\alpha u , \partial^\alpha v )_0 \qquad \hbox{ et }\qquad \| u \|_m = (u,u)_m^{1/2}]\)

\(H^m(\Omega)\) muni du produit scalaire \((.,.)_m\) est un espace de Hilbert.[thr:8]

Si \(\Omega\) est un ouvert de \(\RR^n\) de frontière \(\partial\Omega\) “suffisamment régulière" (par exemple \({\cal C}^1\)), on a l’inclusion : \(H^m(\Omega) \subset {\cal C}^k(\bar{\Omega})\) pour \({ k < m-\frac{n}{2} }\)

En particulier, on voit que pour un intervalle \(I\) de \(\RR\), on a \(H^1(I) \subset {\cal C}^0(\bar{I})\), c’est à dire que, en 1-D, toute fonction \(H^1\) est continue.

L’exemple de \(u(x) = x\, \sin\frac{1}{x}\) pour \(x\in\)0,1]] et \(u(0)=0\) montre que la réciproque est fausse.

L’exemple de \(u(x,y) = | \ln (x^2+y^2) |^k\) pour \(0<k<1/2\) montre qu’en dimension supérieure à 1 il existe des fonctions \(H^1\) discontinues.

1.3.2. Trace d’une fonction

Pour pouvoir faire les intégrations par parties qui seront utiles par exemple pour la formulation variationnelle, il faut pouvoir définir le prolongement (la trace) d’une fonction sur le bord de l’ouvert \(\Omega\).

\(n=1\) (cas 1-D)

on considère un intervalle ouvert \(I=]a,b[\) borné. On a vu que \(H^1(I) \subset {\cal C}^0(\bar{I})\). Donc, pour \(u\in H^1(I)\), \(u\) est continue sur \([a,b\)], et \(u(a)\) et \(u(b)\) sont bien définies.

\(n>1\)

on n’a plus \(H^1(\Omega) \subset {\cal C}^0(\bar{\Omega})\). Comment alors définir la trace ? La démarche est la suivante :

  • On définit l’espace \({\cal C}^1(\bar{\Omega}) = \left\{ \varphi : \Omega \rightarrow \RR \;/\; \exists O \hbox{ ouvert contenant } \bar{\Omega},\; \exists \psi \in {\cal C}^1(O),\; \psi_{|\Omega} = \varphi \right\}\) Autrement dit, \({\cal C}^1(\bar{\Omega})\) est l’espace des fonctions \({\cal C}^1\) sur \(\Omega\), prolongeables par continuité sur \(\partial\Omega\) et dont le gradient est lui-aussi prolongeable par continuité. Il n’y a donc pas de problème pour définir la trace de telles fonctions.

  • On montre que, si \(\Omega\) est un ouvert borné de frontière \(\partial\Omega\) “assez régulière", alors \({\cal C}^1(\bar{\Omega})\) est dense dans \(H^1(\Omega)\).

  • L’application linéaire continue, qui à toute fonction \(u\) de \({\cal C}^1(\bar{\Omega})\) associe sa trace sur \(\partial\Omega\), se prolonge alors en une application linéaire continue de \(H^1(\Omega)\) dans \(L^2(\partial\Omega)\), notée \(\gamma_0\), qu’on appelle application trace. On dit que \(\gamma_0(u)\) est la trace de \(u\) sur \(\partial\Omega\).

Pour une fonction \(u\) de \(H^1(\Omega)\) qui soit en même temps continue sur \(\bar{\Omega}\), on a évidemment \(\gamma_0(u) = u_{|\partial\Omega}\). C’est pourquoi on note souvent par abus simplement \(u_{|\partial\Omega}\) plutôt que \(\gamma_0(u)\).

On peut de façon analogue définir \(\gamma_1\), application trace qui permet de prolonger la définition usuelle de la dérivée normale sur \(\partial\Omega\). Pour \(u\in H^2(\Omega)\), on a \(\partial_i u \in H^1(\Omega)\), \(\forall i=1,\ldots,n\), et on peut donc définir \(\gamma_0(\partial_i u)\). La frontière \(\partial\Omega\) étant “assez régulière" (par exemple, idéalement, de classe \({\cal C}^1\)), on peut définir la normale \(n=\left( \begin{array}{l} n_1 \\ \vdots \\ n_n \end{array} \right)\) en tout point de \(\partial\Omega\). On pose alors \({\gamma_1(u) = \sum_{i=1}^n \gamma_0(\partial_i u) n_i}\). Cette application continue \(\gamma_1\) de \(H^2(\Omega)\) dans \(L^2(\partial\Omega)\) permet donc bien de prolonger la définition usuelle de la dérivée normale. Dans le cas où \(u\) est une fonction de \(H^2(\Omega)\) qui soit en même temps dans \({\cal C}^1(\bar{\Omega})\), la dérivée normale au sens usuel de \(u\) existe, et \(\gamma_1(u)\) lui est évidemment égal. C’est pourquoi on note souvent, par abus, \(\partial_n u\) plutôt que \(\gamma_1(u)\).

1.3.3. Espace \(H^1_0(\Omega)\)

Soit \(\Omega\) ouvert de \(\RR^n\). L’espace \(H^1_0(\Omega)\) est défini comme l’adhérence de \({\cal D}(\Omega)\) pour la norme \(\|.\|_1\) de \(H^1(\Omega)\). (on rappelle que \({\cal D}(\Omega)\) est l’espace des fonctions \({\cal C}^\infty\) sur \(\Omega\) à support compact, encore appelé espace des fonctions tests)

Par construction \(H^1_0(\Omega)\) est un espace complet. C’est un espace de Hilbert pour la norme \(\|.\|_1\)

Si \(n=1\) (cas 1-D)}

on considère un intervalle ouvert \(I=\)a,b[] borné. Alors \[ H^1_0(]a,b[) = \left\{ u \in H^1(]a,b[),\; u(a)=u(b)=0 \right\} \]

Si \(n>1\)

Si \(\Omega\) est un ouvert borné de frontière“assez régulière" (par exemple \({\cal C}^1\) par morceaux), alors \(H^1_0(\Omega) = \ker \gamma_0\). \(H^1_0(\Omega)\) est donc le sous-espace des fonctions de \(H^1(\Omega)\) de trace nulle sur la frontière \(\partial\Omega\).

Pour toute fonction \(u\) de \(H^1(\Omega)\), on peut définir : \[ { |u|1 = \left( \sum{i=1}^n \| \partial_i u \|0^2 \right)^{1/2} = \left( \int\Omega \sum_{i=1}^n \left( \partial_i u \right)^2 dx \right)^{1/2} } \]

Inégalité de Poincaré

Si \(\Omega\) est borné dans au moins une direction, alors il existe une constante \(C(\Omega)\) telle que \[ \forall u \in H^1_0(\Omega), \; \|u\|_0 \le C(\Omega)\; |u|_1. \]

On en déduit que \(|.|_1\) est une norme sur \(H^1_0(\Omega)\), équivalente à la norme \(\|.\|_1\).

Le résultat précédent s’étend au cas où l’on a une condition de Dirichlet nulle seulement sur une partie de \(\partial\Omega\), si \(\Omega\) est connexe.

On suppose que \(\Omega\) est un ouvert borné connexe, de frontière \({\cal C}^1\) par morceaux.

Soit \(V=\left\{ v\in H^1(\Omega),\, v=0 \hbox{ sur }\Gamma_0 \right\}\) où \(\Gamma_0\) est une partie de \(\partial\Omega\) de mesure non-nulle.

Alors il existe une constante \(C(\Omega)\) telle que \(\forall u \in V, \; \|u\|_{0,V} \le C(\Omega)\; |u|_{1,V}\), où \(\|.\|_{0,V}\) et \(|.|_{1,V}\) désignent les norme et semi-norme induites sur \(V\).

On en déduit que \(|.|_{1,V}\) est une norme sur \(V\), équivalente à la norme \(\|.\|_{1,V}\).

1.4. Exercices

  1. Montrer que les fonctions définies par ([eq:fonction-test1]) et ([eq:fonction-test2]) sont bien \({\cal C}^\infty\) à support compact.

  2. Montrer que \({\cal C}^0([a,b\))] est un espace complet pour la norme \(L^\infty\).

  3. Montrer que ce n’est pas le cas pour la norme \(L^1\) (exhiber une suite de Cauchy non convergente dans \({\cal C}^0([a,b\))]).

  4. Démontrer que, lorsqu’elle existe, la dérivée généralisée est unique.

  5. Démontrer que, pour une fonction de classe \({\cal C}^1\), la dérivée généralisée est égale à la dérivée classique.

  6. Soit une fonction de \([a,b\)] vers \(\RR\), formée de deux branches de classe \({\cal C}^1\) sur \([a,c[\) et \(]c,b\)], et discontinue en \(c\). Montrer qu’elle n’admet pas de dérivée généralisée. (il faudrait alors avoir recours à la notion de distribution pour dériver cette fonction).

  7. Montrer que \(|.|_1\) est une norme sur \(H^1_0(\Omega)\), équivalente à la norme \(\|.\|_1\)

Finite Element Method

1. Introduction à la méthode des éléments finis

1.1. Formulation variationnelle

1.1.1. Exemple 1-D

Soit à résoudre le problème (

({\cal P}) \left\{ \begin{array}{l} -u"(x)+c(x)u(x) = f(x),\quad a<x<b,\\ u(a)=u(b)=0 \end{array}\right.

où \(f\) et \(c\) sont des fonctions données continues sur \([a,b\)].

On supposera de plus que la fonction \(c\) est strictement positive sur \([a,b\)].

Un tel problème est appelé problème aux limites.

Une solution classique (ou solution forte) de \(({\cal P})\) est une fonction de \({\cal C}^2([a,b\))] telle que \(u(a)=u(b)=0\) et \(\forall x \in \)a,b[, \; -u"(x)+c(x)u(x) = f(x)].

En faisant le produit scalaire \(L^2(]a,b[)\) de l’équation différentielle avec une fonction-test \(v \in {\cal D}(]a,b[)\) (c’est à dire en intégrant sur \([a,b\)]), on a :

  • \int_a^b u"(x)v(x)\; dx + \int_a^b c(x)u(x)v(x)\; dx = \int_a^b f(x)v(x)\; dx

soit, en intégrant par parties le premier terme :

Formulation faible

\int_a^b u'(x)v'(x)\; dx + \int_a^b c(x)u(x)v(x)\; dx = \int_a^b f(x)v(x)\; dx

car \(v(a)=v(b)=0\) puisque \(v\in{\cal D}(]a,b[)\).

Chaque terme de cette équation a en fait un sens dès lors que \(v\in H^1_0(]a,b[)\).

De plus, \({\cal D}(]a,b[)\) étant dense dans \(H^1_0(]a,b[)\) (§[sec:H10]), cette équation est vérifiée pour tout \(v\in H^1_0(]a,b[)\).

On peut donc définir le nouveau problème :

({\cal Q})\quad \left\{ \begin{array}{l} \hbox{Trouver }u\in H^1_0(]a,b[) \hbox{ tel que }\\ \ds{ \int_a^b u'(x)v'(x)\; dx + \int_a^b c(x)u(x)v(x)\; dx = \int_a^b f(x)v(x)\; dx \quad\forall v\in H^1_0(]a,b[)} \end{array} \right.

Ce problème est la formulation variationnelle (ou formulation faible) du problème \(({\cal P})\). Toute solution de \(({\cal Q})\) est appelée solution faible.

Il est immédiat que toute solution forte de \(({\cal P})\) est aussi une solution faible.

1.1.2. Exemple 2-D

Soit \(\Omega\) ouvert borné de \(\RR^n\).

On veut résoudre le problème

({\cal P})\qquad \left\{ \begin{array}{ll} -\Delta u+u = f & \hbox{ dans }\Omega\\ \partial_n u=0& \hbox{ sur }\partial\Omega\\ \end{array}\right.

[Une solution classique de ce problème est une fonction de \({\cal C}^2(\bar{\Omega})\) vérifiant ([eq:modele-2D]) en tout point de \(\Omega\).

Au passage, on voit que ceci impose que \(f\) soit \({\cal C}^0(\bar{\Omega})\).

Toute solution classique vérifie donc :

\forall v\in {\cal C}^1(\bar{\Omega}) \quad \int_\Omega -\Delta u\; v + \int_\Omega uv = \int_\Omega fv

soit par intégration par parties :

\ds{ \int_\Omega \nabla u \cdot \nabla v + \int_\Omega uv = \int_\Omega fv }

puisque \(\partial_n u = 0 \) sur \(\partial\Omega\).

Or, \(\overline{{\cal C}^1(\bar{\Omega})} = H^1(\Omega)\).

On peut donc définir le nouveau problème :

Formulation variationnelle en 2D

({\cal Q})\quad \left\{ \begin{array}{l} \hbox{Trouver }u\in H^1(\Omega) \hbox{ tel que }\\ \qquad \ds{ \int_\Omega \nabla u \cdot \nabla v + \int_\Omega u v = \int_\Omega f v \quad\forall v\in H^1(\Omega)} \end{array} \right.

C’est la formulation variationnelle de \(({\cal P})\).

On voit aussi que ce problème est défini dès lors que \(f\in L^2(\Omega)\).

1.2. Formulation générale

Les exemples précédents montre que, d’une façon générale, la formulation variationnelle sera obtenue en faisant le produit scalaire \(L^2(\Omega)\) de l’équation avec une fonction \(v\) appartenant à un espace à préciser (c’est à dire en multipliant par \(v\) et en intégrant sur \(\Omega\)), et en intégrant par parties les termes d’ordre les plus élevés en tenant compte des conditions aux limites du problème.

On arrive alors à une formulation du type :

\hbox{Trouver }u\in V \hbox{ tel que } a(u, v) =l(v) \quad\forall v\in V

où \(a(.,.)\) est une forme sur \(V\times V\) (bilinéaire si l’EDP de départ est linéaire) et \(l(.)\) est une forme sur \(V\) (linéaire si les conditions aux limites de l’EDP de départ le sont).

1.3. Théorème de Lax-Milgram

1.3.1. Définitions et théorèmes

On va introduire ici un outil important pour assurer l’existence et l’unicité de solutions à la formulation variationnelle de problèmes aux limites de type elliptique.

Soit \(V\) un espace de Hilbert.

Continuité d’une forme linéaire

Une forme linéaire \(l(u)\) sur \(V\) est continue ssi il existe une constante \(K\) telle que \[ |l(u)| \le K\, \|u\| \quad \forall u \in V \]

Continuité d’une forme bilinéaire

Une forme bilinéaire \(a(u,v)\) sur \(V\times V\) est continue ssi il existe une constante \(M\) telle que

\[ |a(u,v)| \le M\, \|u\| \|v\| \quad \forall (u,v) \in V\times V. \]

Coercivité d’une forme bilinéaire

Une forme bilinéaire \(a(u,v)\) sur \(V\times V\) est coercive (ou V-elliptique) ssi il existe une constante \(\alpha >0\) telle que \[ a(u,u) \ge \alpha \, \|u\|^2, \; \forall u \in V. \]

Lax-Milgram

Soit \(V\) un espace de Hilbert. Soit \(a\) une forme bilinéaire continue coercive sur \(V\). Soit \(l\) une forme linéaire continue sur \(V\). Alors il existe un unique \(u\in V\) tel que \[ a(u,v)=l(v),\; \forall v\in V. \]

De plus, l’application linéaire \(l \rightarrow u\) est continue.

Lax-Milgram

La démonstration générale de ce théorème peut être trouvée par exemple dans (raviart-thomas-1983).

On prend les mêmes hypothèses que pour le théorème de Lax-Milgram, et on suppose de plus que \(a\) est symétrique, c’est à dire que \(a(u,v)=a(v,u)\quad\forall u,v\).

On définit alors la fonctionnelle \(J(v)=\frac{1}{2}\, a(v,v)-l(v)\), et on considère le problème de minimisation :

\[ \hbox{Trouver } u\in V \hbox{ tel que } J(u) = \min_{v\in V} J(v) \]

Alors ce problème admet une solution unique, qui est également la solution du problème variationnel précédent.

La démonstration de ce théorème vient du fait que \(J\) est une fonctionnelle quadratique, et que l’on a \(\nabla J[u\)(v) = a(u,v) - l(v)].

C’est de cette propriété que vient l’utilisation du terme “variationnel", puisqu’elle montre le lien avec le “calcul des variations".

1.3.2. Retour à l’exemple 1-D

En reprenant l’exemple 1-D précédent, on peut poser :

a(u,v) = \int_a^b u'(x)v'(x)\; dx + \int_a^b c(x)u(x)v(x)\; dx

et

l(v) = \int_a^b f(x)v(x)\; dx

\(a\) ainsi définie est une forme bilinéaire symétrique continue coercive sur \(H^1_0(a,b) \times H^1_0(a,b)\), et \(l\) est une forme linéaire continue sur \(H^1_0(a,b)\). Donc le problème ([eq:FV]) admet une solution unique d’après le théorème de Lax-Milgram. Cherchons maintenant à interpréter cette solution \(u\) de ([eq:FV]). Prenons \(v=\varphi \in {\cal D}(]a,b[)\). Alors

\int_a^b u'(x)\varphi'(x)\; dx + \int_a^b c(x)u(x)\varphi(x)\; dx = \int_a^b f(x)\varphi(x)\; dx

soit, en intégrant par parties :

  • \int_a^b u"(x)\varphi(x)\; dx + \int_a^b c(x)u(x)\varphi(x)\; dx = \int_a^b f(x)\varphi(x)\; dx

c’est à dire

(-u"+cu,\varphi)_0 = (f,\varphi)_0\; \forall \varphi \in {\cal D}(]a,b[).

\({\cal D}(]a,b[)\) étant dense dans \(L^2(]a,b[)\), on a donc : \(-u"+cu=f\) dans \(L^2(]a,b[)\).

\(u\) étant dans \(L^2(]a,b[)\), et \(f\) et \(c\) étant dans \({\cal C}^0([a,b])\), donc également dans \(L^2(]a,b[)\), on en déduit que \(u"=cu-f\) est aussi dans \(L^2(]a,b[)\).

Puisque \(u\) est dans \(H^1_0(]a,b[)\) et que \(u"\) est dans \(L^2(]a,b[)\), on en déduit que \(u\) est dans \(H^2(]a,b[)\).

Donc \(u\) est dans \({\cal C}^1([a,b])\) (§[sec:sobolev]).

De ce fait, \(cu-f\), c’est à dire \(u"\), est dans \({\cal C}^0([a,b])\).

Donc \(u'\) est dans\({\cal C}^1([a,b])\), donc \(u\) est dans \({\cal C}^2([a,b])\).

La solution faible \(u\) est donc aussi solution forte du problème de départ.

En résumé :

  • On est parti d’un problème \(({\cal P})\) et on a introduit sa formulation variationnelle \(({\cal Q})\).

  • On a montré l’existence et l’unicité d’une solution faible (en utilisant le théorème de Lax-Milgram). Toute solution forte étant aussi solution faible, ceci prouve qu’il y a au plus une solution forte pour \(({\cal P})\).

  • On a prouvé que cette solution faible est bien une solution forte. Le problème de départ \(({\cal P})\) a donc une solution unique.

L’intérêt de cette démarche est . la formulation variationnelle se prête bien à l’étude de l’existence et de l’unicité de solutions, . on travaille dans des espaces de Hilbert, ce qui va permettre de faire de l’approximation interne.

1.3.3. Équations elliptiques d’ordre 2

Soit \(\Omega\) un ouvert borné de \(\RR^n\), de frontière \(\partial\Omega\) assez régulière. Soient des fonctions \(\alpha_{ij}\) (\(1\le i,j \le n\)) dans \({\cal C}^1(\bar{\Omega})\) et \(\beta\) dans \({\cal C}^0(\bar{\Omega})\).

On considère le problème :

({\cal P})\qquad \left\{ \begin{array}{rl} {\ds -\sum_{i,j=1}^n \partial_i (\alpha_{ij} \, \partial_j\, u) + \beta\, u = f } & \hbox{ dans }\Omega \\ u= 0 & \hbox{ sur }\Gamma_0 \\ {\ds \sum_{i,j=1}^n \alpha_{ij} \, \partial_j u\; n_i = g } & \hbox{ sur }\Gamma_1 % \end{array}\right.

où \(\Gamma_0\) et \(\Gamma_1\) forment une partition de \(\partial\Omega\) (\(\Gamma_0 \cap\Gamma_1 = \emptyset\) et \(\Gamma_0 \cup\Gamma_1 = \partial\Omega\)).

Une solution classique de \(({\cal P})\), sous l’hypothèse que \(f\in{\cal C}^0(\bar{\Omega})\) et \(g\in{\cal C}^0(\Gamma_1)\), sera une fonction de \({\cal C}^2(\bar{\Omega})\) vérifiant l’équation en chaque point de \(\Omega\).La formulation variationnelle de \(({\cal P})\) est obtenue par intégration par parties.

Elle s’écrit :

({\cal Q})\quad \left\{ \begin{array}{l} \hbox{Trouver }u\in V \hbox{ tel que }\\ \qquad \ds{ \int_\Omega \left( \sum_{i,j=1}^n \alpha_{ij} \, \partial_j u\; \, \partial_i v + \beta\, u v \right) = \int_\Omega f v + \int_{\Gamma_1} gv \qquad\forall v\in V} \end{array} \right.

avec \(\ds{ V = \left\{ v \in H^1(\Omega) \; , \; v=0 \hbox{ sur }\Gamma_0 \right\} }\). Cette formulation est en fait définie dès lors que \(\beta\) et les \(\alpha_{ij}\) sont dans \(L^\infty(\Omega)\), \(f\) dans \(L^2(\Omega)\) et \(g\) dans \(L^2(\Gamma_1)\). Posons

\ds{a(u,v) = \int_\Omega \left( \sum_{i,j=1}^n \alpha_{ij} \, \partial_j u\; \partial_i v + \beta\, u v \right) }, \quad \ds{l(v) = \int_\Omega f v + \int_{\Gamma_1} gv. }

Il est immédiat que \(a\) est une forme bilinéaire continue et \(l\) une forme linéaire continue sur \(V\).

Si l’EDP de départ [eq:edp-elliptique] vérifie les deux hypothèses d’ellipticité :

  • il existe \(\alpha >0\) tel que \(\forall \xi=(\xi_1, \ldots , \xi_n)\in\RR^n\), \( {\sum_{i,j=1}^n \alpha_{ij}(x) \, \xi_i \, \xi_j \ge \alpha \, \| \xi \|^2 }\) presque pour tout \(x\in\Omega\)

  • il existe \(\beta_0\) tel que \(\beta(x) \ge \beta_0\) presque pour tout \(x\in\Omega\)

alors \(a\) est coercive :

  • sur \(H^1_0(\Omega)\) dès que \(\ds{\alpha_0 > \frac{-\alpha}{C(\Omega)^2}}\) (et donc en particulier si \(\beta\ge 0\)) où \(C(\Omega)\) est la constante de l’inégalité de Poincaré, voir le théorème [thr:11].

  • sur \(H^1(\Omega)\) si \(\beta > 0\)

Par application du théorème de Lax-Milgram, on a donc existence et unicité d’une solution à la formulation variationnelle \(({\cal Q})\) :

  • si \(\Gamma_0 = \partial\Omega\) (c’est à dire \(\Gamma_1=\emptyset\)) et si \(\ds{\beta > \frac{-\alpha}{C(\Omega)^2}}\)

  • si \(\Gamma_1\ne \emptyset\) et si \(\beta > 0\)

1.4. Approximation interne

1.4.1. Principe général

Soit \(\Omega\) un domaine ouvert de \(\RR^n\) (\(n=1,2\) ou 3 en pratique), de frontière \(\partial\Omega\), et sur lequel on cherche à résoudre une équation aux dérivées partielles, munie de conditions aux limites.

En écrivant la formulation variationnelle, on obtient un problème de la forme

({\cal Q})\qquad \hbox{Trouver } u\in V \hbox{ tel que } a(u,v)=l(v), \quad\forall v\in V

où \(V\) est un espace de Hilbert. Sous réserve que l’équation de départ ait de bonnes propriétés, c’est à dire par exemple qu’on soit dans les hypothèses du théorème de Lax-Milgram, \(({\cal Q})\) admet une solution unique \(u\).

Pour obtenir une approximation numérique de \(u\), on va maintenant remplacer l’espace \(V\) qui est en général de dimension infinie par un sous-espace \(V_h\) de dimension finie, et on va chercher à résoudre le problème approché

\label{eq:6} ({\cal Q}_h)\qquad \hbox{Trouver } u_h\in V_h \hbox{ tel que } a(u_h,v_h)=l(v_h), \quad\forall v_h\in V_h

\(V_h\) étant de dimension finie, c’est un fermé de \(V\). \(V\) étant un espace de Hilbert, \(V_h\) l’est donc aussi. D’où l’existence et l’unicité de \(u_h\), à nouveau par exemple d’après le théorème de Lax-Milgram. L’espace \(V_h\) sera en pratique construit à partir d’un maillage du domaine \(\Omega\), l’indice \(h\) désignant la ``taille typique'' des mailles.

Lorsque l’on construit des maillages de plus en plus fins, la suite de sous-espaces \((V_h)_h\) formera une approximation interne de \(V\), c’est à dire que, pour tout élément \(\varphi\) de \(V\), il existe une suite de \(\varphi_h\in V_h\) telle que \(\|\varphi-\varphi_h\|\longrightarrow 0\) quand \(h\longrightarrow 0\).

Cette méthode d’approximation interne est également appelée méthode de Galerkin.

1.4.2. Interprétation de \(u_h\)

On a \(a(u,v)=l(v), \forall v\in V\), donc en particulier \(a(u,v_h)=l(v_h), \forall v_h\in V_h\), car \(V_h\subset V\).

Par ailleurs, \(a(u_h,v_h)=l(v_h), \forall v_h\in V_h\).

Par différence, on en déduit que \(a(u-u_h,v_h)=0,\quad \forall v_h\in V_h \label{eq:ortho}\)

Dans le cas où \(a(.,.)\) est symétrique, il s’agit d’un produit scalaire sur \(V\). \(u_h\) peut alors être interprétée comme la projection orthogonale de \(u\) sur \(V_h\) au sens de \(a(.,.)\).

1.4.3. Estimation d’erreur

On a :

\begin{array}{ll} a(u-u_h,u-u_h) & = a(u-u_h,u-v_h+v_h-u_h) \quad\forall v_h\in V_h\\ & =a(u-u_h,u-v_h) + a(u-u_h,v_h-u_h) \end{array}

Or \(v_h-u_h \in V_h\). Donc \(a(u-u_h,v_h-u_h)=0\) d’après ([eq:ortho]).

On a donc :

a(u-u_h,u-u_h) = a(u-u_h,u-v_h) \quad\forall v_h\in V_h

\(a\) étant coercive, il existe \(\alpha > 0\) tel que \(a(u-u_h,u-u_h) \ge \alpha \|u-u_h\|^2\), où \(\|.\|\) est une norme sur \(V\).

Par ailleurs, \(a\) étant continue, il existe \(M > 0\) tel que \(a(u-u_h,u-v_h)\le M \|u-u_h\| \, \|u-v_h\|\).

En réinjectant ces deux inégalités de part et d’autre de [eq:estim1] et en simplifiant par \(\|u-u_h\|\) on obtient

\|u-u_h\| \le \frac{M}{\alpha}\; \|u-v_h\| \quad \forall v_h\in V_h

c’est à dire

\|u-u_h\| \le \frac{M}{\alpha}\; d(u,V_h)

où \(d\) est la distance induite par \(\|.\|\).

Cette majoration est appelée lemme de Céa. Elle ramène l’étude de l’erreur d’approximation \(u-u_h\) à l’étude de l’erreur d’interpolation \(d(u,V_h)\).

1.5. Principe général de la méthode des éléments finis

La démarche générale de la méthode des éléments finis est la suivante.

On a une EDP à résoudre sur un domaine \(\Omega\).

On écrit la formulation variationnelle de cette EDP, et on se ramène donc à un problème du type

({\cal Q})\qquad \hbox{Trouver } u\in V \hbox{ tel que } a(u,v)=l(v), \quad\forall v\in V

On va chercher une approximation de \(u\) par approximation interne.

Pour cela, on définit un maillage du domaine \(\Omega\), grâce auquel on va définir un espace d’approximation \(V_h\), s.e.v. de \(V\) de dimension finie \(N_h\) (par exemple \(V_h\) sera l’ensemble des fonctions continues sur \(\Omega\) et affines sur chaque maille).

Le problème approché est alors

({\cal Q}_h)\qquad \hbox{Trouver } u_h\in V_h \hbox{ tel que } a(u_h,v_h)=l(v_h), \quad\forall v_h\in V_h

Soit \((\varphi_1,\ldots,\varphi_{N_h})\) une base de \(V_h\).

En décomposant \(u_h\) sur cette base sous la forme

u_h = \sum_{i=1}^{N_h} \mu_i \; \varphi_i

le problème \(({\cal Q}_h)\) devient

\hbox{Trouver } \mu_1,\ldots,\mu_{N_h} \hbox{ tels que } \sum_{i=1}^{N_h} \mu_i \; a(\varphi_i,v_h)=l(v_h), \quad\forall v_h\in V_h

ou encore par linéarité de \(a\) et \(l\) :

\hbox{Trouver } \mu_1,\ldots,\mu_{N_h} \hbox{ tels que } \sum_{i=1}^{N_h} \mu_i \; a(\varphi_i,\varphi_j)=l(\varphi_j), \quad\forall j=1,\ldots,N_h

c’est à dire résoudre le système linéaire

\left( \begin{array}{ccc} a(\varphi_1,\varphi_1) & \cdots & a(\varphi_{N_h},\varphi_1)\\ \vdots & & \vdots\\ a(\varphi_1,\varphi_{N_h}) & \cdots & a(\varphi_{N_h},\varphi_{N_h})\\ \end{array}\right) \left( \begin{array} \mu_1\\ \vdots\\ \mu_{N_h}\\ \end{array}\right) = \left( \begin{array} l(\varphi_1)\\ \vdots\\ l(\varphi_{N_h})\\ \end{array}\right)

soit

A\mu = b

La matrice \(A\) est a priori pleine.

Toutefois, pour limiter le volume de calculs, on va définir des fonctions de base \(\varphi_i\) dont le support sera petit, c’est à dire que chaque fonction \(\varphi_i\) sera nulle partout sauf sur quelques mailles.

Ainsi les termes \(a(\varphi_i,\varphi_j)\) seront le plus souvent nuls, car correspondant à des fonctions \(\varphi_i\) et \(\varphi_j\) de supports disjoints.

La matrice \(A\) sera donc une matrice creuse, et on ordonnera les \(\varphi_i\) de telle sorte que \(A\) soit à structure bande, avec une largeur de bande la plus faible possible.

A ce niveau, les difficultés majeures en pratique sont de trouver les \(\varphi_i\) et de les manipuler pour les calculs d’intégrales nécessaires à la construction de \(A\).

Sans rentrer pour le moment dans les détails, on peut toutefois indiquer que la plupart de ces difficultés seront levées grâce à trois idées principales :

Le principe d’unisolvance

On s’attachera à trouver des degrés de liberté (ou ddl) tels que la donnée de ces ddl détermine de façon univoque toute fonction de \(V_h\). Il pourra s’agir par exemple des valeurs de la fonction en quelques points. Déterminer une fonction reviendra alors à déterminer ses valeurs sur ces ddl.

Définition des \(\varphi_i\)

On définira les fonctions de base par \(\varphi_i=1\) sur le \(i^{\hbox{\tiny ème}}\) ddl, et \(\varphi_i=0\) sur les autres ddl. La manipulation des \(\varphi_i\) sera alors très simplifiée, et les \(\varphi_i\) auront par ailleurs un support réduit à quelques mailles.

La notion de famille affine d’éléments

Le maillage sera tel que toutes les mailles soient identiques à une transformation affine près. De ce fait, tous les calculs d’intégrales pourront se ramener à des calculs sur une seule maille de référence, par un simple changement de variable.

1.6. Retour à l’exemple 1-D

On reprend le problème 1-D [eq:modele-1D].

On a écrit sa formulation variationnelle [sec:modele-1D> et montré <<sec:modele-1D2] qu’elle admet une solution unique.

On s’intéresse à présent à la construction de l’espace d’approximation \(V_h\).

1.6.1. Construction du maillage

La première étape consiste à construire un maillage de \(\Omega = ]a,b[\) en définissant une subdivision (pas nécessairement régulière) \(a=x_0 < x_1 < \ldots < x_N < x_{N+1}=b\).

Le maillage est donc une collection indexée de (\(=N\)) intervalles \{I_i=[x_{i,1},x_{i,2}]\}_{i=1,...\Nma} et on a

[a,b]=\cup_{i=1}^\Nma [x_{1,i},x_{2,i}] \quad \mbox{et} \quad ]x_{1,i},x_{2,i}[ \cap ]x_{1,j},x_{2,j}[ = \emptyset \quad \mbox{ pour } i\neq j

Les intervalles \(I_i\) sont appelées de mailles ou éléments ou cellules du maillage, on a noté \(\Nma\) le nombre de maillage

Les points \(x_i\) sont appelés les sommets du maillage, on note \(\Nso=N+1\) le nombre de sommets.

On note \(h_i = x_{i+1}-x_i\) et \(h = \max_{1\leq i \leq \Nma} h_i\).

Le maillage est dit uniforme si \(h_i=h\) pour tout \(i=\{1,...,\Nma\}\). Enfin on note \(\calTh=\{I_i\}_{i=\{1,...,\Nma\}}\), \(h\) représentant la finesse globale du maillage.

En 1D on a \(\Nso = \Nma+1\), en dimension supérieure des relations existent entre le nombre de sommets et de mailles en fonction du type de maille, ce sont les relations d’Euler.

1.6.2. Construction de l’espace d’approximation

L’étape suivante est de choisir les fonctions de forme ou fonctions de base sur chaque maillage.

On choisit les fonctions de \(V_h\) telle que leur restriction sur chaque maillage soit un espace polynomial.

Espaces \(\Pk\)

Soit un entier \(k \leq 1\). En dimension 1, on appelle l’espace vectoriel des polynômes à coefficients réels de degré inférieur ou égal à \(k\).

On pose alors

W_h = \{w_h \in L^2(\Omega); \forall i \in \{ 1,…​,\Nma\}, {w_h}_{|I_i} \in \Pk\}

\(W_h\) est un espace de dimension finie égale à \((k+1)*\Nma\) mais il n’est pas inclus dans \(H^1_0(\Omega)\) et ne peut donc pas être utilisé pour l’approximation du problème ([eq:FV]). En effet les fonctions de \(w_h \in W_h\) peuvent être discontinues aux interfaces entre les maillages et un résultat d’analyse fonctionnelle montre que dans ces conditions \(w_h \ni H^1(\Omega)\). Par ailleurs les fonctions de \(W_h\) ne sont pas nécessairement nulles au bord de \(\Omega\).

On pose donc

\label{eq:3} V_h = W_h \cap H^1_0(\Omega).

en d’autres termes, en dimension, on a

\label{} V_h = \left\{ v_h \in {\cal C}^0 (a,b) \; ; \; {v_h}_{|I_i} \in \Pk \hbox{ et } v_h(a)=v_h(b)=0 \right\}

Le problème approché sur \(V_h\) est :

({\cal Q}_h)\qquad \hbox{Trouver } u_h\in V_h \hbox{ tel que } a(u_h,v_h)=l(v_h), \quad\forall v_h\in V_h]

On s’intéresse à présent à des exemples concrets d’espaces d’approximations dans les deux sections suivantes Element fini de Lagrange et Element fini de Lagrange.

1.6.3. Element fini de Lagrange

On introduit les espaces vectoriels suivants:

\Pch{1} = \{ v_h \in C^0(\Omega);\; \forall i \in \{ 1,…​,\Nma\} {v_h}_{|I_i} \in \Pk[1] \}

et

\Pcho{1} = \{ v_h \in \Pch{1};\; v_h(a)=v_h(b)=0 \}

Les éléments de ces espaces sont des fonctions continues et affines par morceaux. Ils sont dérivables par morceaux sur chaque maille et ils sont continus aux interfaces entre les mailles.

On a le résultat d’analyse fonctionnelle suivant:

\(\Pch{1} \subset H^1(\Omega)\) et \(\Pcho{1} \subset H^1_0(\Omega)\).

On introduit la famille de fonctions \(\{\varphi_1,...,\varphi_\Nso\}\) que l’on définit sur chaque maille de la manière suivante, pour tout \(i \in \{2,...,\Nso-1\}\),

\varphi_i(x) = \left\{ \begin{split} \ds{\frac{1}{h_{i-1}} (x-x_{i-1})} & \mbox{ si } x \in I_{i-1}\\ \ds{\frac{1}{h_{i}} (x_{i+1}-x)} & \mbox{ si } x \in I_{i}\\ 0 & \mbox{ sinon}, \end{split} \right.

et

\begin{split}
\varphi_1(x) &= \left\{
  \begin{split}
    \ds{\frac{1}{h_{1}} (x_2-x)} & \mbox{ si } x \in I_{1}\\
    0 & \mbox{ sinon},
  \end{split}
\right.\\
\varphi_\Nso(x) &= \left\{
  \begin{split}
    \ds{\frac{1}{h_{\Nso-1}} (x-x_{\Nso-1})} & \mbox{ si } x \in I_{\Nso-1}\\
    0 & \mbox{ sinon},
  \end{split}
\right.
\end{split}
Les fonctions \((\varphi_i)_{i=1,...,\Nso}\) sont dans \(\Pch{1}\) et \((\varphi_i)_{i=2,...,\Nso-1}\) sont dans \(\Pcho{1}\).
Les fonctions \((\varphi_i)_{i=1,...,\Nso}\) satisfont les relations
\varphi_i(x_j) = \delta_{ij},\quad i,j \in \{1,...,\Nso\},

où \(\delta_{ij}\) désigne le symbole de Kronecker tel que \(\delta_{ij} = 1\) si \(i=j\) et \(\delta_{ij}=0\) si \(i \neq j\).

Les fonctions \(\varphi_i\) sont appelées fonctions chapeau du fait de leur graphe, voir figure [fig:chapeau].

chapeau

  1. La famille \(\{\varphi_1,...,\varphi_\Nso\}\) est une base de \(\Pch{1}\).

  2. La famille \(\{\varphi_2,...,\varphi_{\Nso-1}\}\) est une base de \(\Pcho{1}\).

\(\dim \Pch{1} = \Nso = \Nma+1\) et \(\dim \Pcho{1} = \Nso-2 = \Nma-1\).

On introduit l’opérateur d’interpolation suivant:

\Ich{1} : \Ck{0}(\bar{\Omega}) \ni v \mapsto \sum_{i=1}^\Nso v(x_i)
\varphi_i \in \Pch{1}.

Pour toute fonction \(v \in \Ck{0}(\bar{\Omega})\), \(\Ich{1}{v}\) est la seule fonction continue affine par morceaux prenant les mêmes valeurs que \(v\) aux sommets \(x_i, i=1,...,\Nso\).

\(\Ich{1}{v}\) est appelée l’interpolé de Lagrange de \(v\) de degré \(1\).

En dimension 1, les fonctions de \(H^1(\Omega)\) sont continues, on peut donc voir comme un opérateur de \(H^1(\Omega)\) dans \(H^1(\Omega)\). On montre qu’il est continu et que sa norme \(\|\Ich{1}\|_{\mathcal{L}(H^1(\Omega),H^1(\Omega))}\) est uniformément bornée en \(h\), c’est à dire qu’il existe une constante \(c\), indépendante de \(h\), telle que pour tout \(v \in H^1(\Omega)\)
\|\Ich{1} v \|_{1,\Omega} \leq c \|v\|_{1,\Omega}]

1.6.4. Estimation de l’erreur d’interpolation

Pour tout \(h\), et tout \(v \in H^2(\Omega)\), on a

\[ \|v - \Ich{1} v\|{0,\Omega} \leq h^2 |v|{2,\Omega}\quad \mbox{ et }\quad |v - \Ich{1} v|{1,\Omega} \leq h |v|{2,\Omega} \]

On dit que l’erreur d’interpolation est d’ordre 2 en norme \(L^2\) et d’ordre 1 en semi-norme \(H^1\) et donc en norme \(H^1\).

1.6.5. Element fini de Lagrange

On introduit les espaces vectoriels suivants: stem:[\label{eq:9}

\Pch{k} = \{ v_h \in C^0(\Omega);\; \forall i \in \{ 1,...,\Nma\}, {v_h}_{|I_i} \in \Pk\}

et

\Pcho{k} = \{ v_h \in \Pch;\; v_h(a)=v_h(b)=0 \}

1.6.6. Operateur d’interpolation

On introduit l’opérateur d’interpolation suivant: \(\label{eq:24} \Ich : \Ck{0}(\bar{\Omega}) \ni v \mapsto \sum_{i=1}^\Nno v(x_i) \varphi_i \in \Pch.\) Pour toute fonction \(v \in \Ck{0}(\bar{\Omega})\), \(\Ich v\) est la seule fonction continue polynomial de degré \(k\) par morceaux prenant les mêmes valeurs que \(v\) aux sommets \(x_i, i=1,...,\Nso\). \(\Ich v\) est appelée l’interpolé de Lagrange de \(v\) de degré \(k\).

En dimension 1, les fonctions de \(H^1(\Omega)\) sont continues, on peut donc voir comme un opérateur de \(H^1(\Omega)\) dans \(H^1(\Omega)\). On montre qu’il est continu et que sa norme \(\|\Ich\|_{\mathcal{L}(H^1(\Omega),H^1(\Omega))}\) est uniformément bornée en \(h\), c’est à dire qu’il existe une constante \(c\), indépendante de \(h\) mais dépendante de \(k\), telle que pour tout

v \in H^1(\Omega) \|\Ich{k}{v} \|{1,\Omega} \leq c \|v\|{1,\Omega}

Le résultat suivant permet d’estimer la précision de l’opérateur d’interpolation,

Il existe une constante \(c\), indépendante de \(h\) mais dépendante de \(k\), telle que pour tout \(h\) et pour tout \(v \in H^{k+1}(\Omega)\), on a

\[ \|v - \Ich{k}{v}\|{0,\Omega} + h |v - \Ich{k}{v}|{1,\Omega} \leq c\; h^{k+1}\; |v|_{k+1,\Omega} \] et

\[ \sum_{m=2}^{k+1} h^m \left( \sum_{i=0}^N |v - \Ich{k}{v}|2_{m,I_i}\right){1/2} \leq c\; h^{k+1}\; |v|_{k+1,\Omega}] \]

L’estimation [eq:14] montre que l’erreur d’interpolation est d’ordre \(k+1\) en norme \(\|\cdot\|_{0,\Omega}\) et qu’elle est d’ordre \(k\) en norme \(|\cdot|_{1,\Omega}\). Elle est donc d’ordre \(k\) en norme \(\|\cdot\|_{1,\Omega}\).

1.6.7. Analyse de convergence

Nous nous intéressons à présent à l’analyse de la convergence de \(u_h\) du problème approché de [eq:11] vers la solution \(u\) du problème exact [eq:F] lorsque \(V_h=\Pcho{1}\) ou plus généralement \(V_h=\Pcho{k},\; k\geq 1\).

Estimation en norme \(H^1\)

Il s’agit dans un premier temps d’estimer l’erreur \(u-u_h\) en norme \(H^1\).

Pour cela on part de l’estimation [eq:cea], on a

\begin{align} \|u-u_h\|{1,\Omega} &\leq c\; \inf{v_h \in \Pcho} \|u-v_h\|{1,\Omega}\\ & \leq c\; \|u-\Ich{u}\|{1,\Omega}\\ & \leq c\; h^k |u|_{k+1,\Omega} \end{align}

pourvu que la solution exacte soit suffisamment régulière, c’est à dire \(u \in H^{k+1}(\Omega)\).

On notera que \(\Ich{k}{u} \in \Pcho{k}\) puisque \(u \in H^1_0(\Omega)\) et donc que \(\Ich{u}(a)=\Ich{u}(b)=0\).

On a donc le résultat suivant

Soit un entier \(k\geq 1\).

On suppose que la solution du problème [eq:FV] est dans \(H^{k+1}(\Omega)\).

On note \(u_h\) la solution du problème approché [eq:11] avec l’espace d’approximation \(V_h = \Pcho{k}\).

Alors, il existe une constante \(c\), indépendante de \(h\), telle que

\[ \|u-u_h\|{1,\Omega} \leq c\; h^k |u|{k+1,\Omega} \]

On dit que l’estimation d’erreur [eq:13] est optimale car elle est du même ordre que l’erreur d’interpolation en norme \(H^1\), voir la proposition [prop:2].

1.6.8. Estimation en norme \(L^2\)

Avec les hypothèses de la proposition [prop:1] et en supposant que \(\alpha \in \Ck{1}(\bar{\Omega})\).

Alors, il existe une constante \(c\), indépendante de \(h\), telle que

\[ \|u-u_h\|{0,\Omega} \leq c\; h^{k+1} |u|{k+1,\Omega} \]

On dit que l’estimation d’erreur [eq:16] est optimale car elle est du même ordre que l’erreur d’interpolation en norme \(L^2\), voir la proposition [prop:2].

1.7. Formulation algébrique \(V_h=P_{c,h}^1\)

En décomposant la solution approchée \(u_h\) sur cette base sous la forme \({u_h = \sum_{i=1}^N \mu_i \; \varphi_i}\), on obtient, comme au paragraphe Principe général de la méthode des éléments finis, le système linéaire \(A\mu=b\), avec :

\begin{array}{rcl} A_{ji}=a(\varphi_i,\varphi_j) & = & \ds{\int_a^b \left[ \varphi_i'(x) \varphi_j'(x) + c(x) \varphi_i(x) \varphi_j(x)\right]\; dx }\\ & = & \ds{ \sum_{k=0}^N \int_{x_k}^{x_{k+1}} \left[\varphi_i'(x) \varphi_j'(x) + c(x) \varphi_i(x) \varphi_j(x)\right]\; dx } \end{array}

Le support de \(\varphi_i\) étant réduit à \([x_{i-1},x_{i+1}\)], on en déduit que

\left\{ \begin{array}{lll} a(\varphi_i,\varphi_j) & = & 0 \qquad \hbox{si }|i-j|\ge 2\\ & & \\ a(\varphi_i,\varphi_{i+1}) & = & \ds{ \int_{x_i}^{x_{i+1}} \left[ \varphi_i'(x) \varphi_{i+1}'(x) + c(x) \varphi_i(x) \varphi_{i+1}(x)\right] \; dx}\\ & & \\ a(\varphi_i,\varphi_{i-1}) & = & \ds{ \int_{x_{i-1}}^{x_i} \left[ \varphi_i'(x) \varphi_{i-1}'(x) + c(x) \varphi_i(x) \varphi_{i-1}(x)\right] \; dx}\\ & & \\ a(\varphi_i,\varphi_{i}) & = & \ds{ \int_{x_{i-1}}^{x_{i+1}} \left[ \varphi_i'^2(x) + c(x) \varphi_i^2(x)\right] \; dx}\\ \end{array}\right.

\(A\) est donc tridiagonale.

1.7.1. Exercices

  1. Dans le Définitions et théorèmes, montrer que, dans le cas où \(a\) est symétrique, si \(u\) est solution du problème variationnel, alors elle est solution du problème de minimisation.

  2. Montrer que \(\nabla J[u\)(v) = a(u,v) - l(v)].

  3. Montrer que, si \(a\) est coercive, la matrice \(A\) de ([eq:lin]) est inversible. (C’est donc la démonstration du théorème de Lax-Milgram en dimension finie.)

  4. Pour l’exemple 1-D traité dans ce chapitre, démontrer qu’on est bien dans les hypothèses du théorème de Lax-milgram

  5. Calculer explicitement la matrice \(A\) pour cet exemple.

  6. Pour le problème 2-D du §[sec:modele-2D], montrer que la formulation variationelle ([eq:FV2]) admet une solution unique, qui est aussi solution classique si \(f \in H^2(\Omega)\).

  7. Démontrer les résultats du §[sec:elliptique]

2. Eléments finis de Lagrange

On va présenter dans ce chapitre le type le plus simple et le plus classique d’éléments finis.

2.1. Unisolvance

Soit \(\Sigma=\{ a_1,\ldots,a_N\}\) un ensemble de \(N\) points distincts de \(\RR^n\). Soit \(P\) un espace vectoriel de dimension finie de fonctions de \(\RR^n\) à valeurs dans \(\RR\).On dit que \(\Sigma\) est \(P\)-unisolvant ssi pour tous réels \(\alpha_1,\ldots,\alpha_N\), il existe un unique élément \(p\) de \(P\) tel que \(p(a_i)=\alpha_i,\; i=1,\ldots,N\).

Ceci revient à dire que la fonction :

\begin{array}{rcl} {\cal L} : P & \longrightarrow & \RR^N\\ p & \longrightarrow & (p(a_1),\ldots,p(a_N)) \end{array}

est bijective.

En pratique, on montrera que \(\Sigma\) est \(P\)-unisolvant en vérifiant que \(\opdim P= Card \Sigma\), puis en montrant l’injectivité ou la surjectivité de \({\cal L}\).

L’injectivité de \({\cal L}\) se démontre en établissant que la seule fonction de \(P\) s’annulant sur tous les points de \(\Sigma\) est la fonction nulle.

La surjectivité de \({\cal L}\) se démontre en exhibant une famille \(p_1,\ldots,p_N\) d’éléments de \(P\) tels que \(p_i(a_j)=\delta_{ij}\), c’est à dire un antécédent pour \({\cal L}\) de la base canonique de \(\RR^N\).

En effet, étant donnés des réels \(\alpha_1,\ldots,\alpha_N\), la fonction \({p=\sum_{i=1}^N \alpha_i p_i}\) vérifie alors \(p(a_j)=\alpha_j,\; j=1,\ldots,N\).

2.2. Elément fini de Lagrange

Élément fini de Lagange

Un élément fini de Lagrange est un triplet \((K,\Sigma,P)\) tel que

  • \(K\) est un élément géométrique de \(\RR^n\) (\(n=\) 1, 2 ou 3), compact, connexe, et d’intérieur non vide.

  • \(\Sigma=\{a_1,\ldots,a_N\}\) est un ensemble fini de \(N\) points distincts de \(K\).

  • \(P\) est un espace vectoriel de dimension finie de fonctions réelles définies sur \(K\), et tel que \(\Sigma\) soit \(P\)-unisolvant (donc dim \(P = N\)).

Fonctions de base locales

Soit \((K,\Sigma,P)\) un élément fini de Lagrange. On appelle fonctions de base locales de l’élément les \(N\) fonctions \(p_i\) (\(i=1,\ldots,N\)) de \(P\) telles que

\[ p_i(a_j)=\delta_{ij}\qquad 1\le i,j \le N. \]

On vérifie aisément que \((p_1,\ldots,p_N)\) ainsi définie forme bien une base de \(P\).

On appelle opérateur de interpolation (ou encore \(P-\)interpolation) sur \(\Sigma\) l’opérateur \(\pi_K\) qui, à toute fonction \(v\) définie sur \(K\), associe la fonction \(\pi_K v\) de \(P\) définie par \({\pi_K v = \sum_{i=1}^N v(a_i)\, p_i}\).

\(\pi_K v\) est donc l’unique élément de \(P\) qui prend les mêmes valeurs que \(v\) sur les points de \(\Sigma\).

2.3. Exemples d’éléments finis de Lagrange

2.3.1. Espaces de polynômes

On notera \(\Pk\) l’espace vectoriel des polynômes de degré total inférieur ou égal à \(k\).

  • Sur , \(\Pk=\hbox{Vect} \{ 1,X,\ldots, X^k\}\quad\) et dim \(\Pk=k+1\).

  • Sur \(^2\), \(\Pk=\hbox{Vect} \{ X^i Y^j ; 0\le i+j \le k\}\quad\) et dim \({ \Pk=\frac{(k+1)(k+2)}{2}}\).

  • Sur \(^3\), \(\Pk=\hbox{Vect} \{ X^i Y^j Z^l ; 0\le i+j+l \le k\}\quad\) et dim \({ \Pk=\frac{(k+1)(k+2)(k+3)}{6}}\).

On notera \(\Qk\) l’espace vectoriel des polynômes de degré inférieur ou égal à \(k\) par rapport à chaque variable.

  • Sur \(\RR\), \(\Qk=\Pk\).

  • Sur \(\RR^2\), \(\Qk=\hbox{Vect} \{ X^i Y^j ; 0\le i,j \le k\}\quad\) et dim \({ \Qk=(k+1)^2}\).

  • Sur \(\RR^3\), \(\Qk=\hbox{Vect} \{ X^i Y^j Z^l ; 0\le i,j,l \le k\}\quad\) et dim \({ \Qk=(k+1)^3}\).

2.3.2. Exemples 1-D

Elément \(P_1\)
  • \(K=[a,b\)]

  • \(\Sigma=\{a,b\}\)

  • \(P=\Pk[1\)]

Elément \(P_2\)
  • \(K=[a,b\)]

  • \({\Sigma=\{a,\frac{a+b}{2},b\}}\)

  • \(P=\Pk[2\)]

Elément \(\Pk\)
  • \(K=[a,b\)]

  • \({\Sigma=\{a+i\,\frac{b-a}{m},\quad i=0,\ldots,k\}}\)

  • \(P=\Pk\)

2.3.3. Exemples 2-D triangulaires

Elément \(\Pk[1\)]
  • \(K\)=triangle de sommets \(a_1, a_2, a_3\)

  • \(\Sigma=\{a_1,a_2,a_3\}\)

  • \(P=\Pk[1\)] Les fonctions de base sont définies par \(p_i(a_j)=\delta_{ij}\). Ce sont donc les coordonnées barycentriques : \(p_i=\lambda_i\) (cf annexe [ch:bary]).

Elément \(\Pk[2\)]
  • \(K\)=triangle de sommets \(a_1, a_2, a_3\)

  • \(\Sigma=\{a_1,a_2,a_3, a_{12}, a_{13}, a_{23}\}\), où \({a_{ij}=\frac{a_i+a_j}{2}}\).

  • \(P=\Pk[2\)]

Les fonctions de base sont \(p_i=\lambda_i (2\lambda_i -1)\) et \(p_{ij}=4\lambda_i\lambda_j\). Un exemple de calcul de ces fonctions de base est donné en annexe [ch:bary].

Eléments finis triangulaire \(P_1\), triangulaire \(P_2\) et rectangulaire \(Q_1\)

elements 2D

2.3.4. Exemples 2-D rectangulaires

Elément \(\Qk[1\)]
  • \(K\)=rectangle de sommets \(a_1, a_2, a_3, a_4\), de côtés parallèles aux axes

  • \(\Sigma=\{a_1,a_2,a_3,a_4\}\)

  • \(P=\Qk[1\)] Les fonctions de base sont \({p_i(X,Y)=\frac{(X-x_j)(Y-y_j)}{(x_i-x_j)(y_i-y_j)}}\), où \((x_i,y_i)\) sont les coordonnées de \(a_i\), et où \(a_j\), de coordonnées \((x_j,y_j)\) est le coin opposé à \(a_i\).

2.3.5. Exemples 3-D

Elément tétraèdrique \(\Pk[1\)]
  • \(K\)=tétraèdre de sommets \(a_1, a_2, a_3, a_4\)

  • \(\Sigma=\{a_1,a_2,a_3, a_4\}\)

  • \(P=\Pk[1\)]

Elément tétraèdrique \(\Pk[2\)]
  • \(K\)=tétraèdre de sommets \(a_1, a_2, a_3, a_4\)

  • \({ \Sigma=\{a_i\}_{1\le i\le 4} \cup \{a_{ij}\}_{1\le i < j \le 4} }\)

  • \(P=\Pk[2\)]

Les fonctions de base sont \(p_i=\lambda_i (2\lambda_i -1)\) et \(p_{ij}=4\lambda_i\lambda_j\).

Elément parallélépipèdique \(Q_1\)
  • \(K\)=parallélépipède de sommets \(a_1, \ldots , a_8\) de côtés parallèles aux axes

  • \(\Sigma=\{a_i\}_{1\le i\le 8}\)

  • \(P=\Qk[1\)]

Elément prismatique
  • \(K\)=prisme droit de sommets \(a_1, \ldots , a_6\)

  • \(\Sigma=\{a_i\}_{1\le i\le 6}\)

  • \(P=\{p(X,Y,Z)=(a+bX+cY)+Z(d+eX+fY), \;\; a,b,c,d,e,f \in \RR\}\)

Eléments finis tétraèdriques \(P_1\) et \(P_2\), parallélépipèdique \(Q_1\), et prismatique

elements 3D

2.4. Famille affine d’éléments finis

Équivalence affine

Deux éléments finis \((\hat{K},\hat{\Sigma},\hat{P})\) et \((K,\Sigma,P)\) sont affine-équivalents ssi il existe une fonction affine \(F\) inversible (\(F: \hat{x} \longrightarrow B\hat{x}+b\)) telle que (i) \(K=F(\hat{K})\) (ii) \(a_i=F(\hat{a}_i) \qquad i=1,\ldots,N\) et (iii) \(P=\{ \hat{p}\circ F^{-1} , \quad \hat{p}\in \hat{P} \}\)

Si l’on est dans \(\RR^n\), \(B\) est donc une matrice \(n\times n\) inversible, et \(b\) est un vecteur de \(\RR^n\).

Soient \((\hat{K},\hat{\Sigma},\hat{P})\) et \((K,\Sigma,P)\) deux éléments finis affine-équivalents, via une transformation \(F\).

On note \(\hat{p}_i \; (i=1,\ldots,N)\) les fonctions de base locales de \(\hat{K}\).

Alors les fonctions de base locales de \(K\) sont les \(p_i=\hat{p}_i\circ F^{-1}\).

On appelle famille affine d’éléments finis une famille d’éléments finis tous affine-équivalents à un même élément fini \((\hat{K},\hat{\Sigma},\hat{P})\), appelé élément de référence.

D’un point de vue pratique, le fait de travailler avec une famille affine d’éléments finis permet de ramener tous les calculs d’intégrales à des calculs sur l’élément de référence.

Les éléments de référence sont :

  • En 1-D : le segment \([0,1]\)

  • En 2-D triangulaire : le triangle unité, de sommets \((0,0)\), \((0,1)\) et \((1,0)\).

  • En 2-D rectangulaire : le carré unité \([0,1]\times[0,1]\).

  • En 3-D tétraèdrique : le tétraèdre unité, de sommets \((0,0,0)\), \((1,0,0)\), \((0,1,0)\) et \((0,0,1)\).

  • En 3-D parallélépipèdique : le cube unité \([0,1]\times[0,1]\times[0,1]\).

  • En 3-D prismatique : le prisme unité de sommets \((0,0,0)\), \((0,1,0)\), \((1,0,0)\), et \((0,0,1)\), \((0,1,1)\), \((1,0,1)\).

2.5. Maillages

Nous étendons ici aux dimension 2 et 3 les notions élémentaires de maillage vues en 1D, voir la figure Maillage Feel++.

Maillage Feel++

feelpp mesh

Un maillage est constituée d’une famille d’éléments(ou mailles ou cellules) \(\{K_e\}_{e=1,...,N_e}\) où \(N_e\) est le nombre d’éléments, nous noterons

\[ \calTh = {K_m}_{m=1,…​,N_e} \] avec

\[ h=\max_{1\le e\le N_e} h_{K_e} \] et

\[ h_{K_e} = \diam(K_e) = \max_{x_1,x_2 \in K_e} \|x_1-x_2\|,\, e \in \{1,…​,\Ne\} \]

On travaille par la suite avec des familles de maillage et on les note \(\set{\mathcal{T}_h}_{h > 0}\).

Famille de maillage quasi-uniforme

On dira qu’une famille de maillage \(\set{\mathcal{T}_h}_{h > 0}\) est quasi-uniforme s’il existe une constante \(c\) telle que

\[ \forall h,\ \forall K \in \calTh,\ h_K \geq c h \]

Cela veut dire que les élements sont tous de la même taille pour \(h\) donné.

2.5.1. Transformation géométrique

Un maillage est généré par

  1. un élément de reference noté \(\hat{K}\)

  2. une famille de transformations géométriques mappant \(\hat{K}\) vers les éléments \(K_e, e=1,\ldots,\Ne\) dans le maillage

Nous supposerons que les transformations sont des \(\mathcal{C}^1-\) diffeomorphismes [1].

Pour une cellule \(K \in \mathcal{T}_h\), on note \(T_K\) la transformation géométrique

\[ T_K: \hat{K} \mapsto K \]

Afin de spécifier la transformation géométrique, on considère l’élément fini de Lagrange, noté \((\hat{K},\hat{P}_{\mathrm{geo}}, \hat{\Sigma}_{\mathrm{geo}})\), tel que

  • \(\ngeo = \card{\hat{\Sigma}_{\mathrm{geo}}}\)

  • \(\set{\hat{g}_1,\ldots,\hat{g}_{\ngeo}}\) les noeuds de \(\hat{K}\)

  • \(\set{\hat{\psi}_1,\ldots,\hat{\psi}_{\ngeo}}\) les fonctions de forme

Élément fini géométrique

On dit que \((\hat{K},\hat{P}_{\mathrm{geo}}, \hat{\Sigma}_{\geo})\) est l’élément fini géométrique, \(\set{\hat{g}_1,\ldots,\hat{g}_{\ngeo}}\) sont les noeuds géométriques et \(\set{\hat{\psi}_1,\ldots,\hat{\psi}_{\ngeo}}\) sont les fonctions de formes géométriques

Transformation géométrique associée à un triangle

geomap triangle

Pour chaque \(K \in \mathcal{T}_h\), on a un \(\ngeo\)-uplet \(\set{g^K_1,\ldots,g^K_\ngeo}\).

La transformation géométrique est définie comme suit

T_K: \hat{x} \in \hat{K} \mapsto \sum_{i=1}^\ngeo\ g^K_i \hat{\psi}_i(\hat{x})

et en particulier on a

T_K(\hat{g}_i) = g^K_i, \quad \forall i \in \set{1,\ldots,\ngeo}
On a \(T_K \in [\hat{P}_\geo(\hat{K})\)^d] et que \(\set{g^K_1,\ldots,g^K_\ngeo}\) sont les noeuds géométriques de \(K\).

\(T_K\) est un \(\mathcal{C}^1\)-diffeomorphism donc la numérotation des noeuds \(\set{g^K_1,\ldots,g^K_\ngeo}\) doit être compatible avec les noeuds de l’élément finit géométrique.

La numérotation locale des entités géométriques dans doit être consistente avec la numérotation locale des générateurs de maillage. voir \(\triangleright\) format de fichier Gmsh pour une numérotation locale.

Un cas particulier est la transformation géométrique affine.

Maillage affine

Quand toutes les transformations géométriques \(\set{T_K}_{K \in \mathcal{T}_h}\) sont affines, cela veut dire que pour tout \(K \in \mathcal{T}_h\), il existe un vecteur \(b_K \in \RR^d\) et une matrice \(J_K \in \RR^{d\times d}\) tels que

\[ T_K : \hat{x} \in \hat{K} \mapsto b_K + J_K \hat{x} \in K \] On dit que le maillage est affine.

Si l’élément fini géométrique est \((\hat{K},\poly{P}_1,\Sigma_\ngeo)\) alors les éléments \(K\) sont soit des triangles soit des tétrahèdres.
Mesh<Simplex<d,1> > ou Mesh<Simplex<d> > est le type pour les maillages affines formés de simplexes dans \(\RR^d\). 1 indique l’ordre de l’élément fini géométrique et est la valeur par défaut.

2.5.2. Quelques calculs avec la transformation géométrique

Gradient, Inverse et Jacobien

On note \(\xi\) un ensemble de \(n\) points dans \(\hat{K}\) et on note \(\nabla T_K(\xi)\) le gradient de \(T_K\) aux points \(\xi\)

\nabla T_K( \xi )\ =\ \sum_{i=0}^{\ngeo}\ g^K_i\ \nabla \psi_i (\xi)

et \(B_K(\xi) = \nabla T_K^{-1}(\xi)\) l’inverse \(\xi\) et finalement \(J_K(\xi)\) le jacobien de \(T_K\) en \(\xi\)

J_K(\xi)\ =\ |\det( \nabla T_K(\xi) )|

Dérivation dans l’élément de référence

Afin de dériver un polynome dans l’élément réel \(K\), grâce à la transformation géométrique et la règle de différentiation des fonctions composées, nous dérivons seulement dans l’élément de référence \(\hat{K}\).

Soit \(f: K \mapsto \RR\) et \(\hat{f}: \hat{K} \mapsto \RR\) telle que \(\hat{f} = f \circ T_K\)

\nabla f\ =\  \hat{\nabla} \underbrace{\hat{f}(\xi)}_{f \circ T_K(\xi)} B_K(\xi)

en 2D, on a, en notant \(x=T_K(\xi)\),

\nabla f(x) =
\begin{pmatrix}
  \frac{\hat{\partial} \hat{f} (\xi)}{\partial \xi_1} & \frac{\hat{\partial} \hat{f} (\xi)}{\partial \xi_2}
\end{pmatrix}
\begin{pmatrix}
  B_{K_{11}}(\xi) & B_{K_{12}}(\xi)\\
  B_{K_{21}}(\xi) & B_{K_{22}}(\xi)\\
\end{pmatrix}
Intégration dans l’élément de référence

De manière similaire, au lieu de calculer les intégrales sur l’élément réel \(K\), nous appliquons un changement de variables et calculons les intégrales sur l’élément de réference \(\hat{K}\).

Soit \(f: K \mapsto \RR\) et \(\hat{f}: \hat{K} \mapsto \RR\) telle que \(\hat{f} = f \circ T_K\), et \({\mathbf{F}}: K \mapsto \RR^d\) et \({\hat{\mathbf{F}}}: \hat{K} \mapsto \RR^d\) telle que \(\hat{{\mathbf{F}}} = {\mathbf{F}} \circ T_K\),

On a alors les relations suivantes:

\int_{K} \ f\ dx\ =\ \int_{\hat{K}} f(T_K(\xi) ) J_K( \xi )\ d \xi \ =\ \int_{\hat{K}} \hat{f}(\xi) J_K( \xi )\ d \xi]

\int_{K}\ \nabla f\ dx\ =\ \int_{\hat{K}} \Big(\hat{\nabla} \underbrace{\hat{f}(\xi)}_{f \circ T_K(\xi)} B_K(\xi)\Big) J_K( \xi )\ d \xi

\int_{\partial K} f( x ) dx = \int_{\partial \hat{K}} \hat{f}(\xi) \| B_K(\xi) {\hat{\mathbf}}(\xi) \| J_K( \xi ) d \xi

\int_{\partial K}\ {\mathbf{F}}( x )\ \cdot\ {\mathbf}(x) dx = \int_{\partial \hat{K}} {\hat{\mathbf{F}}}( \xi )\ \cdot \Big(B_K(\xi) {\hat{\mathbf}}(\xi) \Big) J_K( \xi )\ d \xi

où \({\mathbf{n}}(x)\) est la normale extérieure unitaire à \(\partial K\) évaluée en \(x \in \partial K\), et \({\hat{\mathbf{n}}}(\xi)\) la normale unitaire extérieure à \(\hat{K}\) évaluée en \(\xi \in \partial \hat{K}\).

Feel++ effectue automatiquement pour vous les changements de variables dans les intégrales.
Maillage conforme

Un maillage est dit conforme si l’intersection de deux éléments est soit vide, un sommet, une arête ou une face.

On ne manipule que des maillages conformes dans le cours mais Feel++ peut traiter des maillages non-conformes par exemple dans le contexte de méthode de décomposition de domaines.

2.6. Espaces élément fini de Lagrange

Soit \(\mathcal{T}_h\) un maillage généré par \((\hat{K}, \hat{P}_{\mathrm{geo}}, \hat{\Sigma}_{\mathrm{geo}})\), une cellule \(K \in \mathcal{T}_h\) est alors l’image de \(\hat{K}\) par la transformation géométrique \(T_K\) défini par [eq:23].

L’objectif à présent est de générer la famille d’éléments finis de Lagrange grâce à l’élément fini de référence \((\hat{K},\hat{P}, \hat{\Sigma})\)

\{(K,P_K,\Sigma_K)\}_{K \in \mathcal{T}_h}

On note \(\{\hat{x}_1,...,\hat{x}_{\nf}\}\) les noeuds de l’élément fini.

On note \(\{\hat{\Psi}_1,...,\hat{\Psi}_{\nf}\}\) les fonctions de forme élément fini.

Soit \(K \in \mathcal{T}_h\) et \(P_K=\{\hat{p} \circ T_K^{-1}; \hat{p} \in \hat{P}\}\).

Pour tout \(i \in \{1,...,\nf\}\), on pose \(x_{K,i} = T_K(\hat{x}_i)\)

On note \(\Sigma_K\) l’ensemble des degrés de liberté associé à \(\{x_{K,1},...,x_{K,\nf}\}\)

Alors \((K,P_K,\Sigma_K)\) est un élément fini de Lagrange.

Les fonctions de forme sont définies de la façon suivante

\[ \Psi_{K,i} = \hat{\Psi}_i \circ T_K^{-1}, \quad i=1,…​,\nf \] et l’opérateur d’interpolation local comme

\[ \Ilag{K}: v \in \mathcal^0(K) \mapsto \sum_{i=1}^\nf\ v(x_{K,i})\ \Psi_{K,i}\ \in P_K \] Une propriété importante de \(\Ilag{K}\) est que

\[ \forall v \in \mathcal^0(K),\quad \Ilag{K}( v \circ T_K ) = \Ilag{K}( v ) \circ T_K. \]

Soit \(T_K\) une transformation affine.

Soit \(\mathbb{P}_k \subset \hat{P}\) et \(k+1 > \frac{d}{2}\).

Soit \(h_K\) le diamètre de \(K\) et \(\rho_K\) le diamètre de la plus grande boule inscrite dans \(K\) et \(\omega_K = \frac{h_K}{\rho_K}\).

Alors il existe une constante \(c\) independente de \(K\) telle que \(\forall v \in H^{k+1}(K)\) et pour tout \(m \in \{0,...,k+1\}\),

\[ |v - \Ilag{K}(v)|{m,K} \leq c h^{k+1-m} \omega_K^m |v|{k+1,K} \]

  • \(\omega_K\) devrait être aussi proche de \(1\) que possible

  • La deuxième hypothèse technique permet d’assurer que \(H^{k+1}(K) \subset C^0(K)\).

  • on obtient des résultats similaires si \(v\) n’est pas suffisamment régulière

Espace \(H^1\)-conforme

Un espace vectoriel \(V_h\) de fonctions définies sur un domaine \(\Omega_h\) est dit être \(H^1\)-conforme si \(V_h \subset H^1(\Omega_h)\).

Afin de construire un tel espace on introduit tout d’abord

W_h = \{v_h \in L^2(\Omega_h); \forall K \in \mathcal{T}_h, v_h|_K \in P_K\}

mais ce n’est pas suffisant: les fonctions \(W_h\) peuvent avoir des sauts entre les éléments du maillage. Nous avons donc besoin d’assurer la continuité de ces fonctions

V_h = W_h \cap C^0(\Omega_h) = \{ v_h \in W_h; \forall F \in \mathcal{F}^i_h. \jump{v_h}_F = 0\}

Concernant l’implémentation, nous avons besoin de d’indentifier les degrés de liberté communs entre les éléments quand nous construisons la tables des degrés de liberté.

Voici deux exemples d’espace \(H^1\)-conforme

Espace \(H^1\)-conforme sur des simplexes (e.g. triangle)
P^k_{c,h} = \{ v_h \in C^0(\Omega_h); \forall K \in \mathcal{T}_h, v_h
    \circ T_K \in \mathbb{P}_k\}
Espace \(H^1\)-conforme sur des hypercubes (e.g quadrilatère)
Q^k_{c,h} = \{ v_h \in C^0(\Omega_h); \forall K \in \mathcal{T}_h, v_h
\circ T_K \in \mathbb{Q}_k\}

2.6.1. Projections orthogonales

On note

\begin{aligned} \Pi^{0,k}_{c,h} : Lˆ2(\Omega) \rightarrow P^k_{c,h}\\ \Pi^{1,k}_{c,h} : H^1(\Omega) \rightarrow P^k_{c,h} \end{aligned}

les projections orthogonales associées respectivement aux produits scalaires \((u,v)_{0,\Omega} = \int_\Omega u v\) and \((u,v)_{1,\Omega} = \int_\Omega u v + \int_\Omega \nabla u \cdot \nabla v\)

On a pour \(l=1,...,k\) et si \(v \in H^{l+1}(\Omega)\)

\begin{array}{rl} \|v - \Pi^{0,k}_{c,h}(v)\|{0,\Omega} & \leq c h^{l+1} |v|{k+1,\Omega} \\ \|v - \Pi^{1,k}_{c,h}(v)\|{1,\Omega} & \leq c h^{l} |v|{k+1,\Omega} \end{array}

Pour calculer \(\Pi^{0,k}_{c,h}\) et \(\Pi^{1,k}_{c,h}\) on a besoin de résoudre les problèmes

Projection \(l^2\)

Soit \(v\) une fonction de \(L^2\), calculer \(\Pi^{0,k}_{c,h}(v) \in P^k_{c,h}\) tel que \(\forall v_h \in P^k_{c,h}\) on a

\[ (\Pi^{0,k}_{c,h}(v), v_h){0,\Omega} = (v, v_h){0,\Omega} \]

Projection \(H^1\)

Soit \(v\) une fonction de \(H^1\), calculer \(\Pi^{1,k}_{c,h}(v) \in P^k_{c,h}\) tel que \(\forall v_h \in P^k_{c,h}\) on a

\[ (\Pi^{1,k}_{c,h}(v), v_h){1,\Omega} = (v, v_h){1,\Omega} \]

2.6.2. Interpolant de Lagrange

Notons \(V_h\) un espace \(H^1\) conforme, \(\{\Psi_i\}_{1,...,N}\) une base nodale de \(V_h\) et \(\{x_1,...,x_N\}\) les noeuds associés alors l’interpolant de Lagrange est défini par

\[ \Ilag{h}: v \in C^0(\Omega_h) \mapsto \sum_{i=1}^N v(x_{i}) \Psi_i \in V_h \]

La restriction de l’interpolant de Lagrange à une cellule \(K\) coincide avec l’interpolant de Lagrange appliqué à la fonction dans la cellule \(K\):

\[ \Ilag{h}(v)|_K = \Ilag{h}(v|_K) \]

Supposons \(\{\mathcal{T}_h\}_{h>0}\) une famille de maillage quasi-uniforme et conformes, \(\mathbb{P}_k \subset \hat{P}\) et \(k+1 > \frac{d}{2}\).

Alors il existe une constante \(c\) telle que pour tout \(h\) et \(v \in H^{k+1}(\Omega_h)\)

\[ \|v - \Ilag{h}(v)\|{0,\Omega_h} + h |v - \Ilag{h}(v)|{1,\Omega_h} \leq c h^{k+1} |v|_{k+1,\Omega_h} \]

Nous allons vérifier sur un exemple ce théorème.

Nous considérons pour cela \(\alpha\) un réel et \(\mathbf{x}=(x_1,...,x_d) \in \mathbb{R}^d\) un point de \(\Omega = [0,1]^d, d=1,2,3\) et \(v\) la fonction définie par

\begin{array}{rl} v : \Omega &\rightarrow \mathbb{R}\\ \mathbf{x} &\rightarrow ( \mathbf{x} \cdot \mathbf{x} )^{\alpha/2}\ \Pi_{i=1}^d( 1-x_i^2) \end{array}

Nous construisons l’interpolant de Lagrange de \(v\) dans \(P^k_{c,h}\) avec \(k=1,...,5\) et \(d=1,2,3\) et étudions l’erreur d’interpolation \(L^2\) et \(H^1\) du théorème [eq:62] en échelle log-log.

Nous devons obtenir des droites de pentes \(k\) (resp. \(k+1\)) pour la norme \(L^2\) (resp. \(H^1\).)

2.7. Approximation iso-paramétrique

Quand le domaine est courbe, si nous désirons obtenir des propriétés de convergence optimale nous avons besoin de discretiser le bord du domaine avec suffisamment de précision

Notons \((\hat{K},\hat{P}_{\mathrm{geo}}, \hat{\Sigma}_{\geo})\) l’élément fini géométrique et \((\hat{K},\hat{P}, \hat{\Sigma})\) l’élément fini de référence pour \(V_h\).

Si \(\hat{P}_{\mathrm{geo}} = \hat{P}\) l’approximation est dite iso-parametrique.

Si \(\hat{P}_{\mathrm{geo}} \subset \hat{P}\) l’approximation est dite sub-parametrique.

Si \(\hat{P} \subset\hat{P}_{\mathrm{geo}}\) l’approximation est dite sur-parametrique.

Gmsh (gmsh) un mailleur libre permet de générer des maillage d’ordre élevé jusqu’à l’ordre 5 en 2D et 4 en 3D

Supposons que \(\{\mathcal{T}_h\}_{h>0}\) une famille de maillage quasi-uniformes et conformes, \(\mathbb{P}_k \subset \hat{P}\) et \(k+1 > \frac{d}{2}\) et \(k_{\mathrm{geo}} = k\).

Alors il existe une constante \(c\) telle que pour tout \(h\) et \(v \in H^{k+1}(\Omega_h)\)

\[ \|v - \Ilag{h}(v)\|{0,\Omega_h} + h |v - \Ilag{h}(v)|{1,\Omega_h} \leq c h^{k+1} |v|_{k+1,\Omega_h} \]

Les résultats sont identiques à ceux du theorème [thr:15].

Nous allons vérifier sur un exemple à l’erreur d’approximation de l’élément géométrique. Considérons les cercles unité généré par une transformation affine, noté \(\Omega^1_h\), et d’ordre 2, noté \(\Omega^2_h\), et calculons leur aire respective.

Construisons une famille de maillage \(\{\calTh\}_{h >0}\), par exemple \(h\in \{0.4, 0.2, 0.1, 0.05\}\) et calculons l’erreur entre le calcul exact de l’aire \(\pi\) et le calcul numérique \(\int_{\Omega^1_h} 1\) et \(\int_{\Omega^2_h} 1\) respectivement.

Le listing suivant présente le code C++ pour effectuer cela

Calcul de l’aire d’un cercle
Unresolved directive in fem/ch-ef-lagrange.adoc - include::../../../codes/isoparam/circle.cpp[tag=isoparam,indent=0]

La table Convergence de l’approximation \(P_1\) et \(P_2\) présente les erreurs d’approximation et la figure [fig:circle] présente les courbes de convergence en échelle log-log ainsi que les pentes associées à ces courbes.

On s’attend d’après le théorème [eq:62] appliqué à des pentes à l’élément fini géométrique.

Table 1. Convergence de l’approximation \(P_1\) et \(P_2\)
Unresolved directive in fem/ch-ef-lagrange.adoc - include::../../../codes/isoparam/circle.csv[]

La table toto présente les résultats de convergence.

On observe un phénomène de super-convergence pour le cas \(\Omega^2_h\), on obtient un ordre de convergence \(4\) et nous devrions obtenir \(3\).

2.8. Feel++

En Feel++, un maillage Mesh est décomposé en un ensemble d’éléments décomposés en sous entités (volume,face,arête,point):

  • faces(decomposés en sous entités),

  • arêtes(decomposés en sous entités) et

  • points.

et à chaque élément \(K\) est associé une transformation géométrique \(T_K\).

L’ordre polynomial de la transformation géométrique est donné par le second argument template

Type d’un object maillage en Feel++
Mesh<Simplex<d, k_geo>> (1)
Mesh<Hypercube<d, k_geo>> (2)
1 Un maillage de simplexes en dimension \(d\) avec transformation géométrique d’ordre \(k_{\mathrm{geo}}\).
2 Un maillage d’hypercubes en dimension \(d\) avec transformation géométrique d’ordre \(k_{\mathrm{geo}}\).

Ain de parcourir les éléments et faces du maillage, Feel++ fournit des fonctions renvoyant des itérateurs (début et fin) sur ces ensembles

elements(mesh)

retourne 2 itérateurs sur l’ensemble des éléments du maillage

markedelements(mesh,<int>)

et markedelements(mesh,<string>) retourne 2 itérateurs sur les éléments marqués par l’entier <int> et la chaîne des caractères <string> respectivement, cela correspondra typiquement à des propriétés de matériau

boundaryfaces(mesh)

retourne 2 itérateurs sur les faces au bord du maillage

markedfaces(mesh,<int>) et markedelements(mesh,<string>)

retourne 2 itérateurs sur les faces marquées par l’entier <int> et la chaîne des caractères <string> respectivement, ca correspondra typiquement aux conditions aux limites

link pending to Feel++ user manual mesh iterators section.

L’espace d’approximation \(V_h\) \(H^1\) conforme (espaces de functions continues sur \(\Omega\) polynomiales par morceaux de degré \(\leq k\)) est défini comme suit

FunctionSpace<Mesh<Simplex<d, k_geo> >, bases<Lagrange<k> > > V_h; (1)
FunctionSpace<Mesh<Hypercube<d, k_geo> >, bases<Lagrange<k> > > V_h; (2)
1 \(\Pch{k}\)
2 \(\Qch{k}\)

2.9. Du problème global aux éléments locaux

On va maintenant faire le lien entre la résolution d’un problème par méthode d’éléments finis et les notions qui viennent d’être introduites.

Soit une EDP à résoudre sur un domaine \(\Omega\), et \(V\) l’espace de Hilbert dans lequel on cherche une solution de la formulation variationnelle du problème.

On réalise un maillage de \(\Omega\) par une famille affine de \(N_e\) éléments finis \((K_i,\Sigma_i,P_i)_{i=1,\ldots,N_e}\).

Par unisolvance, la solution approchée \(u_h\) sera entièrement définie sur chaque élément \((K_i,\Sigma_i,P_i)\) par ses valeurs sur les points de \(\Sigma_i\), qu’on appellera les noeuds du maillage.

Il est à noter qu’un noeud sera en général commun à plusieurs éléments adjacents.

Le nombre total de noeuds \(N_h\) est donc inférieur à \(N_e\times\hbox{Card} \Sigma_i\), et on a dim \(V_h = N_h\).

Notons \(a_1,\ldots,a_{N_h}\) les noeuds du maillage.

Le problème approché se ramène donc à la détermination des valeurs de \(u_h\) aux points \(a_i\): ce sont les degrés de liberté du problème approché.

On va construire une base de \(V_h\) en associant à chaque ddl \(a_i\) un vecteur de la base. On définit ainsi les fonctions de base globales \(\varphi_i\) (\(i=1,\ldots,N_h\)) par

Fonctions de base globales

{\varphi_i}_{|K_j} \in P_j, \quad j=1,\ldots,N_e \mbox{ et } \varphi_i(a_j)=\delta_{ij}, 1\le i,j \le N_h

L’espace d’approximation interne est donc alors :

V_h = \hbox{Vect }\left\{\varphi_1,\ldots,\varphi_{N_h}\right\}

Exemple de fonction de base globale (élément triangulaire \(P_1\))

fonction globale

Il est facile de remarquer qu’une telle fonction \(\varphi_i\) est nulle partout, sauf sur les éléments dont \(a_i\) est un noeud.

En effet, si \(a_i\) n’appartient pas à un élément \(K\), \(\varphi_i\) est nulle sur tous les noeuds de \(K\), et donc sur \(K\) tout entier par unisolvance.

De plus, sur un élément \(K\) dont \(a_i\) est un noeud, \(\varphi_i\) vaut 1 sur \(a_i\) et 0 sur les autres noeuds de \(K\).

Donc \(\varphi_i\, _{|K}\)est une fonction de base locale de \(K\).

la fonction de base globale \(\varphi_i\) est donc construite comme la réunion des fonctions de base locales sur les éléments du maillage dont \(a_i\) est un noeud.

C’est à ce niveau que se situe le lien entre les définitions locales introduites au Elément fini de Lagrange et le problème global approché à résoudre.

Par ailleurs, ceci implique que tous les calculs à effectuer sur les fonctions de base globales peuvent se ramener à des calculs sur les fonctions de base locales, et donc simplement à des calculs sur l’élément de référence (car on a maillé le domaine avec une famille d’éléments finis affine-équivalents).

Maillage non conforme

conforme

Ce type de définition des fonctions de base n’est possible que si le maillage est conforme, c’est à dire si l’intersection entre deux éléments est soit vide, soit réduite à un sommet ou une arête en dimension 2 (ou à un sommet, une arête ou une face en dimension 3).

On interdit ainsi les situations du type de celle de la figure.

2.10. Exercices

  1. Calculer les fonctions de base locales des éléments finis de Lagrange introduits dans ce chapitre.

  2. Donner l’allure des fonctions de base globales correspondantes. Sont-elles continues ? dérivables ?

  3. Pour les éléments finis de Lagrange introduits dans ce chapitre, écrire le changement de variable affine entre élément quelconque et élément de référence.

Convergence de la méthode des éléments finis

1. Introduction

On suppose ici que l’on résout un problème sur un domaine \(\Omega\in\RR^n\) de façon approchée par méthode d’éléments finis.

Le but de ce chapitre est de fournir une estimation de l’erreur \(\|u-u_h\|_m\) où \(\|.\|_m\) désigne la norme \(H^m\).

La régularité de \(u\) et de \(u_h\) (et donc les valeurs possibles pour \(m\)) dépendant évidemment du problème continu et du type d’éléments finis choisis pour sa résolution, on exposera ici la démarche de façon générale, en supposant les fonctions suffisamment régulières par rapport à la valeur de \(m\).

En pratique, on aura le plus souvent \(m=0, 1, 2.\)

On notera \({\cal T}_h\) le maillage de \(\Omega\) considéré.

On supposera ici le domaine \(\Omega\) polygonal, ce qui permet de recouvrir exactement \(\Omega\) par le maillage.

Si ce n’est pas le cas, les calculs qui suivent doivent être modifiés pour tenir compte de l’écart entre le domaine couvert par le maillage et le domaine réel.

Les différentes étapes du calcul seront, de façon assez schématique, les suivantes :

L’erreur d’approximation est bornée par l’erreur d’interpolation

\(\|u-u_h\|_m \le C \|u-\pi_h u\|_m\)

L’erreur d’approximation est bornée par l’erreur d’interpolation

\(\|u-u_h\|_m \le C \|u-\pi_h u\|_m\)

On se ramène à des majorations locales sur chaque élément

\(\|u-\pi_h u\|^2_{m} = \sum_{K\in{\cal T}_h} \|u-\pi_h u\|_{m,K}^2\)

On se ramène à l’élément de référence

\({\|u-\pi_h u\|_{m,K} \le C(K) \|\hat{u}-\hat{\pi}\hat{u}\|_{m,\hat{K}} }\)

Majoration sur l’élément de référence

\(\|\hat{u}-\hat{\pi} \hat{u}\|_{m,\hat{K}} \le \hat{C} |\hat{u}|_{k+1,\hat{K}}\)

Assemblage des majorations locales

\({\|u-\pi_h u\|_m \le C' h^{k+1-m} |u|_{k+1} }\)

2. Calcul de majoration d’erreur

2.1. Etape 1: majoration par l’erreur d’interpolation

L’équation ([eq:cea]) établie dans la section Estimation d’erreur indique que

\|u-u_h\|_m \le \frac{M}{\alpha}\; \|u-v_h\|_m \quad \forall v_h\in V_h

Cette majoration est aussi appelée Lemme de Céa

On peut l’appliquer dans le cas particulier où \(v_h=\pi_h u\), ce qui donne

\|u-u_h\|_m \le \frac{M}{\alpha}\; \|u-\pi_h u\|_m

2.2. Etape 2: Décomposition sur les éléments

On a, avec des notations évidentes :

\[\begin{array}{lll} {\|u-\pi_h u\|_m^2 }& = &{\sum_{K\in{\cal T}_h} \|u-\pi_h u\|_{m,K}^2 }\\ & = & {\sum_{K\in{\cal T}_h} \sum_{l=0}^m |u-\pi_h u|_{l,K}^2} \end{array}\]

Le calcul est donc ramené à un calcul sur chaque élément, pour toutes les semi-normes \(|.|_{l,K},\; l=0,\ldots,m\).

2.3. Etape 3: Passage à l’élément de référence

Soit \(K\) un élément quelconque de \({\cal T}_h\), et \(\hat{K}\) l’élément de référence.

Soit \(F\) la transformation affine de \(\hat{K}\) vers \(K\) :

\[ \forall v\in H^l(K), \qquad |\hat{v}|{l,\hat{K}} \le C(l,n)\; \|B\|^l_2 \; |\hbox{det} B|^{-1/2} \, |v|{l,K} \]

Nous avons le corollaire suivant du théorème [thr1]

On a de même :

\[ \forall v\in H^l(K), \qquad |v|{l,K} \le C(l,n)\; \|B{-1}\|l_2 \; |\mathrm{det} B|^{1/2} \; |\hat{v}|{l,\hat{K}} \]

2.3.1. Preuve du théorème [thr1]

Ce théorème se montre comme suit. Il s’agit là en fait d’un simple résultat de changement de variable dans une intégrale.

Soit \(v\) une fonction \(l\) fois différentiable au point \(x\).

On note \(D^l v(x)\) sa dérivée \(l^{\text{\tiny ème}}\) au sens de Fréchet au point \(x\).

Il s’agit donc d’une forme \(l\)-linéaire symétrique sur \(\RR^n\).

On notera \(D^l v(x).(\xi_1,\ldots,\xi_l)\) sa valeur pour \(l\) vecteurs \(\xi_i\in\RR^n\) (\(1\le i \le l\)).

Reprenons les notations de la section Les espaces \(H^m\).

\(\alpha=(\alpha_1,\ldots,\alpha_n)\) désigne un multi-entier, et on note \(|\alpha|=\alpha_1+\cdots+\alpha_n\). On a alors :

\[|v|_{l,K}^2 = \int_{x\in K} \sum_{|\alpha|=l}\left\|\partial^{|\alpha|}v (x)\right\|^2 dx\]

et

\[\partial^{|\alpha|}v (x) = D^{|\alpha|}v(x).(\underbrace{e_1,\ldots,e_1}_{\alpha_1 \text{{\tiny fois}}}, \ldots , \underbrace{e_n,\ldots,e_n}_{\alpha_n \text{{\tiny fois}}})\]

où \((e_1,\ldots,e_n)\) désigne la base canonique de \(\RR^n\).

Alors, en posant :

\[\|D^l v(x)\| = \sup_{\xi_1,\ldots,\xi_l \in \left(\RR^{\ast}\right)^n} \; \frac{D^l v(x).(\xi_1,\ldots,\xi_l)}{|\xi_1| \ldots |\xi_l|}\qquad ,\]

on déduit qu’il existe des constantes \(\gamma_1\) et \(\gamma_2\) dépendant uniquement de \(n\) et \(l\) (donc en particulier indépendantes de \(v\)) telles que

\[\gamma_1 \; |v|_{l,K} \le \left( \int_{x\in K} \| D^l v(x)\|^2 \, dx\right)^{1/2} \le \gamma_2 \; |v|_{l,K}\]

Par ailleurs, si l’on utilise le changement de variable \(x=F(\hat{x})=B\hat{x}+b\) dans \(D^l v(x)\), il vient :

\forall \xi_1,\ldots,\xi_l \in \RR^n,\qquad D^l \hat{v}(\hat{x}).(\xi_1,\ldots,\xi_l) = D^l v(x).(B\xi_1,\ldots,B\xi_l)

d’où

\|D^l \hat{v}(\hat{x})\| \le \|B\|^l\; \|D^l v(x)\|

Or, \(D^l v(x) = D^l v(F(\hat{x}))\). Donc

\int_{\hat{x}\in \hat{K}}\|D^l \hat{v}(\hat{x})\|^2 \; d\hat{x} \le \|B\|^{2l}\; \int_{\hat{x}\in \hat{K}} \|D^l v(F(\hat{x}))\|^2 \; d\hat{x} = \|B\|^{2l}\; |\hbox{det }B|^{-1} \int_{x\in K} \|D^l v(x)\|^2 \; dx

En minorant et majorant ([eq:maj2]) grâce à la (majoration), on obtient :

\gamma^2_1 \; |\hat{v}|^2_{l,\hat{K}} \le \|B\|^{2l}\; |\hbox{det }B|^{-1} \gamma^2_2 \; |v|^2_{l,K}

d’où le résultat de (majoration du théorème) ce qui conclut la preuve de [thr1] \(\blacksquare\)

2.3.2. Estimation de \(\|B\|\)

Soit \(h_K\) le diamètre de \(K\), c’est à dire le maximum des distances euclidiennes entre deux points de \(K\).

Soit \(\rho_K\) la rondeur de \(K\), c’est à dire le diamètre maximum des sphères incluses dans \(K\).

On a :

\|B\| = \sup_{x\ne 0} \frac{\|Bx\|}{\|x\|} = \sup_{\|x\|=\hat{\rho}} \frac{\|Bx\|}{\hat{\rho}}

Soit \(x\) un vecteur de \(\RR^n\) tel que \(\|x\|=\hat{\rho}\).

Par définition de \(\hat{\rho}\), il existe deux points \(\hat{y}\) et \(\hat{z}\) de \(\hat{K}\) tels que \(x=\hat{y}-\hat{z}\).

Alors \(Bx=B\hat{y}-B\hat{z}=F(\hat{y})-F(\hat{z})=y-z\) avec \(y\) et \(z\) appartenant à \(K\).

Par définition de \(h_K\), \(\|y-z\| \le h_K\). Donc \(\|Bx\| \le h_K\).

En reportant dans la définition de \(\|B\|\), on obtient donc :

\|B\| \le \frac{h_K}{\hat{\rho}}

Et on a évidemment de même :

\|B^{-1}\| \le \frac{\hat{h}}{\rho_K}

2.4. Etape 4: Majoration sur l’élément de référence

Le résultat principal est le suivant :

Soient \(l\) et \(k\) deux entiers tels que \(0\le l \le k+1\). Si \(\hat{\pi} \in {\cal L}(H^{k+1}(\hat{K}),H^l(\hat{K}))\) laisse \(P_k(\hat{K})\) invariant (c’est à dire vérifie \(\forall \hat{p}\in P_k(\hat{K}), \hat{\pi}\hat{p}=\hat{p}\)), alors

\[ \exists C(\hat{K},\hat{\pi}) ,\; \forall \hat{v} \in H^{k+1}(\hat{K}), \; |\hat{v}-\hat{\pi}\hat{v}|{l,\hat{K}} \le C |\hat{v}|{k+1,\hat{K}} \]

2.4.1. Preuve du théorème [thr2]

On montre ce résultat comme suit:

\(\hat{\pi} \in {\cal L}(H^{k+1}(\hat{K}),H^l(\hat{K}))\), et donc \(I-\hat{\pi} \in {\cal L}(H^{k+1}(\hat{K}),H^l(\hat{K}))\) car \(l\le k+1\).

Et donc

\[|\hat{v}-\hat{\pi}\hat{v}|_{l,\hat{K}} \le \|I-\hat{\pi}\|_{\mathcal{L}(H^{k+1}(\hat{K}),H^l(\hat{K}))}\; \|\hat{v}\|_{k+1,\hat{K}}\]

On utilise maintenant l’invariance de \(P_k(\hat{K})\):

On aura donc démontré le théorème si l’on montre que

\[\exists C,\; \forall \hat{v}\in H^{k+1}(\hat{K}) \; \inf_{\hat{p}\in P_k(\hat{K})} \|\hat{v}+\hat{p}\|_{k+1,\hat{K}} \le C |\hat{v}|_{k+1,\hat{K}}\]

Soit \((f_i)_{i=0,\ldots,k}\) une base du dual de \(P_k(\hat{K})\).

D’après le théorème d’Hahn-Banach, il existe des formes linéaires continues sur \(H^{k+1}(\hat{K})\), que l’on notera encore \(f_i\), et qui prolongent les \(f_i\).

En particulier, si \(\hat{p}\in P_k(\hat{K})\) vérifie \(f_i(\hat{p})=0,\, (i=0,\ldots,k)\), alors \(\hat{p}=0\).

Nous allons montrer que

\[\exists C, \, \forall \hat{v}\in H^{k+1}(\hat{K}), \; \|\hat{v}\|_{k+1,\hat{K}} \le C \left\{ |\hat{v}|_{k+1,\hat{K}} + \sum_{i=0}^k |f_i(\hat{v})| \right\}\]
On aura le résultat souhaité en appliquant ([eqref2]) à \(\hat{v}+\hat{q}\), avec \(\hat{q}\) tel que \(f_i(\hat{q})=f_i(-\hat{v})\).

On montre la relation ([eqref2]) par l’absurde comme suit:

Si [eqref2] n’est pas vraie, alors il existe une suite de fonctions \(\hat{v}_n\) de \(H^{k+1}(\hat{K})\) telles que :

\[ \|\hat{v}_n\|_{k+1,\hat{K}} =1, \;\; |\hat{v}_n|_{k+1,\hat{K}} \longrightarrow 0,\; \hbox{ et } \forall i \; f_i(\hat{v}_n)\longrightarrow 0\]

Par complétude de \(H^{k+1}(\hat{K})\), on extrait une sous-suite convergente vers \(\hat{v} \in H^{k+1}(\hat{K})\).

Mais \(|\hat{v}_n|_{k+1,\hat{K}} \longrightarrow 0\).

Donc \(\hat{v} \in P_k(\hat{K})\) et \(f_i(\hat{v})=0\).

D’où une contradiction. \(\blacksquare\).

2.5. Etape 5: Assemblage des majorations locales

2.5.1. Majoration sur un élément quelconque

En rassemblant les résultats précédents, on peut établir une majoration sur un élément quelconque \(K\) du maillage.

On a :

\[\begin{array}{rclr} |v-\pi_K v|_{l,K} & \le & C(l,n)\; \|B^{-1}\|^l\; |\hbox{det }B|^{1/2} \; |\hat{v}-\hat{\pi}\hat{v}|_{l,\hat{K}}&\hbox{d'après (<<eq:majref2>>)} \\ & \le & C(l,n)\; \|B^{-1}\|^l\; |\hbox{det }B|^{1/2} \; C(\hat{K},\hat{\pi})\; |\hat{v}|_{k+1,\hat{K}} &\hbox{d'après (<<eq:majref0>>)}\\ & \le & C(l,n)\; \|B^{-1}\|^l\; |\hbox{det }B|^{1/2} \; C(\hat{K},\hat{\pi})\; C(k+1,n) \; \|B\|^{k+1} |\hbox{det }B|^{-1/2}\; |v|_{k+1,K} & \hbox{d'après (<<eq:majref>>)}\\ & \le & C(l,n)\; \frac{\hat{h}^l}{\rho_K^l} \; \; C(\hat{K},\hat{\pi})\; C(k+1,n) \; \frac{h_K^{k+1}}{\hat{\rho}^{k+1}} \; |v|_{k+1,K} & \hbox{d'après (<<eq:kk1>>) et (<<eq:kk2>>)}\\ \end{array}\]

D’où finalement :

\[|v-\pi_K v|_{l,K} \le \hat{C}(\hat{\pi},\hat{K},l,k,n)\; \frac{h_K^{k+1}}{\rho_K^l} \; |v|_{k+1,K}\]
Il est important de remarquer à ce niveau que \(\hat{C}\) est indépendant de \(K\).

2.5.2. Assemblage des résultats locaux

On va maintenant reprendre la majoration ([eqmajloc]) pour tous les éléments du maillage et toutes les valeurs de \(l=0,\ldots,m\).

On va définir deux quantités représentatives du maillage :

  • \(h\quad\) tel que \(h_K \le h, \; \forall K\in {\cal T}_h\qquad\) (diamètre maximum des éléments)

  • \(\sigma\quad\) tel que \({\frac{h_K}{\rho_K}} \le \sigma, \; \forall K\in {\cal T}_h\qquad\) (caractérise l’aplatissement des éléments)

On a

\[\begin{array}{rcl} \|v-\pi_K v\|^2_{m,K} & = & \sum_{l=0}^m |v-\pi_K v|^2_{l,K} \\ & \le & \sum_{l=0}^m \hat{C}^2(\hat{\pi},\hat{K},l,k,n)\; \left(\frac{h_K^{k+1}}{\rho_K^l}\right)^2 \; |v|^2_{k+1,K}\qquad\hbox{d'apr\`es (\ref{eqmajloc})}\\ & \le & \sum_{l=0}^m \hat{C}^2(\hat{\pi},\hat{K},l,k,n)\; \left\{\left(\frac{h_K}{\rho_K}\right)^l\; h_K^{m-l}\; h_K^{k+1-m}\right\}^2 \; |v|^2_{k+1,K}\\ & \le & \left\{ \sum_{l=0}^m \hat{C}^2(\hat{\pi},\hat{K},l,k,n)\; \sigma^{2l} h^{2m-2l} \right\} \; \left[ h^{k+1-m}\; |v|_{k+1,K} \right]^2 \end{array}\]

Le terme entre accolades ne tend ni vers 0 ni vers l’infini quand \(h\) tend vers 0.

D’où :

\[\|v-\pi_K v\|_{m,K} \le \hat{C}'(\hat{\pi},\hat{K},l,k,n,\sigma,h)\; h^{k+1-m}\; |v|_{k+1,K}\]

En sommant ensuite sur tous les éléments du maillage :

\[\begin{array}{rcl} \|v-\pi_h v\|^2_{m} & = & \sum_{K\in {\cal T}_h} \|v-\pi_K v\|^2_{m,K} \\ & \le & \sum_{K\in {\cal T}_h} \left[ \hat{C}'(\hat{\pi},\hat{K},l,k,n,\sigma,h)\; h^{k+1-m} \; |v|_{k+1,K} \right]^2 \end{array}\]

On obtient finalement :

\[\|v-\pi_h v\|_{m} \le C(\mathcal{T}_h,m,k,n) \; h^{k+1-m}\; |v|_{k+1}\]

2.6. Résultat final

En reportant ([eqmajfinal]) dans ([eq:cea2]), on obtient le résultat final classique de majoration d’erreur :

\[\|u-u_h\|_{m} \le \mathcal{C} \; h^{k+1-m}\; |u|_{k+1}\]

3. Quelques commentaires

Une utilisation fréquente de ([eqmajfinal2]) a lieu dans le cas \(m=1\). Alors si l’espace de polynômes \(P_k(\hat{K})\subset H^1(\hat{K})\) (ce qui est toujours le cas) et si \(\hat{\pi}\) est bien défini sur \(H^{k+1}(\hat{K})\), on a :

\[\mbox{si }u\in H^{k+1}(\Omega),\quad \|u-u_h\|_1 \le \mathcal{C} \; h^k \; |u|_{k+1}\]
Si le domaine \(\Omega\) n’est pas polygonal, la majoration précédente n’est plus valable. On peut alors établir d’autres majorations du même type – se référer par exemple à (raviartthomast1983). De même, si les calculs d’intégrales ne sont pas faits exactement mais à l’aide d’une intégration numérique, une erreur supplémentaire doit être prise en compte, qui conduit à une nouvelle majoration d’erreur – voir là-aussi par exemple (raviartthomast1983).

Approximation de problèmes coercifs

Dans ce chapitre, on s’intéresse à l’approximation de problèmes coercifs:

1. Le Laplacien

On s’instéresse dans cette section à l’approximation élément fini conforme du problème suivant:

Formulation forte du Laplacien

On cherche \(u\) telle que

\[ \begin{split} -\Delta u &= f \mbox{ dans } \Omega\\ u &= 0 \mbox{ sur } \partial \Omega \end{split} \]

1.1. Cadre Mathematique

On suppose que \(f \in L^2(\Omega)\).

La formulation faible du problème Formulation forte du Laplacien est la suivante:

Formulation faible pour des conditions de Dirichlet homogènes

On cherche \(u \in H^1_0(\Omega)\) telle que \[\label{eq:65} \int_\Omega \nabla u \cdot \nabla v = \int_\Omega f\, v,\quad \forall v \in H^1_0(\Omega) \]

1.2. Problème bien posé

Introduisons

  • \(V = H^1_0(\Omega)\) doté de la norme \(\|\cdot\|_{1,\Omega}\) telle que \(\|v\|_{1,\Omega} = (\|v\|^2_{0,\Omega} + \|\nabla v\|^2_{0,\Omega})^{1/2}\)

  • la forme bilinéaire \(a \in \mathcal{L}(V \times V, \R)\) telle que \(a(u,v) = \int_\Omega \nabla u \cdot \nabla v \)

  • la forme linéaire \(\ell \in \mathcal{L}(V, \R)\) telle que \(l(v) = \int_\Omega f \nabla v \)

Le problème Formulation faible pour des conditions de Dirichlet homogènes s’écrit sous forme abstraite

On cherche \(u \in V\) telle que \[ a(u, v) = \ell(v), \quad \forall v \in V \]

L’espace \(V\) est un espace de Hilbert et les formes \(a\) et \(\ell\) sont continues sur \(V\times V\) et \(V\) respectivement.

Il ne reste plus qu’à vérifier si le problème est bien posé (existence d’une solution unique).

Pour cela on utilise démontre la coercivité de la forme bilinéaire \(a\) sur \(V \equiv H^1_0(\Omega)\).

Ceci se fait grâce au lemme suivant:

Soit \(\Omega\) un ouvert borné de \(\R\). Il existe une constante \(c_\Omega\) (dépendente de \(\Omega\) telle que \[ \forall v \in H^1_0(\Omega),\quad \|v\|{0,\Omega} \le c\Omega \|\nabla v\|_{0,\Omega} \]

\(c_\Omega\) est homogène à une longeur et peut être interprétée comme une mesure caractéristique de \(\Omega\).

Grâce à l’inégalité de Poincaré, on a le résultat suivant

La forme bilinéaire \(a\) du problème Formulation faible pour des conditions de Dirichlet homogènes est coercive

On note tout d’abord que par l’inégalité de Poincaré et la définition de \(\|\cdot\|_{1,\Omega}\)

\[ \|v\|^2_{1,\Omega} \le (1 + c^2_\Omega) \|\nabla v\|^2_{0,\Omega} \] On en déduit que

\[ \forall v \in H^1_0(\Omega),\quad a(v,v) = \|\nabla v\|^2_{0,\Omega} \ge \frac{1}{1+c^2_\Omega} \|v\|^2_{1,\Omega} \]

Le Lemme de Lax-Milgram Lax-Milgram permet alors de conclure sur l’existence d’une solution unique pour le problème Formulation faible pour des conditions de Dirichlet homogènes.

1.3. Approximation conforme

On utilise une approximation conforme par éléments finis de Lagrange.

On considère \(\Omega\) un polygone ou polyhèdre régulier de \(\RR^2\) ou \(\RR^3\) respectivement et un maillage \(\calTh = \{K_e\}_{e=1...\Ne}\) de \(\Omega\).

On considère un élément fini de référence \((\hat{K},\hat{P},\hat{\Sigma})\) tel que \(\Pk \subset \hat{P}\) et \(k+1 > \frac{d}{2}\), voir le théorème [thr:16].

On note

\[ L^k_{c,h} = \{ v_h \in C^0(\bar{\Omega}); \forall K \in \mathcal{T}_h,\ v_h \circ T_K \in \hat{P}\} \] où \(T_K\) est la tranformation géométrique de \(\hat{K}\) dans \(K\).

  • Si on utilise \(\hat{P}=\Pk{k}\) on a \(L^k_{c,h} = P^k_{c,h}\).

  • Si on utilise \(\hat{P}=\Qk{k}\) on a \(L^k_{c,h} = Q^k_{c,h}\).

Afin de construire un espace d’approximation conforme \(V_h \subset V =H^1_0(\Omega)\) on prend \[ \label{eq:71} V_h = L^k_{c,h} \cap H^1_0(\Omega) \] c’est à dire que les fonctions de \(V_h\) satisfont les conditions aux limites outre le fait d’être dans \(L^k_{c,h}\).

Le problème discret s’écrit alors

Trouver \(u_h \in V_h\) telle que \[ a(u_h,v_h) = \ell(v_h),\, \forall v_h \in V_h \]

qui est bien posé (existence et unicité de \(u_h\)) car \(a\) est coercive sur \(V\) et que \(V_h \subset V\).

On a le résultat suivant:

On suppose que \(u\) solution de Formulation faible pour des conditions de Dirichlet homogènes est dans \(H^{k+1}(\Omega) \cap H^1_0(\Omega)\), alors il existe une constante \(c_1\) telle que pour tout \(h\)

\[ \|u-u_v\|{1,\Omega} \le c_1 h^k |u|{k+1,\Omega} \] et il existe une constante \(c_2\) telle que pour tout \(h\)

\[ \|u-u_v\|{0,\Omega} \le c_2 h^{k+1} |u|{k+1,\Omega} \]

La preuve de ([eq:73]) est obtenue grâce au lemme de Cea ([eq:cea]) et du théorème d’interpolation ([eq:62]).

La preuve de ([eq:74]) est obtenue grâce au lemme d’Aubin-Nitsche qui permet d’affirmer qu’il existe une constante \(c\) telle que

\[ \|u-u_h\|{0,\Omega} \leq c h |u-u_h|{1,\Omega} \] et donc que ([eq:74]) se déduit de ([eq:74]).

1.4. Implémentation avec Feel++

Avec Feel++, on ne construit pas explicitement l’espace \(V_h\) mais \(L^k_{c,h}\).

Le traitement des conditions aux limites de Dirichlet du problème ([eq:64]) peut être effectué de diverses façons, nous en verrons une.

Commencons par le maillage, dans un premier temps nous définissons le type du maillage contenant soit des éléments de type simplexe (segment,triangle, tetrahèdre) ou de type hypercube (segment, quadrangle, hexahèdre).

  // un maillage de simplexe dans $\R$ telle que la transformation
  // géométrique $T_K,\ K \in \calTh$  soit affine
  typedef Mesh<Simplex<d,1> > mesh_type;

  // un maillage d'hypercube dans $\R$ telle que la transformation
  // géométrique $T_K,\ K \in \calTh$  soit affine en chacune des variables
  // typedef Mesh<Hypercube<d,1> > mesh\_type;

  // generate the mesh associated to the unit square $[0,1]^2$ using triangles
  auto mesh = unitSquare();
Le mot clé auto permet de faire de l’inférence de type, pour plus de détails consultez la page C++11 de Wikipedia.

Ensuite nous pouvons définir l’espace \(L^k_{c,h}\),

  // Vh est une structure de donnée allouée dynamiquement
  auto Vh = Pch<1>( mesh );
  // u est un élément de Vh
  auto u = Vh->element();
  // u est un autre élément de Vh
  auto u = Vh->element();

À présent, nous définissons les formes bilinéaires \(a\) et \(\ell\) qui sont respectivement des formes bilinéaires et linéaires.

  auto a = form2( _test=Vh, _trial=Vh ); (1)
  a = integrate( _range=elements(mesh), _expr=gradt(u)*trans(grad(v)) ); (2)

  auto l = form1( _test=Vh ); (3)
  l = integrate( _range=elements(mesh), _expr=f*id(v) ); (4)
1 \(a \in \mathcal{L}(V_h \times V_h,\ \RR)\)
2 \(a = \sum_{e=1...\Ne} \int_{K_e} \nabla \varphi_j \cdot \nabla \varphi_i,\quad i,j=1...,\dim{V_h}\)
3 \(\ell \in \mathcal{L}(V_h,\ \RR)\)
4 \(\ell = \sum_{e=1...\Ne} \int_{K_e} f \varphi_i,\quad i=1...,\dim{V_h}\)

Afin de traiter les conditions aux limites de Dirichlet homogènes, on peut utiliser le mot-clé on qui permet de les imposer de manière forte.

  a += on(_range=boundaryfaces(mesh), _element=u, _rhs=l, _expr=constant(0.) );
Le mot-clé constant permet de transformer une type numérique ( entier, flottant) en une expression utilisable par le langage de Feel++. Notez également l’opération += qui permet de rajouter le traitement des conditions de Dirichlet tout en gardant les contributions précédentes. L’opération = aurait d’abord remis à \(0\) les entrées de la matrice associée à \(a\).

Enfin nous pouvons résoudre le problème [prob:6]

  a.solve( _rhs=l, _solution=u );

Le listing complet

1.5. Conditions aux limites

1.5.1. Conditions aux limites de Dirichlet non homogène

On suppose toujours \(f \in L^2(\Omega)\) et on se donne une fonction \(g \in C^{0,1}(\partial \Omega)\)

[\(g\) est Lipschitzienne sur \(\partial \Omega\)].

On s’intéresse au problème suivant:

On cherche \(u : \Omega \rightarrow \RR\) telle que

\[ \begin{split} -\Delta u &= f \mbox{ dans } \Omega\\ u &= g \mbox{ sur } \partial \Omega \end{split} \]

L’hypothèse \(g \in C^{0,1}(\partial \Omega)\) permet d’affirmer qu’il existe \(u_g \in H^1(\Omega)\) telle que \(u_{g_{|\partial \Omega}} = g\).

On se ramène au problème avec conditions de Dirichlet homogène en faisant le change d’inconnue \(u_0=u-u_g\) et on s’intéresse au problème suivant:

On cherche \(u_0 \in H^1_0(\Omega)\) telle que

\[ a(u_0,v) = \ell(v) - a(u_g,v),\quad \forall v \in H^1_0(\Omega) \]

Ce problème est bien posé d’après Lax-Milgram, voir section précédente.

On suppose que \(u\) solution de [prob:8] est dans \(H^{k+1}(\Omega) \cap H^1_0(\Omega)\), alors il existe une constante \(c_1\) telle que pour tout \(h\)

\[ \|u-u_v\|{1,\Omega} \le c_1 h^k |u|{k+1,\Omega} \] et il existe une constante \(c_2\) telle que pour tout \(h\)

\[ \|u-u_v\|{0,\Omega} \le c_2 h^{k+1} |u|{k+1,\Omega} \]

Avec Feel++, les conditions Dirichlet non-homogènes sont traitées par exemple avec le mot-clé on.

Conditions de Dirichlet non homogènes
  auto g = sin(2*pi*Px() ); (1)
  (2)
  ...
  a += on( _range=boundaryfaces(mesh), _expr=g ); (3)
1 définition de la fonction, p.ex \(g=sin(2 \pi x)\)
2 définition de \(a\)
3 ajout des conditions de Dirichlet non-homogènes
Il n’y a pas besoin de rajouter le terme \(a(u_g,v)\) au second membre \(\ell(v)\), il est pris en compte automatiquement par on.

Voici le listing complet de l’exemple du laplacien avec conditions de Dirichlet non-homogène

1.6. Condition aux limites de Neumann

Étant donnés un réel \(\mu\) strictement positif, \(f \in L^2(\Omega)\) et \(g \in L^2(\partial \Omega)\), on s’intéresse au problème suivant:

On cherche \(u : \Omega \rightarrow \RR\) telle que

\[ \begin{split} -\Delta u + \mu u &= f, \mbox{ dans } \Omega\\ \partial_\Next u &= g, \mbox{ sur } \partial\Omega \end{split} \]

où \(\partial_\Next u = \nabla u \cdot \Next = \sum_{i=1}^d n_i \partial_i u\) dénote la dérivée normale de \(u\) avec \(\Next=(n_1,...,n_d) \in \RR^d\) la normale extérieure unitaire en un point du bord de \(\Omega\).

La formulation faible s’écrit

On cherche \(u \in H^1(\Omega)\) telle que

\[ a( u, v ) = \ell(v),\ \forall v \in H^1(\Omega) \] avec

\[ a( u, v ) = \int_\Omega \nabla u \cdot \nabla v + \mu u v \] et

\[ \ell( v ) = \int_\Omega f v + \int_{\partial\Omega} g v \]

On a

\[ a(v, v) = \int_\Omega \nabla v \cdot \nabla v + \mu v v \ge \min(1,\mu) \int_\Omega \nabla v \cdot \nabla v + v v = \min(1,\mu) \|v\|_{1,\Omega}, \quad \forall v \in H^1(\Omega) \] ce qui nous permet d’affirmer que \(a\) est coercive sur \(H^1(\Omega)\) et que le problème [prob:13] est bien posé grâce à Lax-Milgram.

On peut approcher le problème [prob:13] par des éléments finis de Lagrange.

On utilise la formulation développée dans la section Approximation conforme

On cherche \(u_h \in P^k_{c,h}\) tel que \[ a_h(u_h,w_h)=\ell(w_h), \quad \forall w_h \in P^k_{c,h} \]

Par rapport aux conditions de Dirichet, les conditions de Neumann sont directement (naturellement) traitées par la formulation faible. Elles ne sont pas directement imposées dans l’espace et les fonctions tests peuvent prendre des valeurs non-nulles au bord. Ces conditions sont traitées de manière approchée et non pas exacte.

L’analyse du problème [prob:13_1] est identique par le théorème [thr:17] aux mêmes estimations que dans le cas Dirichlet homogène.

1.6.1. Cas \(\lambda=0\)

Le cas \(\lambda=0\) présente quelques difficultés techniques.

On a

On cherche \(u : \Omega \rightarrow \RR\) telle que

\[ \begin{split} -\Delta u &= f, \mbox{ dans } \Omega\\ \partial_\Next u &= g, \mbox{ sur } \partial\Omega \end{split} \]

On observe ici qu’une condition nécessaire d’existence de solution est que
\[\int_\Omega f + \int_{\partial \Omega} g = 0\]

de par le théorème de la divergence.

la deuxième observation est que la solution de [prob:13_2] est connue à une constante additive près: si \(u\) est solution et \(c\in \RR\) alors \(u+c\) est également solution.

Il convient alors de chercher la solution dans l’espace fonctionel suivant

H^1_*(\Omega) = \{ v \in H^1(\Omega):\quad \int_\Omega v = 0 \}

La formulation faible du problème [prob:13_2] s’écrit alors

On cherche \(u\in H^1_*(\Omega)\), telle que \[ a(u,v) = \ell(v) \quad \forall v\in H^1_*(\Omega) \] avec \(a(u,v)=\int_\Omega \nabla u \cdot \nabla v\).

Le caractère bien posé du problème [prob:13_3] résulte de la coercivité de \(a\) sur \(H^1_*(\Omega)\) qui elle-même résulte du Lemme suivante

Lemme d’inégalité de Poincaré-Wirtinger

Soit \(\Omega\) un ouvert borné de \(\RR^d\), il existe une constante \(C_\Omega\) telle que \[ \|v\| \leq C_\Omega \|\nabla v\|{0,\Omega},\quad \forall v \in H^1*(\Omega) \]

Le problème [prob:13_3] peut être approché par des éléments finis de Lagrange ce qui conduit au problème discret suivant

On cherche \(u_h \in P^k_{c,h}\) telle que a_h(u_h,v_h)=\ell(v_h) \quad \forall v_h\in P^k_{c,h}

L’espace d’approximation \(P^k{c,h}\) n’est pas conforme dans \(H^1_*(\Omega)\)

le problème [prob:13_4] est singulier, il a une infinité de solution. L’une d’entre elle peut être approchée par une méthode itérative de type gradient conjugué. Notons \(u^*_h\) cette solution alors la solution à moyenne nulle est

\[u_h=u^*_h-\frac{1}{|\Omega|}\int_\Omega u^*_h\]

Il s’agit donc d’effectuer a posteriori du calcul un post-processing pour se ramener à la solution à moyenne nulle. Dans Feel++ on utilisera la fonction mean

// ...
// u est la solution du problème discret
a.solve(_rhs=l,_solution=u);

// calcul de la valeur moyenne de u
double meanu = mean(_range=elements(mesh),_expr=idv(u))(0,0);

// calcul de la solution à moyenne nulle u = u - meanu
u.add(-meanu);

1.7. Conditions aux limites de Robin

Étant donnés un réel \(\mu\) strictement positif, \(f \in L^2(\Omega)\) et \(g \in L^2(\partial \Omega)\), on s’intéresse au problème suivant:

On cherche \(u : \Omega \rightarrow \RR\) telle que

\[ \begin{split} -\Delta u &= f, \mbox{ dans } \Omega\\ \mu u + \partial_\Next u &= g, \mbox{ sur } \partial\Omega. \end{split} \]

La formulation faible s’écrit

On cherche \(u \in H^1(\Omega)\) telle que \[\label{eq:84} a( u, v ) = \ell(v),\ \forall v \in H^1(\Omega) \] avec

\[ a( u, v ) = \int_\Omega \nabla u \cdot \nabla v + \int_{\partial \Omega} \mu u v \] et

\[ \ell( v ) = \int_\Omega f v + \int_{\partial\Omega} g v \]

On a

\[ \begin{split} a(v, v) & = \int_\Omega \nabla v \cdot \nabla v + \int_{\partial\Omega} \mu v v \\ & \geq \min(1,\mu)\left( \int_\Omega \nabla v \cdot \nabla v
\int_{\partial\Omega} v v\right) \\ &\geq \min(1,\mu) \|v\|_{1,\Omega}, \quad \forall v \in H^1(\Omega) \end{split} \]

La forme bilinéaire \(a\) est donc coercive et le problème [prob:15] est bien posé grâce à Lax-Milgram.

On considère le problème discret suivant

On cherche \(u_h \in P^k_{c,h}(\Omega)\) telle que

\[ a_h(u_h,v_h) = \ell(v_h)\quad \forall v_h \in P^k_{c,h} \]

Le problème est bien posé (\(P^k_{c,h} \subset H^1(\Omega)\)).

Comme pour le problème avec conditions de Neumann, les fonctions tests peuvent prendre des valeurs non nulles au bord. Les conditions de Robin(ou Fourier) ne sont satisfaites que de manière approchée et non pas exacte.

La convergence de \(u_h\) est donnée par le théorème [thr:17].

Considérons \(\Omega=[0,1\)^2] et les données \(\mu=0.01\), \(f=1\) et \(g=0\).

2. Advection-diffusion-réaction avec diffusion dominante

On s’intéresse au problème suivant:

On cherche \(u : \Omega \rightarrow \RR\) telle que

\[ \begin{split} -\nabla \cdot ( \mathbf{\alpha} \nabla u ) + \mathbf{\beta} \cdot \nabla u + \mu u &= f \mbox{ dans } \Omega\\ \mathcal{B} u &= 0 \mbox{ sur } \partial \Omega\\ \end{split} \]

où \(\mathcal{B}\) est l’opérateur prenant en compte les conditions aux limites. La variation sur les conditions aux limites de la section [sec:cond-aux-limit] s’appliquent.

  • \(\mathcal{B}u=u\) condition de Dirichlet

  • \(\mathcal{B}u=(\mathbf{\alpha}\nabla{u}) \cdot n\) condition de Neumann

  • \(\mathcal{B}u=\gamma u + (\mathbf{\alpha} \nabla u )\cdot n\) condition de Robin

Les différents termes de l’opérateur ont une interprétation physique

  • \(-\nabla \cdot ( \mathbf{\alpha} \nabla u )\) est un terme de diffusion,

  • \(\mathbf{\beta} \cdot \nabla u\) est un terme d’advection,

  • \(\mu u\) est un terme de réaction (production ou destruction).

Ce type d’équation est très fréquente en ingéniérie, biologie ou encore finance.

Transfert de chaleur

\(u\) is the temperature, \(\alpha=\kappa \mathcal{I}\) où \(\kappa\) est la conductivité thermique, \(\beta\) est le champ d’écoulement, \(\mu = 0\) et \(f\) est chaleur apportée exterieurement par unité de volume.

Advection Diffusion

\(u\) est la concentration d’un soluté transportée dans un champ d’écoulement \(\beta\). La matrice \(\alpha\) modélise la diffusivité du soluté soit de la diffusion moléculaire soit du mélange turbulent du fluide transporteur. La production ou destruction par réaction chimique is pris en compte par le terme linéaire \(\mu u\). Le second membre \(f\) modélise les sources ou puits.

On suppose que \(\mathbf{\alpha} \in [L^{\infty}(\Omega)]^{d\times d}\), \(\mathbf{\beta} \in [L^{\infty}(\Omega)]^d\), \(\nabla \cdot \mathbf{\beta} \in L^{\infty}(\Omega)\) et \(\mu \in L^\infty(\Omega)\).

On suppose que l’opérateur \(\mathcal{L}\) tel que \(\mathcal{L} u = -\nabla \cdot ( \alpha \nabla u ) + \beta \cdot \nabla u + \mu u\) est elliptique au sens suivant:

L’opérateur \(\mathcal{L}\) est elliptique si il existe une constante \(\alpha_0\) telle que presque pour tout \(x\in\Omega\)

\[ \forall \xi=(\xi_1, \ldots , \xi_n)\in\RR^n,\quad { \sum_{i,j=1}^n \alpha_{ij}(x) \, \xi_i \, \xi_j \ge \alpha_0 \, \| \xi \|^2 } \]

Le Laplacien est dans la catégorie des opérateurs elliptiques, il correspond à \(\mathbf{\beta} = 0\), \(\mu = 0\) et \(\mathbf{\alpha}=\mathcal{I}_d\) avec \(\mathcal{I}_d\) la matrice identité de \(\RR^{d\times d}\).

2.1. Condition de Dirichlet homogène

Nous imposons \(u=0\) sur \(\partial \Omega\).

Nous multiplions \(\mathcal{L} u = f\) par une fonction test (suffisamment régulière) s’annulant au bord et on intègre sur \(\Omega\). En intégrant par parties (formule de Green) nous avons

\[\int_\Omega \nabla \cdot (\mathbf{\alpha} \nabla u ) v = \int_\Omega (\mathbf{\alpha} \nabla u ) \cdot \nabla v -\int_{\partial \Omega} ((\sigma \nabla u)\cdot n) v \nabla u) v + \mu u v\]

ce qui donne

\[\int_\Omega (\mathbf{\alpha} \nabla u ) \cdot \nabla v + (\mathbf{\beta} \cdot \nabla u) v + \mu u v = \int_\Omega f v\]

Afin que les intégrales sur \(\Omega\) aient un sens, nous demandons par exemple que \(u\) et \(v\) aient la régularité suivante:

\[u\in H^1(\Omega)\quad\text{et}\quad v\in H^1(\Omega).\]

Comme \(u\in H^1(\Omega)\), le théorème [toto] implique que \(u\) ait une trace sur le bord. Comme \(u_{|\partial \Omega} = 0\), on cherche en fait \(u\) dans \(H^1_0(\Omega)\). Les fonctions tests sont également prises dans \(H^1_0(\Omega)\) ce qui donne la formulation faible suivante

On cherche \(u \in H^1_0(\Omega)\) telle que

\[ a(u,v) = \ell(v) \quad \forall v \in H^1_0(\Omega) \] avec

\[ a(u,v)=\int_\Omega (\mathbf{\alpha} \nabla u ) \cdot \nabla v + (\mathbf{\beta} \cdot \nabla u) v + \mu u v \] et

\[ \ell(v) = \int_\Omega f v \]

Nous avons le résultat suivant .Une solution du problème faible est solution du problème fort

Si \(u\) résout [eq:89], alors \(\mathcal{L}u=f\) presque partout dans \(\Omega\) et \(u=0\) presque partout sur \(\partial \Omega\).

Soit \(\phi\in \mathcal{D}(\Omega)\) et \(u\) une solution de [prob:12]. On a \[ \begin{split} \left\langle-\nabla \cdot (\alpha \nabla u), \phi \right\rangle_{\mathcal{D}',\mathcal{D}} &= \left\langle (\alpha \nabla u), \nabla \phi \right\rangle_{\mathcal{D}',\mathcal{D}} = \int_\Omega \alpha \nabla u \cdot \nabla \phi \\ & = \int_\Omega (f - \beta \cdot \nabla u - \mu u ) \phi, \end{split} \] ce qui donne \(\langle\mathcal{L}u,\phi\rangle_{\mathcal{D}',\mathcal{D}} = \int_\Omega f \phi\). Or \(\mathcal{D}(\Omega)\) est dense dans \(L^2(\Omega)\), nous avons donc que \(\mathcal{L}u=f\) dans \(L^2(\Omega)\). C’est à dire que \(\mathcal{L}u=f\) presque partout dans \(\Omega\). Enfin \(u=0\) presque partout sur \(\partial \Omega\) par définition de \(H^1_0(\Omega)\) d’après le théorème [toto] \(\blacksquare\)

2.2. Conditions de Dirichlet Non homogène

Nous imposons \(u=0\) sur \(\partial \Omega\) où \(g : \partial \Omega \rightarrow \RR\) est une fonction donnée.

Nous supposons que \(g\) soit suffisamment régulière telle qu’il existe un relèvement \(u_g\) de \(g\) dans \(H^1(\Omega)\), c’est à dire qu’il existe une fonction \(u_g\) de \(H^1(\Omega)\) telle que \(u_g=g\) sur \(\partial \Omega\).

On cherche \(u \in H^1(\Omega)\) telle que \[ u=u_g+\phi,\quad \phi \in H^1_0(\Omega) \] où \(\phi\) est solution de

\[ a(\phi,v) = \ell(v) - a(u_g,v) \quad \forall v \in H^1_0(\Omega) \] avec

\[ a(u,v)=\int_\Omega (\mathbf{\alpha} \nabla u ) \cdot \nabla v + (\mathbf{\beta} \cdot \nabla u) v + \mu u v \] et

\[ \ell(v) = \int_\Omega f v \]

Soit \(g\in H^{\frac{1}{2}}(\partial \Omega)\), si \(u\) résout [prob:12_1], alors \(\mathcal{L}u=f\) presque partout dans \(\Omega\) et \(u=g\) presque partout sur \(\partial \Omega\).

La preuve est similaire à celle de [prob:12].

Lorsque l’opérateur \(\mathcal{L}\) est le Laplacien, [prob:12_1] est appelé un problème de Poisson.

2.3. Conditions de Neumann

Soit une fonction \(g:\partial \Omega \rightarrow \RR\), nous voulons imposer \(\mathcal{B}u=(\alpha \nabla u)\cdot n = g\) sur \(\partial \Omega\).

Dans le cas où \(\alpha = \mathcal{I}\), la condition de Neumann correspond à spécifier la dérivée normale de \(u\) car \(\nabla u \cdot n = \partial_n u = \frac{\partial u}{\partial n}\).

Nous procédons de la même façon que précédemment et en utilisant la condition de Neumann dans l’intégrale de surface [eq:adr-green], on obtient la formulation faible suivante:

On cherche \(u \in H^1(\Omega)\) telle que

\[ a(u,v) = \ell(v) \quad \forall v \in H^1(\Omega) \] avec

\[ a(u,v)=\int_\Omega (\mathbf{\alpha} \nabla u ) \cdot \nabla v + (\mathbf{\beta} \cdot \nabla u) v + \mu u v \] et

\[ \ell(v) = \int_\Omega f v + \int_{\partial \Omega} g v \]

Soit \(g\in L^2(\partial \Omega)\), si \(u\) résout [prob:12_2], alors \(\mathcal{L}u=f\) presque partout dans \(\Omega\) et \(\alpha \nabla u \cdot n = g\) presque partout sur \(\partial \Omega\).

En prenant les fonctions tests dans \(\mathcal{D}(\Omega)\) , on a immediatement que \(\mathcal{L}u = f\) presque partout dans \(\Omega\). Nous avons donc que \(-\nabla \cdot (\alpha \nabla u ) \in L^2(\Omega)\). Le corollaire [B59] implique que \(\alpha \nabla u \cdot n \in H^{\frac{1}{2}}(\partial \Omega)' = H^{-\frac{1}{2}}(\partial \Omega)\) du fait que \[ \forall \phi \in H^{\frac{1}{2}}(\partial \Omega), \quad \left\langle\alpha \nabla u \cdot n, \phi\right\rangle_{H{-\frac{1}{2}},H{\frac{1}{2}}} = \int_{\Omega} -\nabla \cdot (\alpha \cdot \nabla u) u_\phi + \int_\Omega \alpha \nabla u \cdot \nabla u_\phi \] où \(u_\phi \in H^1(\Omega)\) est un relèvement de \(\phi\) dans \(H^1(\Omega)\). On a alors que [prob12_2] donne \[ \forall \phi \in H^{\frac{1}{2}}(\partial \Omega), \quad \left\langle\alpha \nabla u \cdot n, \phi\right\rangle_{H{-\frac{1}{2}},H{\frac{1}{2}}} = \int_{\partial \Omega} g \phi \] montrant que \(\alpha \nabla u \cdot n = g\) dans \(H^{-\frac{1}{2}}(\partial \Omega)\) et donc, par conséquent, dans \(L^2(\Omega)\) du fait que \(g\) soit dans \(L^2(\Omega)\).\(\blacksquare\)

2.5. Conditions de Robin

documentation en cours

Soient deux fonctions \(g,\gamma : \partial \Omega \rightarrow \RR\), nous voulons imposer \(\gamma u + (\alpha \nabla u)\cdot n = g\) sur \(\partial \Omega\). En utilisant la relation sur l’intégrale de surface [eq:adr-green], nous avons la formulation faible suivante

On cherche \(u \in H^1(\Omega)\) telle que

\[ a(u,v) + \int_{\partial \Omega} \gamma u v = \ell(v) \quad \forall v \in H^1(\Omega) \] avec

\[ a(u,v)=\int_\Omega (\mathbf{\alpha} \nabla u ) \cdot \nabla v + (\mathbf{\beta} \cdot \nabla u) v + \mu u v \] et

\[ \ell(v) = \int_\Omega f v + \int_{\partial \Omega} g v \]

Soit \(g\in L^2(\partial \Omega)\) et \(\gamma \in L^\infty(\partial \Omega\). Si \(u\) résout [prob:12_4], alors \(\mathcal{L}u=f\) presque partout dans \(\Omega\) et \(\gamma u + \alpha \nabla u \cdot n = g\) presque partout sur \(\partial \Omega\).

On procède de manière similaire aux preuves précédentes. \(\blacksquare\)

Nous récapitulons dans la table suivante les différentes formulations

Table 2. Formulation faible pour différentes conditions aux limites
Problème Espace Forme bilinéaire Forme linéaire

Dirichlet homogène

\(H^1_0(\Omega)\)

\(a(u,v)\)

\(\int_\Omega f v\)

Neumann

\(H^1(\Omega)\)

\(a(u,v)\)

\(\int_\Omega f v + \int_{\partial \Omega} g v\)

Dirichlet Neumann

\(H^1_{\partial \Omega_D}(\Omega)\)

\(a(u,v)\)

\(\int_\Omega f v + \int_{\partial \Omega_N} g v\)

Robin (Fourier)

\(H^1(\Omega)\)

\(a(u,v) + \int_{\partial \Omega} \gamma uv\)

\(\int_\Omega f v + \int_{\partial \Omega} g v\)

Résumé

Les problèmes considérés précédent se mettent tous sous la forme suivante, sauf le problème de Dirichlet non homogène

\[\left\{\begin{array}{l} \text{ Chercher } u \in V \text{ telle que}\\ a(u,v)=f(v), \quad \forall v \in V \end{array}\right.\]

où \(V\) est un espace de Hilbert satisfaisant

\[H^1_0(\Omega) \subset V \subset H^1(\Omega)\]

Dans le cas de conditions de Dirichlet non-homogènes, nous avons \(u\in H^1(\Omega)\) et \(u=u_g+\phi\) où \(u_g\) est un relèvement de la donnée au bord \(g\) et \(\phi\) résout le problème générique avec des conditions de Dirichlet homogènes.

Conditions aux limites essentielles et naturelles

Il est important d’observer le traitement différent entre les conditions de Dirichlet et Neumann ou Robin conditions.

Les conditions de Dirichlet sont imposées explicitement dans l’espace fonctionnel où la solution est recherchée, et les fonctions de test disparaissent (i.e. \(v=0\)) sur la partie correspondante de la frontière. Pour cette raison, les conditions de Dirichlet sont souvent appelées conditions aux limites essentielles .

Les conditions de Neumann et Robin ne sont pas imposées par le cadre fonctionnel mais par la formulation faible elle-même. Le fait que les fonctions test ont des degrés de liberté sur la partie correspondante de la frontière est suffisant pour faire respecter les conditions limites en question. Pour cette raison, ces conditions sont souvent appelées conditions aux limites naturelles.

2.6. Coercivité

Cette section est en cours d’écriture

Soient \(f\in L^2(\Omega), stem:[\mathbf{\alpha} \in [L^{\infty}(\Omega)]^{d\times d}\), \(\mathbf{\beta} \in [L^{\infty}(\Omega)]^d\), \(\nabla \cdot \mathbf{\beta} \in L^{\infty}(\Omega)\) et \(\mu \in L^\infty(\Omega)\).

On note \(p = \essinf_{x \in \Omega} (\mu -\frac{1}{2} \nabla \cdot \mathbf{\beta})\) et \(c_\Omega\) est la constante de l’inégalité de Poincaré [B23].

(i) Les problèmes avec conditions de Dirichlet homogènes [prob12] et non-homogènes [prob12_2] sont bien posés si

\[ \alpha_0 > - \min( 0, \frac{\gamma}{c_\Omega} ) \]

(ii) Le problème avec condition de Neumann [prob12_1] est bien posé si

\[ p > 0\quad\text{et}\quad \essinf_{x \in \Omega}(\beta \cdot n ) \geq 0 \]

(iii) Le problème avec condition de Dirichlet-Neumann [prob12_3] est bien posé si [eq:92] est vérifiée et

\[ \mathrm{meas}(\partial \Omega_D) > 0\quad\text{et}\quad \partial \Omega^- = \{ x\in \partial \Omega ; (\beta \cdot n)(x) < 0\} \subset \partial \Omega_D \]

(iv) On pose \(q = \essinf_{x \in \Omega} (\gamma +\frac{1}{2} \mathbf{\beta}\cdot n)\). Le problème avec conditions de Robin [prob12_4] non-homogènes [prob12_2] est posé si

\[ p \geq 0,\quad q \geq 0,\quad\text{et}\quad pq \neq 0. \]

Nous prouvons (i) et (iv).

Preuve de (i) En utilisant l’ellipticité de \(\mathcal{L}\) et l’identité suivante(Formule de Green) \[ \int_\Omega u(\beta \cdot \nabla u) = -\frac{1}{2} \int_\Omega (\nabla \cdot \beta) u^2 + \frac{1}{2} \int_{\partial \Omega} (\beta \cdot n) u^2 \] alors on a \[ \forall u \in H^1_0(\Omega),\quad a(u,u) \geq \alpha_0|u|^2_{1,\Omega} + p \|u\|^2_{0,\Omega}. \] En posant \(c=\min(0,\frac{p}{c_\Omega}\) et en utilisant l’inégalité de Poincaré [B53] on a \[ \forall u \in H^1_0(\Omega),\quad a(u,u) \geq \left(\alpha_0+\frac{c_\Omega}\right)|u|^2_{1,\Omega} \geq \sigma \|u\|^2_{1,\Omega}. \] avec \(\sigma=\frac{c_\Omega(c_\Omega\alpha_0+c)}{1+c^2_\Omega}\).

Cela montre que \(a\) est coercive sur \(H^1_0(\Omega)\).

Le caractère bien posé résulte alors du Lemme de Lax-Milgram pour les conditions de Dicirichlet homogènes et de la proposition [TODO] pour les conditions de Dirichlet non-homogènes.\(\square\)

Preuve de (iv) Notons \(a(u,v)=\int_\Omega \alpha \nabla u\cdot \nabla u + (\beta \cdot \nabla u) v + \mu u v + \int_{\partial \Omega} \gamma u v\). Nous avons l’inégalité suivante: \[ \forall u \in H^1(\Omega),\quad a(u,u) \geq \alpha_0|u|^2_{1,\Omega} + p \|u\|^2_{0,\Omega} + q\|u\|^2_{0,\partial \Omega}. \]

Si \(p>0\) et \(q\geq 0\), la forme bilinéaire est coercive sur \(H^1(\Omega)\) avec comme constante stem:[\sigma = \min(\alpha_0,p).

Si \(p\geq 0\) et \(q > 0\), la forme bilinéaire est coercive sur \(H^1(\Omega)\) grâce au Lemme [TODO].

Dans les deux cas, le caractère bien posé est obtenu par le Lemme de Lax-Milgram.\(\square\) Ceci termine les preuves de (i) et (iv). \(\blacksquare\)

La coercivité est garantie si \(\alpha_0\) est suffisamment grand c’est à dire que si la diffusion est dominante.

2.7. Approximation conforme dans \(H^1\)

Cette section est en cours d’écriture

L’approximation élément fini est similaire à celle du Laplacian, de plus les variantes sur les conditions aux limites s’appliquent également: condition de Dirichlet non homogène, de Neumann ou de Robin.

Soit \(\Omega\) un polyèdre dans \(\RR^d\), soit \(\{\mathcal{T}_h\}\) une famille de maillages de \(\Omega\), et soit \(\{\hat{K}, \hat{P}, \hat{E}\}\) un élément fini Lagrange de référence du degré \(k \geq 1\).

Soit \(P^k_{c,h}\) l’espace d’approximation conformé \(H^1\) défini par

\[P^k_{c,h}=\{v \in C^0(\bar{\Omega}); \quad \forall K \in \mathcal{T}_h v_{|K} \in \mathbb{P}^k (K) \}\]

Pour obtenir une approximation conforme dans \(V\), nous rajoutons les conditions aux limites, i.e, nous avons

\[V_h = P^k_{c,h} \cap V\]

Dans le cas de conditions de Dirichlet homogène cela donne

\[V_h=\{ v_h \in P^k_{c,h} ;\ v_h = 0 \mbox{ sur } \partial \Omega\}\]

Dans le cas Neumann et Robin, nous avons \(Vh=P^k_{c,h}\). Enfin dans le cas Dirichlet-Neumann, nous avons

\[V_h=\{v_h \in P^k_c; \ v_h = 0 \mbox{ sur } \partial \Omega_D\}\]

Le problème discret s’ecrit

Trouver \(u_h\) dans \(V_h\) telle que \[ a(u_h,v_h)=\ell(v_h), \quad \forall v_h \in V_h \]

Nous cherchons à présent à estimer l’erreur \(u-u_h\) en norme \(H^1\) et \(L^2\).

Estimation \(H^1\)

Soit \(\Omega\) un polyèdre dans \(\RR^d\), soit \(\{\mathcal{T}_h\}\) une famille de maillages de \(\Omega\), et soit \(\{\hat{K}, \hat{P}, \hat{E}\}\) un élément fini Lagrange de référence du degré \(k \geq 1\). Nous avons \(\lim_{h\rightarrow 0} \|u-u_h\|_{1,\Omega} = 0\) et pour \(\frac{d}{2} < s < k+1, il existe une constante stem:[c\) telle que \[ \forall h, \quad \|u-u_h\| \leq c h^{s-1} |u|_{s,\Omega} \]

Estimation \(L^2\)

Soient les hypothèse du théorème précédent, auxquelles nous ajoutons des hypothèses initial a des propriétés régularisantes Nous avons \(\lim_{h\rightarrow 0} \|u-u_h\|_{1,\Omega} = 0\) et pour \(\frac{d}{2} < s < k+1, il existe une constante stem:[c\) telle que \[ \forall h, \quad \|u-u_h\|{0,\Omega} \leq c h |u-u_h|{1,\Omega} \]

nous n’avons pas défini ce que sont ces propriétés régularisantes. Nous supposerons qu’elles sont vérifiées.
Exemple du Laplacien avec conditions de Dirichlet homogène en P1

\forall h, \quad \|u-u\|{0,\Omega} + h \|u-u_h\|{1,\Omega} \leq c h^2 \|f\|_{0,\Omega}

Convection dominated flows

Let \(\Omega \subset \mathbb{R}^d, d=1,2,3\), consider this equation defined in \(\Omega\), find \(u\) a scalar field (e.g. temperature or concentration) such that

\[ \underbrace{\frac{\partial u}{\partial t}}_{time} \underbrace{- \nabla \cdot ( \epsilon \nabla u )}_{diffusion} + \underbrace{\mathbf{ \beta} \cdot \nabla u}_{advection/convection} + \underbrace{\mu u}_{reaction} = \underbrace{f}_{source},\quad \nabla \cdot \mathbf{ \beta} = 0\]

with the given data

From now on we consider only steady problems (\(\frac{\partial u}{\partial t}=0\))

1. Applications

  • Transport of temperature

  • Transport of pollutants

  • Transport of chemical species (\(O_2, CO_2,...\))

These equations are often coupled with fluid flow equations such as incompressible Navier-Stokes equations (air, blood, water…​) or darcy equations (porous media) to obtain the velocity field \(\mathbf{ \beta}\). In this lesson we consider \(\mathbf{ \beta}\) given.

2. Boundary conditions

Boundary conditions can be of 2 types on parts of the boundary or the whole boundary of the domain \(\Omega\):

  • Dirichlet conditions : \(\label{eq:2} u|_{{\partial \Omega}_D} = g\)

  • Neumann conditions :

    • Outflow \(\label{eq:4} ( - \epsilon \nabla u + \mathbf{ \beta} u ) \cdot \mathbf{ n}|_{{\partial \Omega}_N} = (\mathbf{ \beta} u) \cdot \mathbf{ n}\)

    • (Heat) Flux \(\label{eq:5} ( - \epsilon \nabla u ) \cdot \mathbf{ n}|_{{\partial \Omega}_N} = Q \quad (e.g. \mathbf{\beta} = \mathbf{0})\)

3. Variational Formulation

The variational problem reads, Find \(u \in V\) such that for all \(v\ in V\)

\[ a(u,v) \equiv \int_\Omega \epsilon \nabla u \cdot \nabla v + (\mathbf{ \beta} \cdot \nabla u ) v + \underbrace{\int_{\partial \Omega} (- \epsilon \nabla u ) \cdot \mathbf{ n}\ v}_{\mbox{apply boundary conditions}} = \int_\Omega f v \equiv \ell (v)\]

where

  • \(V\) is an Hilbert space

  • we used the identity below for the integration by part

\[\nabla \cdot (\mathbf{ c} w ) = \mathbf{ c} \cdot \nabla w + \underbrace{w \nabla \cdot \mathbf{ c}}_{=0}\]
  • Suppose that \(-\frac{1}{2} \nabla \cdot ( \mathbf{ \beta} ) + \mu \geq 0\) almost everywhere in \(\Omega\) then we can show that \(a\) is coercive provide \(\epsilon \geq \epsilon_0 > 0\)

\[a(v,v) \geq \alpha ||v||_{1,\Omega}, \quad \alpha = \frac{\epsilon_0}{1+C_\Omega^2}\]
  • \(a\) is continuous, there exists \(M\) a constant such that

\[|a(u,v)| \leq < M ||u||_{H^1{(\Omega)}} ||v||_{H^1{(\Omega)}}, \quad M = ||\mu||_{0,\Omega} + ||\epsilon||_{\infty,\Omega} + ||\mathbf{\beta}||_{\infty,\Omega}\]

4. Discrete formulation

Let \(V_h \subset V \equiv H^1(\Omega)\) a discrete space, consider the standard Galerkin approximation on \(V_h\) for . The problem reads

Find \(u_h \in V_h\) such that for all \(v_h \in V_h\) we have:

\[ \int_\Omega \epsilon \nabla u_h \cdot \nabla v_h + (\mathbf{ \beta} \cdot \nabla u_h ) v_h + \underbrace{\int_{\partial \Omega} (- \epsilon \nabla u_h ) \cdot \mathbf{ n}\ v_h}_{\mbox{apply boundary conditions}} = \int_\Omega f v_h \]

We can show that

\[||u_h||_{1,\Omega} \leq \frac{1}{\alpha} ||f||_{0,\Omega}, \quad ||\nabla u_h||_{1,\Omega} \leq \frac{C_\Omega}{\epsilon_0} ||f||_{0,\Omega},\]

which means that the Galerkin error inequality (best approximation) gives

\[||u-u_h||_V \leq \frac{M}{\alpha} \mathrm{inf}_{v_h \in V_h} ||u-v_h||_V\]

which means that given \(M\) and \(\alpha\), the estimate \(\epsilon_0 \rightarrow 0\) In such case, the standard Galerkin method can yield to inaccurate solutions unless:

  • we use a very small \(h\) (cost!!)

  • we use a stabilisation method (loss of \(\frac{1}{2}\) in convergence rate)

5. Stabilisation methods for convection dominated flows

We now introduce

  • \(\mathrm{Pe}=\frac{|\mathbf{ \beta}|L}{2 \epsilon}\) the global number. \(L\) is the characteristic size of the domain.

  • the local Péclet number \(\mathrm{Pe}_h=\frac{|\mathbf{\beta}|h}{2 \epsilon}\)

The dominating convection and inacurrate behavior occurs when, on at least some cells, \(\mathrm{Pe} > 1\) or \(\mathrm{Pe}_h > 1\).

A remedy is to use a sufficiently small \(h\) same to ensure \(\mathrm{Pe}_h < 1\). For example if \(|\mathbf{beta}| = 1\) and \(\epsilon=5e-4\), we should take \(h=1e-4\).
Another remedy is to use a different approximation space for the unknown \(u_h\) and the test functions \(v_h\). We talk about Petrov-Galerkin approximation.

Find \(u_h \in V_h\) such that

\[ a_h(u_h,v_h) = \ell_h(v_h),\quad \forall v_h \in V_h \] with

\[ a_h(u_h,v_h) = a(u_h,v_h) + b_h(u_h,v_h),\quad \ell_h(v_h) = \ell(v_h) + g_h(v_h),\quad \].

The purpose of \(b_h\) and \(g_h\) is to eliminate(or reduce) the numerical oscillations produced by the standard Galerkin method: they are called stabilisation methods and depend on \(h\). The term stabilisation method is not exact. Indeed the Galerkin method is already stable (i.e. continuity). Here it is to be understood as the aim of reducing (or elimination) numerical oscillations when \(\mathrm{Pe} > 1\).

Without doing anything wiggles occur.

There are remedies so called stabilisation techniques, here some some examples:

  • Artificial diffusion (streamline diffusion) (SDFEM)

  • Galerkin Least Squares method (GaLS)

  • Streamline Upwind Petrov Galerkin (SUPG)

  • Continuous Interior Penalty methods (CIP)

5.1. Artificial diffusion (or streamline diffusion) (SDFEM)

Method The method consists in adding an

\[\epsilon_h =\epsilon(1+\phi(\mathrm{Pe}))\]

with \(\phi(\mathrm{Pe}) \rightarrow 0\) as \(h \rightarrow 0\), e.g. \(\phi(\mathrm{Pe}) = \mathrm{Pe}-1+B(2*\mathrm{Pe})\) where \(B\) is the so-called Bernoulli function \(B(t) = \frac{t}{e^t-1}\) if \(t > 0\) and \(B(0) = 1\) (also exponential fitting scheme)

\[ b_h(u_h,v_h) = \int_\Omega \epsilon \Phi(\mathrm{Pe}) \nabla u_h \cdot \nabla v_h, \quad g_h(v_h) = 0\]

for a given \(\epsilon\) and for \(h\) tending to \(0\), we have for \(u \in H^{r+1}(\Omega)\)

\[ ||u-u_h||{1,\Omega} \leq C_1 \Big[ h^r||u||{r+1,\Omega} + \phi(\mathrm{Pe})||u||_{1,\Omega}\Big] \] and for a given \(h\) and \(\epsilon\) tending to 0,

\[ ||u-u_h||{1,\Omega} \leq C_1 \Big[ h^{r-1}||u||{r+1,\Omega} + ||u||_{1,\Omega}\Big] \] If \(\phi(\mathrm{Pe})=\frac{|\mathbf{ \beta}|h}{2 \epsilon}\), the convergence is linear, with the exponential fitting scheme it is quadratic if \(r \geq 2\).

5.2. GaLS and SUPG

First we decompose our operators into a symmetric (\(<Lu,v> = <u,Lv>\) and skew symmetric (\(<L u, v> = -<u,L v>\)) contributions, we start with

\[ L u = -\epsilon \Delta u + \nabla \cdot (\mathbf{ \beta} u ) + \mu u\]
\[L u = \underbrace{-\epsilon \Delta u + \Big[ \mu + \frac{1}{2} \nabla \cdot \mathbf{ \beta} \Big] u}_{L_S u} + \underbrace{\frac{1}{2}\Big[ \nabla \cdot ( \mathbf{ \beta} u) + \mathbf{ \beta} \cdot \nabla u \Big]}_{L_{SS} u}\]
Consistent schemes

We say that a method is consistent when adding a term to a problem such as:

Find \(u_h \in V_h\) such that

\[ a(u_h,v_h) + \mathcal{L}_h(u_h,f;v_h) = (f,v_h), \quad \forall v_h \in V_h\] the term added statisfies

\[ \mathcal{L}_h(u,f;v_h) = 0, \forall v_h \in V_h \]

5.2.1. Choice for consistent methods

A possible choice for \(\mathcal{L}_h\) is the following

\[ \mathcal{L}_h(u_h,f;v_h) = \mathcal{L}^{(\rho)}_h(u_h,f;v_h) = \sum_{K \in \mathcal{T}_h} \delta (L u_h - f, \mathcal{S}^{(\rho)}_K(v_h))_{0,\Omega}\]

where

  • \((\cdot,\cdot)_{0,\Omega}\) is the \(L^2\) scalar product

  • \(\rho\) and \(\delta\) are parameters

and we have set

\[\mathcal{S}^{(\rho)}_K(v_h) = \frac{h_K}{|\mathbf{\beta}|}\Big[ L_{SS} v_h + \rho L_S v_h\Big]\]
Galerkin Least-Square

if \(\rho = 1\) we have the Galerkin Least Square method (GaLS)

\[\mathcal{S}^{(\rho)}_K(v_h) = \frac{h_K}{|\mathbf{ \beta}|}\Big[ L v_h\Big]\]
Streamline Upwind Petrov-Galerkin

if \(\rho = 0\) we have the Streamline Upwind Petrov-Galerkin (SUPG)

\[\mathcal{S}^{(0)}_K(v_h) = \frac{h_K}{|\mathbf{ \beta}|}\Big[ L_{SS} v_h\Big]\]
Douglas and Wang

if \(\rho = -1\) we have the Douglas and Wang (DW)

\[\mathcal{S}^{(-1)}_K(v_h) = \frac{h_K}{|\mathbf{ \beta}|}\Big[ (L_{SS} -L_S )v_h\Big]\]

We define the \(\rho\) Norm

\[||v||_{(\rho)} = \Big\{\epsilon ||\nabla u||^2_{0,\Omega} + ||\sqrt{\gamma} v||^2_{0,\Omega} + \sum_{K \in \mathcal{T_h}} \delta \Big( (L_{SS}+\rho L_S )v, \mathcal{S}^{(\rho)}_K(v) \Big)_{0,\Omega} \Big\}^{1/2}\]

where \(\gamma\) is a positive constant such that \(-\frac{1}{2} \nabla \cdot \mathbf{\beta} + \mu \geq \gamma > 0\)

We have the following result

if \(u \in H^{k+1}(\Omega)\), then the following error estimates hold:

\[ {\|u-u_h\|{(\rho)}} \leq C {h^{k+1/2}} |u|{k+1,\Omega} \]

GaLS

In practice for GaLS (\(\rho = 1\)) we take \(\delta\) such that

\[\delta(h_K,\epsilon) \frac{h_K}{|\mathbf{ \beta}|} = \Big( \frac{1}{h_K} + \frac{\epsilon}{h^2_K} \Big)^{-1}\]

and we can prove the following estimates if \(u\in H^{k+1}(\Omega)\),

\[\forall \epsilon \quad {\|u-u_h\|_{0,\Omega}} \leq c {h^{k+1/2}} \|u\|_{k+1,\Omega}\]
\[\forall \epsilon \geq c h \quad {\|u-u_h\|_{1,\Omega}} \leq c {h^{k}} \|u\|_{k+1,\Omega}\]

and finally if the family \(\{\mathcal{T}_h\}_{h > 0}\) is quasi-uniform and \(\epsilon \leq c h \), then

\[\| \beta \cdot \nabla (u -u_h) \|_{0,\Omega} \leq c h^k \|u \|_{k+1,\Omega}\]

5.3. Continuous Interior Penalty

In the continuous interior penalty we add the following term

\[\sum_{F \in \Gamma_\mathrm{int} } \int_{F} \eta\ h_F^2\ |\mathbf{ \beta} \cdot \mathbf{n}|\ \jump{\nabla u} \jump{\nabla v}\]

where

  • \(\Gamma_\mathrm{int}\) is the set of internal faces

  • the \(\mathrm{Pe}>>1\) (typically it is applied to all internal faces)

  • we have

\jump{\nabla u} = \nabla u \cdot \mathbf{n}|_1 + \nabla u \cdot \mathbf{n}|_2

is the so called jump of \(\nabla u\)(scalar valued) across the face.

In the case of scalar valued functions

\[ \jump{u} = u \mathbf{n}|_1 + u \mathbf{n}|_2\]

Choice for \(\eta\) \(\eta\) can be taken in the range \([1e-2;1e-1\)]. A typical value is \(\eta=2.5e-2\). A similar error estimate \(O(h^{r+1/2})\) holds for CIP.

Example CIP

// define the stabilisation coefficient expression
auto stab_coeff = (eta*abs(idv(beta))*abs(trans(N())*idv(beta)))*vf::pow(hFace(),2.0));

// assemble the stabilisation operator
form2( Xh, Xh, M ) +=
 integrate( internalfaces(Xh->mesh()), // faces of the mesh
            stab_coeff*(trans(jumpt(gradt(u)))*jump(grad(v))));

6. Élasticité Linéaire

On s’intéresse dans cette section à l’approximation par éléments finis de problèmes de mécanique des milieux continus en 3D.

Soit \(\mathbf{f}: \Omega \rightarrow \mathbb{R}^3\) la charge extérieur s’appliquant au domaine \(\Omega\).

On note \(\disp: \Omega \rightarrow \mathbb{R}^3\) le déplacement de la structure induit par cette charge \(\mathbf{f}\).

En supposant que les déformations soient petites pour être modélisées dans le cadre de l’elasticité linéaire, on a la relation suivante à l’équilibre:

\[ \nabla \cdot \stresst + \mathbf{f} = \mathbf{0} \mbox{ dans } \Omega \] où \(\stresst: \Omega \rightarrow \RR^{d\times d}\) est le tenseur des contraintes défini par la relation

\[ \stresst = \lambda \tr(\deformt) \Id + 2 \mu \deformt \] et \(\deformt : \Omega \rightarrow \RR^{d\times d}\) est le tenseur des déformations défini par

\[ \deformt = \frac{1}{2} \left( \nabla \disp + \nabla \disp^T \right), \] \(\lambda\) et \(\mu\) les coefficients de Lamé et \(\Id\) la matrice identité de \(\RR^{d\times d}\).

On a alors

\[ \stresst = \lambda( \nabla \cdot \disp ) \Id + \mu( \nabla \disp + \nabla \disp^T) \]

6.1. Coefficients de Lamé

Ils sont des coefficients phénoménologiques contraints par les relations suivants:

  • \(\mu >0\)

  • \(\lambda + \frac{2}{3} \mu \ge 0\)

Dans ce qui suit, on supposera que \(\lambda \ge 0\) et ces coefficients constants.

D’un point de vue pratique, ces coefficients sont obtenus par les module d’Young \(E\) et coefficient de Poisson \(\nu\) tels que

\[ \lambda = \frac{E \nu}{( 1+\nu )*( 1-2 \nu )} , \quad \mu =\frac{E}{2 ( 1+\nu )} \]

7. Approximation de problèmes mixtes

7.1. Model Problems

We consider now model problems as systems of PDEs where several functions are unknowns and which don’t play the same roles mathematically and physically.

Stokes

\[ \left\{\begin{array}[c]{rl} -\Delta u + \nabla p & = f\ \mbox{ in } \Omega\\ \nabla \cdot u & = 0\ \mbox{ in } \Omega \end{array}\right. \] where \(u: \Omega \mapsto \RR^d\) is a velocity and \(p: \Omega \mapsto \RR\) is a pressure.

Darcy

\[ \left\{\begin{array}[c]{rl} \sigma + \nabla u & = f\ \mbox{ in } \Omega\\ \nabla \cdot \sigma & = g\ \mbox{ in } \Omega \end{array}\right. \] where \(\sigma: \Omega \mapsto \RR^d\) is a velocity and \(u: \Omega \mapsto \RR\) is a hydraulic charge(pressure).

7.2. Applications

We shall focus on Stokes, but the abstract setting of the next section is the same for Stokes and Darcy.

Stokes and incompressible Navier-Stokes for Newtonian fluids

The Stokes model is the basis for fluid mechanics models and is a simplication of the Navier-Stokes equations where the viscous effects/terms are much bigger than the convective ones

\[ \left\{\begin{array}[c]{rl} \rho( \frac{\partial u}{\partial t} + u \cdot \nabla u) - \nu \Delta u + \nabla p & = f\ \mbox{ in } \Omega\\ \nabla \cdot u & = 0\ \mbox{ in } \Omega \end{array}\right. \] The first equation results from the conservation of momentum and the second from the conservation of mass.

The well-posedness of these problems results from a so-called which is not automatically transfered at the discrete level.

In practice In order to ensure that the finite element approximation is well-posed, we will need to choose approximation spaces that satisfy a compatibility condition that ensures that a discrete inf-sup condition is satisfied.

7.3. Saddle point problems

7.3.1. Abstract Continuous Setting

Denote

  • \(X\) and \(M\) two Hilbert spaces.[2]

  • two linear forms \(f \in X'=\mathcal{L}(X, \RR)\) and \(g \in M'=\mathcal{L}(M, \RR)\)

  • \(a \in \mathcal{L}(X\times X, \RR)\) and \(b \in \mathcal{L}(X\times M, \RR)\) two bilinear forms

We are interested in the following abstract problem:

Look for \((u,p) \in X \times M\) such that

\[ \left\{ \begin{array}[c]{rl} a(u,v) + b(v,p) & = f(v), \quad \forall v \in X\\ b(u,q) & = g(q), \quad \forall q \in M \end{array} \right. \]

Definition of a saddle point problem

If the bilinear form \(a\) is symmetric and positive on \(X\times X\), we say that [prob:chmixte:1] is a saddle point problem.

The structure of the problem is as follows

  • the space of solution is the same of the test space

  • the unknown \(p\) does not appear in the second equation

  • the unknown functions \(u\) and \(p\) are coupled via the same bilinear form \(b\) is the first and second equation.

The next question is :

7.3.2. Well posed problem

Reformulation

Let’s rewrite Problem [prob:chmixte:1].

Denote \(V=X\times M\) and introduce \(c \in \mathcal{L}(V\times V, \RR)\) such that

\[ cu,p),(v,q = a(u,v)+b(v,p)+b(u,q) \] and \(h\in \mathcal{L}(V,\RR)\) such that

\[ h(v,q) = f(v)+g(q) \] then problem [prob:chmixte:1] reads

Look for \((u,p) \in V\) such that

\[ \begin{array}[c]{rl} cu,p), (v,q & = h(v,q), \quad \forall (v,q) \in V \end{array} \]

We suppose that \(a\) is coercive over \(X\), the [prob:chmixte:2] is well-posed if and only if the bilinear form \(b\) satisfies the following inf-sup condition:

there exists \(\beta > 0\) such that

\[ \inf_{q \in M} \sup_{v \in X} \frac{b(v,q)}{||v||_X ||q||_M} \geq \beta \]

Lax-Milgram provides only a sufficient condition for well-posedness
The form \(c\) in [prob:chmixte:2] does not satisfy Lax-Milgram.

Let’s introduce the so-called Lagrangian \(l \in \mathcal{L}(X\times M, \RR)\) defined by

\[ l(v,q) = \frac{1}{2} a(v,v) + b(v,q) - f(v) - g(q) \]

We say that the point \((u,p)\in X\times M\) is a saddle point of \(l\) if

\[ \forall (v,q) \in X\times M, \quad l(u,q) \leq l(u,p) \leq l(v,p) \]

Under the hypothesys oF [thr:chmixte:1], the Lagrangian \(l\) defined by has a unique saddle point. Moreover this saddle point is the unique solution of problem [prob:chmixte:1].

7.4. Finite element approximation

7.4.1. Abstract Discrete Problem

We now turn to the approximation of the problem [prob:chmixte:1] by a standard Galerkin method in a conforming way.

Denote the two spaces \(X_h \subset X\) and \(M_h \subset M\), we consider the following problem:

Look for \((u_h,p_h) \in X_h \times M_h\) such that

\[ \left\{ \begin{array}[c]{rl} a(u_h,v_h) + b(v_h,p_h) & = f(v_h), \quad \forall v_h \in X_h\\ b(u_h,q_h) & = g(q_h), \quad \forall q_h \in M_h \end{array} \right. \]

We suppose that \(a\) is coercive over \(X\) and that \(X_h \subset X\) and \(M_h \subset M\).

Then the [prob:chmixte:3] is well-posed if and only if the following discrete inf-sup condition is satisfied:

there exists \(\beta_h > 0\) such that

\[ \inf_{q_h \in M_h} \sup_{v_h \in X_h} \frac{b(v_h,q_h)}{||v_h||{X_h} ||q_h||{M_h}} \geq \beta_h \]

The compatibility condition problem [prob:chmixte:3], to be well posed, requires that the spaces \(X_h\) and \(M_h\) satisfy the condition.

This is known as the Babuska-Brezzi (BB) or Ladyhenskaya-Babuska-Brezzi (LBB).

Regarding error analysis, we have the following lemma

Thanks to the Lemma of Céa applied to Saddle-Point Problems, the unique solution \((u,p)\) of problem [prob:chmixte:3] satisfies

\[ \begin{array}[c]{rl} ||u-u_h||X & \leq c{1h} \inf_{v_h \in X_h} ||u-v_h||X + c{2} \inf_{q_h \in M_h} ||q-q_h||M\\ ||p-p_h||_X & \leq c{3h} \inf_{v_h \in X_h} ||u-v_h||X + c{4h} \inf_{q_h \in M_h} ||q-q_h||_M \end{array} \] where

  • \(c_{1h} = (1+\frac{||a||_{X,X}}{\alpha})(1+\frac{||b||_{X,M}}{\beta_h})\) with \(\alpha\) the coercivity constant of \(a\) over X.

  • \(c_{2} = \frac{||b||_{X,M}}{\alpha}\)

  • \(c_{3h} = c_{1h} \frac{||a||_{X,X}}{\beta_h}\), \(c_{4h} = 1+ \frac{||b||_{X,M}}{\beta_h}+\frac{||a||_{X,X}}{\beta_h}\)

The constants \(c_{1h}, c_{3h}, c_{4h}\) are as large as \(\beta_h\) is small.

7.4.2. Linear system associated

The discretisation process leads to a linear system.

We denote

  • \(N_u = \dim {X_h}\)

  • \(N_p = \dim {M_h}\)

  • \(\{\phi_i\}_{i=1,...,N_u}\) a basis of \(X_h\)

  • \(\{\psi_k\}_{k=1,...,N_p}\) a basis of \(M_h\)

  • for all \(u_h = \sum_{i=1}^{N_u} u_i \phi_i\), we associate \(U \in \R{N_u}\), \(U=(u_1,\ldots,u_{N_u})^T\), the component vector of \(u_h\) is \(\{\phi_i\}_{i=1,\ldots,N_u}\)

  • for all \(p_h = \sum_{k=1}^{N_p} u_k \psi_k\), we associate \(P \in \R{N_p}\), \(P=(p_1,\ldots,p_{N_p})^T\), the component vector of \(p_h\) is \(\{\psi_k\}_{k=1,\ldots,N_p}\)

The matricial form of problem [prob:chmixte:3] reads

\[ \begin{bmatrix} \mathcal{A} & \mathcal{B}^T\\ \mathcal{B} & 0 \end{bmatrix} \begin{bmatrix} U \\ P \end{bmatrix} = \begin{bmatrix} F\\ G \end{bmatrix} \]

where the matrix \(\mathcal{A} \in \R{N_u,N_u}\) and \(\mathcal{B} \in \R{N_p,N_u}\) have the coefficients

\[ \mathcal{A}_{ij} = a(\phi_j,\phi_i), \quad \mathcal{B}_{ki} = b(\phi_i,\psi_k) \]

and the vectors \(\mathcal{F} \in \R{N_u}\) and \(\mathcal{G} \in \R{N_p}\) have the coefficients

  • \(F_i=f(\phi_i)\)

  • \(G_k=g(\psi_k)\)

  1. Since \(a\) is symmetric and coercive, \(\mathcal{A}\) is symmetric positive definite

  2. The matrix of the system is symmetric but not positive

  3. The inf-sup condition  is equivalent to the fact that \(\mathcal{B}\) is of maximum rank, i.e. \(\ker(\mathcal{B}^T) = \{0 \}\).

  4. From theorem [thr:chmixte:2], the matrix of the system  is invertible

When the inf-sup is not satisfied

The counter examples when the inf-sup condition  is not satisfied(e.g. \(\mathcal{B}\) is not maximum rank ) occur usually in two cases:

Locking

\(\dim {M_h} > \dim {X_h}\): the space of pressure is too large for the matrix \(\mathcal{B}\) to be maximum rank. In that case \(\mathcal{B}\) is injective (\(\ker(\mathcal{B}) = \{0\})\). We call this *locking*.

Spurious modes

there exists a vector \(Q^* \neq 0\) in \(\ker(\mathcal{B}^T)\). The discrete field\(q^*_h\) in \(M_h\), \(q^*_h=\sum_{k=1}^{N_p} Q^*_k \psi_k\), associated is called a *spurious mode*. \(q^*_H\) is such that

\[ b(v_h,q^*_h)=0. \]

We now introduce the Uzawa matrix as follows

The matrix

\[ \mathcal{U} = \mathcal{B} \mathcal{A}^{-1} \mathcal{B}^T \] is called the Uzawa matrix. It is symmetric positive definite from the properties of \(\mathcal{A}\), \(\mathcal{B}\)

Applications

The Uzawa matrix occurs when eliminating the velocity in system  and get a linear system on \(P\):

\[ \mathcal{U} P = \mathcal{B} \mathcal{A}^{-1} F - G \] then one application is to solve by solving iteratively and compute the velocity afterwards.

7.5. Mixed finite element for Stokes

7.5.1. Variational formulation

We start with the Well-posedness at the continuous level

  • We consider the model problem  with homogeneous Dirichlet condition on velocity \(u = 0\) on \(\partial \Omega\)

  • We suppose the \(f \in [L^2(\Omega)\)^d] and \(g \in L^2(\Omega)\) with \[ \int_\Omega g = 0 \] Introduce

    \[ L^2_0(\Omega) = \Big\{ q \in L^2(\Omega): \int_\Omega q = 0 \Big\} \]

The condition comes from the divergence theorem applied to the divergence equation and the fact that \(u=0\) on the boundary

\[ \int_\Omega g = \int_\Omega \nabla \cdot u = \int_{\partial \Omega} u \cdot n = 0 \] This is a necessary condition for the existence of a solution \((u,p)\) for the Stokes equations with these boundary conditions.

We turn now to the variational formulation.

The Stokes problem reads

Look for \((u,p) \in [H^1_0(\Omega)\)^d \times L^2_0(\Omega)] such that

\[ \left\{ \begin{array}[c]{rl} \int_\Omega \nabla u : \nabla v -\int_\Omega p \nabla \cdot v & = \int_\Omega f \cdot v, \quad \forall v \in [H1_0(\Omega)]d\\ \int_\Omega q \nabla \cdot u & = - \int_\Omega g q, \quad \forall q \in L^2_0(\Omega) \end{array} \right. \]

We retrieve the problem [prob:chmixte:1] with \(X=[H^1_0(\Omega)\)^d] and \(M=L^2_0(\Omega)\) and

\[ \begin{array}[c]{rlrl} a(u,v) &= \int_\Omega \nabla u : \nabla v,& \quad b(v,p) &= -\int_\Omega p \nabla \cdot v,\\ \quad f(v) &= \int_\Omega f \cdot v,& \quad g(q) &= - \int_\Omega g q \end{array} \]

Pressure up to a constant
The pressure is known up to a constant, that’s why we look for them in \(L^2_0(\Omega)\) to ensure uniqueness.

7.5.2. Finite element approximation

Denote \(X_h \subset [H^1_0(\Omega)\)^d] and \(M_h \subset L^2_0(\Omega)\)

Look for \((u_h,p_h) \in X_h \times M_h\) such that

\[ \left\{ \begin{array}[c]{rl} \int_\Omega \nabla u_h : \nabla v_h + \int_\Omega p_h \nabla \cdot v_h & = \int_\Omega f \cdot v_h, \quad \forall v_h \in X_h\\ \int_\Omega q_h \nabla \cdot u_h & = -\int_\Omega g q_h, \quad \forall q_h \in M_h \end{array} \right. \]

This problem, thanks to theorem [thr:chmixte:2] is well-posed if and only if \(X_h\) and \(M_h\) are such that there exists \(\beta_h > 0\)

\[ \inf_{q_h \in M_h} \sup_{v_h \in X_h} \frac{\int_\Omega q_h \nabla \cdot v_h}{||v_h||{X_h} ||q_h||{M_h}} \geq \beta_h \]

7.5.3. Some counter examples: bad finite element for Stokes

In this section, we present two classical bad finite element approximations.

Finite element \(\poly{P}_1/\poly{P}_0\): locking

Thanks to the Euler relations, we have

\[ \begin{array}[c]{rl} N_{\mathrm{cells}} - N_{\mathrm{edges}} + N_{vertices} &= 1-I\\ N^\partial_{\mathrm{vertices}} - N^\partial_{\mathrm{edges}} &= 0 \end{array} \]

where \(I\) is the number of holes in \(\Omega\).

We have that \(\dim {M_h} = N_{\mathrm{cells}}\),\(\dim {X_h} = 2 N^i_{\mathrm{vertices}}\) and so

\[ \dim {M_h} - \dim {X_h} = N_{\mathrm{cells}} - 2 N^i_{\mathrm{vertices}} = N^\partial_{\mathrm{edges}} - 2 > 0 \]

so \(M_h\) is too rich for the condition and we have \(\ker(\mathcal{B}) = \{0\}\) such that the only discrete \(u_h^*\), with components \(U^*\), satisfying \(\mathcal{B} U^*\) is the null field, \(U^*=0\).

Finite element \(\poly{Q}_1/\poly{P}_0\): spurious mode

We can construct in that case a function \(q_h^*\) on a uniform grid which is equal alternatively -1, +1 (chessboard) in the cells of the mesh, then

\[ \forall v_h \in [Q1_{c,h}]d, \quad \int_\Omega q^*_h \nabla \cdot v_h = 0 \] and thus the associated \(X_h\), \(M_h\) do not satisfy the condition.

Finite element \(\poly{P}_1/\poly{P}_1\): spurious mode

We can construct in that case a function \(q_h^*\) on a uniform grid which is equal alternatively -1, 0, +1 at the vertices of the mesh, then

\[ \forall v_h \in [P1_{c,h}]d, \quad \int_\Omega q^*_h \nabla \cdot v_h = 0 \] and thus the associated \(X_h\), \(M_h\) do not satisfy the condition.

7.5.4. Mini-Element

The problem with the \(\poly{P}_1/\poly{P}_1\) mixed finite element is that the velocity is not rich enough.

A cure is to add a function \(v_h^*\) in the velocity approximation space to ensure that

\[ \int_\Omega q^_h \nabla \cdot v_h^ \neq 0 \] where \(q_h^*\) is the spurious mode.

To do that we add the bubble function to the \(\poly{P}_1\) velocity space.

Recall the construction of finite elements on a reference convex \(\hat{K}\). We say that \(\hat{b}: \hat{K} \mapsto \RR\) is a bubble function if:

  • \(\hat{b} \in H^1_0(\hat{K})\)

  • \(0 \leq \hat{b}(\hat{x}) \leq 1, \quad \forall \hat{x} \in \hat{K}\)

  • \(\hat{b}(\hat{C}) = 1, \quad \mbox{where} \hat{C}\) is the barycenter of \(\hat{K}\)

Example

The function

\[ \hat{b} = (d+1)^{d+1} \Pi_{i=0}^d\ \hat{\lambda}_i \] where \((\hat{\lambda}_0, \ldots, \hat{\lambda}_d)\) denote the barycentric coordinates on \(\hat{K}\)

Denote now \(\hat{b}\) a bubble fonction on \(\hat{K}\), we set

\[ \hat{P} = [\poly{P}_1(\hat{K}) \oplus \mathrm{span} (\hat{b})]^d, \] and introduce

X_h &=& \Big\{ v_h \in [C^0(\bar{\Omega})]^d : \forall K \in \mathcal{T}_h, v_h
\circ T_K \in \hat{P}; v_{h_|{\partial \Omega}} = 0 \Big\}\\
M_h &=& P^1_{c,h}

The spaces \(X_h\) and \(M_h \cap L^2_0(\Omega)\) satisfy the compatibility condition  uniformly in \(h\).

Suppose that \((u,p)\), solution of [prob:chmixte:1], is smooth enough, ie. \(u \in [H^2(\Omega)\)^d \cap [H1_0(\Omega)]d] and \(p\in H^1(\Omega) \cap L^2_0(\Omega)\).

Then there exists a constant \(c\) such that for all \(h >0\)

\[ \| u- u_h \|{1,\Omega} + \|p-p_h\|{0,\Omega} \leq c h (\|u\|{2,\Omega} + \|p\|{1,\Omega}) \] and if the Stokes problem is stabilizing then

\[ \|u-u_h\|{0,\Omega} \leq c h^2 ( \|u\|{2,\Omega} +\|p\|_{1,\Omega}). \]

Stabilizing Stokes problem

We say that the Stokes problem is stabilizing if there exists a constant \(c_S\) such that for all \(f \in [L^2(\Omega)\)^d], the unique solution \((u,p)\) of with \(g=0\) is such that:

\[ \|u\|{2,\Omega} + \|p\|{1,\Omega} \leq c_S \|f\|_{0,\Omega} \] A sufficient condition for stabilizing Stokes problem is that the \(\Omega\) is a polygonal convex in 2D or of class \(C^1\) in \(\RR^d, d=2,3\).

7.5.5. Taylor-Hood Element

The mini-element solved the compatibility condition problem, but the error estimation in equation is not optimal in the sense that

  1. the pressure space is sufficiently rich to enable a \(h^2\) convergence in the pressure error,

  2. but the velocity space is not rich enough to ensure a \(h^2\) convergence in the velocity error.

The idea of the Taylor-Hood element is to enrich even more the velocity space to ensure optimal convergence in \(h\).

Here we will take \([\poly{P}_2\)^d] for the velocity and \(\poly{P}_1\) for the pressure.

Introduce \[\begin{aligned} \label{eq:chmixte:39} X_h &=& [P2_{c,h}]d\\ M_h &=& P^1_{c,h} \end{aligned} \]

The spaces \(X_h\) and \(M_h \cap L^2_0(\Omega)\) satisfy the compatibility condition  uniformly in \(h\).

Suppose that \((u,p)\), solution of problem [prob:chmixte:1], is smooth enough, ie. \(u \in [H^3(\Omega)\)^d \cap [H1_0(\Omega)]d] and \(p\in H^2(\Omega) \cap L^2_0(\Omega)\).

Then there exists a constant \(c\) such that for all \(h >0\)

\[ \| u- u_h \|{1,\Omega} + \|p-p_h\|{0,\Omega} \leq c h^2 (\|u\|{3,\Omega} + \|p\|{2,\Omega}) \] and if the Stokes problem is stabilizing then

\[ \|u-u_h\|{0,\Omega} \leq c h^3 ( \|u\|{3,\Omega} +\|p\|_{2,\Omega}). \]

Generalized Taylor-Hood element

We consider the mixed finite elements \(\poly{P}_k/\poly{P}_{k-1}\) and \(\poly{Q}_k/\poly{Q}_{k-1}\) which allows to approximate the velocity and pressure respectively with, on Simplices \[\begin{aligned} \label{eq:chmixte:42} X_h &=& [P{k}_{c,h}]d\\ M_h &=& P^{k-1}_{c,h} \end{aligned}\]] On Hypercubes \(\[\begin{aligned} \label{eq:chmixte:43} X_h &=& [Q^{k}_{c,h}\)^d\\ M_h &=& Q^{k-1}_{c,h} \end{aligned} \] We then have

\[ \|u-u_h\|{0,\Omega} + h ( \| u- u_h \|{1,\Omega} + \|p-p_h\|{0,\Omega} ) \leq c h^{k+1} (\|u\|{k+1,\Omega} +\|p\|_{k,\Omega}) \]

There are other stable discretization spaces

  • Discrete inf-sup condition: dictates the choice of spaces

  • Inf-sup stables spaces:

    • \(\mathbb Q_k\)-\(\mathbb Q_{k-2}\), \(\mathbb Q_k\)-\(\mathbb Q^{disc}_{k-2}\)

    • \(\mathbb P_k\)-\(\mathbb P_{k-1}\), \(\mathbb P_k\)-\(\mathbb P_{k-2}\), \(\mathbb P_k\)-\(\mathbb P^{disc}_{k-2}\)

    • Discrete inf-sup constant independent of \(h\), but dependent on \(k\)

Numerical validation: Test case

We consider the Kovasznay solution of the steady Stokes equations.

The exact solution reads as follows

\[\begin{array}{r c l} \mathbf{u}(x,y) & = & \left(1 - e^{\lambda x } \cos (2 \pi y), \frac{\lambda}{2 \pi} e^{\lambda x } \sin (2 \pi y)\right)^T \\ p(x,y) & = & -\frac{e^{2 \lambda x}}{2} \\ \lambda & = & \frac{1}{2 \nu} - \sqrt{\frac{1}{4\nu^2} + 4\pi^2}. \end{array}\]

The domain is defined as \(\domain = (-0.5,1) \times (-0.5,1.5)\) and \(\nu = 0.035\).

The forcing term for the momentum equation is obtained from the solution and is

\[ \mathbf{f} = \left( e^{\lambda x} \left( \left( \lambda^2 - 4\pi^2 \right) \nu \cos (2\pi y) - \lambda e^{\lambda x} \right), e^{\lambda x} \nu \sin (2 \pi y) (-\lambda^2 + 4 \pi^2) \right)^T\]

Dirichlet boundary conditions are derived from the exact solution.

Solving Linear Systems

documentation pending

Appendix A: Bibliography

  1. raviart-thomas-1983

  2. gmsh

  3. raviartthomast1983


1. la transformation et son inverse sont \(\mathcal{C}^1\) et bijectives
2. An euclidian space which is complete for the norm induced by the scalar product