Convergence, erreur et maillage

Objectifs du chapitre

À la fin de ce chapitre, vous serez capable de:

• énoncer l'estimation a priori $\Vert u-u_h{\Vert }_E\le C{h}^p|u{|}_{p+1}$, démontrer l'estimation d'interpolation qui la fonde pour les éléments linéaires en dimension 1, et dire ce que chaque facteur — le pas $h$, le degré $p$, la régularité $|u{|}_{p+1}$ — y contrôle;
• mesurer un ordre de convergence sur une suite de maillages, avec ou sans solution exacte, et estimer l'erreur du maillage le plus fin par l'extrapolation de Richardson;
• comparer le raffinement h et le raffinement p à nombre de degrés de liberté égal, et reconnaître la convergence exponentielle du second sur une solution régulière;
• repérer les singularités d'un modèle — coin rentrant, pointe de fissure, force ponctuelle en dimension 2 —, en calculer l'exposant, expliquer pourquoi l'ordre chute pour tous les degrés et comment un maillage gradué le rétablit;
• appliquer le test de la pièce à une routine d'élément, et construire puis justifier un estimateur a posteriori — par le résidu ou par la reconstruction de Zienkiewicz et Zhu (estimateur «ZZ») — qui pilote un maillage adaptatif;
• distinguer vérification et validation, et dérouler une liste de contrôle avant de faire confiance à un résultat d'éléments finis.

Les chapitres précédents ont mesuré beaucoup d'erreurs: la table du problème modèle aux chapitres 1 et 4, celle des éléments quadratiques au chapitre 7, la valeur au centre de la plaque carrée au chapitre 8. Ils ont aussi pris des engagements: démontrer l'ordre 1 en énergie des éléments linéaires (chapitres 3 et 4), expliquer la perte d'ordre aux singularités (chapitres 2, 4 et 7), tirer une estimation d'erreur des sauts de l'effort normal (chapitre 4). Ce chapitre tient ces engagements. Il répond à la question que tout calcul pose à la fin: de combien le résultat est-il faux, et que faire pour qu'il le soit moins?

Estimations a priori

Une estimation est dite a priori quand elle borne l'erreur avant tout calcul, à partir de propriétés de la solution exacte — que l'on ne connaît pas. Elle ne donne donc pas un nombre, mais une loi: comment l'erreur décroît quand on raffine, et de quoi elle dépend. Le lemme de Céa (théorème 3.3) en a fait la moitié du travail: la solution de Galerkin est, en énergie, la meilleure approximation de $u$ dans $V_h$. Il reste à savoir à quelle distance de $u$ se trouve l'espace $V_h$; pour le majorer, il suffit d'exhiber un élément de $V_h$ proche de . Le candidat naturel est l'interpolant.

Mesurer la régularité

Une estimation d'erreur fait intervenir les dérivées d'ordre supérieur de $u$. Il faut les mesurer dans les mêmes termes que l'énergie, c'est-à-dire en moyenne quadratique.

Pour le problème modèle, $u=\sin \pi x$ et $u^{\prime \prime }=-{\pi }^2\sin \pi x$, d'où $|u{|}_2={\pi }^2/\sqrt{2}=6,9789$. Plus généralement, pour , on a : . Cette remarque, banale en dimension 1, cesse d'être vraie en dimension 2 dès que le domaine a un coin rentrant, et c'est toute l'histoire des singularités.

L'interpolant et son erreur

Le théorème suivant est la démonstration exigée de ce chapitre. Il ne fait intervenir que l'intégration par parties et l'inégalité de Cauchy–Schwarz.

Démonstration. Étape 1: une inégalité de Poincaré sur l'élément. Soit $v\in H^1(a,b)$ nulle aux deux extrémités. Pour $x$ dans la moitié gauche, $v(x)={\int }_a^x{v}^{\prime }$, et Cauchy–Schwarz donne, comme dans la démonstration du théorème 2.3, . Pour dans la moitié droite, on part de l'autre extrémité: et . En intégrant sur chaque moitié,

Donc $\Vert v{\Vert }_{L^2(e)}\le \frac{h_e}{2}\Vert v^{\prime }{\Vert }_{L^2(e)}$.

Étape 2: l'erreur sur la pente. Sur l'élément, $I_{h}u$ est affine, donc $w^{\prime \prime }=u^{\prime \prime }$, et $w$ s'annule en $a$ et en $b$ par construction. Une intégration par parties, dont les termes de bord $w{w}^{\prime }$ s'annulent, donne

par Cauchy–Schwarz puis l'étape 1 appliquée à $v=w$. Si $w^{\prime }=0$ sur l'élément, il n'y a rien à démontrer; sinon on divise par $\Vert w^{\prime }{\Vert }_{L^2(e)}$, ce qui donne la première inégalité (11.2). La seconde s'obtient en appliquant encore l'étape 1: .

Étape 3: la somme. Les carrés s'additionnent élément par élément: $|u-I_{h}u{|}_1^2={\sum }_e\Vert w^{\prime }{\Vert }_{L^2(e)}^2\le {\sum }_e\frac{h_e^2}{4}\Vert u^{\prime \prime }{\Vert }_{L^2(e)}^2\le \frac{h^2}{4}|u{|}_2^2$, et de même en norme avec .

La démonstration dit où se loge l'erreur: chaque élément contribue selon sa longueur et selon la courbure de $u$ sur lui. Un élément long dans une région où $u$ est presque affine coûte peu; un élément court dans une région très courbée peut coûter beaucoup. Toute la stratégie de maillage du chapitre est dans ce produit $h_e\Vert u^{\prime \prime }{\Vert }_{L^2(e)}$.

La constante $\frac{1}{2}$ n'est pas la meilleure. L'étape 1 est l'inégalité de Poincaré sur un intervalle avec deux extrémités nulles, dont le chapitre 2 a donné la constante optimale $h_e/\pi$, sans la démontrer; avec elle, la même démonstration donne $\Vert w^{\prime }{\Vert }_{L^2(e)}\le \frac{h_e}{\pi }\Vert u^{\prime \prime }{\Vert }_{L^2(e)}$, et cette constante est atteinte: pour sur l'élément, et les deux membres sont égaux. Pour une fonction et un maillage fin, enfin, est presque constante sur chaque élément, l'erreur de pente y vaut à peu près autour du milieu , et son carré intégré vaut : la constante est .

Du lemme de Céa à l'ordre de convergence

Il suffit maintenant d'enchaîner. Le lemme de Céa (théorème 3.3) dit que $u_h$ fait au moins aussi bien, en énergie, que n'importe quel élément de $V_h$, et en particulier que $I_{h}u$. Pour la barre de rigidité $EA\le E{A}_{\max }$, la norme d'énergie vérifie par (2.16). D'où, avec le théorème 11.1,

Pour le problème modèle, $EA=1$ et $|u{|}_2=\Vert f{\Vert }_{L^2}$: l'erreur en énergie est majorée par $\frac{h}{2}\Vert f{\Vert }_{L^2}$, une quantité que l'on peut calculer avant toute résolution. L'ordre 1 de la table du cours est démontré, comme les chapitres 3 et 4 l'avaient annoncé. Remarquez la structure de l'argument, qui est celle de toute la théorie: (Céa, qui ne dépend que de l'équation) fois (l'interpolation, qui ne dépend que de l'espace).

Le même schéma vaut pour tous les éléments de Lagrange, en toute dimension.

Esquisse de démonstration. Par Céa, il suffit d'estimer l'erreur d'interpolation, et on le fait élément par élément en se ramenant à l'élément de référence. Plaçons-nous en dimension 1, avec $x=x_e+\frac{h}{2}(\xi +1)$ et $\hat{u}(\xi )=u(x(\xi ))$. Sur l'élément de référence , l'interpolant de degré reproduit exactement les polynômes de degré (théorème 7.1); l'opérateur «erreur d'interpolation» s'annule donc sur eux, et un résultat d'analyse fonctionnelle — le , admis ici — en déduit , avec une constante qui ne dépend que de . Reste à revenir à l'élément réel. Chaque dérivation en apporte un facteur et l'élément de longueur un facteur :

On somme les carrés sur les éléments. En dimension 2, le même calcul passe par la matrice jacobienne de la transformation isoparamétrique, et c'est là qu'intervient l'hypothèse sur les angles: un triangle aplati a un jacobien presque nul, et la constante explose. L'estimation (11.6) s'obtient par un argument de dualité dû à Aubin et Nitsche, que l'exercice 11.5 fait en entier pour les éléments linéaires en dimension 1. $\square$

Trois lectures de (11.5), qui organisent la suite du chapitre.

• L'exposant est le degré. Diviser $h$ par 2 divise l'erreur en énergie par $2^p$: par 2 pour les éléments linéaires, par 4 pour les quadratiques, par 8 pour les cubiques. C'est le raffinement h. Augmenter $p$ change l'exposant lui-même: c'est le raffinement p.
• La régularité est une hypothèse. Si $u$ n'est pas dans $H^{p+1}$ — parce que le domaine a un coin rentrant, parce que la charge est ponctuelle —, l'estimation ne s'applique pas, et l'ordre observé sera plus bas que $p$, quel que soit $p$.

La superconvergence

L'estimation (11.5) est une borne globale, en moyenne quadratique. En certains points, l'approximation fait mieux que ce qu'elle promet. Le cours en a déjà rencontré trois exemples, qu'il est temps de rassembler sous leur nom de superconvergence:

• les nœuds en dimension 1: pour $-(EA{u}^{\prime }{)}^{\prime }=q$ avec charge intégrée exactement, les valeurs nodales des éléments linéaires sont exactes (théorème 4.3), et celles des extrémités d'éléments quadratiques aussi (chapitre 7). C'est une propriété de la dimension 1, que la fonction de Green rend possible et qui disparaît en dimension 2;
• les dérivées au milieu des éléments linéaires: au chapitre 4, l'erreur sur $u_h^{\prime }$ au milieu des éléments est d'ordre 2 ($1,2600\cdot 10^{-3}$ pour ), contre l'ordre 1 aux nœuds ();

Ces points privilégiés ne sont pas des curiosités: ils fournissent une approximation du gradient meilleure que celle de $u_h^{\prime }$ elle-même. La comparer à $u_h^{\prime }$ donnera une estimation de l'erreur — c'est l'idée de Zienkiewicz et Zhu, plus loin dans ce chapitre.

Mesurer un ordre de convergence

Une estimation a priori est un énoncé sur une famille de maillages. Pour la confronter à un calcul, on résout le même problème sur plusieurs maillages emboîtés et l'on mesure la pente.

La question 1.8 a programmé (11.7). Appliquée aux tables du cours, elle donne:

Les ordres en énergie approchent $p$ par en dessous, exactement comme (11.5) le prédit. (Pour $p=3$, les erreurs, que le chapitre 7 ne donnait que par leurs rapports, valent $2,6332\cdot 10^{-2}$, $3,3650\cdot 10^{-3}$, $4,2295\cdot 10^{-4}$, et pour à 32, recalculées ici par le même programme.)

Le pas $h$ n'est pas toujours la bonne abscisse. Pour comparer des éléments de degrés différents, ce qui compte est le coût, dont le nombre de degrés de liberté est la mesure la plus simple. La figure 11.1 porte donc l'erreur en fonction de $N$. En dimension 1, $N\approx pn$ et les pentes valent $-p$, mais sur les maillages grossiers elles sont nettement plus faibles: entre $n=2$ et $n=4$, les éléments linéaires passent de 1 à 3 inconnues pour une erreur divisée par 1,94 seulement — une pente de 0,60 en $N$.

En dimension 2: la plaque carrée

Le chapitre 8 a fixé la valeur au centre de la plaque carrée, $-\Delta u=1$ sur $(0,1{)}^2$ avec $u=0$ au bord, calculée avec des triangles linéaires sur une grille $N\times N$. Prolongeons la table de deux maillages, avec le même programme et les mêmes règles (diagonales du coin inférieur gauche au coin supérieur droit, charge $A/3$ par nœud):

L'ordre tend vers 2, celui de (11.6) avec $p=1$ — une valeur en un point n'est pas une norme $L^2$, et l'estimation (11.6) ne s'y applique pas telle quelle, mais elle se comporte ici de la même façon. Notez combien le régime asymptotique est lent à s'installer: l'ordre observé n'est que 1,73 entre $N=2$ et $N=4$. Et le coût: en dimension 2, diviser $h$ par 2 multiplie le nombre d'inconnues par 4 environ, si bien qu'une erreur d'ordre 2 en $h$ ne décroît que comme $N^{-1}$.

Sans solution exacte: l'extrapolation de Richardson

Dans un calcul réel, la colonne «erreur» n'existe pas. Trois maillages suffisent pourtant à mesurer l'ordre et à estimer l'erreur. Si $u_k\approx u^{\ast }+C{h}_k^p$ sur trois maillages de pas $h$, $h/r$, , alors et : leur rapport vaut , . On en tire l'ordre, puis la limite:

C'est l'extrapolation de Richardson, que l'Analyse numérique a appliquée à l'intégration numérique (chapitre 8, théorème 8.9), où elle transforme la méthode des trapèzes en celle de Simpson; ici, on l'applique à des maillages. Le second terme de (11.8) est l'estimation de l'erreur du maillage le plus fin.

Raffinement h et raffinement p

La figure 11.1 met en regard les deux façons de réduire l'erreur que distingue l'estimation (11.5). Le raffinement h garde le degré et divise les éléments: l'erreur décroît comme une puissance de $N$, de pente $-p/d$. Le raffinement p garde le maillage et monte le degré: pour une solution analytique comme $\sin \pi x$, l'erreur décroît plus vite que toute puissance de $N$. Les erreurs de la courbe p de la figure 11.1, avec deux éléments et $p=1$ à 8, sont

pour 1, 3, 5, …, 15 inconnues. Le rapport entre deux degrés successifs vaut 4,90, puis 7,49, 10,06, 12,63, 15,19, 17,74 et 20,30: il croît sans cesse, alors qu'une loi en $N^{-q}$ donnerait des rapports qui tendent vers 1. C'est la signature d'une convergence exponentielle, en ${\mathrm{e}}^{-bN}$ ou mieux. À 15 inconnues, le raffinement p est à $3,8\cdot 10^{-8}$, les éléments quadratiques à $1,27\cdot 10^{-2}$ et les linéaires à : six ordres de grandeur d'écart au même coût nominal.

Ce coût nominal n'est pas tout. Un élément de degré élevé a une matrice élémentaire pleine de taille $(p+1{)}^2$ en dimension 1 — bien plus en dimension 2 et 3 —, il demande des quadratures plus riches (chapitre 7), et la matrice assemblée est plus large de bande. Mais l'écart est tel qu'il n'y a pas de discussion pour une solution régulière. Toute la question est là: le raffinement p exploite la régularité, et il s'effondre là où elle manque. Le remède, quand elle manque en quelques points connus, est la méthode hp: des éléments de degré élevé là où la solution est lisse, de petits éléments là où elle ne l'est pas. Babuška et ses collaborateurs ont montré que, sur des maillages géométriquement raffinés vers une singularité, elle retrouve une convergence exponentielle.

Les singularités

Les chapitres 2, 4 et 7 ont prévenu: tout ce qui précède suppose $u\in H^{p+1}$, et cette hypothèse tombe en des points que le modèle désigne lui-même.

D'où viennent-elles

Plaçons-nous au sommet d'un coin du domaine, d'angle intérieur $\omega$, avec $u=0$ sur ses deux côtés, et cherchons les solutions de $\Delta u=0$ en coordonnées polaires $(r,\theta )$ centrées au sommet, $0\le \theta \le \omega$. Les fonctions

sont harmoniques (exercice 11.3) et s'annulent sur les deux côtés. Près du sommet, la solution d'un problème de Poisson s'écrit comme une partie régulière plus une combinaison de ces termes, et c'est le premier, $k=1$, qui commande.

L'exposant $\alpha =\pi /\omega$ ne dépend que de l'angle. Pour un coin saillant, $\omega <\pi$, il est supérieur à 1 et le gradient reste borné: un carré ne pose aucun problème. Pour un coin rentrant, $\omega >\pi$, il passe sous 1:

• le domaine en L, $\omega =3\pi /2$: $\alpha =2/3$, et $u\approx c{r}^{2/3}\sin (2\theta /3)$;
• la fissure, $\omega =2\pi$ (les deux lèvres sont deux côtés confondus): . En élasticité plane, le déplacement se comporte comme et les contraintes comme près de la pointe; le coefficient de ce terme est le de la mécanique de la rupture.

Trois autres sources sont fréquentes en pratique. Un changement de type de condition aux limites le long d'un bord droit — encastré d'un côté d'un point, libre de l'autre — produit lui aussi un exposant $1/2$ (exercice 11.3). Un coin d'interface entre deux matériaux de rigidités différentes en produit un qui dépend du contraste. Et une force ponctuelle en dimension 2 produit une singularité d'un autre genre, plus forte, que nous allons regarder de près.

Une fonction de $H^1$ qui n'est pas bornée

Le chapitre 2 a démontré qu'en dimension 1 toute fonction de $H^1$ est continue (théorème 2.1), et il a annoncé que cela cesse d'être vrai dans le plan. Le chapitre 8 en a donné l'exemple classique, que nous reprenons ici pour en tirer les conséquences sur la convergence. Sur le disque de rayon $\frac{1}{2}$, posons

Elle tend vers $+\infty$ quand $r\to 0$, lentement. Son gradient est radial, de module $|v^{\prime }(r)|=\frac{1}{r\ln (1/r)}$, et en coordonnées polaires, avec le changement de variable $t=\ln (1/r)$,

L'énergie est finie, et la fonction n'est pas bornée. En dimension 2, une fonction de $H^1$ n'a donc pas de valeur en un point, et la forme linéaire $\ell (v)=Pv(x_0)$ d'une force concentrée n'est pas continue sur $H^1$: le théorème de Lax–Milgram ne s'applique pas, et il n'y a pas de solution d'énergie finie.

On le voit sur la solution. La flèche d'une membrane infinie de tension $T$ sous une force $P$ est $u=-\frac{P}{2\pi T}\ln r$ (à une constante près), dont l'énergie dans la couronne $\epsilon <r<R$ vaut — infinie quand . Le déplacement la force est infini. Un calcul par éléments finis donne pourtant un nombre, et c'est ce nombre qu'il faut savoir lire.

Le domaine en L

Le domaine en L, $(-1,1{)}^2$ privé du carré $[0,1]\times [-1,0]$, est l'exemple type d'une singularité de coin. Pour mesurer l'erreur proprement, prenons comme solution exacte la fonction singulière elle-même,

harmonique, nulle sur les deux côtés du coin rentrant, et imposée comme condition de Dirichlet (interpolée aux nœuds) sur le reste du bord. Son énergie vaut $a(u,u)=1,83623$, soit $\Vert u{\Vert }_E=1,35507$. Les maillages sont ceux de la plaque carrée — $N$ cellules par unité de longueur, coupées en deux triangles —, soit uniformes, soit gradués vers le coin par la transformation radiale $P\mapsto P{\rho }^{\mu -1}$, où et : elle laisse le bord extérieur et les deux côtés du coin en place et resserre les éléments près de l'origine (figure 11.2). L'erreur en énergie est calculée exactement par , où les deux premiers termes se ramènent, étant harmonique, à des intégrales sur le bord extérieur, loin de la singularité. Résultats:

Sur maillage uniforme, le rapport tend vers $2^{2/3}=1,587$, pas vers 2: l'erreur décroît comme $h^{2/3}$, c'est-à-dire comme $N^{-1/3}$, au lieu de l'ordre 1 promis par (11.5). Sur maillage gradué, le rapport tend vers 2, l'ordre 1 est rétabli, et l'erreur de $2,0\cdot 10^{-2}$ s'obtient avec 2 945 inconnues au lieu de 48 641, .

Le résultat général est le suivant; nous l'admettons.

La condition (11.10) se comprend par un bilan. Sur un élément de taille $h_K$ situé à la distance $r$ du coin, les dérivées d'ordre $p+1$ de $r^{\alpha }$ valent environ $r^{\alpha -p-1}$, et la contribution de l'élément à l'erreur au carré est de l'ordre de . Avec la transformation radiale, , et le nombre d'éléments à une distance comprise entre et est de l'ordre de . En multipliant, l'erreur au carré totale est de l'ordre de , intégrale finie si et seulement si : c'est (11.10). Pour le domaine en L et des éléments linéaires, il faut ; nous avons pris . Pour des éléments quadratiques, il faudrait .

Le premier point du théorème est la raison de l'avertissement du chapitre 7: monter en degré ne rapporte rien près d'une singularité si l'on ne raffine pas en même temps. Des éléments quadratiques sur le maillage uniforme du domaine en L convergeraient, eux aussi, comme $h^{2/3}$ — avec une meilleure constante, mais la même pente.

Le test de la pièce

L'estimation (11.5) suppose une méthode conforme et des intégrales élémentaires correctes. Un élément qui viole l'une de ces hypothèses — un élément non conforme, un élément sous-intégré, ou simplement une routine d'élément boguée — peut-il converger quand même? Bruce Irons a proposé dans les années 1960 un critère pratique, le test de la pièce (patch test), que le chapitre 10 a défini pour tous les éléments isoparamétriques (définition 10.5), en annonçant que ce chapitre-ci en ferait une condition nécessaire de convergence. Rappelons-le, sous la forme qui sert en vérification de code.

On prend une petite pièce d'éléments, volontairement irrégulière — des éléments de tailles et de formes différentes, au moins un nœud intérieur —, et l'on impose sur son bord les valeurs d'une solution exacte affine (en élasticité: un champ de déplacement à déformation constante), sans charge intérieure. Puisque l'espace des éléments contient les fonctions affines (théorème 7.1), la solution de Galerkin doit reproduire cette solution exactement: aux nœuds intérieurs, et dans le gradient (la déformation) de chaque élément. Un élément qui échoue ne peut pas converger vers la bonne solution quand on raffine, parce qu'à l'échelle d'un élément assez petit toute solution régulière ressemble à une solution affine: le test de la pièce est une condition nécessaire de convergence.

Une routine écrite pour des maillages uniformes. Voici le défaut le plus courant qu'il attrape. Une routine d'élément de barre calcule sa rigidité avec le pas global $h=1/n$ au lieu de la longueur $h_e$ de l'élément — elle a été écrite et testée sur des maillages uniformes, où les deux coïncident. Appliquons le test de la pièce: trois éléments sur les nœuds $0$; $0,2$; $0,5$; $1$, avec $u(0)=0$ et , sans charge. La solution exacte est , de déformation 1 partout. La routine fautive assemble sur les deux nœuds intérieurs, avec le second membre , et trouve et au lieu de et ; les déformations des trois éléments valent , et au lieu de 1. Le test échoue. Sur une pièce de trois éléments égaux, il aurait réussi: c'est pourquoi la pièce doit être irrégulière.

Une pièce de triangles. Dans le carré unité, plaçons un nœud intérieur en $(0,4;0,55)$, relié aux quatre coins par quatre triangles, et imposons aux coins les valeurs de $u=1+2x+3y$. L'équation du nœud intérieur, assemblée à partir des matrices élémentaires du chapitre 8, a pour coefficient diagonal $1625/396$, et elle donne exactement $u_P=69/20=3,45=1+2\cdot 0,4+3\cdot 0,55$; le gradient vaut dans chacun des quatre triangles. Le triangle linéaire passe le test, comme il se doit pour un élément conforme intégré exactement.

Le test n'est pas suffisant: un élément peut le passer et rester inutilisable parce qu'il est instable — l'intégration réduite et ses modes parasites (chapitres 7 et 10) en sont l'exemple. Il n'est pas non plus réservé aux éléments exotiques: c'est le premier test à faire passer à tout code d'éléments finis que l'on écrit ou que l'on modifie, sur une pièce distordue, et il détecte la plupart des erreurs de matrice élémentaire, de jacobien et d'assemblage. Strang et Fix ont analysé son lien avec la convergence des méthodes qui commettent des «crimes variationnels» — non-conformité, quadrature inexacte —, et c'est sous cette forme qu'il entre dans la théorie.

L'estimation a posteriori

Les estimations a priori disent comment l'erreur varie; elles ne disent pas combien elle vaut pour le maillage que l'on vient de calculer. L'extrapolation de Richardson le dit, mais elle demande trois calculs emboîtés et un régime asymptotique. Une estimation a posteriori utilise la solution calculée $u_h$ elle-même — et les données — pour estimer l'erreur, élément par élément. C'est ce qui permet ensuite de raffiner là où il le faut.

L'indicateur de résidu

Reprenons la démonstration du théorème 11.1 du point de vue de ce que l'on connaît. Pour le problème modèle $-u^{\prime \prime }=f$, avec des éléments linéaires et la charge intégrée exactement, les valeurs nodales sont exactes (théorème 4.3), donc l'erreur $e=u-u_h$ s'annule aux deux extrémités de chaque élément. Sur l'élément, $u_h^{\prime \prime }=0$, donc : la dérivée seconde de l'erreur est , qui est une donnée. C'est le de l'équation: ce que laisse de .

Démonstration. Sur l'élément $e$, $w=u-u_h$ s'annule aux extrémités et $w^{\prime \prime }=-f$. Par intégration par parties, comme à l'étape 2 du théorème 11.1,

par Cauchy–Schwarz puis l'inégalité de Poincaré sur l'élément avec sa constante optimale $h_e/\pi$, admise au chapitre 2. On divise par $\Vert w^{\prime }{\Vert }_{L^2(e)}$ et l'on somme les carrés. $\square$

Avec la constante $h_e/2$, démontrée au théorème 11.1, la borne reste garantie, mais plus pessimiste d'un facteur $\pi /2$. Sur le problème modèle et des maillages uniformes, ${\sum }_e{\eta }_e^2=(h/\pi {)}^2\Vert f{\Vert }_{L^2}^2$, d'où : l'indice d'efficacité vaut pour , puis , , et pour , et il tend vers , le rapport des constantes et de l'exemple 11.1.

Le théorème 11.4 est propre à la dimension 1, où les nœuds sont exacts. En dimension 2, l'erreur ne s'annule plus sur les côtés des éléments, et l'indicateur de résidu prend une forme à deux termes:

où le premier terme est le résidu intérieur (pour des triangles linéaires, $\Delta u_h=0$ et il ne reste que $f$) et le second les sauts du flux normal à travers les côtés $E$ de l'élément. Ce second terme est l'exact analogue du saut de l'effort normal $N_h$ entre deux éléments de barre, que le chapitre 4 avait désigné comme une «erreur visible»: la solution exacte a un flux continu, $u_h$ non. Babuška et Rheinboldt ont établi en 1978 que cet estimateur est fiable et efficace, avec des constantes qui ne dépendent que de la forme des éléments — mais qui ne valent pas 1, et que l'on ne connaît pas exactement. La garantie de (11.13) est un luxe de la dimension 1.

La reconstruction de Zienkiewicz et Zhu

L'autre famille d'estimateurs part de la superconvergence. Puisque la pente $u_h^{\prime }$ est particulièrement précise au milieu des éléments linéaires, construisons à partir de ces valeurs une pente reconstruite $G$, continue et affine par morceaux, et comparons-la à $u_h^{\prime }$. En chaque nœud intérieur $x_i$, est la valeur en de la droite qui passe par les pentes des deux éléments voisins, portées en leurs milieux — sur un maillage uniforme, la simple moyenne des deux pentes —; aux deux nœuds du bord, prend la pente de l'unique élément. C'est la (1987), et l'indicateur qu'on en tire, que nous noterons ZZ, est

En un nœud d'un maillage uniforme, $G-u_h^{\prime }$ vaut, de part et d'autre, plus ou moins la moitié du saut de la pente: l'estimateur ZZ est un estimateur de sauts, déguisé (exercice 11.4). Sur le problème modèle, son indice d'efficacité vaut 1,1943 pour $n=2$, puis 1,0909, 1,0267, 1,0070 et 1,0018 pour $n=32$: il est asymptotiquement exact, parce que la pente reconstruite converge plus vite que $u_h^{\prime }$. En élasticité, la même idée reconstruit les contraintes aux nœuds à partir de leurs valeurs aux points de Gauss. Ne la confondez pas avec la simple du post-traitement (chapitre 9), qui produit des isovaleurs continues sans rien estimer: la reconstruction de Zienkiewicz et Zhu s'en sert pour mesurer l'écart entre le champ lissé et le champ brut, et c'est cet écart qui estime l'erreur. Beaucoup de logiciels proposent cet estimateur, sous le nom de ses auteurs ou sous celui d'«erreur en contrainte».

Sa faiblesse est l'envers de sa force: il repose sur la superconvergence, donc sur un maillage déjà assez fin pour que $u_h$ décrive la solution. Sur un maillage grossier qui rate un phénomène — une couche que trop peu d'éléments traversent —, $G$ et $u_h^{\prime }$ se trompent ensemble, et l'estimateur ZZ peut sous-estimer l'erreur d'un facteur 3. L'exemple 11.4 le montrera.

Le maillage adaptatif

Un estimateur qui dit où est l'erreur permet de ne raffiner que là. Toute méthode adaptative répète la même boucle:

1. résoudre le problème sur le maillage courant;
2. estimer l'erreur par un indicateur ${\eta }_e$ sur chaque élément;
3. marquer les éléments à raffiner;
4. raffiner les éléments marqués, et recommencer — jusqu'à ce que l'estimateur passe sous la tolérance voulue.

Le marquage le plus étudié est celui de Dörfler (1996): on trie les éléments par indicateurs décroissants et l'on marque les premiers, en nombre minimal, jusqu'à ce que leurs ${\eta }_e^2$ représentent une fraction $\theta$ du total ${\eta }^2$. Avec $\theta$ proche de 1, on raffine presque partout; avec $\theta$ petit, on ne raffine que le pire. (Une remarque de notation: dans ce chapitre, $\theta$ désigne trois choses, chaque fois selon l'usage du domaine — l'angle polaire de (11.9), l'indice d'efficacité de (11.11) et la fraction de Dörfler —, et le chapitre 12 l'emploiera encore pour le paramètre du -schéma. De même, l'exposant de gradation de (11.10) n'a rien à voir avec la variable auxiliaire du chapitre 12.) En dimension 1, raffiner un élément, c'est le couper en deux (). En dimension 2, il faut couper les triangles marqués sur le côté d'un voisin non raffiné, ce qui oblige à couper aussi certains voisins: c'est l'objet des algorithmes de bissection par le plus récent sommet ou de raffinement rouge-vert, que nous ne détaillons pas.

Le programme

Prenons le problème $-u^{\prime \prime }=f$ sur $(0,1)$, $u(0)=u(1)=0$, avec une solution fabriquée qui présente une couche raide: $u(x)=\arctan (100(x-0,4))$ moins la corde qui la ramène à zéro aux deux bouts. On en tire par dérivation, et l'on connaît donc la solution exacte — c'est le principe des solutions fabriquées, sur lequel la dernière section reviendra. La couche, en , a une largeur de l'ordre de : la solution y monte de plus de 3 unités, et elle est presque affine partout ailleurs.

Le programme réutilise resoudre(K, F) du chapitre 1, recopiée à l'identique; comme dans tout le cours, les listes commencent à 0 en Python, et x[i] contient le nœud $x_i$ du texte.

from math import atan, pi, sqrt


def resoudre(K, F):
    """Resout K d = F par elimination de Gauss avec pivot partiel."""
    n = len(F)
    M = [ligne[:] for ligne in K]
    c = list(F)
    for k in range(n - 1):
        p = max(range(k, n), key=lambda i: abs(M[i][k]))

La solution fabriquée, sa dérivée et l'intégrale de $f^2$ sur un élément, dont la primitive se vérifie par dérivation:

ALPHA, X0 = 100.0, 0.4
A0, A1 = atan(-ALPHA * X0), atan(ALPHA * (1 - X0))


def u(x):
    return atan(ALPHA * (x - X0)) - ((1 - x) * A0 + x * A1)


def du(x):
    return

La résolution et l'estimation. La charge est intégrée exactement, par parties: sur un élément $[a,b]$, $\int f{N}_1=u^{\prime }(a)-\frac{u(b)-u(a)}{h}$ et , ce qui garantit l'exactitude nodale dont le théorème 11.4 a besoin.

def resoudre_et_estimer(x):
    """Elements lineaires sur les noeuds x; renvoie u_h aux noeuds et les eta_e^2."""
    n = len(x) - 1
    K = [[0.0] * (n + 1) for _ in range(n + 1)]
    F = [0.0] * (n + 1)
    for e in range(n):
        noeuds 

Avant d'écrire la boucle, il faut le marquage. C'est l'objet de la question suivante.

La boucle, enfin, avec $\theta =0,7$ et la fonction marquer corrigée, en partant de quatre éléments égaux:

x = [0.0, 0.25, 0.5, 0.75, 1.0]
for k in range(14):
    uh, eta2 = resoudre_et_estimer(x)
    print(f"{k:2d}  {len(x) - 1:4d} elements  eta = {sqrt(sum(eta2)):.4e}")

 0     4 elements  eta = 7.0524e+01
 1     5 elements  eta = 3.5262e+01
 2     6 elements  eta = 1.7812e+01
 3     7 elements  eta = 9.2731e+00
 4     9 elements  eta = 5.3076e+00
 5    12 elements  eta = 3.1316e+00
 6    16 elements  eta = 1.9507e+00
 7    21 elements  eta = 1.3312e+00
 8    28 elements  eta = 9.1075e-01
 9    39 elements  eta = 6.2071e-01
10    55 elements  eta = 4.1976e-01
11    81 elements  eta = 2.8697e-01
12   112 elements  eta = 1.9745e-01
13   166 elements  eta = 1.3529e-01

L'explorateur suivant refait la boucle sous vos yeux, et la compare au raffinement uniforme.

Trois expériences. D'abord, restez en adaptatif et avancez de l'étape 0 à l'étape 6: les éléments marqués, en couleur dans la bande des indicateurs, sont tous dans la couche ou à côté; le nombre de degrés de liberté passe de 3 à 15 et l'erreur de 90,8 % à 13,0 %. Ensuite, passez en uniforme à la même étape: 255 degrés de liberté et 8,2 % d'erreur — dix-sept fois plus d'inconnues pour un gain modeste —, et une bande d'indicateurs dont presque tous les éléments, loin de la couche, sont inutiles. Enfin, regardez les deux indices d'efficacité en parcourant les étapes: celui du résidu ne descend jamais sous 1, celui de ZZ y descend sur les premiers maillages, dans les deux stratégies.

Vérification et validation

Le chapitre 1 a nommé quatre sources d'erreur et promis de distinguer deux activités que l'on confond souvent. Les voici.

Vérifier le code. Les outils sont ceux du cours. Le test de la pièce contrôle chaque élément. Les solutions exactes du problème modèle, de la barre étagée, du treillis à trois barres et de la console contrôlent l'assemblage et les conditions aux limites. Et la méthode des solutions fabriquées — choisir une fonction $u$, en déduire la charge $f$ et les conditions aux limites, puis vérifier que le code retrouve $u$ avec l'ordre de convergence théorique — contrôle le tout, sur des géométries et des coefficients quelconques. La couche raide de l'exemple 11.4 était une solution fabriquée. Le critère n'est pas que l'erreur soit petite, mais que l'ordre observé soit le bon: une erreur de signe dans un terme d'ordre inférieur peut laisser l'erreur petite sur un maillage et détruire l'ordre, et c'est l'ordre qui la trahit.

Vérifier le calcul. Sur l'ouvrage étudié, il n'y a pas de solution exacte. On dispose de l'extrapolation de Richardson sur trois maillages, des estimateurs a posteriori, et de leur version la plus simple: raffiner et comparer. Le chapitre 1 l'a présentée comme la pratique de l'ingénieur, en prévenant qu'elle peut tromper. Elle trompe précisément aux singularités: sous une force ponctuelle, le déplacement change de 0,11 à chaque raffinement sans jamais se stabiliser; au coin rentrant, la contrainte augmente de 26 % à chaque division du pas. Une grandeur qui ne converge pas n'a pas d'erreur de discrétisation finie — et aucune comparaison de deux maillages ne le révèle si l'on ne regarde pas l'ordre.

Valider le modèle. Aucun raffinement ne valide un modèle: une solution vérifiée à $10^{-6}$ près d'un modèle faux reste fausse. La validation compare des prédictions à des mesures, en tenant compte des incertitudes des deux côtés — sur les données du modèle (modules, charges, conditions d'appui) et sur les mesures. Elle n'a de sens qu'après la vérification: tant que l'erreur de discrétisation n'est pas connue, un accord ou un désaccord avec l'expérience ne prouve rien, puisqu'il peut venir du maillage. Des normes professionnelles, comme les guides de vérification et de validation de l'ASME pour la mécanique des solides numérique, codifient cet ordre des opérations.

Faire confiance à un résultat

Les chapitres de ce cours ont accumulé des contrôles. Rassemblons-les, dans l'ordre où on les fait.

Aucun de ces contrôles n'est coûteux, et chacun attrape une famille d'erreurs que les autres laissent passer. Le dixième est le seul qui sorte de l'ordinateur — et c'est le seul qui réponde à la question que pose le maître d'ouvrage.

Synthèse

• Le lemme de Céa ramène la convergence à une question d'approximation: pour les éléments linéaires en dimension 1, l'interpolant vérifie $|u-I_{h}u{|}_1\le \frac{h}{2}|u{|}_2$ (théorème 11.1), d'où pour le problème modèle; la constante optimale est , la constante asymptotique , retrouvée à près sur la table du cours.

Exercices

Vous pouvez afficher le corrigé directement sous chaque énoncé après avoir cherché la solution.

Références

• Strang, G. et Fix, G., An Analysis of the Finite Element Method, Prentice Hall / Wellesley-Cambridge (estimations d'interpolation, argument de dualité, singularités de coin, test de la pièce et «crimes variationnels»).
• Ern, A. et Guermond, J.-L., Theory and Practice of Finite Elements, Springer (théorie a priori: lemme de Bramble–Hilbert, régularité des maillages, estimations $H^1$ et $L^2$).
• Zienkiewicz, O. C., Taylor, R. L. et Zhu, J. Z., The Finite Element Method: Its Basis and Fundamentals, Butterworth-Heinemann (convergence, test de la pièce, estimation d'erreur par reconstruction et raffinement adaptatif).
• Quarteroni, A., Numerical Models for Differential Problems, Springer (estimations a priori et a posteriori, estimateurs de résidu, maillages adaptatifs).
• Hughes, T. J. R., The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover (convergence des éléments, test de la pièce, superconvergence).
• Cook, R. D., Malkus, D. S., Plesha, M. E. et Witt, R. J., Concepts and Applications of Finite Element Analysis, Wiley (modélisation, singularités en pratique, contrôles et fiabilité des résultats).