Calcul des Matrices Élémentaires

$\newcommand{\Cb}{\mathbb{C}}$ $\newcommand{\Nb}{\mathbb{N}}$ $\newcommand{\Pb}{\mathbb{P}}$ $\newcommand{\Rb}{\mathbb{R}}$ $\newcommand{\PS}[2]{\left(#1,#2\right)}$ $\newcommand{\PSV}[2]{\PS{#1}{#2}_V}$ $\newcommand{\PSL}[2]{\PS{#1}{#2}_{L^2(\Omega)}}$ $\newcommand{\PSH}[2]{\PS{#1}{#2}_{H^1(\Omega)}}$ $\newcommand{\norm}[1]{\left\|#1\right\|}$ $\newcommand{\normV}[1]{\left\|#1\right\|_{V}}$ $\newcommand{\normH}[1]{\left\|#1\right\|_{H^1(\Omega)}}$ $\newcommand{\normL}[1]{\left\|#1\right\|_{L^2(\Omega)}}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\ee}{\mathbf{e}}$ $\newcommand{\nn}{\mathbf{n}}$ $\newcommand{\ssb}{\mathbf{s}}$ $\newcommand{\sK}{\mathbf{s}^K}$ $\newcommand{\xx}{\mathbf{x}}$ $\newcommand{\yy}{\mathbf{y}}$ $\newcommand{\zz}{\mathbf{z}}$ $\newcommand{\Ccal}{\mathcal{C}}$ $\newcommand{\Cscr}{\mathscr{C}}$ $\newcommand{\Sscr}{\mathscr{S}}$ $\newcommand{\Tscr}{\mathscr{T}}$ $\newcommand{\omegai}{\omega_i}$ $\newcommand{\dsp}{\displaystyle}$ $\newcommand{\diff}{{\rm d}}$ $\newcommand{\conj}[1]{\overline{#1}}$ $\newcommand{\dn}{\partial_\nn}$ $\newcommand{\supp}{\mathrm{supp}}$ $\newcommand{\restrict}{\mathclose{}|\mathopen{}}$ $\newcommand{\enstq}[2]{\left\{#1 \mathrel{}\middle|\mathrel{}#2\right\}}$ $\newcommand{\Image}{\mathrm{Im}}$ $\newcommand{\Ker}{\mathrm{Ker}}$ $\newcommand{\dxi}{\partial_{x_i}}$ $\newcommand{\di}{\partial_{i}}$ $\newcommand{\dj}{\partial_{j}}$ $\newcommand{\dxj}{\partial x_{j}}$ $\newcommand{\Ho}{H^1(\Omega)}$ $\newcommand{\Lo}{L^2(\Omega)}$ $\newcommand{\Cinfc}{\Cscr^{\infty}_c}$ $\newcommand{\CinfcO}{\Cinfc(\Omega)}$ $\newcommand{\hme}[1]{#1_h}$ $\newcommand{\vh}{v_h}$ $\newcommand{\Vh}{V_h}$ $\newcommand{\uh}{u_h}$ $\newcommand{\Nh}{N_h}$ $\newcommand{\ui}{u_i}$ $\newcommand{\uj}{u_j}$ $\newcommand{\Sscrh}{\hme{\Sscr}}$ $\newcommand{\deltaij}{\delta_{i,j}}$ $\newcommand{\Kp}{K_p}$ $\newcommand{\Kq}{K_q}$ $\newcommand{\Kl}{K_\ell}$ $\newcommand{\Pun}{\Pb_1}$ $\newcommand{\Punw}{\Pun(\omega)}$ $\newcommand{\grandO}[1]{O\left(#1\right)}$ $\newcommand{\qh}{\widehat{\qq}}$ $\newcommand{\sh}{\widehat{\ssb}}$ $\newcommand{\phih}{\widehat{\phi}}$ $\newcommand{\Meh}{\widehat{M}^e}$ $\newcommand{\varphih}{\widehat{\varphi}}$ $\newcommand{\psih}{\widehat{\psi}}$ $\newcommand{\gh}{\widehat{g}}$ $\newcommand{\Kh}{\widehat{K}}$ $\newcommand{\sumit}[1]{\ssb_{#1}}$ $\newcommand{\sumitK}[2]{\ssb_{#2}^{#1}}$ $\newcommand{\tri}[1]{K_{#1}}$ $\newcommand{\loctoglob}{\mathrm{Loc2Glob}}$ $\newcommand{\aK}[1]{a_{#1}}$ $\newcommand{\ellK}[1]{\ell_{#1}}$ $\newcommand{\Ns}{N_s}$ $\newcommand{\mphi}[1]{\varphi_{#1}}$ $\newcommand{\mphiK}[2]{\mphi{#2}^{#1}}$ $\newcommand{\xK}[2]{x^{#1}_{#2}}$ $\newcommand{\yK}[2]{y^{#1}_{#2}}$ $\newcommand{\TK}[1]{T^{\tri{#1}}}$ $\newcommand{\BK}[1]{B_{\tri{p}}}$ $\newcommand{\JK}[1]{J_{\tri{#1}}}$ $\newcommand{\Me}[1]{M^e_{{#1}}}$ $\newcommand{\Mep}{\Me{p}}$ $\newcommand{\De}[1]{D^e_{{#1}}}$ $\newcommand{\Dep}{\De{p}}$

Matrices de Masse et de Rigidité

Nous devons maintenant calculer les contributions élémentaires, c’est à dire les quantités $\aK{p}(\mphiK{p}{j}, \mphiK{p}{i})$ et $\ellK{p}(\mphiK{p}{i})$. Pour simplifier, les sommets d’un triangle $\tri{p}$ seront notés $(\sumitK{p}{1}, \sumitK{p}{2}, \sumitK{p}{3})$ et ordonnés dans le sens trigonométrique. Nous noterons $\sumitK{p}{i}=(\xK{p}{i},\yK{p}{i})$ et $\mphiK{p}{j}$ une fonction de forme du triangle $\tri{p}$ sans la supposer linéaire. D’autre part, nous scindons les contributions élémentaires en deux parties :

  • La matrice de masse, notée $M$ de coefficient $M(I,J)$ donné par $$ M(I,J) = \int_{K}\mphi{J}(\xx)\overline{\mphi{I}(\xx)}\diff\xx. $$
  • Matrice de rigidité, notée $D$ de coefficient $D(I,J)$ donné par $$ D(I,J)= \int_{K}\nabla\mphi{J}(\xx)\overline{\nabla\mphi{I}(\xx)}\diff\xx. $$

Deux remarques :

  • La matrice $M$ représente l’opérateur identité dans la base des éléments finis (ce n’est pas la matrice diagonale remplie de 1…). En effet, si “l’EDP” était juste $u = f$, alors on aurait $a(u,v) = \int_{\Omega}u\overline{v}$ et sa matrice $\Pun$ serait la matrice de masse.
  • Dans la littérature, la matrice de rigidité est souvent notée $K$, mais nous l’appelons $D$ pour éviter toute confusion avec les triangles, nommés $K$ également.

Matrice de masse élémentaire

Nous nous focalisons sur la matrice de masse, le principe est similaire pour la matrice $D$ et est détaillé juste après.

Pour construire la matrice $M$, nous avons vu qu’il était préférable de parcourir les triangles plutôt que les sommets, autrement dit, plutôt que de calculer $M(I,J)$ directement, mieux vaut calculer, pour tout triangle $p$, la contribution élémentaire $\Mep(i,j)$ pour $i,j = 1,2,3$, définie par : \begin{equation} \label{eq:matelem} \Mep(i,j) = \int_{\tri{p}} \mphiK{p}{j} \overline{\mphiK{p}{i}}. \end{equation} Chaque contribution élémentaire $\Mep(i,j)$ est ensuite ajoutée à $M(I,J)$, avec $I=\loctoglob(p,i)$ et $J=\loctoglob(p,j)$, de sorte qu’une fois toutes les contributions calculées, nous retrouvons bien la matrice de masse. Nous nous focalisons ainsi maintenant sur le calcul de \eqref{eq:matelem}.

La quantité $\Mep$ est une matrice $3\times 3$ appelée matrice de masse élémentaire.

Triangle de référence

Plaçons nous dans un triangle “simple” $\Kh$, appelé triangle de référence, souvent choisi comme étant le triangle rectangle de sommets $\sh_{1}=(0,0)$, $\sh_{2}=(1,0)$ et $\sh_{3}=(0,1)$, ordonnés dans le sens trigonométrique. Pour différencier ce triangle d’un triangle du maillage, nous lui adjoignons un repère $(\xi,\eta)$ dit repère paramétrique.

Triangle de référence $\Kh$ et son repère paramétrique $(\xi,\eta)$.

Plutôt que d’indicer par $p$, nous notons $\varphih_i \in \Pun(\Kh)$ les trois fonctions de forme associés aux sommets $\sh_{i}$, pour $i=1,2,3$, définies par $$ \varphih_j(\sh_{i}) = \delta_{ij}. $$ Dans ce triangle $\Kh$, la contribution élémentaire $\Meh(i,j)$ pour $i,j = 1,2,3$ est donnée par $$ \Meh(i,j) = \int_{\Kh}\varphih_j\overline{\varphih_i}. $$ Ces fonctions $\varphih_j$ étant des polynômes de degré un, nous pouvons donc calculer analytiquement leurs coefficients : $$ \left\{ \begin{array}{l} \varphih_1(\xi,\eta) = 1-\xi-\eta\\
\varphih_2(\xi,\eta) = \xi\\
\varphih_3(\xi,\eta) = \eta\\
\end{array} \right. $$ Par suite, nous en déduisons les valeurs de $\Meh(i,j)$ pour $i,j=1,2,3$ :

Lemma.

Pour $i,j=1,2,3$, avec $i=j$ : $$ \int_{\Kh} \varphih_j\overline{\varphih_j} \diff(x,y) =\frac{1}{12}, $$ et pour $j\neq i$: $$ \int_{\Kh} \varphih_j\overline{\varphih_i} \diff(x,y) =\frac{1}{24}. $$

Proof.

Prenons tout d’abord le cas $i=j=2$, soit $\overline{\varphih_i} = \varphih_j = \varphih_2(\xi,\eta) = \xi$. Dans ce cas : $$ \int_{\Kh} \xi^2 \diff (\xi,\eta) = \int_0^1\int_0^{1-\xi} \xi^2 \diff\eta\diff\xi = \int_0^1(1-\xi)\xi^2\diff\xi = \left[\frac{\xi^3}{3} - \frac{\xi^4}{4}\right]_0^1=\frac{1}{3}-\frac{1}{4} = \frac{1}{12}. $$ Les calculs sont similaires pour $j=1$ et $j=3$.

Prenons maintenant $i\neq j$, par exemple $i=3$ et $j=2$ : $$ \int_{\Kh} \xi\eta \diff (\xi,\eta) = \int_0^1\left(\int_0^{1-\xi} \eta \diff\eta\right)\xi\diff\xi = \frac{1}{2}\int_0^1(1-\xi)^2\xi\diff\xi
= \frac{1}{2}\left[ \frac{1}{2} - \frac{2}{3} +\frac{1}{4}\right] =\frac{1}{24}. $$ Les calculs sont similaires avec $i=1$ et $j=2,3$.

La matrice de masse élémentaire $\Meh$ du triangle de référence $\Kh$ est donc donnée par $$ \Meh = \frac{1}{24}\left( \begin{array}{c c c} 2 & 1 & 1\\
1 & 2 & 1\\
1 & 1 & 2 \end{array} \right). $$

Triangle quelconque

Soit un triangle $\tri{p}$ du maillage et supposons que nous disposions d’une transformation bijective et linéaire $\TK{p}$ permetteant de transformer le triangle de référence $\Kh$ en $\tri{p}$ avec en plus $\TK{p}(\sh_i) = \sumitK{p}{i}$. Cette fonction $\TK{p}$ transforme les coordonnées paramétriques $(\xi,\eta)$ en coordonnées physiques $(x,y)$ avec $(x,y)=\TK{p}(\xi,\eta)\in\tri{p}$, et “conserve l’ordre des sommets”.

Passage du triangle de référence $\Kh$ vers un triangle $\tri{p}$ par la transformation $\TK{p}$.

Changement de coordonnées

Nous avons alors $\varphi_j^p(x,y) = \varphi_j^p(\TK{p}(\xi,\eta))$ avec :

  1. $\varphi_j^p\circ\TK{p}\in\Pun(\Kh)$
  2. $\varphi_j^p\circ\TK{p}(\sh_i) = \delta_{ij}$

En d’autres termes, nous avons $\varphi_j^p\circ\TK{p} = \varphih_j$ et donc la suite d’égalités suivantes : $$ \varphi_j^p(x,y) = \varphi_j^p(\TK{p}(\xi,\eta)) = \varphih_j(\xi,\eta). $$ En notant $\JK{p}$ la matrice Jacobienne de $\TK{p}$, alors la quantité $\Mep(i,j)$ peut alors s’écrire, par changement de variables : $$ \begin{array}{r c l} \Mep(i,j) &=& \dsp\int_{\tri{p}}\mphiK{p}{j}(x,y)\overline{\mphiK{p}{i}(x,y)} \diff(x,y)\\
&=&\dsp \det(\JK{p})\underbrace{\int_{\Kh}\varphih_{j}(\xi,\eta)\overline{\varphih_{i}(\xi,\eta)}\diff(\xi,\eta)}_{\text{Déjà calculé !}}\\
\end{array} $$

Pour calculer la matrice élémentaire d’un triangle $\tri{p}$ quelconque, nous avons maintenant uniquement besoin du déterminant de la Jacobienne : $\det(\JK{p})$.

Expression et Jacobienne de la transformation

La transformation que nous cherchons, $\TK{p}$, est linéaire et “conserve” les sommets et leur ordre. Pour obtenir son expression, nous construisons des fonctions d’interpolation géométrique, $(\psih_i)_{1\leq i \leq 3}$, linéaires sur $\Kh$ et telles que : $$ \forall i,j=1,2,3, \quad \psih_i(\sh_{j}) = \deltaij. $$ La transformation aura alors pour expression : $$ \begin{array}{r c c l} \TK{p}\colon & \Kh & \to & \tri{p}\\
& (\xi,\eta) & \mapsto & \TK{p}(\xi,\eta) = (x,y) = \psih_{1}(\xi,\eta) \sumitK{p}{1} + \psih_{2}(\xi,\eta) \sumitK{p}{2} + \psih_{3}(\xi,\eta) \sumitK{p}{3}. \end{array} $$

Les fonctions d’interpolation géométrique $\varphih_j$ sont ici identiques aux fonctions de forme $\varphih_j$, c’est pourquoi nous parlons d’éléments finis isoparamétriques. Cependant, il faut bien se rappeler que ce n’est pas obligatoire et que le choix des fonctions $\psih_j$ et $\varphih_j$ est découplé. En particulier, pour obtenir des éléments courbes, nous pourrions choisir par exemple des fonctions $\psih_j$ quadratiques.

Comme $\psih_j = \varphih_j$ pour tout $j=1,2,3$, nous disposons de leur expression analytique : $$ \left\{ \begin{array}{l} \psih_{1}(\xi,\eta) = 1 - \xi - \eta\\
\psih_{2}(\xi,\eta) = \xi\\
\psih_{3}(\xi,\eta) = \eta\\
\end{array} \right. $$

La matrice Jacobienne de la transformation est alors donnée par $$ \JK{p} = \left( \begin{array}{c c} \dsp\frac{\partial x}{\partial \xi} &\dsp \frac{\partial x}{\partial \eta} \\
\dsp\frac{\partial y}{\partial \xi} &\dsp \frac{\partial y}{\partial \eta} \end{array} \right) = \left( \begin{array}{c c} \xK{p}{2} - \xK{p}{1} & \xK{p}{3} - \xK{p}{1}\\
\yK{p}{2} - \yK{p}{1} & \yK{p}{3} - \yK{p}{1} \end{array} \right) $$ Le déterminant de cette matrice vaut $$ \det(\JK{p}) = (\xK{p}{2}-\xK{p}{1})(\yK{p}{3}-\yK{p}{1}) - (\xK{p}{3}-\xK{p}{1})(\yK{p}{2}-\yK{p}{1}). $$ En valeur absolue, cela correspond à 2 fois l’aire du triangle $\tri{p}$. Le déterminant est donc non nul puisque le triangle n’est pas dégénéré et la transformation $\TK{p}$ est donc bien inversible.

Expression finale de la matrice élémentaire

Lemma.

Pour $i,j=1,2,3$ : $$ \Mep(i,j) = \begin{cases} \dsp\frac{\abs{\tri{p}}}{6} & \text{ si } i = j,\\
\dsp \frac{\abs{\tri{p}}}{12} & \text{ sinon.} \end{cases} $$ Mise sous forme matricielle : $$ \Mep = \frac{\abs{\tri{p}}}{12} \left( \begin{array}{c c c} 2 & 1 & 1\\
1 & 2 & 1 \\
1 & 1 & 2 \end{array} \right). $$

Matrice de rigidité élémentaire

Nous appliquons la même procédure pour la matrice de rigidité, autrement dit, nous calculons les matrices de rigidité élémentaire $\Dep$ définie par $$ \Dep(i,j) = \int_{\tri{p}}\nabla \mphiK{p}{j}\cdot \overline{\nabla\mphiK{p}{i}}. $$

Nous pouvons d’ors et déjà remarqué que, en éléments finis $\Pun$, les quantités $\mphiK{p}{j}$ sont constantes sur chaque triangle et peuvent alors être sorties de l’intégrales : $$ \Dep(i,j) = \nabla \mphiK{p}{j}\cdot \overline{\nabla\mphiK{p}{i}}\int_{\tri{p}}\diff\xx = \abs{\tri{p}}\left[\nabla \mphiK{p}{j}\cdot \overline{\nabla\mphiK{p}{i}}\right] $$ Les fonctions $\mphiK{p}{j}$ peuvent être calculées analytiquement ainsi que leur gradient. Nous étudions tout de même le cas plus abstrait d’éléments finis non forcément $\Pun$.

Triangle de référence

Bien que nous puissions obtenir une expression analytique dans le cas d’un triangle quelconque, nous nous en tenons ici au triangle de référence. Notons que nous disposons des expressions analytiques des gradients des fonctions de forme $\varphih_j$ : $$ \nabla_{\xi,\eta}\varphih_{1} = \left( \begin{array}{l} -1\\
-1 \end{array} \right) , \quad \nabla_{\xi,\eta}\varphih_{2} = \left( \begin{array}{l} 1\\
0 \end{array} \right), \quad \nabla_{\xi,\eta}\varphih_{3} = \left( \begin{array}{l} 0\\
1 \end{array} \right). $$ La matrice de rigidité élémentaire $\hat{D}^e$ du triangle de référence $\Kh$ est alors donnée par : $$ \hat{D}^e(i,j) = \int_{\Kh}\nabla \varphih_{j}\cdot \overline{\nabla\varphih_{i}}. $$

Lemma.

Dans le triangle de référence, la matrice de rigidité élémentaire $\hat{D}^e$ vaut $$ \hat{D}^e = \frac{1}{2} \left( \begin{array}{l l c} 2 & -1 & -1 \\
-1 & 1 & 0 \\
-1 & 0 & 1 \end{array} \right) $$

Proof.

Clairement, la matrice est symétrique. Nous avons : $$\hat{D}_{1,1} = \int_{\Kh}\nabla\varphih_1\cdot\overline{\nabla\varphih_1} \diff (\xi,\eta) = \int_{\Kh} (-1,-1)\left(\begin{array}{l}-1\\\ -1\end{array}\right) \diff (\xi,\eta) = 2 \int_{\Kh} \diff(\xi,\eta) = 1. $$ $$\hat{D}_{2,2} = \int_{\Kh}\nabla\varphih_2\cdot\overline{\nabla\varphih_2} \diff (\xi,\eta) = \int_{\Kh} (1,0)\left(\begin{array}{l}1\\\ 0\end{array}\right) \diff (\xi,\eta) = \int_{\Kh} \diff(\xi,\eta) = \frac{1}{2} =\hat{D}_{3,3}. $$ $$\hat{D}_{1,2} = \int_{\Kh}\nabla\varphih_1\cdot\overline{\nabla\varphih_2} \diff (\xi,\eta) = \int_{\Kh} (-1,-1)\left(\begin{array}{l}1\\\ 0\end{array}\right) \diff (\xi,\eta) = -\int_{\Kh} \diff(\xi,\eta) = -\frac{1}{2}. $$ $$\hat{D}_{1,3} = \int_{\Kh}\nabla\varphih_1\cdot\overline{\nabla\varphih_3} \diff (\xi,\eta) = \int_{\Kh} (-1,-1)\left(\begin{array}{l}0\\\ 1\end{array}\right) \diff (\xi,\eta) = -\int_{\Kh} \diff(\xi,\eta) = -\frac{1}{2}. $$ $$\hat{D}_{2,3} = \int_{\Kh}\nabla\varphih_2\cdot\overline{\nabla\varphih_3} \diff (\xi,\eta) = \int_{\Kh} (1,0)\left(\begin{array}{l}0\\\ 1\end{array}\right) \diff (\xi,\eta) = 0. $$

Triangle quelconque

Pour calculer les dérivées partielles selon $x$ et $y$ de $\varphih_{j}$, nous utilisons la dérivée de fonction composée : $$ \left( \begin{array}{c} \dsp \frac{\partial \mphiK{p}{j}}{\partial x}\\
\dsp \frac{\partial \mphiK{p}{j}}{\partial y} \end{array} \right) = \left( \begin{array}{c c} \dsp \frac{\partial \xi}{\partial x} & \dsp \frac{\partial \eta}{\partial x}\\
\dsp \frac{\partial \xi}{\partial y} & \dsp \frac{\partial \eta}{\partial y} \end{array} \right) \left( \begin{array}{c} \dsp \frac{\partial \varphih_{j}}{\partial \xi}\\
\dsp \frac{\partial \varphih_{j}}{\partial \eta} \end{array} \right) $$ En notant $\BK{p}$ la matrice de passage, nous avons $$ \nabla_{x,y}\mphiK{p}{j}(x,y) = \BK{p}\nabla_{\xi,\eta}\varphih_{j}(\xi,\eta). $$ Pour obtenir la matrice $\BK{p}$, nous utilisons la Jacobienne $\JK{p}$, en remarquant que : $$ \left( \begin{array}{c} \dsp \frac{\partial \varphih_{j}}{\partial \xi}\\
\dsp \frac{\partial \varphih_{j}}{\partial \eta} \end{array} \right) = \left( \begin{array}{c c} \dsp \frac{\partial x}{\partial \xi} & \dsp \frac{\partial y}{\partial \xi}\\
\dsp \frac{\partial y}{\partial \eta} & \dsp \frac{\partial y}{\partial \eta} \end{array} \right) \left( \begin{array}{c} \dsp \frac{\partial \mphiK{p}{j}}{\partial x}\\
\dsp \frac{\partial \mphiK{p}{j}}{\partial y} \end{array} \right). $$ Autrement dit, nous avons $$ \nabla_{\xi,\eta}\varphih_{j}(\xi,\eta) = (\JK{p})^T\nabla_{x,y}\mphiK{p}{j}(x,y). $$ Nous en déduisons que $\BK{p} = (\JK{p}^T)^{-1}$, en particulier, dans le cas d’une transformation linéaire de triangle, nous obtenons : $$ \BK{p} = \frac{1}{\det(\JK{p})} \left( \begin{array}{c c} \yK{p}{3}-\yK{p}{1} & \yK{p}{1}-\yK{p}{2}\\
\xK{p}{1}-\xK{p}{3} & \xK{p}{2}-\xK{p}{1} \end{array} \right). $$ Au final, comme $X\cdot Y = X^TY$, nous obtenons \begin{equation}\label{eq:intRigidite} \int_{\tri{p}} (\nabla\mphiK{p}{j})^T\overline{\nabla\mphiK{p}{i}} \diff(x,y) = \det(\JK{p})\int_{\Kh} (\nabla\varphih_{j})^T (\BK{p}^T \overline{\BK{p}})\overline{\nabla\varphih_{i}} \diff (\xi,\eta). \end{equation} La matrice $\BK{p}$ étant réelle, nous pouvons supprimer la conjugaison portant sur $\BK{p}$.

Pour les éléments finis $\Pun$, la matrice $\BK{p}$ ne dépend pas de $(\xi,\eta)$, de même des gradients $\nabla\varphih_{j}$ : tout peut être sorti de l’intégrale, et comme $\det(\JK{p}) = 2\abs{\tri{p}}$ et $\abs{\Kh}= \frac{1}{2}$ : $$ \int_{\tri{p}} \nabla\mphiK{p}{j}\cdot\overline{\nabla\mphiK{p}{i}} \diff(x,y) = \abs{\tri{p}}(\nabla\varphih_{j})^T (\BK{p}^T \overline{\BK{p}})\overline{\nabla\varphih_{i}}. $$

Quadratures

Sur un triangle

En général, nous ne pouvons pas calculer analytiquement les intégrales. Par exemple, pour les termes du membre de droite, nous devons calculer : $$ \int_{\tri{p}}f(\xx)\overline{\mphiK{p}{i}(\xx)}\diff \xx. $$ Sauf pour certaines fonctions $f$ particulières, nous ne disposerons certainement pas de formules explicites pour ce terme. En pratique, nous passons à l’éléments de référence et nous approchons l’intégrale à l’aide d’une formule de quadrature : $$ \begin{array}{r l} \dsp \int_{\tri{p}}f(\xx)\overline{\mphiK{p}{i}(\xx)}\diff \xx &= \dsp \det(\JK{p})\int_{\Kh}f(\xx(\xi,\eta))\overline{\varphih_i(\xi,\eta)}\diff (\xi,\eta) \\
& \dsp \simeq \det(\JK{p})\sum_{m=1}^M\omega_m f(\xx(\xi_m,\eta_m))\overline{\varphih(\xi_m,\eta_m)}. \end{array} $$ Les points $(\xi_m,\eta_m)$ sont appelés points de quadrature (parfois Points de Gauss, même si la règle de quadrature n’est pas de Gauss) et les quantités $\omega_m\in\Rb$ les poids associés. Notons que le point $\xx_m = \xx(\xi_m,\eta_m)$ s’obtient par l’expression vue précédemment : $$ \xx_m = \sum_{i=1}^3\sumitK{p}{i}\psih_i(\xi_m,\eta_m). $$ Nous présentons ici deux règles de quadrature pour l’intégration sur $\Kh$ d’une fonction $g$ quelconque (intégrale) : $\int_{\Kh}\gh(\xx)\diff\xx$. La première règle est exacte pour des polynômes de degré 1, la deuxième pour des polynômes de degré 2 (règles de Hammer) :

$\xi_m$ $\eta_m$ $\omega_m$ Degré de précision
13 13 16 1
16 16 16 2
46 16 16
16 46 16

Sur une arrête

Voici quelques formules de quadrature sur un segment $[\sumitK{p}{1}, \sumitK{p}{2}]$ avec le degré de précision, i.e la formule est exacte si $g$ est un polynôme de degré égal ou inférieur. Nous notons $\abs{\sigma} = \norm{\sumitK{p}{1} - \sumitK{p}{2}}$ la taille du segment et $\sumitK{p}{12} = \frac{\sumitK{p}{1} + \sumitK{p}{2}}{2}$ le milieu du segment:

Nom Degré de
précision
$\dsp \int_{\sumitK{p}{1}}^{\sumitK{p}{2}}g(x)\diff x\simeq \ldots$
Point du milieu 1 $\dsp g(\sumitK{p}{12})$
Trapèze 1 $\dsp\frac{\abs{\sigma}}{2}\left(g(\sumitK{p}{1}) + g(\sumitK{p}{2})\right)$
13 Simpson 2 $\dsp\frac{\abs{\sigma}}{6}\left(g(\sumitK{p}{1}) + 4g(\sumitK{p}{12}) + g(\sumitK{p}{2})\right)$

Les formules de quadrature ont évidemment un impact sur la qualité de l’approximation, toutefois, elles jouent un rôle relativement mineur par rapport aux autres approximations (et l’on peut choisir plus de points d’intégration !).