La méthode de Galerkin

Objectifs du chapitre

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

• remplacer l'espace $V$ de la formulation faible par un sous-espace de dimension finie $V_h$, et transformer le problème faible en un système linéaire $K\mathbf{d}=\mathbf{F}$ dont vous savez écrire chaque coefficient;
• énoncer les méthodes de Ritz (minimiser l'énergie sur $V_h$) et de Galerkin (tester l'équation par les fonctions de $V_h$), et démontrer qu'elles coïncident quand la forme $a$ est symétrique;
• calculer à la main, en fractions exactes, une approximation de Ritz à une, deux ou trois fonctions polynomiales globales, et en mesurer l'erreur en énergie;
• démontrer l'orthogonalité de Galerkin et le lemme de Céa: la solution de Galerkin est la meilleure approximation de $u$ dans $V_h$ au sens de la norme d'énergie;
• démontrer que la matrice de rigidité est symétrique définie positive, expliquer pourquoi des fonctions de base locales la rendent creuse, et compter ses coefficients non nuls;
• calculer le conditionnement de la matrice des fonctions chapeau, montrer qu'il croît comme $h^{-2}$, et le comparer à celui d'une base polynomiale globale.

Le chapitre 2 a établi la formulation faible et ses propriétés: l'espace $V$, la forme bilinéaire $a(u,v)$, la forme linéaire $\ell (v)$, la norme d'énergie et le principe du minimum de l'énergie potentielle. Ce chapitre fait le pas qui rend tout cela calculable: on renonce à chercher $u$ dans un espace de dimension infinie, et l'on choisit la meilleure fonction dans un espace de dimension finie. Toute la méthode des éléments finis est un cas particulier de ce qui suit; le chapitre 4 n'aura plus qu'à choisir l'espace et à organiser le calcul.

Du problème faible au système linéaire

Rappelons le problème du chapitre 2 sous sa forme abstraite. On se donne un espace de Hilbert $V$ — pour le problème modèle, $V=H_0^1(0,1)$ —, une forme bilinéaire $a$ et une forme linéaire $\ell$ sur $V$, et l'on cherche

Pour le problème modèle $-u^{\prime \prime }=f$ sur $(0,1)$, $u(0)=u(1)=0$, les deux formes valent

Nous utiliserons les deux propriétés établies au chapitre 2, avec les constantes $M$ et $\alpha$ de la définition 2.7: la continuité, $|a(u,v)|\le M\Vert u{\Vert }_V\Vert v{\Vert }_V$, et la coercivité, $a(v,v)\ge \alpha \Vert v{\Vert }_V^2$ avec , pour tous . Lorsque est de plus symétrique, , elle est un produit scalaire sur , et la est la norme associée, . Pour le problème modèle, .

Le problème (3.1) porte sur une infinité d'inconnues — une fonction entière — et une infinité d'équations — une par fonction test. L'idée de Galerkin est de réduire les deux infinités à la même dimension finie.

L'indice $h$ rappelle que, pour les éléments finis, $V_h$ dépendra d'un maillage de pas $h$; mais rien dans ce chapitre n'exige un maillage. $V_h$ peut être engendré par des polynômes sur tout le domaine, par des sinus, ou par des fonctions chapeau: c'est précisément ce choix que nous allons comparer.

Le système linéaire

Choisissons une base ${\phi }_1,\dots ,{\phi }_N$ de $V_h$. Toute fonction de $V_h$ s'écrit de façon unique comme combinaison de ces fonctions; en particulier

et les $N$ coefficients $d_j$, rangés dans le vecteur $\mathbf{d}$, sont les inconnues. Pour les fonctions chapeau du chapitre 1, ce sont les valeurs nodales; pour une base polynomiale, ce sont des coefficients sans signification ponctuelle.

Il suffit d'écrire (3.2) pour les $N$ fonctions de base. En effet, si $a(u_h,{\phi }_i)=\ell ({\phi }_i)$ pour $i=1,\dots ,N$, alors pour quelconque, la linéarité de et de en leur second argument donne . En reportant (3.3) dans ces équations et en utilisant la linéarité en le premier argument,

C'est un système linéaire $N\times N$:

La matrice $K$ est la matrice de rigidité et $\mathbf{F}$ le vecteur des charges, comme au chapitre 1. L'ordre des indices dans $K_{ij}=a({\phi }_j,{\phi }_i)$ n'est pas une coquetterie: la ligne est l'équation testée par , la colonne la contribution de l'inconnue . Il n'importe pas tant que est symétrique, et il importe dès qu'elle ne l'est plus (exercice 3.4).

Résumons la démarche, car elle est celle de tous les chapitres suivants: (1) choisir $V_h$ et sa base; (2) calculer les intégrales $a({\phi }_j,{\phi }_i)$ et $\ell ({\phi }_i)$; résoudre ; reconstruire par (3.3) et en déduire les grandeurs utiles. Le seul choix que fait l'ingénieur est le premier. Tout le reste en découle, y compris — nous le verrons — la qualité du résultat, la structure de et le coût de la résolution.

La méthode de Ritz: minimiser l'énergie sur un sous-espace

Le chapitre 2 a démontré que, lorsque $a$ est symétrique et coercive, la solution de (3.1) est l'unique minimum de l'énergie potentielle totale

Pour la barre, $\frac{1}{2}a(v,v)$ est l'énergie de déformation et $\ell (v)$ le travail des charges. Cette caractérisation suggère une méthode d'approximation d'une autre nature que (3.2): au lieu de minimiser $\Pi$ sur tout $V$, ce qui est impossible, on la minimise sur le sous-espace $V_h$.

Walther Ritz, physicien né à Sion en 1878, publie cette méthode en 1909, l'année de sa mort, et l'applique notamment aux vibrations d'une plaque carrée libre; Boris Galerkin, ingénieur à Saint-Pétersbourg, publie en 1915 la formulation par l'équation testée. Pour les problèmes symétriques, les deux méthodes donnent la même approximation.

Démonstration. Soient $u_h,w\in V_h$ et $t$ réel. En développant par bilinéarité, et en regroupant $a(u_h,w)$ et $a(w,u_h)$ grâce à la symétrie,

Galerkin entraîne Ritz. Si $u_h$ vérifie (3.2), le terme en $t$ est nul pour tout $w\in V_h$. Toute $v_h\in V_h$ s'écrit avec ; avec , (3.7) donne , avec égalité seulement si , c'est-à-dire par coercivité. Donc est l'unique minimum.

Ritz entraîne Galerkin. Si $u_h$ minimise $\Pi$ sur $V_h$, alors pour tout $w\in V_h$ la fonction $g(t)=\Pi (u_h+tw)$, un polynôme du second degré en , atteint son minimum en ; sa dérivée y est nulle: . C'est (3.2).

Coordonnées. Pour $v_h={\sum }_j{y}_j{\phi }_j$, $a(v_h,v_h)={\sum }_{i,j}{y}_{i}a({\phi }_j,{\phi }_i)y_j={\mathbf{y}}^{\mathsf{T}}K\mathbf{y}$ et . Enfin, en prenant dans (3.2), , d'où , qui est (3.6).

La méthode de Ritz a deux mérites pédagogiques. Elle dit ce que l'on optimise — l'énergie —, ce qui permet de comparer deux approximations sans connaître $u$: celle qui a l'énergie la plus basse est la meilleure. Et l'identité du chapitre 2, $\Vert u-v{\Vert }_E^2=2(\Pi (v)-\Pi (u))$, montre que minimiser l'énergie, c'est minimiser l'erreur en énergie. Le lemme de Céa, plus bas, le redira sans passer par $\Pi$.

La méthode de Galerkin, elle, va plus loin. Elle n'a pas besoin d'une énergie: (3.2) a un sens pour une forme $a$ non symétrique, comme celle d'un problème de convection–diffusion (exercice 3.4), où $\Pi$ n'a plus de minimum qui caractérise $u$. C'est pourquoi on parle aujourd'hui de méthode de Galerkin dans tous les cas, et de Ritz quand on veut insister sur l'énergie.

Bases polynomiales globales

Commençons par l'espace le plus naturel: des polynômes définis sur tout l'intervalle. Pour être dans $V=H_0^1(0,1)$, une fonction d'essai doit s'annuler en $0$ et en $1$; c'est le cas de toute fonction de la forme $x(1-x)p(x)$ avec $p$ polynôme. On prend la base

qui engendre les polynômes de degré au plus $N+1$ nuls aux deux extrémités. La première fonction est $x(1-x)$, la parabole de l'explorateur du chapitre 2; avec les deux premières, $V_h$ contient toutes les fonctions $x(1-x)(\beta +\gamma x)$, soit, quand , les fonctions avec .

Les coefficients de $K$ sont des intégrales de polynômes, donc des nombres rationnels. Avec ${\phi }_k^{\prime }=k{x}^{k-1}-(k+1)x^k$ et ,

Pour $N=3$:

Tous les coefficients sont non nuls: chaque fonction de base «voit» toutes les autres, puisque toutes sont non nulles sur tout l'intervalle (figure 3.1, en haut). Nous y reviendrons.

Le calcul exact se programme en quelques lignes avec le module fractions de la bibliothèque standard, qui fait l'arithmétique des rationnels sans arrondi. La fonction resoudre(K, F) est celle du chapitre 1, inchangée; elle fonctionne telle quelle avec des fractions, puisqu'elle n'utilise que les quatre opérations et abs. Comme dans tout le cours, les listes Python commencent à 0 alors que le texte numérote les fonctions de base à partir de 1: ici d[0] contient $d_1$, et la boucle sur range(1, m + 1) parcourt les indices $k=1,\dots ,m$ du texte, avec $m=N$.

from fractions import Fraction


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

def ritz_polynomial(m):
    """Ritz pour -u'' = 12 x^2, base phi_k = x^k (1 - x), k = 1..m (indices 1..m)."""
    K = [[Fraction(i * j, i + j - 1) - Fraction(i * (j + 1) + j * (i + 1), i + j)
          + Fraction((i + 1) * (j + 1), i + j + 1) 

1 ['9/5'] 36/175 0.45356
2 ['4/5', '2'] 1/175 0.07559
3 ['1', '1', '1'] 0 0.00000

Le problème modèle avec la charge sinusoïdale

Avec $f={\pi }^2\sin \pi x$, la solution exacte $u=\sin \pi x$ n'est pas un polynôme: aucun $V_h$ polynomial ne la contient, et l'erreur ne s'annule jamais. Les charges $F_k={\pi }^2{\int }_0^1{x}^k(1-x)\sin (\pi x)\mathrm{d}x$ s'obtiennent par la récurrence sur , avec et (deux intégrations par parties). Avec une seule fonction, et , la valeur de Ritz du chapitre 2.

La précision spectaculaire des polynômes tient à la régularité de la solution: $\sin \pi x$ est analytique, et ses approximations polynomiales convergent plus vite que toute puissance de $1/N$. Elle fond dès que la solution l'est moins. Avec la charge $f=1$ sur $(0,\frac{1}{2})$ et $f=0$ sur , la solution exacte, à gauche et à droite, a une dérivée seconde qui saute en . Le même calcul, en fractions exactes, donne pour les polynômes les erreurs , , et pour , et pour les chapeaux , , et . Les polynômes restent devant, mais leur erreur ne diminue plus que d'un facteur 2 à 4 quand double: la convergence est devenue , comme celle des chapeaux. Une solution singulière — un coin rentrant, une charge ponctuelle — aggrave encore les choses, et c'est le cas courant en ingénierie (chapitre 11).

Bases locales: les fonctions chapeau

Si les polynômes globaux sont si précis, pourquoi la méthode des éléments finis ne les utilise-t-elle pas? Le tableau de l'exemple 3.3 contient la première raison: le conditionnement. Il y en a trois autres.

La matrice est pleine. Toutes les fonctions de (3.8) sont non nulles sur tout $(0,1)$, donc toutes les intégrales $a({\phi }_j,{\phi }_i)$ sont a priori non nulles: il faut stocker $N^2$ coefficients et la résolution coûte de l'ordre de $\frac{2}{3}{N}^3$ opérations. Avec inconnues, ce qui est modeste en dimension 3, il faudrait stocker nombres.

La géométrie. En dimension 1, l'intervalle est simple. En dimension 2 ou 3, il faut des fonctions qui s'annulent sur le bord ${\Gamma }_D$ d'un domaine quelconque — une plaque percée, une pièce mécanique —, et il n'existe pas de recette générale pour construire une base polynomiale globale qui le fasse.

Les données par morceaux. Un matériau qui change à une interface, une charge appliquée sur une partie seulement, une section qui saute: autant de lieux où la solution perd sa régularité, et où une base globale, très raide, propage l'erreur à tout le domaine.

Les fonctions chapeau du chapitre 1 répondent aux trois difficultés à la fois. Sur un maillage $0=x_0<x_1<\cdots <x_n=1$, la fonction chapeau , , est affine sur chaque élément, vaut 1 au nœud et 0 à tous les autres nœuds. Elle est continue et dérivable par morceaux, donc dans (chapitre 2), et nulle en et en : l'espace

des fonctions continues, affines par morceaux et nulles aux extrémités, est un sous-espace de $V=H_0^1(0,1)$, de dimension $N=n-1$. La méthode de Galerkin dans cet espace est exactement la méthode du chapitre 1: les coefficients $d_j$ sont les valeurs nodales , et sur un maillage uniforme de pas , (3.4) donne

Ce qui était au chapitre 1 une recette — minimiser l'énergie sur les lignes brisées — est maintenant un cas particulier de la définition 3.1, et tout ce que nous allons démontrer pour Galerkin vaut pour elle.

La figure 3.1 met les deux bases côte à côte. Les fonctions chapeau ne sont pas «meilleures» une à une: avec quatre fonctions, elles approchent $\sin \pi x$ plus de quarante fois moins bien que quatre polynômes. Leur force est collective. On peut en ajouter autant qu'on veut sans dégrader la base; on peut les concentrer là où la solution varie vite; on peut les construire sur n'importe quelle géométrie à partir d'un maillage; et chacune n'interagit qu'avec ses voisines, ce qui rend $K$ creuse. L'explorateur suivant permet de comparer les deux bases sur deux charges.

L'explorateur résout les deux systèmes à chaque mouvement du curseur et mesure les erreurs par l'identité de Pythagore du théorème 3.2, sans intégrer l'écart. Trois lectures. Avec la charge sinusoïdale, la courbe polynomiale se confond avec la solution exacte dès $N=3$, alors que la ligne brisée bleue reste visiblement anguleuse; les paliers de l'erreur polynomiale (de $N=1$ à 2, de 3 à 4, de 5 à 6) sont ceux de l'exemple 3.3. Avec la charge sur la moitié gauche, les paliers changent de place — la partie symétrique de la solution, due à la moyenne $\frac{1}{2}$ de la charge, est captée exactement dès $N=1$, et seules les fonctions qui apportent une partie antisymétrique comptent —, et l'avance des polynômes fond: à , elle n'est plus que d'un facteur 5. , enfin, ne dépend pas de la charge: il vaut pour six polynômes et 19,20 pour six chapeaux.

Orthogonalité de Galerkin et meilleure approximation

Nous arrivons au résultat central du chapitre. Il dit ce que la méthode de Galerkin calcule, indépendamment de la base choisie, et il ramène la question «la méthode converge-t-elle?» à une question d'approximation pure.

Démonstration. Comme $V_h\subset V$, toute $v_h\in V_h$ est une fonction test admissible de (3.1): $a(u,v_h)=\ell (v_h)$. Par (3.2), . En soustrayant et par linéarité en le premier argument, : c'est (3.11).

Pour (3.12), écrivons $u=u_h+e$ et développons: $a(u,u)=a(u_h,u_h)+2a(e,u_h)+a(e,e)$, où l'on a utilisé la symétrie. Le terme croisé est nul par (3.11), puisque . Enfin par (3.6).

L'hypothèse $V_h\subset V$ — la conformité — est la seule qui serve, et elle sert de façon essentielle: c'est elle qui autorise à tester l'équation continue par $v_h$. Les méthodes non conformes, où $V_h\not\subset V$, perdent (3.11) et demandent une autre analyse.

L'identité (3.12) a une conséquence pratique immédiate: on peut mesurer l'erreur en énergie sans calculer l'erreur, dès que l'on connaît $a(u,u)$ — ce qui est le cas pour une solution de référence, et ce qui a servi dans tous les exemples précédents. Elle montre aussi que l'énergie de l'approximation est toujours inférieure à celle de la solution exacte: $a(u_h,u_h)\le a(u,u)$. Pour une structure, $a(u_h,u_h)={\mathbf{F}}^{\mathsf{T}}\mathbf{d}$ est le travail des charges, et sous une charge donnée, un travail plus petit veut dire des déplacements plus petits en moyenne: . C'est une propriété générale des approximations conformes en déplacements, que le chapitre 9 retrouvera avec le triangle à déformation constante.

Le lemme de Céa

L'orthogonalité (3.11) a une lecture géométrique simple quand $a$ est symétrique: $a$ est alors un produit scalaire, et (3.11) dit que l'erreur $u-u_h$ est orthogonale au sous-espace $V_h$ pour ce produit scalaire. Autrement dit, $u_h$ est la de sur (figure 3.2). Et dans un espace euclidien, la projection orthogonale est le point du sous-espace le plus proche.

Démonstration. Soit $v_h\in V_h$. Écrivons $u-v_h=(u-u_h)+(u_h-v_h)$, où . En développant par bilinéarité et symétrie,

le terme croisé étant nul par (3.11). Donc $\Vert u-v_h{\Vert }_E\ge \Vert u-u_h{\Vert }_E$ pour tout , avec égalité si et seulement si : le minimum est atteint en , et seulement là. Si , on peut prendre , et le minimum vaut 0: , donc par coercivité.

Le théorème répond à la question posée à l'exemple 3.1: avec $f=1$ et la base $x(1-x)$, ou avec $f=6x$ et la famille $x(1-x)(1+cx)$, l'approximation est exacte parce que la solution est dans $V_h$. Il dit aussi ce que la méthode promet : rien sur l'erreur en un point, rien sur l'erreur maximale. L'exemple 3.2 l'a montré: en passant de une à deux fonctions, l'erreur en énergie a été divisée par 6 et l'erreur au milieu n'a pas bougé.

Lorsque $a$ n'est pas symétrique, il n'y a plus de produit scalaire ni de projection orthogonale, mais il reste une quasi-optimalité.

Démonstration. Soit $v_h\in V_h$. Par coercivité, puis par l'orthogonalité (3.11) appliquée à $u_h-v_h\in V_h$, puis par continuité,

Si $u-u_h=0$, il n'y a rien à démontrer; sinon, on divise par $\alpha \Vert u-u_h{\Vert }_V$, et l'on prend le minimum sur $v_h$. Dans le cas symétrique, coercivité et continuité donnent , et (3.13) entraîne

Jean Céa a publié ce résultat en 1964, dans sa thèse sur l'approximation variationnelle des problèmes aux limites. Sa portée est considérable, parce qu'il sépare deux questions. La méthode de Galerkin est-elle stable? Oui, avec une constante $M/\alpha$ qui ne dépend que du problème continu, ni de $V_h$, ni de $h$. Converge-t-elle, et à quelle vitesse? Exactement aussi vite que $V_h$ est capable d'approcher $u$ — une question de théorie de l'approximation, qui ne fait plus intervenir l'équation. Pour les fonctions chapeau, il suffira de majorer l'écart entre $u$ et la ligne brisée qui l'interpole aux nœuds: c'est ce que fera le chapitre 11 pour obtenir l'ordre 1 en énergie de la table du cours, et l'exercice 3.5 montre que, pour en dimension 1, cette ligne brisée est même exactement .

La matrice de rigidité est symétrique définie positive

Le système (3.4) a-t-il une solution, et une seule? Pour les exemples ci-dessus, nous l'avons résolu sans nous poser la question. La réponse tient en une ligne de calcul, et elle a des conséquences pratiques sur le choix du solveur.

Démonstration. 1. $K_{ij}=a({\phi }_j,{\phi }_i)=a({\phi }_i,{\phi }_j)=K_{ji}$.

2. Soit $\mathbf{y}\ne 0$ et $v_h={\sum }_j{y}_j{\phi }_j$. Comme les sont linéairement indépendantes, . Par bilinéarité,

Si $K\mathbf{y}=0$, alors ${\mathbf{y}}^{\mathsf{T}}K\mathbf{y}=0$, donc $\mathbf{y}=0$: le noyau de $K$ est réduit à zéro et $K$ est inversible. $\square$

La démonstration mérite d'être relue pour ce qu'elle dit physiquement: ${\mathbf{y}}^{\mathsf{T}}K\mathbf{y}$ est deux fois l'énergie de déformation du champ $v_h$ dont les degrés de liberté sont $\mathbf{y}$. Une matrice de rigidité est définie positive parce que tout déplacement non nul déforme la structure et lui coûte de l'énergie. Pour le problème modèle, $a(v,v)={\int }_0^1(v^{\prime }{)}^2$ ne s'annule que si , c'est-à-dire si est constante, et la seule constante qui vérifie est la fonction nulle: c'est la condition de Dirichlet qui rend coercive.

Retirons cette condition — une barre libre aux deux bouts, chargée par des forces qui s'équilibrent. La fonction constante $v=1$ est alors un déplacement d'ensemble, un mode rigide, qui ne coûte aucune énergie: $a(1,1)=0$. Si l'on garde toutes les fonctions chapeau, y compris celles des deux nœuds extrêmes, leur somme vaut 1 sur tout l'intervalle, et le vecteur $\mathbf{y}=(1,\dots ,1{)}^{\mathsf{T}}$ vérifie $K\mathbf{y}=0$: chaque ligne de somme à zéro, ce que la question 1.5 faisait vérifier. La matrice est seulement positive, et singulière.

Les conséquences pratiques du théorème 3.5 sont directes. Une matrice symétrique définie positive se factorise par Cholesky, $K=L{L}^{\mathsf{T}}$, sans pivotage et avec deux fois moins d'opérations que l'élimination de Gauss (Algèbre linéaire, chapitre 11, pour la définie positivité; Analyse numérique, chapitre 3, pour la factorisation). Les grands codes résolvent aussi ces systèmes par la méthode itérative du gradient conjugué, qui ne s'applique qu'aux matrices symétriques définies positives et que ce cours ne développe pas. Inversement, un $K$ non symétrique, ou un pivot négatif rencontré pendant une factorisation de Cholesky, sont des signaux d'alarme: une erreur dans les matrices élémentaires, une condition aux limites mal imposée, ou un matériau à module négatif.

Pour une forme non symétrique et coercive, la partie 2 du théorème reste vraie: $K$ n'est plus symétrique, mais ${\mathbf{y}}^{\mathsf{T}}K\mathbf{y}>0$ suffit à la rendre inversible. L'exercice 3.4 en donne un exemple.

Fonctions locales, matrice creuse

La deuxième propriété structurale de $K$ ne dépend pas de la forme $a$ mais de la base.

Pour une forme intégrale comme $a(u,v)={\int }_0^1{u}^{\prime }{v}^{\prime }$, le coefficient $K_{ij}=\int {\phi }_j^{\prime }{\phi }_i^{\prime }$ est nul dès que les supports de et ne se chevauchent pas, ou ne se touchent qu'en un point, puisque l'intégrande est alors nul presque partout. La fonction chapeau a pour support ; deux fonctions chapeau ont des supports qui se chevauchent sur un élément entier si et seulement si leurs nœuds sont voisins. Donc

et $K$ est tridiagonale. Elle a $N$ coefficients diagonaux et $2(N-1)$ coefficients voisins, soit $3N-2$ coefficients non nuls sur $N^2$: 22 sur 64 pour $N=8$ (figure 3.3), 298 sur 10 000 pour $N=100$, et une proportion qui tend vers zéro comme . C'est la propriété qui rend la méthode : le stockage et le coût de la résolution — l'algorithme de Thomas (, chapitre 3), cité au chapitre 1, en opérations — croissent comme , et non comme et .

Le même raisonnement vaut en toute dimension: $K_{ij}\ne 0$ seulement si les nœuds $i$ et $j$ appartiennent à un même élément. Sur le maillage en triangles de la plaque carrée du chapitre 8, où chaque carré est coupé selon la diagonale qui monte vers la droite, un nœud intérieur a six voisins, et chaque ligne de $K$ compte sept positions dans la structure creuse, dont cinq coefficients non nuls seulement, quel que soit le nombre de nœuds: les couplages le long des diagonales sont exactement nuls, et l'on retrouve le schéma à cinq points (chapitre 8; exemple A.6). La matrice n'est plus tridiagonale, mais elle reste creuse, et sa largeur de bande dépend de la numérotation des nœuds: numéroter une plaque rectangulaire ligne par ligne dans le sens de sa plus petite dimension donne une bande plus étroite que dans l'autre sens. Les logiciels renumérotent automatiquement les nœuds avant la factorisation pour cette raison.

La base globale (3.8), à l'inverse, donne une matrice pleine quel que soit $N$. Pour un problème symétrique en dimension 1, on peut faire mieux en orthogonalisant la base — des polynômes de Legendre intégrés donnent une matrice de rigidité diagonale pour le problème modèle —, et c'est l'idée des méthodes spectrales et des éléments de degré élevé (méthode $p$, chapitre 11). Mais cela suppose une géométrie simple, et ne règle pas la question des solutions peu régulières.

Le conditionnement croît comme $h^{-2}$

La troisième propriété de $K$ est quantitative: son conditionnement. L'Analyse numérique (chapitre 4) a montré qu'en résolvant $K\mathbf{d}=\mathbf{F}$ en virgule flottante, l'erreur relative sur $\mathbf{d}$ peut atteindre $\kappa (K)$ fois l'erreur relative commise sur les données, où $\kappa (K)=\Vert K\Vert \Vert K^{-1}\Vert$. Pour une matrice symétrique définie positive et la norme euclidienne, , le rapport de ses valeurs propres extrêmes. Pour la matrice des fonctions chapeau, ces valeurs propres se calculent exactement.

Démonstration. La formule $\sin (\theta -\beta )+\sin (\theta +\beta )=2\sin \theta \cos \beta$, avec $\theta =k\pi x_i$ et , donne pour

Pour $i=1$ et $i=n-1$, il manque dans $K{\mathbf{s}}_k$ les termes en $x_0=0$ et ; mais et , de sorte que la formule vaut aussi. C'est le calcul de l'exercice 1.5 avec au lieu de . Les vecteurs sont non nuls (leur composante , , l'est pour ) et associés à des valeurs propres distinctes, puisque est strictement croissant sur et que parcourt cet intervalle: ce sont toutes les valeurs propres. La plus petite est , la plus grande , avec , donc . Le quotient donne (3.16), et quand donne l'équivalent.

Les deux valeurs propres extrêmes ont une lecture physique. Le quotient de Rayleigh ${\mathbf{y}}^{\mathsf{T}}K\mathbf{y}/{\mathbf{y}}^{\mathsf{T}}\mathbf{y}$ est l'énergie d'un champ rapportée à la somme des carrés de ses valeurs nodales. Le mode le plus lisse, $\sin \pi x$, a une énergie faible par rapport à ses $n-1$ valeurs nodales: ${\lambda }_1\approx {\pi }^{2}h$. Le mode le plus , qui alterne de signe d'un nœud au suivant, a une pente de l'ordre de sur chaque élément: . Le rapport de ces deux énergies est le conditionnement, et il croît comme parce que la forme contient deux dérivées.

Le programme suivant vérifie le théorème: il construit $K$, évalue le quotient de Rayleigh sur ${\mathbf{s}}_1$ et ${\mathbf{s}}_{n-1}$, et imprime ${\lambda }_{\min }$, ${\lambda }_{\max }$, et .

from math import pi, sin


def rayleigh(K, v):
    """Quotient de Rayleigh v.Kv / v.v."""
    Kv = [sum(K[i][j] * v[j] for j in range(len(v))) for i in range(len(v))]
    return sum(a * b for a, b in zip(v, Kv)) / sum(a * a 

   4  2.34315    13.6569        5.83  0.36428
   8  1.21793    30.7821       25.27  0.39491
  16  0.61487    63.3851      103.09  0.40268
  32  0.30818   127.6918      414.35  0.40463
  64  0.15418   255.8458     1659.38  0.40512
 128  0.07710   511.9229     6639.52  0.40524

La dernière colonne tend vers $4/{\pi }^2=0,40528$: le conditionnement est multiplié par 4 chaque fois que $h$ est divisé par 2 — exactement par $4,019$ entre $n=16$ et $n=32$, par $4,005$ entre $n=32$ et . Pour , , la valeur de l'exemple 3.3. La matrice des différences finies (1.11) étant , elle a le même conditionnement, et l'on retrouve la valeur pour citée au chapitre 1.

Que faut-il en penser? Une croissance en $h^{-2}$ n'est pas une catastrophe, et c'est même le comportement normal de toute discrétisation d'un problème du second ordre, en dimension 1, 2 ou 3. Le chapitre 1 a montré que l'erreur d'arrondi qui en résulte, de l'ordre de ${\kappa }_2{\epsilon }_{\text{mach}}$ dans le pire des cas, ne rattrape l'erreur de discrétisation qu'au-delà de $10^4$ éléments en dimension 1, une finesse qu'aucun calcul pratique n'atteint. Le conditionnement compte davantage pour les , dont le nombre d'itérations croît avec , et c'est pourquoi les grands codes utilisent des préconditionneurs. Il compte aussi quand des coexistent dans un modèle — un élément d'acier collé à un élément de caoutchouc, un ressort d'appui de raideur énorme pour imiter un encastrement (la méthode de pénalisation du chapitre 4) — car est alors multiplié par le rapport des raideurs.

La base polynomiale globale, elle, a un conditionnement qui croît exponentiellement avec $N$: $1,687\cdot 10^6$ pour six fonctions, $1,04\cdot 10^9$ pour huit (figure 3.3). Avec seize fonctions, on serait au-delà de ce que la double précision peut résoudre avec un seul chiffre exact. La raison est visible sur la figure 3.1: les fonctions $x^k(1-x)$ se ressemblent de plus en plus quand grandit, et une base presque liée donne une matrice presque singulière.

Synthèse

• La méthode de Galerkin cherche $u_h$ dans un sous-espace $V_h\subset V$ de dimension finie $N$ et teste l'équation faible par les fonctions de ce même espace; avec une base ${\phi }_1,\dots ,{\phi }_N$, elle conduit au système , , .

Exercices

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

Références

• Strang, G. et Fix, G. J., An Analysis of the Finite Element Method, Prentice Hall / Wellesley-Cambridge, chap. 1 (méthode de Ritz–Galerkin, projection en énergie, conditionnement de la matrice de rigidité).
• Ern, A. et Guermond, J.-L., Theory and Practice of Finite Elements, Springer (approximation de Galerkin, lemme de Céa, propriétés des matrices de rigidité).
• Quarteroni, A., Numerical Models for Differential Problems, Springer (méthode de Galerkin, quasi-optimalité, conditionnement et méthodes spectrales).
• Hughes, T. J. R., The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover, chap. 1 (Galerkin sur le problème modèle, propriété de meilleure approximation, exactitude nodale en dimension 1).
• Zienkiewicz, O. C., Taylor, R. L. et Zhu, J. Z., The Finite Element Method: Its Basis and Fundamentals, Butterworth-Heinemann, chap. 3 (résidus pondérés, Galerkin et principes variationnels).
• Ritz, W., «Über eine neue Methode zur Lösung gewisser Variationsprobleme der mathematischen Physik», Journal für die reine und angewandte Mathematik, 1909; Céa, J., «Approximation variationnelle des problèmes aux limites», Annales de l'Institut Fourier, 1964.