5.4 Éléments finis quadrangulaire $\mathcal{Q}^{1}$

La seconde façon de mailler un domaine $\Omega$ est d'utiliser des quadrangles. Bien qu'ayant des possibilités de maillage moins générale que les triangles, les maillages quadrangulaires sont cependant très utilisés, en particulier en mécanique des solides, car ils peuvent être plus précis. L'autre intérêt des approximations quadrangulaires est que ce sont des produits tensoriels d'approximation 1D.

Le maillage quadrangulaire est décrit par une structure de données G identique à celle du maillage triangulaire du paragraphe 6.3.1, mais avec G.ddl=4 et une table de connection G.Tbc à 4 éléments ( de dimension 4*ne). Un exemple de maillage du cylindre est donné sur la figure 6.24.

Figure: Maillage $\mathcal{Q}^{1}$ du cylindre

L'élément $\mathcal{Q}^{1}$ est l'élément quadrangulaire le plus simple, qui possède comme degrés de liberté les valeurs aux 4 sommets du quadrangle.

5.4.1 Élément de référence et fonctions de forme

L'approximation sur un élément $\mathcal{Q}^{1}$ est définie à l'aide d'une transformation vers un élément de référence, qui est le carré $[-1,+1]*[-1,+1]$.

Figure: transformation $\mathcal{T}_{k}: (x,y)\Leftrightarrow(\xi,\eta)$ vers l'élément de référence $\mathcal{Q}^{1}$

Cette transformation $\mathcal{T}^{k}$ est une transformation bilinéaire, qui transforme les 4 sommets $\{ S_{q}\}_{q=1,4}$ de $e_{k}$ vers les 4 sommets $\{\hat{S}_{q}\}_{q=1,4}$ de l'élément de référence $\hat{e}$.

Sur cet élément de référence, l'approximation $f^{h}$ d'une fonction $f$ est une fonction bilinéaire (i.e. le produit d'un polynôme de degré 1 en $\xi$ par un polynôme de degré 1 en $\eta$):

$$f_{\vert e_{k}}^{h}(\xi,\eta)=a\xi\eta+b\xi+c\eta+d$$

Les coefficients $\{ a,b,c,d\}$ dépendent linéairement des valeurs $\{ F_{n_{q}}\}_{q=1,4}$ de $f^{h}$ aux 4 sommets de l'élément (on a noté $\{ n_{q}\}_{q=1,4}$ les numéros des 4 sommets). Cette approximation s'écrit donc sous la forme:

$$f_{\vert e_{k}}^{h}(\xi,\eta)=F_{n_{1}}N_{1}(\xi,\eta)+F_{n_... ...{2}(\xi,\eta)+F_{n_{3}}N_{3}(\xi,\eta)+F_{n_{4}}N_{4}(\xi,\eta)$$ (5.39)

Les fonctions $\{ N_{q}(\xi,\eta)\}_{q=1,4}$ sont les fonctions de forme de l'élément $\mathcal{Q}^{1}$: ce sont des fonctions bilinéaires qui vérifient:

$$N_{q}(\hat{S}_{p})=\delta_{p,q} \mbox{ pour } p,q=1,4$$

Ce sont les polynômes de Lagrange de degré 1, dont l'expression est la suivante:

$$\begin{eqnarray*} N_{1}(\xi,\eta)=\frac{1}{4}(1-\xi)(1-\eta)\\ N_{2}(\xi,\eta)=... ...}{4}(1+\xi)(1+\eta)\\ N_{4}(\xi,\eta)=\frac{1}{4}(1-\xi)(1+\eta)\end{eqnarray*}$$

Ces fonctions de forme $\mathcal{Q}^{1}$ sont tracées sur la figure 6.26. Ce sont des surfaces réglées et non plus des plans comme pour les fonctions de forme $\mathcal{P}^{1}$.

Figure: fonctions de forme $\mathcal{Q}^{1}$

Ces fonctions de formes vont nous permettre aussi de déterminer la transformation $\mathcal{T}^{k}$ puisque les coordonnées d'un point $(x,y)$ de l'élément $e_{k}$ dépend linéairement des coordonnées $\{ x_{n_{q}},y_{n_{q}}\}$ des sommets $\{ S_{q}\}$ de l'élément et est une fonction bilinéaire de $(\xi,\eta)$

$$x=x_{n_{1}}N_{1}(\xi,\eta)+x_{n_{2}}N_{2}(\xi,\eta)+x_{n_{3}}N_{3}(\xi,\eta)+x_{n_{4}}N_{4}(\xi,\eta)$$ (5.40)

$$y=y_{n_{1}}N_{1}(\xi,\eta)+y_{n_{2}}N_{2}(\xi,\eta)+y_{n_{3}}N_{3}(\xi,\eta)+y_{n_{4}}N_{4}(\xi,\eta)$$

La matrice jacobienne $\mathbf{J}_{k}$= $\frac{D(x,y)}{D(\xi,\eta)}$ de cette transformation se calcule alors simplement:

$$\mathbf{J}_{k}=\frac{1}{4}\left[\begin{array}{cc} (x_{n_{2}}... ..._{n_{1}})(1-\xi)+(y_{n_{3}}-y_{n_{2}})(1+\xi)\end{array}\right]$$ (5.41)

On constate que contrairement à l'élément $\mathcal{P}^{1}$, la matrice jacobienne de cette transformation n'est pas constante. Dans le cas d'un quadrangle rectangle de cotés $h_{1}$ et $h_{2}$ alignés avec les axes, cette matrice jacobienne $\mathbf{J}_{k}$ se simplifie:

$$\mathbf{J}_{k}=\frac{1}{2}\left[\begin{array}{cc} h_{1} & 0\\ 0 & h_{2}\end{array}\right]$$

Exercice: montrer que si les cotés du quadrangles sont parallèles, alors la matrice jacobienne est constante. Calculer $\mathbf{J}_{k}$ dans ce cas.

Le programme Matlab JacobienQ1 ci dessous calcule la matrice Jacobienne en un point $(\xi,\eta)$ d'un élément $\mathcal{Q}^{1}$.

programme matlab 6.4.1: calcul de la matrice jacobienne : JacobienQ1.m

function [J]=JacobienQ1(G,k,xi,eta)
% calcul le jacobien de l'elt k au pt xi,eta
n=G.Tbc(k,:); % numero des nds
X=G.X(n,1);
Y=G.X(n,2);
J= 1/4*[(X(2)-X(1))*(1-eta)+(X(3)-X(4))*(1+eta),...
        (X(4)-X(1))*(1-xi)+ (X(3)-X(2))*(1+xi);...
        (Y(2)-Y(1))*(1-eta)+(Y(3)-Y(4))*(1+eta),...
        (Y(4)-Y(1))*(1-xi)+ (Y(3)-Y(2))*(1+xi)];
return;

5.4.2 Intégration numérique

La transformation vers l'élément de référence est pour un quadrangle quelconque plus complexe que pour les éléments $\mathcal{P}^{1}$. Le calcul des intégrales sur l'élément de référence, qui fait intervenir le déterminant du Jacobien, voir la jacobienne inverse, ne se fait plus en générale analytiquement, mais utilise une intégration numérique.

L'intégration numérique sur l'élément de référence $\mathcal{Q}^{1}$ est cependant plus simple que l'intégration numérique sur un élément $\mathcal{P}^{1}$, puisque les bornes d'intégration en $(\xi,\eta)$ sont fixes. Cette intégration est en faite un produit tensoriel d'intégrations numériques en 1D.

5.4.2.1 intégration numérique 1D

Sur le segment de référence $[-1,+1]$, l'intégration numérique d'une fonction $f(\xi)$ utilisant $ng$ points de Gauss s'écrit:

$$\int_{-1}^{+1}f(\xi) d\xi\simeq \sum_{i=1}^{ng}w_{i} f(\xi_{i})$$

Les paramètres de cette formule de quadrature sont tels que l'intégration soit exacte pour les polynômes de degré le plus élevé possible. Cette formule ayant $2ng$ paramètres $\{ w_{i},\xi_{i}\}$, on va pouvoir imposer $2ng$ conditions: i.e intégrer exactement tous les monômes $\xi^{m}$ pour $m=0,2ng-1$. La formule d'intégration avec $ng$ points de Gauss va permettre d'intégrer exactement tous les polynômes de degré $\le2ng-1$.

Par exemple avec 2 points de Gauss on peut intégré exactement des polynômes de degré $\le3$. Cette formule s'écrit:

$$\int_{-1}^{+1}f(\xi) d\xi\simeq w_{1} f(\xi_{1})+w_{2} f(\xi_{2})$$ (5.42)

En écrivant que cette formule (6.42) est exacte pour $f=1,\xi,\xi^{2},\xi^{3}$ on obtiens les 4 relations permettant de déterminer les 4 coefficients $\{ w_{1},\xi_{1},w_{2},\xi_{2}\}$. On peut aussi noter que cette formule doit être symétrique, et donc que l'on a forcément: $\xi_{1}=-\xi_{2}$ et $w_{1}=w_{2}$. La condition pour $f=1$ donne immédiatement la valeur $w_{1}=1$, et la condition pour $f=\xi^{2}$ la valeur $\xi_{1}=\frac{\sqrt{3}}{3}$. D'où la formule de quadrature à 2 points:

$$\int_{-1}^{+1}f(\xi) d\xi\simeq f(-\frac{\sqrt{3}}{3})+f(\frac{\sqrt{3}}{3})$$

Les positions $\xi_{i}$ des points de Gauss sont les racines de polynômes caractéristiques: les polynômes de Legendre.On a cependant pas d'expressions analytiques générales pour ces racines, et il faut les déterminer pour chaque valeur de $ng$.

A l'aide de Maple, on a déterminer analytiquement ces valeurs jusqu'à $ng=5$ en utilisant les propriétés de symétrie de ces points de Gauss. Le programme Maple 6.4.2 ci dessous détermine la formule de quadrature de Gauss pour $ng=3.$

programme Maple 6.4.2: calcul des points de Gauss-Legendre

> restart:
# Formule de Gauss en 1D 3pts
> f:='f':
> Ie:=int(f(xi),xi=-1..1);Ih:=w1*f(xi1)+w2*f(xi2)+w3*f(xi3);
# Conditions de symetrie
> w2:=w1; xi2:=-xi1; xi3:=0; Ih;
# 6 conditions
> f:=xi->1; Ie=Ih; eq1:=%:
> w3:=solve(eq1,w3);
> f:=xi->xi; Ie=Ih; eq1:=%:
> f:=xi->xi^2; Ie=Ih; eq1:=%:
> w1:=solve(eq1,w1);
> f:=xi->xi^3; Ie=Ih; 
> f:=xi->xi^4; Ie=Ih;eq1:=%:
> solve(eq1,xi1);
> xi1:=sqrt(3/5);
> f:=xi->xi^5; Ie=Ih;
> f:='f':Ih;

Exercice: en s'inspirant du programme précédent, déterminer la formule de quadrature de Gauss pour $ng=4.$

Le tableau 6.4 donne les expressions analytiques, ainsi que les valeurs numériques approchées (que l'on comparera avec les valeurs numériques données classiquement dans les livres), pour les formules de quadrature de Gauss d'ordre 1 à 9 (i.e. pour $ng=1,5$).

Tableau 6.4: Intégration numérique 1D

ordre

formule

poids

position

paramêtres

valeurs 1

5.4.2.2 intégration numérique en 2D

Ayant les formules de quadrature de Gauss en 1D, on en déduit les formules de quadrature en 2D en intégrant successivement suivant $\xi$ et $\eta$:

$$\int_{-1}^{+1}\int_{-1}^{+1}f(\xi,\eta) d\xi d\eta \approx \sum_{i=1}^{ng}w_{i}\int_{-1}^{+1}f(\xi_{i},\eta) d\eta$$

$$\approx \sum_{i=1}^{ng}w_{i}\sum_{j=1}^{ng}w_{j}f(\xi_{i},\eta_{j})$$

$$\approx \sum_{i=1}^{ng}\sum_{j=1}^{ng}w_{i}w_{j}f(\xi_{i},\eta_{j})$$ (5.43)

Cette formule de quadrature (6.43) est bien le produit tensoriel de formule 1D.

Remarque: on peut évidement choisir un nombre de points de Gauss différent suivant chaque direction.

5.4.3 Matrices élémentaires

Par rapport à la formulation élément finis $\mathcal{P}^{1}$, le seul changement provient du calcul des matrices élémentaires, puisque l'assemblage et le calcul des intégrales de bords restent identiques.

5.4.3.1 matrice élémentaire de rigidité

L'expression de la matrice élémentaire de rigidité $\mathbf{K}^{k}$ s'écrit:

$$\mathbf{K}_{pq}^{k}=\int_{-1}^{+1}\int_{-1}^{+1}\left(\left[... ...y}\right]\right)\det(\mathbf{J}_{k}) d\eta d\xi p,q=1,4$$

Pour calculer ces intégrales, nous utiliserons une intégration numérique avec 2 points de Gauss dans chaque direction. La fonction Matlab MatriceRigiditeQ1 6.4.3 implémente le calcul de la matrice de rigidité pour un élément $e_{k}$ quelconque.

programme matlab 6.4.3: Matrice de rigidité $\mathcal{Q}^1$: MatriceRigiditeQ1.m

function [Ke]=MatriceRigiditeQ1(G,k)
% calcul de la matrice de masse de l'element k (Q1)
% gradient des fonction de forme
dN1=inline('[-(1-eta)/4,-(1-xi)/4]','xi','eta');
dN2=inline('[ (1-eta)/4,-(1+xi)/4]','xi','eta');
dN3=inline('[ (1+eta)/4, (1+xi)/4]','xi','eta');
dN4=inline('[-(1+eta)/4, (1-xi)/4]','xi','eta');
% pts de Gauss 2pts
Xg=[-sqrt(3)/3,sqrt(3)/3];
Wg=[1,1];
ng=size(Wg,2);
% integration numerique
Ke=zeros(4,4);
for i=1:ng
    for j=1:ng
        J=JacobienQ1(G,k,Xg(i),Xg(j));
        J1=inv(J);
        detJ=det(J);
        dN=[dN1(Xg(i),Xg(j)); dN2(Xg(i),Xg(j)); ...
            dN3(Xg(i),Xg(j)); dN4(Xg(i),Xg(j))];
        % gradient dans l'espace physique
        dP=dN(:,:)*J1;
        Ke=Ke+Wg(i)*Wg(j)*det(J)*(dP(:,1)*dP(:,1)'+dP(:,2)*dP(:,2)');
    end;
end;
return;

Remarque: on peut vérifier numériquement que la précision d'intégration est suffisante, puisque que pour des quadrangles pas trop déformés, on obtiens la solution exacte.

Dans le cas d'un triangle rectangle de cotés $h_{1}$ et $h_{2}$ alignés suivant les axes $x$ et $y$, la matrice jacobienne est constante et on peut calculer analytiquement la matrice de rigidité:

$$\mathbf{K}^{k}=\frac{1}{6h_{1}h_{2}}\left[\begin{array}{cccc... ...h_{1}^{2}-2h_{2}^{2} & 2(h_{1}^{2}+h_{2}^{2})\end{array}\right]$$

On vérifie que la matrice est symétrique et que la somme des lignes et des colonnes est nulle.

Exercice: démontrer la relation précédente pour $\mathbf{K}^{k}$

5.4.3.2 Matrice élémentaire de masse

L'expression de la matrice élémentaire de masse $\mathbf{M}^{k}$ s'écrit:

$$\mathbf{M}_{pq}^{k}=\int_{-1}^{+1}\int_{-1}^{+1}N_{q}(\xi,\eta)N_{p}(\xi,\eta)\det(\mathbf{J}_{k}) d\eta d\xi p,q=1,4$$

Pour calculer ces intégrales nous utiliserons une intégration numérique de Gauss avec 2 points de Gauss dans chaque direction, bien que l'on puisse les calculer analytiquement. En effet le déterminant du Jacobien $\det(\mathbf{J}_{k})$ est un polynôme de degré 1 en $\xi$ et $\eta,$ et donc $\mathbf{M}_{pq}^{k}$ est une intégrale de polynômes de degré 3 en $\xi$ et $\eta$. Mais l'expression obtenue reste cependant complexe et n'est pas généralisable au cas où le coefficient $\alpha$ devant la matrice de masse est variable.

La fonction Matlab MatriceMasseQ1 6.4.3 implémente le calcul de la matrice de masse pour un élément $e_{k}$ quelconque.

programme matlab 6.4.3: Matrice de masse $\mathcal{Q}^1$: MatriceMasseQ1.m

function [Me]=MatriceMasseQ1(G,k)
% calcul de la matrice de masse de l'element k (Q1)
% fonction de forme
N1=inline('(1-xi)/2*(1-eta)/2','xi','eta');
N2=inline('(1+xi)/2*(1-eta)/2','xi','eta');
N3=inline('(1+xi)/2*(1+eta)/2','xi','eta');
N4=inline('(1-xi)/2*(1+eta)/2','xi','eta');
% pts de Gauss 2pts
Xg=[-sqrt(3)/3,sqrt(3)/3];
Wg=[1,1];
%
ng=size(Wg,2);
% integration numerique
Me=zeros(4,4);
for i=1:ng
    for j=1:ng
      J=JacobienQ1(G,k,Xg(i),Xg(j));
      detJ=det(J);
      N=[N1(Xg(i),Xg(j)),N2(Xg(i),Xg(j)),N3(Xg(i),Xg(j)),N4(Xg(i),Xg(j))];
      Me=Me+Wg(i)*Wg(j)*det(J)*N'*N;
    end;
end;
return;

Dans le cas d'un triangle rectangle de cotés $h_{1}$ et $h_{2}$ alignés suivant les axes $x$ et $y$, le Jacobien de la transformation est constant et on peut calculer simplement la matrice de masse:

$$\mathbf{M}_{pq}^{k}=\frac{h_{1}h_{2}}{36}\left[\begin{array}... ...& 4 & 2 & 1\\ 1 & 2 & 4 & 2\\ 2 & 1 & 2 & 4\end{array}\right]$$

Exercice: montrer que l'on obtiens la même expression de la matrice de masse pour un rectangle non aligné avec les axes.

5.4.3.3 second membre élémentaire

Comme dans le cas $\mathcal{P}^{1}$, le second membre élémentaire s'écrit comme le produit de la matrice de masse par le vecteur $F^{k}$ élémentaire des valeurs nodales de $f$:

$$\mathbf{B}^{k}=\mathbf{M}^{k}\mathbf{F}^{k} \mbox{ avec } \mathbf{F}^{k}=[f_{n_{1}},f_{n_{2}},f_{n_{3}},f_{n_{4}}]$$

La fonction Matlab SmbElementQ1 6.4.3 est alors une simple adaptation de la fonction SmbElement

programme matlab 6.4.3: Second membre $\mathcal{Q}^1$: SmbElementQ1.m

function [Be]=SmbElementQ1(G,F,k)
% calcul du second membre elmentaire (F valeurs nodales de f)
n=G.Tbc(k,:); % numero des sommets de l'element k
Fk=F(n);
Mk=MatriceMasseQ1(G,k);
Be=Mk*Fk;
return;

5.4.4 Étude de la précision

Nous allons étudier le même cas qu'au paragraphe 6.3, mais avec des éléments $\mathcal{Q}^{1}$. Le maillage est obtenu par transformation “quadrilatère courbe” et utilise les programmes Matlab de maillage décrit en annexe.

Le programme 6.4.4 ci dessous donne la fonction Matlab cylindremeshQ1 utilisée pour le maillage de $\Omega$, paramétré par le nombre de noeuds $n_{1}$et $n_{2}$ sur les cotés AB et BC.

programme matlab 6.4.4: Maillage éléments finis 2D $\mathcal{Q}^1$: cylindremeshQ1.m

function G=cylindremeshQ1(n1,n2)
% maillage du cylindre en Q1
addpath ../Mailleur
% construction de la geometrie
R1=1; R2=3;
A=[-R2,0]; B=[-R1,0]; C=[0,R1]; D=[0,R2];
AB=[linspace(A(1),B(1),n1)',linspace(A(2),B(2),n1)'];
theta=linspace(pi,pi/2,n2)';
BC=[R1*cos(theta),R1*sin(theta)];
CD=[linspace(C(1),D(1),n1)',linspace(C(2),D(2),n1)'];
theta=linspace(pi/2,pi,n2)';
DA=[R2*cos(theta),R2*sin(theta)];
G=Quacou(AB,BC,CD,DA,[1,1,3,2],2);
return;

Le script Matlab 6.4.4 qui enchaıne la suite des opérations pour calculer la solution approchée du problème (6.37) avec des éléments $\mathcal{Q}^{1}$ est donné ci-dessous.

programme matlab 6.4.4: Résolution par éléments finis 2D $\mathcal{Q}^1$

% maillage du cylindre en Q1
G=cylindremeshQ1(20,20);
% second membre
F=zeros(G.nn,1);
% assemblage
[A,B]=AssemblageQ1(G,1.0,0.0,F);
% conditions aux limites
Psie=inline('y-y./(x.^2+y.^2)','x','y');
[A,B]=Climite(G,A,B,0.0,0.0,Psie);
% resolution
[LA,UA]=lu(A);
Psi=UA\(LA\B);
% calcul de l'erreur
err=ErreurQ1(G,Psie,Psi)

Le calcul a été effectué sur un maillage ayant le même nombre de noeuds ( $n_{1}=n_{2}=20$) que celui utilisé au paragraphe 6.3.2pour les éléments $\mathcal{P}^{1}$. La solution obtenue est tracée sur la figure (), que l'on peut comparer avec la solution $\mathcal{P}^{1}$ sur un maillage équivalent (figure 6.21). On constate que les deux solutions sont quasiment identiques (le maillage étant relativement fin pour ce type de calcul).

Figure 6.27: solution $\psi ^{h}$ par éléments finis $\mathcal{Q}^{1}$ avec un maillage équivalent à 6.19

Pour étudier la précision, nous avons fait varier la taille caractéristique $h$ des éléments du maillage et nous avons calculer l'erreur d'approximation $\left\Vert \psi_{ex}-\psi^{h}\right\Vert$. Comme en $\mathcal{P}^{1}$, nous avons calculer l'erreur élémentaire sur un élément à l'aide d'une intégration numérique utilisant 3 points de Gauss dans chaque direction.

Remarque: pour le calcul de l'erreur nous avons utilisé plus de points de Gauss $(3*3)$ que pour la calcul de la solution approchée ($2*2)$, de façon à minimiser l'erreur numérique d'intégration par rapport à l'erreur d'approximation.

La fonction Matlab ErreurQ1 (programme 6.4.4) correspondante qui calcul l'erreur globale (6.38),est donnée ci dessous.

programme matlab 6.4.4: calcul de l'erreur d'intégration $\mathcal{Q}^1$: ErreurQ1.m

function [err]=ErreurQ1(G,f,F)
% calcul de l'erreur en norme L2
% pts de Gauss
% pts de Gauss 3pts
Xg=[-sqrt(15)/5,0,sqrt(15)/5];
Wg=[5/9,8/9,5/9];
ng=size(Wg,2);
% fct de forme
N1=inline('(1-xi)/2*(1-eta)/2','xi','eta');
N2=inline('(1+xi)/2*(1-eta)/2','xi','eta');
N3=inline('(1+xi)/2*(1+eta)/2','xi','eta');
N4=inline('(1-xi)/2*(1+eta)/2','xi','eta');
err=0.0;
for k=1:G.ne
    % calcul de l'erreur / elt
    n=G.Tbc(k,:); % numero des sommets de l'element k
    for i=1:ng
      for j=1:ng
        J=JacobienQ1(G,k,Xg(i),Xg(j));
        detJ=det(J);
        % calcul des coordonnees des pts de Gauss
        N=[N1(Xg(i),Xg(j)),N2(Xg(i),Xg(j)),...
           N3(Xg(i),Xg(j)),N4(Xg(i),Xg(j))];
        X=N*G.X(n,:);
        % erreur sur l'elt k
        errk=(f(X(:,1),X(:,2))-N*F(n));
        err=err+Wg(i)*Wg(j)*detJ*errk^2;
      end;
    end;
end;
err=sqrt(err);
return;

Figure: Norme de l'erreur d'approximation $\mathcal{Q}^{1}$ en fonction de $n\propto h^{-1}$

Nous avons donc construit une série de maillage de $\Omega,$ avec une taille $h\rightarrow0$ et un rapport d'aspect $R$ constant en fixant un nombre de noeuds identiques sur les cotés lors de la génération du maillage: i.e. $n_{1}=n_{2}=n$. Dans ce cas la taille des éléments $h$ est proportionnelle à $\frac{1}{n}$. On a donc tracé sur la figure 6.28, la norme de l'erreur en fonction de $n$ et on vérifie que l'on obtiens une décroissance de l'erreur en $n^{-2}$, soit en $\theta(h^{2})$, ce qui montre que l'approximation par éléments finis $\mathcal{Q}^{1}$ est d'ordre 2. En comparant la courbe obtenue avec celle de la figure 6.23, on constate que dans ce cas l'erreur est un peu plus faible avec les éléments $\mathcal{Q}^{1}$ qu'avec les éléments $\mathcal{P}^{1}$. C'est une constatation assez générale dans le cas de solutions régulières alignées avec le maillage.