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

\( \newcommand{\Feel}{\textsc{Feel++}} \newcommand{\eg}{\emph{e.g.}} \newcommand{\ie}{\emph{i.e.}} \newcommand{\ds}{\displaystyle} \newcommand{\cf}{{\it cf \hspace*{1 mm}}} \newcommand{\bdm}{\begin{displaymath}} \newcommand{\edm}{\end{displaymath}} \newcommand{\bi}{\begin{itemize}} \newcommand{\ei}{\end{itemize}} \newcommand{\be}{\begin{equation}} \newcommand{\ee}{\end{equation}} \newcommand{\diver}{\hbox{div}} \newcommand{\hun}{\hbox{H}^1(\Omega)} \newcommand{\czero}{\hbox{C}^0(\bar{\Omega})} \newcommand{\dpa}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\NN}{\hbox{I\hspace*{-2 mm} N}} \newcommand{\CC}{\mathbf{C}} \newcommand{\RR}{{\mathbb{R}}} \newcommand{\Rn}{\RR$^n$} \newcommand{\ZZ}{{\bf Z}} \newcommand{\dd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\saut}{\vspace*{5 mm}\\} \newcommand{\Th}{\mathcal{T}_h} \newcommand{\Thi}{\mathcal{T}_h^i} \newcommand{\Thb}{\mathcal{T}_h^b} \newcommand{\Rh}{\mathcal{R}_h} % simplicial conformal refinement \newcommand{\Sh}{\mathcal{S}_h} % stencil \newcommand{\Nh}{\mathcal{N}_h} % neighbours of neighbours \newcommand{\Fh}{\mathcal{F}_h} \newcommand{\Fhi}{\mathcal{F}_h^i} \newcommand{\Fhb}{\mathcal{F}_h^b} \newcommand{\jump}[1]{[\![ #1 ]\!]} \newcommand{\Nma}{{{N_{\mathrm{ma}}}}} \newcommand{\Nso}{{{N_{\mathrm{so}}}}} \newcommand{\Nno}{{{N_{\mathrm{no}}}}} \newcommand{\Ne}{{{N_{\mathrm{e}}}}} \newcommand{\Ilag}[1]{{{\mathcal{I}^{\mathrm{lag}}_{#1}}}} \newcommand{\diam}{{{\operatorname{diam}}}} \renewcommand{\max}{{{\operatorname{max}}}} \newcommand{\essinf}{{{\operatorname{ess}\, \operatorname{inf}}}} \newcommand{\calTh}{{{\mathcal{T}_h}}} \newcommand{\Ck}[1]{{{\mathcal{C}^#1}}} \newcommand{\Ck}[1]{{{\mathcal{C}^#1}}} \newcommand{\Pk}[1]{{{{\mathbb{P}_#1}}}} \newcommand{\Qk}[1]{{{{\mathbb{Q}_#1}}}} \newcommand{\Pch}[1]{{{{P^#1_{c,h}}}}} \newcommand{\Pcho}[1]{{{{P^#1_{c,h,0}}}}} \newcommand{\Qch}[1]{{{{Q^#1_{c,h}}}}} \newcommand{\Ich}[2]{{{{\mathcal{I}^#1_{c,h}}#2}}} \newcommand{\set}[1]{\left\{#1\right\}} \newcommand{\Next}{{{\mathrm{n}}}} \newcommand{\nf}{{{n_f}}} \newcommand{\ngeo}{{{n_{\mathrm{geo}}}}} \newcommand{\geo}{{{\mathrm{geo}}}} \renewcommand{\dim}[1]{{{\operatorname{dim}(#1)}}} \newcommand{\card}[1]{{{\operatorname{card}(#1)}}} \newcommand{\poly}[1]{{{\mathbb{#1}}}} \newcommand{\R}[1][d]{{\mathbb{R}^{#1}}} \newcommand{\tr}{{{\operatorname{tr}}}} \newcommand{\Id}{\mathcal{I}} \newcommand{\disp}[1][u]{{\mathbf{#1}}} \newcommand{\stresst}[1][(\mathbf{u})]{{\mathbf{\sigma#1}}} \newcommand{\deformt}[1][(\mathbf{u})]{{\mathbf{\varepsilon#1}}} \newcommand{\domain}{\Omega} \newcommand{\prect}[1]{\left\(#1\right\)} %%% Local Variables: %%% mode: latex %%% TeX-master: "mef-intro" %%% End: \)

1.1. Quelques rappels

1.1.1. Normes et produits scalaires

Soit \(E\) un espace vectoriel.

Definition
\(\|.\|\) : \(E \rightarrow \RR\) est une norme sur \(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| \]

Definition
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 \| }\)

Example
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.
Definition
Un espace vectoriel muni d’une norme est appelé espace normé.
Definition
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

Definition
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\]
Theorem
Toute suite convergente est de Cauchy.
La réciproque de Theorem est fausse.
Definition
Un espace vectoriel est complet ssi toute suite de Cauchy y est convergente.
Definition
Un espace normé complet est un espace de Banach.
Definition
Un espace préhilbertien complet est un espace de Hilbert.
Definition
Un espace de Hilbert de dimension finie est appelé espace euclidien.

1.1.3. Espaces fonctionnels

Definition
Un espace fonctionnel est un espace vectoriel dont les éléments sont des fonctions.
Example
\({\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\).

Example
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)\) :

Definition
\(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.
Theorem
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

Theorem
\(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\).

Definition
On appelle support de \(\bf \varphi\) l’adhérence de \(\{ x \in \Omega / \varphi(x) \ne 0 \}\).
Example
Pour \(\Omega = ]-1,1\[\), et \(\varphi\) la fonction constante égale à 1, \(\hbox{Supp}\, \varphi = [-1,1\)].
Definition
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.
Example
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)\).

Theorem
\(\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.

Definition
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' \]
Example
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.
Definition
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\).

Theorem
Quand elle existe, la dérivée généralisée est unique.
Theorem
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\)

Definition
\({ 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.

Definition
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\}\] </div>+

\(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\} \]

Theorem
\(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 stem:[(.,.)_0] le produit scalaire stem:[L^2]. On notera stem:[\|.\|_1] la norme associée à stem:[(.,.)_1]. </div>+

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}]\)

Theorem
\(H^m(\Omega)\) muni du produit scalaire \((.,.)_m\) est un espace de Hilbert.[thr:8]
Theorem
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} }\)
Example
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)\)

Definition
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)
Theorem
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\).

Definition
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} } \]
Theorem
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.

Definition
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 :

\[ \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 :

\[ ({\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.

Definition
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 \]
Definition
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. \]
Definition
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. \]
Theorem
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.
Proof
La démonstration générale de ce théorème peut être trouvée par exemple dans cite:[raviart-thomas-1983].
Theorem
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.
Proof
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 (equation) 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 equation 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 equation 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 equation.

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 \] ()
Definition
Les intervalles \(I_i\) sont appelées de mailles ou éléments ou cellules du maillage, on a noté \(\Nma\) le nombre de maillage
Definition
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.

Definition
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:

Theorem
\(\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

Proposition
. La famille \(\{\varphi_1,...,\varphi_\Nso\}\) est une base de \(\Pch{1}\). . La famille \(\{\varphi_2,...,\varphi_{\Nso-1}\}\) est une base de \(\Pcho{1}\).
Corollary
\(\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

Proposition
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,

Proposition
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 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 equation 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 equation, 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

Proposition
Soit un entier \(k\geq 1\). On suppose que la solution du problème equation est dans \(H^{k+1}(\Omega)\). On note \(u_h\) la solution du problème approché equation 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 est optimale car elle est du même ordre que l’erreur d’interpolation en norme \(H^1\), voir la proposition Proposition.

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

Proposition
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 est optimale car elle est du même ordre que l’erreur d’interpolation en norme \(L^2\), voir la proposition Proposition.

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

Definition
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.

Proof
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

Definition
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\)).
Definition
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\).

Definition
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

Definition
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\).
Property
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}\).
Definition
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

Definition
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}\).

Definition
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].

Definition
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

Definition
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.

Definition
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 \(\hatglossterm::[\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.
Definition
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.

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} \] ()
Notation
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.
Proposition
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{C}^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{C}^0(K),\quad \Ilag{K}( v \circ T_K ) = \Ilag{K}( v ) \circ T_K. \]
Theorem
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

Definition
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

\[ 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\} \] ()
\[ 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

Problem
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} \]
Problem
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

Definition
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) \]

Theorem
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 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\).

Definition
Si \(\hat{P}_{\mathrm{geo}} = \hat{P}\) l’approximation est dite iso-parametrique.
Definition
Si \(\hat{P}_{\mathrm{geo}} \subset \hat{P}\) l’approximation est dite sub-parametrique.
Definition
Si \(\hat{P} \subset\hat{P}_{\mathrm{geo}}\) l’approximation est dite sur-parametrique.
Gmsh cite:[gmsh] un mailleur libre permet de générer des maillage d’ordre élevé jusqu’à l’ordre 5 en 2D et 4 en 3D
Theorem
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 Theorem.

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
// genere le maillage pass:[\(\Omega^1_h\)] associe une transformation affine
auto mesh1 = unitCircle<1>();
auto area1 = integrate( _range=elements(mesh1),
                        _expr=cst(1.)).evaluate().norm();
std::cout << "|pi-surface<1>| = " << math::abs(pi-area1)
          << "\n";

// genere le maillage pass:[\(\Omega^2_h\)] associe une transformation d'ordre 2
auto mesh2 = unitCircle<2>();
auto area2 = integrate( _range=elements(mesh2),
                        _expr=cst(1.)).evaluate().norm();
std::cout << "|pi-surface<2>| = " << math::abs(pi-area2)
          << "\n";

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 appliqué à des pentes à l’élément fini géométrique.

Table 1. Convergence de l’approximation \(P_1\) et \(P_2\)
h error \(P_1\) error \(P_2\)

0.4000

1.41593e-01

4.87932e-04

0.2000

2.62996e-02

1.65709e-05

0.1000

5.73876e-03

7.86831e-07

0.0500

1.34419e-03

4.31428e-08

0.0250

3.36079e-04

2.69659e-09

La table Convergence de l’approximation \(P_1\) et \(P_2\) 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

\[ {\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 (equation) é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

Theorem
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 Theorem

Corollary
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 Theorem

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 \textglossterm::[\tiny fois]}, \ldots , \underbrace{e_n,\ldots,e_n}_{\alpha_n \textglossterm::[\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 (equation) grâce à la ([eq), 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 ([eq) ce qui conclut la preuve de Theorem \(\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 :

Theorem
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 Theorem

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 (equation), 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 à cite:[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 cite:[raviartthomast1983].

Approximation de problèmes coercifs

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

  • le Laplacien et des variantes sur les conditions aux limites et le l’équation elle-même;

  • l’élasticité linéaire qui permet de modéliser de petites déformations d’un milieu continu déformable.

1. Le Laplacien

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

Problem
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 Problem est la suivante:

Problem
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 Problem s’écrit sous forme abstraite

Problem
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:

Lemma
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

Proposition
La forme bilinéaire \(a\) du problème Problem est coercive
Proof
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 Theorem permet alors de conclure sur l’existence d’une solution unique pour le problème Problem.

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 Theorem.

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

Problem
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:

Theorem
On suppose que \(u\) solution de Problem 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} \]
Proof
La preuve de ([eq) est obtenue grâce au lemme de Cea (equation) et du théorème d’interpolation ([eq). La preuve de ([eq) 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) se déduit de ([eq).

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) 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 pass:[\(\R\)] telle que la transformation
  // géométrique pass:[\(T_K,\ K \in \calTh\)]  soit affine
  typedef Mesh<Simplex<d,1> > mesh_type;

  // un maillage d'hypercube dans pass:[\(\R\)] telle que la transformation
  // géométrique pass:[\(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 pass:[\([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 Problem

  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:

Problem
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:

Problem
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.

Theorem
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:

Problem
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

Problem
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 Problem est bien posé grâce à Lax-Milgram.

On peut approcher le problème Problem par des éléments finis de Lagrange.

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

Problem
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 Problem est identique par le théorème Theorem 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

Problem
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 Problem 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 Problem s’écrit alors

Problem
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 Problem résulte de la coercivité de \(a\) sur \(H^1_*(\Omega)\) qui elle-même résulte du Lemme suivante

Lemma
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 Problem peut être approché par des éléments finis de Lagrange ce qui conduit au problème discret suivant

Problem
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 Problem 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:

Problem
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

Problem
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 Problem est bien posé grâce à Lax-Milgram.

On considère le problème discret suivant

Problem
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 Theorem.

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:

Problem
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:

Definition
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

Problem
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

Proposition
Si \(u\) résout [eq, alors \(\mathcal{L}u=f\) presque partout dans \(\Omega\) et \(u=0\) presque partout sur \(\partial \Omega\).
Proof
Soit \(\phi\in \mathcal{D}(\Omega)\) et \(u\) une solution de Problem. 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\).

Problem
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 \]
Proposition
Soit \(g\in H^{\frac{1}{2}}(\partial \Omega)\), si \(u\) résout Problem, alors \(\mathcal{L}u=f\) presque partout dans \(\Omega\) et \(u=g\) presque partout sur \(\partial \Omega\).
Proof
La preuve est similaire à celle de Problem.
Lorsque l’opérateur \(\mathcal{L}\) est le Laplacien, Problem 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, on obtient la formulation faible suivante:

Problem
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 \]
Proposition
Soit \(g\in L^2(\partial \Omega)\), si \(u\) résout Problem, alors \(\mathcal{L}u=f\) presque partout dans \(\Omega\) et \(\alpha \nabla u \cdot n = g\) presque partout sur \(\partial \Omega\).
Proof
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, nous avons la formulation faible suivante

Problem
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 \]
Proposition
Soit \(g\in L^2(\partial \Omega)\) et \(\gamma \in L^\infty(\partial \Omega\). Si \(u\) résout Problem, alors \(\mathcal{L}u=f\) presque partout dans \(\Omega\) et \(\gamma u + \alpha \nabla u \cdot n = g\) presque partout sur \(\partial \Omega\).
Proof
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 [eqpbadr] 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
Theorem
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 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. \]
Proof
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}{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 \(\sigma = \min(\alpha_0,p). Si stem:[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

Problem
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\).

Theorem
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} \]
Theorem
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

  • \(\epsilon\) is the diffusion coefficient (matrix \(d \times d\)), \(\epsilon(\mathbf{x}) > 0\)

  • \(\mathbf{\beta}\) velocity and \(\nabla \cdot \mathbf{ \beta} \in L^2(\Omega)\)

  • \(\mu > 0\) reaction (or absorption) coefficient.

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

Problem
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.
Problem
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\]
Theorem
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:

Problem
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 [[eq:stab21]] \[ \mathcal{L}_h(u,f;v_h) = 0, \forall v_h \in V_h \] </div>+

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

Theorem
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:

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
Definiton
If the bilinear form \(a\) is symmetric and positive on \(X\times X\), we say that Problem 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 Problem.

Denote \(V=X\times M\) and introduce \(c \in \mathcal{L}(V\times V, \RR)\) such that

\[ c((u,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 Problem reads

Problem
Look for \((u,p) \in V\) such that \[ \begin{array}[c]{rl} c((u,p), (v,q)) & = h(v,q), \quad \forall (v,q) \in V \end{array} \]
Theorem
We suppose that \(a\) is coercive over \(X\), the Problem 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 Problem 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) \]

Definition
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) \]
Proposition
Under the hypothesys oF Theorem, the Lagrangian \(l\) defined by has a unique saddle point. Moreover this saddle point is the unique solution of problem Problem.

7.4. Finite element approximation

7.4.1. Abstract Discrete Problem

We now turn to the approximation of the problem Problem 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:

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. \]
Theorem
We suppose that \(a\) is coercive over \(X\) and that \(X_h \subset X\) and \(M_h \subset M\). Then the Problem 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 Problem, 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

Lemma
Thanks to the Lemma of Céa applied to Saddle-Point Problems, the unique solution \((u,p)\) of problem Problem 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 Problem 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

Definition
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

Problem
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 [H^1_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 Problem 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)\)

Problem
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 Theorem 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 [Q^1_{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 [P^1_{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.

Definition
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

Aligned
X_h &=& \Big\{ v_h \in [C0(\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}
Lemma
The spaces \(X_h\) and \(M_h \cap L^2_0(\Omega)\) satisfy the compatibility condition  uniformly in \(h\).
Theorem
Suppose that \((u,p)\), solution of Problem, 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}). \]
Definition
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 &=& [P^2_{c,h}]^d\\ M_h &=& P^1_{c,h} \end{aligned} \]

Lemma
The spaces \(X_h\) and \(M_h \cap L^2_0(\Omega)\) satisfy the compatibility condition  uniformly in \(h\).
Theorem
Suppose that \((u,p)\), solution of problem Problem, 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 latexmath:[\[\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.

8. Weakly imposed Dirichlet boundary conditions: Nitsche method

8.1. Weak treatment

In order to treat the boundary conditions uniformly (i.e. the same way as Neumann and Robin Conditions), we wish to treat the Dirichlet BC (e.g. \(u=g\)) weakly.

Initial Idea add the penalisation term

\[\int_{\partial \Omega} \mu( u - g)\]

where \(\mu\) is a constant. But this is not enough: this is not consistent with the initial formulation.

One can use the Nitsche method to implement weak Dirichlet conditions and follows the next steps:

  • write the equations in conservative form (i.e. identify the flux);

  • add the terms to ensure consistency (i.e the flux on the boundary);

  • symmetrize to ensure adjoint consistency;

  • add a penalisation term with factor \(\gamma (u-g)/h\) that ensures that the solution will be set to the proper value at the boundary;

8.2. Penalisation parameter

8.2.1. Choosing \(\gamma\)

\(\gamma\) must be chosen such that the coercivity(or inf-sup) property is satisfied. It is difficult to do in general. We increase \(\gamma\) until the weak Dirichlet boundary conditions are properly satisfied, e.g. start with \(\gamma=1\), typical values are between 1 and 10.

The choice of \(\gamma\) is a problem specially when \(h\) is small.

8.2.2. Advantages, disadvantages

We compare the advantages and disadvantages of strong and weak Dirichlet boundary conditions treatment.

We start with the weak conditions

  • advantage uniform(weak) treatment of all boundary conditions type

  • advantage if boundary condition is independant of time, the terms are assembled once for all

  • disadvantage Introduction of the penalisation parameter \(\gamma\) that needs to be tweaked

Strong treatment: Advantages

  • advantage Enforce exactely the boundary conditions

  • disadvantage Need to modify the matrix once assembled to reflect that the Dirichlet degree of freedom are actually known. Then even if the boundary condition is independant of time, at every time step if there are terms depending on time that need reassembly (e.g. convection) the strong treatment needs to be reapplied.

  • disadvantage it can be expensive to apply depending on the type of sparse matrix used, for example using CSR format setting rows to 0 except on the diagonal to 1 is not expensive but one must do that also for the columns associated with each Dirichlet degree of freedom and that is expensive.

8.3. Example: Laplacian

We look for \(u\) such that

-\nabla\cdot( k \nabla u )= f\ \mbox{in} \Omega,\quad u=g|_{\partial \Omega}

\[\begin{gathered} \label{eq:51} \int_\Omega k \nabla u \cdot \nabla v + \int_{\partial \Omega} \underbrace{-k \frac{\partial u}{\partial n}v}_{\text{integration by part}} \underbrace{-k \frac{\partial v}{\partial n} u}_{\text{adjoint consistency: symetrisation}} \\ + \underbrace{\frac{\gamma}{h} u v}_{\text{penalisation: enforce Dirichlet condition}} =\\ \int_\Omega f v + \int_{\partial \Omega} (\underbrace{-k \frac{\partial v}{\partial n} }_{\text{adjoint consistency}} + \underbrace{\frac{\gamma}{h} v) g}_{\text{penalisation: enforce Dirichlet condition}} \end{gathered}\]

8.3.1. Implementation

// contribution to bilinear form (left hand side)
form2( _test=Xh, _trial=Xh ) +=
integrate( boundaryfaces(mesh),
           // integration by part
           -(gradt(u)*N())*id(v)
           // adjoint consistency
           -(grad(v)*N())*idt(u)
           // penalisation
           +gamma*id(v)*idt(u)/hFace());
// contribution to linear form (right hand side)
form1( _test=Xh ) +=
integrate( boundaryfaces(mesh),
           // adjoint consistency
           -(grad(v)*N())*g
           // penalisation
           +gamma*id(v)*g/hFace());

8.4. Example: Convection-Diffusion

Convection Diffusion Consider now the following problem, find \(u\) such that

\[\nabla \cdot ( -\epsilon\nabla u + \mathbf{c} u ) = f,\quad u = g|_{\partial \Omega},\quad \nabla \cdot \mathbf{c} = 0\]

the flux vector field is \(\mathbf{F}=-\nabla u\).

Note that here the condition, \(\nabla \cdot \mathbf{c} = 0\) was crucial to expand \(\nabla \cdot (\mathbf{c} u )\) into \(\mathbf{c} \cdot \nabla u\) since

\[\nabla \cdot (\mathbf{c} u ) = \mathbf{c} \cdot \nabla u + \underbrace{u \nabla \cdot \mathbf{c}}_{=0}\]

Weak formulation for convection diffusion Multiplying by any test function \(v\) and integration by part of ([eq:2]) gives

\[\int_\Omega \epsilon \nabla u \cdot \nabla v + (\mathbf{c} \cdot \nabla u)v + \int_{\partial \Omega} (\mathbf{F}\cdot \mathbf{n}) v = \int_\Omega f v\]

where \(\mathbf{n}\) is the outward unit normal to \(\partial \Omega\).

We now introduce the penalisation term that will ensure that \(u \rightarrow g\) as \(h \rightarrow 0\) on \(\partial \Omega\). ([eq:4]) reads now

\[\int_\Omega \epsilon \nabla u \cdot \nabla v + (\mathbf{c} \cdot \nabla u)v + \int_{\partial \Omega} (\mathbf{F}\cdot \mathbf{n}) v + {\frac{\gamma}{h} u v} = \int_\Omega f v + {\int_{\partial \Omega} \frac{\gamma}{h} g v}\]

Finally we can incorporate the symetrisation

\[\begin{gathered} \int_\Omega \epsilon \nabla u \cdot \nabla v + (\mathbf{c} \cdot \nabla u)v + \int_{\partial \Omega} ((-\epsilon \nabla u)\cdot \mathbf{n}) v+ {(-\epsilon\nabla v\cdot \mathbf{n} + {\mathbf{c}}\cdot \mathbf{n} v ) u} + \frac{\gamma}{h} u v = \\ \int_\Omega f v + \int_{\partial \Omega} {(-\epsilon\nabla v\cdot \mathbf{n} + {\mathbf{c}}\cdot \mathbf{n} v) g}+ \frac{\gamma}{h} g v \end{gathered}\]

8.4.1. Implementation

// bilinear form (left hand side)
form2( _trial=Xh, _test=Xh ) +=
integrate( boundaryfaces(mesh),
  // integration by part
  -(pass:[\(\epsilon\)] gradt(u)*N())*id(v) + (idt(u)*trans(idv(c))*N())*id(v)
  // adjoint consistency
  -(pass:[\(\epsilon\)] grad(v)*N())*idt(u) + (idt(u)*trans(idv(c))*N())*id(v)
  // penalisation
  +gamma*id(v)*idt(u)/hFace());
// linear form (right hand side)
form1( _test=Xh ) +=
integrate( boundaryfaces(mesh),
  // adjoint consistency
  -(pass:[\(\epsilon\)] grad(v)*N())*g
  + g*trans(idv(c))*N())*id(v)
  // penalisation
  +gamma*id(v)*g/hFace());

Solving Linear Systems

documentation pending

Appendix A: Bibliography

bibliography::[]


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