Éléments d'ordre supérieur et intégration numérique

Objectifs du chapitre

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

construire les fonctions de forme de Lagrange de degré 1, 2 et 3 sur l'élément de référence $[- 1, 1]$ , et en vérifier les propriétés — valeur 1 en leur nœud et 0 aux autres, partition de l'unité, reproduction des polynômes de degré $p$ ;
écrire la transformation isoparamétrique $x (ξ)$ d'un élément, calculer son jacobien $J = d x / d ξ$ , et en déduire les dérivées $d N_{a} / d x$ et l'élément de longueur ; reconnaître un élément dont le jacobien s'annule;

Le chapitre 4 a bâti l'élément linéaire, et le chapitre 6 l'élément de poutre d'Hermite, en intégrant tout à la main. Deux choses changent ici. On monte en degré, ce qui améliore la précision à nombre d'inconnues égal; et l'on cesse d'intégrer à la main, ce qui est la condition pour traiter des coefficients variables, des éléments déformés et, au chapitre 10, des éléments bidimensionnels de forme quelconque. Les deux outils qui le permettent — l'élément de référence et la quadrature de Gauss — sont ceux de tout code d'éléments finis.

Pourquoi monter en degré

Rappelons la table fixée du cours pour les éléments linéaires sur le problème modèle $- u^{''} = π^{2} sin (π x)$ , $u (0) = u (1) = 0$ , de solution $u (x) = sin (π x)$ et de norme d'énergie :

$n$	$∥ u - u_{h} ∥_{E}$	rapport	$∥ u - u_{h} ∥_{L^{2}}$

L'erreur en énergie est divisée par 2 quand le maillage est divisé par 2: c'est l'ordre 1. Pour gagner un chiffre significatif sur les pentes — c'est-à-dire sur les déformations et les contraintes —, il faut dix fois plus d'éléments. Avec 32 éléments, l'erreur relative en énergie vaut encore $0, 06295/2, 2214 = 2, 8$ %.

La raison est géométrique. Sur un élément, une fonction affine ne peut pas suivre la courbure de $u$ : la dérivée de $u_{h}$ est constante par élément, et elle approche $u^{'}$ comme un escalier approche une pente. Le chapitre 11 démontrera l'estimation $∥ u - u_{h} ∥_{E} \leq C h^{p} ∣ u ∣_{p + 1}$ , où est le degré des polynômes utilisés sur chaque élément. Elle dit que l'exposant de est . Il y a donc deux manières de réduire l'erreur: diminuer (on raffine le maillage, c'est le ) ou augmenter (on enrichit les éléments, c'est le ). Ce chapitre construit les outils du second, et mesure ce qu'il rapporte.

Construire un élément de degré 2 ou 3 demande plus de nœuds par élément, des fonctions de forme plus riches et des intégrales plus lourdes. Pour ne pas refaire ce travail sur chaque élément du maillage, on le fait une fois, sur un élément fixe, et l'on transporte le résultat.

L'élément de référence et la transformation isoparamétrique

Un seul élément pour tous

Le chapitre 4 a déjà introduit l'élément de référence (définition 4.2) pour l'élément linéaire; on le reprend ici pour en faire l'outil général. Sur un élément $e = [x_{1}^{e}, x_{2}^{e}]$ de longueur $h_{e}$ , les fonctions de forme linéaires dépendent de $x_{1}^{e}$ et de . Toutes ces fonctions sont pourtant , vue à une autre échelle et à un autre endroit. On les écrit donc une fois pour toutes sur l'intervalle fixe , de coordonnée , puis on les transporte sur chaque élément par un changement de variable.

Dans ce chapitre, les nœuds locaux sont numérotés de gauche à droite, $a = 1, \dots, p + 1$ , et équidistants sur l'élément de référence: $ξ_{a} = - 1 + 2 (a - 1) / p$ . En Python, où les listes commencent à 0, le nœud local $a$ est ; et le nœud global qui porte la valeur du texte, avec au bout d'un maillage de éléments de degré , est l'indice de la liste. Chaque listing le rappelle dans sa docstring.

Pour l'élément linéaire, les deux fonctions de forme sont

N_{1} (ξ) = \frac{1 - ξ}{2}, N_{2} (ξ) = \frac{1 + ξ}{2},

et la transformation (7.1) est affine:

x (ξ) = x_{1}^{e} \frac{1 - ξ}{2} + x_{2}^{e} \frac{1 + ξ}{2} = \frac{x _{1}^{e} + x _{2}^{e}}{2} + \frac{h _{e}}{2} ξ, J = \frac{h _{e}}{2} . (7.2)

Le centre de $[- 1, 1]$ va au milieu de l'élément, et la longueur 2 devient $h_{e}$ : le jacobien est le facteur d'échelle des longueurs, constant pour une transformation affine. C'est exactement le changement de variable que le chapitre 8 d'Analyse numérique utilise pour transporter une formule de quadrature de $[- 1, 1]$ sur $[a, b]$ .

Dériver et intégrer sur l'élément de référence

Toutes les intégrales élémentaires se ramènent à l'élément de référence par deux règles. La première est la dérivation composée: puisque $N_{a}$ est une fonction de $ξ$ et que $ξ$ est une fonction de $x$ ,

\frac{d N _{a}}{d x} = \frac{d N _{a}}{d ξ} \frac{d ξ}{d x} = \frac{1}{J} \frac{d N _{a}}{d ξ}, d x = J d ξ . (7.3)

La seconde est le changement de variable dans l'intégrale, qui exige que $x (ξ)$ soit une bijection croissante de $[- 1, 1]$ sur l'élément, c'est-à-dire $J > 0$ partout. Avec les deux, la matrice de rigidité élémentaire de la barre, $k_{ab}^{e} = \int_{e} E A N_{a}^{'} N_{b}^{'} d x$ , devient

k_{ab}^{e} = \int_{- 1}^{1} E A (x (ξ)) \frac{d N _{a}}{d ξ} \frac{d N _{b}}{d ξ} \frac{1}{J ( ξ )} d ξ, f_{a}^{e} = \int_{- 1}^{1} q (x (ξ)) N_{a} (ξ) J (ξ) d ξ . (7.4)

Le facteur $1/ J$ de la rigidité est le produit de deux facteurs $1/ J$ (un par dérivée) et d'un facteur $J$ (l'élément de longueur). Pour l'élément linéaire avec $E A$ constant, $d N_{1} / d ξ = - \frac{1}{2}$ , et :

k_{11}^{e} = \int_{- 1}^{1} E A \cdot \frac{1}{4} \cdot \frac{2}{h _{e}} d ξ = \frac{E A}{h _{e}}, k_{12}^{e} = - \frac{E A}{h _{e}},

ce qui redonne la matrice $\frac{E A}{h _{e}} (1 - 1 - 1 1)$ que le chapitre 4 a déjà obtenue sur l'élément de référence (définition 4.2); la nouveauté est que la même démarche servira telle quelle pour tout degré et toute transformation.

Les éléments de Lagrange quadratiques et cubiques

Les fonctions de forme de degré p

Sur l'élément de référence, on veut $p + 1$ fonctions polynomiales de degré $p$ telles que la valeur de $u_{h}$ au nœud $a$ soit le coefficient $u_{a}^{e}$ lui-même. C'est la question de l'interpolation de Lagrange (Analyse numérique, chapitre 6), et elle a une réponse explicite.

Pour $p = 2$ , avec les nœuds $- 1$ , $0$ , $1$ :

N_{1} (ξ) = \frac{ξ ( ξ - 1 )}{2}, N_{2} (ξ) = 1 - ξ^{2}, N_{3} (ξ) = \frac{ξ ( ξ + 1 )}{2} . (7.5)

Pour $p = 3$ , avec les nœuds $- 1$ , $- \frac{1}{3}$ , $\frac{1}{3}$ , :

N_{1} (ξ) N_{3} (ξ) = - \frac{9}{16} (ξ + \frac{1}{3}) (ξ - \frac{1}{3}) (ξ - 1), = - \frac{27}{16} (ξ + 1) (ξ + \frac{1}{3}) (ξ - 1), N_{2} (ξ) N_{4} (ξ) = \frac{27}{16} (ξ + 1) (ξ - \frac{1}{3}) (ξ - 1), = \frac{9}{16} (ξ + 1) (ξ + \frac{1}{3}) (ξ - \frac{1}{3}) . (7.6)

On vérifie (7.5) en trois lignes: $N_{2} (0) = 1$ , $N_{2} (\pm 1) = 0$ , et de même pour les deux autres. Les coefficients de (7.6) s'obtiennent en évaluant le produit au nœud propre: pour $N_{2}$ , , dont l'inverse est .

Les fonctions intérieures ( $N_{2}$ pour $p = 2$ , $N_{2}$ et $N_{3}$ pour $p = 3$ ) s'annulent aux deux extrémités de l'élément: elles ne «débordent» jamais sur l'élément voisin. Les fonctions , elles, se recollent avec celles de l'élément voisin au nœud commun, exactement comme les deux moitiés d'une fonction chapeau. La fonction assemblée est donc sur — elle est dans , ce qui est la condition du chapitre 2 —, mais sa dérivée saute en général aux nœuds d'extrémité.

Démonstration. Le polynôme $Q = P - \sum_{a} P (ξ_{a}) N_{a}$ est de degré au plus $p$ et s'annule aux $p + 1$ nœuds distincts , puisque . Un polynôme non nul de degré au plus a au plus racines, donc . Les cas et donnent les deux identités.

Pour la dernière affirmation, soit $u (x) = α + β x$ et $u_{a} = u (x_{a}^{e})$ ses valeurs nodales. Alors, par (7.1) et la partition de l'unité,

a \sum u_{a} N_{a} (ξ) = α a \sum N_{a} (ξ) + β a \sum x_{a}^{e} N_{a} (ξ) = α + β x (ξ) = u (x (ξ)) . □

La reproduction des fonctions affines n'est pas un luxe: un élément qui ne représenterait pas exactement un déplacement d'ensemble (une constante) ou une déformation uniforme (une fonction affine) ne pourrait pas converger. C'est le sens du test de la pièce (patch test), que les chapitres 10 et 11 introduiront. Remarquez aussi ce que le théorème ne dit pas: si les nœuds ne sont pas équidistants sur l'élément réel, l'interpolant d'une fonction quadratique de $x$ n'est plus exact, car $ξ$ n'est plus une fonction affine de $x$ . Nous le retrouverons avec le jacobien.

Question 7.1

On maille $[0, 1]$ avec $n$ éléments quadratiques et l'on impose $u (0) = u (1) = 0$ . Combien le système réduit a-t-il d'inconnues?

La matrice de rigidité de l'élément quadratique

Pour l'élément cubique à nœuds équidistants, le même calcul — plus long, à faire faire par la machine en fractions exactes — donne

k^{e} = \frac{E A}{40 h} 148 - 189 54 - 13 - 189 432 - 297 54 54 - 297 432 - 189 - 13 54 - 189 148, f^{e} = \frac{q h}{8} 1331,

les poids de la règle des trois huitièmes, cette fois. L'exercice 7.3 l'utilise.

Question 7.2

Un élément quadratique a ses nœuds en $x_{1} = 2, 0$ , $x_{2} = 2, 3$ et $x_{3} = 2, 6$ (mètres). Que vaut la dérivée $d N_{2} / d x$ au point $ξ = 1/2$ de l'élément de référence, en m⁻¹?

L'intégration numérique de Gauss–Legendre

Pourquoi intégrer numériquement

Les intégrales (7.4) se sont faites à la main parce que tout y était polynomial: $E A$ constant, $q$ constant, $J$ constant. Il suffit d'une des trois complications suivantes pour que ce ne soit plus vrai, ou plus raisonnable:

un coefficient variable, $E A (x)$ ou $q (x)$ — une section qui varie, une charge de forme quelconque: l'intégrande contient alors une fonction donnée, que l'on ne sait pas toujours intégrer;
une transformation non affine: dès que le nœud milieu d'un élément quadratique n'est plus au milieu, $J (ξ)$ dépend de $ξ$ , et l'intégrande de la rigidité, qui contient $1/ J$ , devient une fraction rationnelle;
le nombre d'intégrales: un code qui offre dix types d'éléments ne peut pas contenir dix tables de matrices écrites à la main, surtout en dimension 2 et 3, où les éléments isoparamétriques du chapitre 10 ont des jacobiens variables par construction.

Tous les codes d'éléments finis remplacent donc les intégrales élémentaires par une quadrature: une somme pondérée de valeurs de l'intégrande en quelques points de l'élément de référence. Le choix de ces points est le sujet de cette section.

Sur un élément réel, la quadrature s'applique après le changement de variable (7.3):

\int_{e} G (x) d x = \int_{- 1}^{1} G (x (ξ)) J (ξ) d ξ \approx i = 1 \sum n w_{i} G (x (ξ_{i})) J (ξ_{i}) . (7.10)

Le facteur $J (ξ_{i})$ est le plus souvent oublié des trois; la question 7.3 vous fera réparer une routine qui l'oublie.

Les formules de Gauss–Legendre

Les formules de Newton–Cotes — trapèze, Simpson — fixent des points équidistants et n'optimisent que les poids. Avec $n$ points et $n$ poids, il y a $2 n$ paramètres libres, et l'on peut espérer intégrer exactement les $2 n$ monômes $1, ξ, \dots, ξ^{2 n - 1}$ . C'est ce que réalisent les formules de Gauss–Legendre, dont le chapitre 8 d'Analyse numérique a donné la construction complète. Leurs points sont les racines du polynôme de Legendre $L_{n}$ , et leurs poids sont tous strictement positifs:

$n$	points $ξ_{i}$	poids $w_{i}$	degré d'exactitude
1	$0$	$2$	1
2	$\pm 1/ 3 = \pm 0, 5773502692$

La formule à un point est la règle du point milieu. Celle à deux points fait, avec deux évaluations, ce que Simpson fait avec trois.

Démonstration pour $n = 2$ , en entier. Les points sont $ξ_{1, 2} = \mp 1/ 3$ et les poids $w_{1} = w_{2} = 1$ . La formule et l'intégrale sont toutes deux en : il suffit donc de vérifier l'égalité sur la base des polynômes de degré au plus 3.

$g = 1$ : $\int_{- 1}^{1} 1 = 2$ , et $1 + 1 = 2$ .
$g = ξ$ et : les deux intégrales sont nulles par imparité, et les deux sommes aussi, pour impair, puisque les points sont symétriques et les poids égaux.

La formule est donc exacte au degré 3. Elle ne l'est pas au degré 4: pour $g = ξ^{4}$ , $\int_{- 1}^{1} ξ^{4} = \frac{2}{5}$ tandis que . Plus précisément, le polynôme , de degré 4, s'annule aux deux points, donc la formule rend 0, alors que son intégrale est strictement positive (, non identiquement nul).

D'où viennent ces points? Si l'on impose l'exactitude sur $1, ξ, ξ^{2}, ξ^{3}$ à une formule symétrique $w g (- a) + w g (a)$ , les monômes impairs sont automatiquement satisfaits, $g = 1$ donne et donne : le système impose et . C'est l'unique formule symétrique à deux points de degré 3.

Démonstration pour $n$ quelconque. Elle repose sur deux propriétés établies au chapitre 8 d'Analyse numérique (définition 8.6 et théorème 8.10): le polynôme de Legendre $L_{n}$ , de degré $n$ , est orthogonal sur $[- 1, 1]$ à tout polynôme de degré inférieur à $n$ , et il a $n$ racines simples $ξ_{1}, \dots, ξ_{n}$ dans ; les poids sont , où est le polynôme de Lagrange de degré associé à , de sorte que la formule est exacte pour les polynômes de degré au plus . Soit alors de degré au plus . La division euclidienne par donne , avec et de degré au plus . D'une part

\int_{- 1}^{1} g = \int_{- 1}^{1} Q L_{n} + \int_{- 1}^{1} R = 0 + \int_{- 1}^{1} R,

par orthogonalité. D'autre part, $L_{n} (ξ_{i}) = 0$ entraîne $g (ξ_{i}) = R (ξ_{i})$ , donc , puisque est de degré au plus . Les deux membres sont égaux. Enfin, pour , de degré , la formule donne 0 et l'intégrale est strictement positive; le même polynôme construit sur les points d'une formule à points montre qu'aucune n'atteint le degré .

Le programme suivant vérifie le théorème monôme par monôme. Il imprime, pour $n = 1, 2, 3$ et $d = 0, \dots, 7$ , l'erreur de la formule sur $\int_{- 1}^{1} ξ^{d} d ξ$ , en écrivant 0 quand elle est au niveau de l'arrondi:

from math import sqrt

GAUSS = {
    1: ([0.0], [2.0]),
    2: ([-1 / sqrt(3), 1 / sqrt(3)], [1.0, 1.0]),
    3: ([-sqrt(3 / 5), 0.0, sqrt(3 / 5)], [5 /

1   0     0   7e-01   0   4e-01   0   3e-01   0  
2   0     0     0     0   2e-01   0   2e-01   0  
3   0     0     0     0     0     0   5e-02   0

La première erreur non nulle apparaît au degré $2 n$ : 2, 4, 6. Les degrés impairs sont toujours intégrés exactement, mais c'est la symétrie qui travaille, pas le théorème: un polynôme impair a une intégrale nulle et une somme de Gauss nulle. C'est pourquoi l'explorateur suivant prend une intégrande sans symétrie.

Explorateur 7.1 · La quadrature de Gauss et son degré d'exactitude

Choisissez le nombre $n$ de points de Gauss et le degré $d$ de l'intégrande $g (ξ) = ((1 + ξ) /2)^{d}$ . La valeur de Gauss est l'aire des rectangles: chacun a pour largeur un poids $w_{i}$ et pour hauteur $g (ξ_{i})$ . Tant que $d \leq 2 n - 1$ , cette aire est exactement celle sous la courbe; au degré $2 n$ , l'erreur apparaît d'un coup.

Nombre de points n2

Degré d de l'intégrande4

Intégrale exacte 2/(d + 1)

0,400000

Gauss à n points

0,388889

Erreur

0,0111

Degré d'exactitude 2n − 1

3 (pas exacte)

L'explorateur dessine la valeur de Gauss comme l'aire de $n$ rectangles juxtaposés, de largeurs $w_{1}, \dots, w_{n}$ (qui somment à 2) et de hauteurs $g (ξ_{i})$ ; chaque point de Gauss tombe à l'intérieur de son propre rectangle, ce qu'on vérifie numériquement pour $n \leq 5$ et qui est un résultat classique pour toutes les formules de Gauss. Trois lectures. , pour , les rectangles ne coïncident pas avec la courbe — ils la débordent par endroits et la manquent ailleurs —, et pourtant leur aire totale est exactement l'intégrale: c'est la compensation que le théorème 7.2 garantit, et qu'aucun dessin ne rend évidente. , au degré , l'erreur apparaît d'un coup: avec et , la valeur exacte est et Gauss donne , soit une erreur de . , à fixé, ajouter un point divise l'erreur par un facteur important: pour , elle vaut avec deux points et avec trois.

Question 7.3

La fonction integre_gauss(g, a, b, n) doit approcher l'intégrale de g sur [a, b] par la formule de Gauss à n points, transportée par x = (a + b)/2 + (b - a)/2 ξ. Elle oublie le jacobien de ce changement de variable. Réparez-la. Le programme doit alors afficher les approximations à deux points de l'intégrale de x³ sur [0, 3], qui est exacte (81/4), et de x⁴ sur [0, 3], qui ne l'est pas (la valeur exacte est 48,6).

gauss.py

Combien de points faut-il?

Compter les degrés

La règle pratique découle directement du théorème 7.2: on compte le degré de l'intégrande sur l'élément de référence, et l'on prend le plus petit $n$ tel que $2 n - 1$ l'atteigne. Pour un élément de degré $p$ dont la transformation est affine ( $J$ constant):

la rigidité $\int E A N_{a}^{'} N_{b}^{'}$ avec $E A$ constant a une intégrande de degré $2 (p - 1) = 2 p - 2$ ; il faut , soit points. Un point pour l'élément linéaire, deux pour le quadratique, trois pour le cubique;

Pour l'élément de poutre d'Hermite du chapitre 6, les fonctions de forme sont cubiques et l'énergie de flexion fait intervenir leurs dérivées secondes, de degré 1: l'intégrande de la rigidité est de degré 2, et deux points de Gauss suffisent. L'exemple 7.4 le vérifie sur la console du cours.

Question 7.4

Un élément cubique de Lagrange (quatre nœuds équidistants) a une rigidité $E A$ qui varie linéairement sur l'élément. Combien de points de Gauss faut-il au minimum pour intégrer exactement sa matrice de rigidité?

Ce que coûte une intégration trop pauvre

Que se passe-t-il si l'on intègre la rigidité de l'élément quadratique avec un seul point, en $ξ = 0$ ? Les dérivées y valent $- \frac{1}{2}$ , $0$ et $\frac{1}{2}$ , et la formule (7.10) avec et donne, pour ,

k_{r \overset{e}{ˊ} duite}^{e} = 2 \cdot \frac{h}{2} \cdot \frac{4}{h ^{2}} \frac{1}{4} 0 - \frac{1}{4} 000 - \frac{1}{4} 0 \frac{1}{4} = \frac{1}{h} 10 - 1 000 - 1 01 .

Le nœud milieu a disparu de la matrice: sa ligne et sa colonne sont nulles. Le déplacement $d^{e} = (0, 1, 0)$ — le nœud milieu bouge, les extrémités non — est une déformation réelle de l'élément, de dérivée $- 2 ξ \cdot \frac{2}{h}$ , dont l'énergie exacte vaut $\frac{16}{3 h} > 0$ ; mais sa dérivée s'annule précisément au seul point de Gauss, et l'intégration réduite lui attribue une énergie . C'est un (): un au sens de la définition 5.3, mais créé par la sous-intégration, et non par un mode rigide ou un mécanisme réel — une déformation que l'élément ne «voit» pas.

Le mode ne disparaît pas à l'assemblage, parce que le nœud milieu n'appartient qu'à un seul élément: aucun voisin ne vient le retenir. Sur le problème modèle avec deux éléments, la matrice réduite des trois inconnues intérieures vaut $diag (0, 4, 0)$ : elle est singulière, et la fonction resoudre s'arrête sur une division par zéro au premier pivot. Le même phénomène frappe l'élément cubique intégré à deux points: sa matrice $4 \times 4$ n'a plus que le rang 2, au lieu du rang 3 qu'exige un seul mode rigide.

Le critère qui se dégage est un comptage: le nombre de points d'intégration, multiplié par le nombre de composantes de déformation évaluées en chaque point, doit être au moins égal au nombre de degrés de liberté de l'élément diminué du nombre de ses modes rigides. Ici, une seule déformation par point; l'élément quadratique a 3 ddl et 1 mode rigide, donc il lui faut au moins 2 points; le cubique, 4 ddl et 1 mode rigide, au moins 3.

Pourquoi, alors, l'intégration réduite existe-t-elle? En dimension 1 et pour la barre, elle n'a aucun intérêt. Mais pour la poutre de Timoshenko, évoquée au chapitre 6, et pour les éléments bidimensionnels en flexion du chapitre 10, l'intégration complète rend certains éléments beaucoup trop raides — c'est le verrouillage (locking) — et sous-intégrer la partie fautive de l'énergie est l'un des remèdes classiques. Il se paie précisément par des modes parasites, dits modes en sablier (hourglass modes) en 2D, qu'il faut alors contrôler. Le chapitre 10 en traitera; retenez ici que l'intégration réduite est un choix délibéré, à justifier, jamais un moyen d'économiser des points.

Les coefficients variables: une barre conique

Les éléments distordus et le signe du jacobien

Quand le nœud milieu d'un élément quadratique n'est pas au milieu, la transformation (7.1) n'est plus affine. Avec les dérivées de l'exemple 7.1,

J (ξ) = \frac{x _{3}^{e} - x _{1}^{e}}{2} + ξ (x_{1}^{e} - 2 x_{2}^{e} + x_{3}^{e}), (7.11)

une fonction affine de $ξ$ , qui se réduit à $h /2$ quand $x_{2}^{e}$ est au milieu.

Démonstration. $J$ est affine en $ξ$ , donc il est strictement positif sur $[- 1, 1]$ si et seulement s'il l'est aux deux bornes. En $ξ = - 1$ , (7.11) donne $J (- 1) = \frac{h}{2} - x_{1}^{e} + 2 x_{2}^{e} - x_{3}^{e} = 2 x_{2}^{e} - \frac{3 x _{1}^{e} + x _{3}^{e}}{2}$ , qui est positif si et seulement si . En , , positif si et seulement si .

Si le nœud milieu sort de la moitié centrale, $J$ change de signe dans l'élément: la transformation n'est plus injective, l'élément «se replie» sur lui-même, et les intégrales (7.4), qui supposent $J > 0$ , n'ont plus de sens. Le cas limite est instructif. Avec $x_{1}^{e} = 0$ , $x_{3}^{e} = h$ et , le calcul donne , d'où : toute fonction quadratique en devient une fonction de , et sa dérivée en se comporte comme près du nœud . Un élément mal formé fabrique une singularité que le problème n'a pas. La mécanique de la rupture exploite d'ailleurs exactement ce mécanisme, à dessein: le champ de contraintes y est singulier en à la pointe d'une fissure, et un élément dont le nœud milieu est placé au quart le reproduit (l'élément «au quart», ).

Le jacobien doit rester positif: un élément mal formé ruine silencieusement le résultat

Un code d'éléments finis calcule $J$ en chaque point de Gauss et s'arrête, ou devrait s'arrêter, si $J \leq 0$ . Mais un jacobien positif et très petit passe ce contrôle et dégrade la solution sans le moindre message: l'intégrande $1/ J$ de la rigidité devient très grande près du point où $J$ est petit, la quadrature de Gauss l'intègre mal, et les dérivées — donc les contraintes — sont faussées dans tout l'élément. En dimension 1, cela n'arrive qu'à un nœud milieu mal placé; en dimension 2 et 3, un sommet déplacé, un angle trop ouvert ou une arête courbe trop prononcée suffisent, et le chapitre 10 le montrera sur des quadrilatères. Regardez la qualité du maillage avant de regarder les résultats.

L'effet sur la quadrature se mesure. Pour l'élément $x_{1}^{e} = 0$ , $x_{2}^{e} = 0, 3$ , $x_{3}^{e} = 1$ (), le jacobien vaut , de en à en . Le coefficient se calcule exactement — c'est le problème guidé 7.1 — et vaut . Gauss donne avec 2 points, avec 3, avec 4 et avec 6: l'erreur relative passe de 22 % à 0,1 %, mais nombre de points n'est exact, car l'intégrande n'est pas un polynôme. Le champ affine , lui, a toujours l'énergie exacte (théorème 7.1 et (7.4): son intégrande d'énergie vaut , un polynôme de degré 1). Le champ quadratique , en revanche, n'est plus reproduit: ses valeurs nodales donnent, en , c'est-à-dire en , la valeur au lieu de .

Question 7.5

La fonction rigidite_quadratique(x, EA, npg) doit calculer la matrice 3 × 3 d'un élément quadratique isoparamétrique de nœuds x[0], x[1], x[2], par Gauss à npg points, selon (7.4). Elle calcule bien le jacobien J et les dérivées B = dN/dx, mais elle oublie l'élément de longueur dx = J dξ dans la somme. Réparez-la. Pour l'élément [0, 1] avec son nœud milieu en 0,5 et EA = 3, elle doit afficher la matrice (7.7), c'est-à-dire (7, −8, 1; −8, 16, −8; 1, −8, 7).

rigidite_quadratique.py

Le problème modèle avec des éléments quadratiques

Le programme

Il est temps de mettre ensemble les fonctions de forme, le jacobien, la quadrature et l'assemblage. Le programme ci-dessous résout le problème modèle avec $n$ éléments de Lagrange de degré $p$ quelconque. Il réutilise la fonction resoudre(K, F) du chapitre 1, à recopier telle quelle en tête du fichier, et commence par calculer les points et les poids de Gauss pour un nombre de points quelconque, par la méthode de Newton sur la récurrence des polynômes de Legendre:

from math import pi, sin, cos


def gauss_legendre(n):
    """Points et poids de Gauss-Legendre a n points (Newton sur le polynome de Legendre)."""
    points, poids = [], []
    for i in range(n):
        x = cos(pi * (i + 0.75) / (n + 0.5))           # estimation initiale de la racine
        for _ in range(50):
            p0, p1 = 1.0, x

Viennent ensuite les fonctions de forme de Lagrange et leurs dérivées, calculées par la définition 7.2 (la dérivée d'un produit s'accumule facteur par facteur), puis l'élément et l'assemblage. L'élément $e$ (compté à partir de 0) porte les nœuds globaux $p e, p e + 1, \dots, p e + p$ ; la liste renvoyée contient les $p n + 1$ valeurs nodales, conditions aux limites comprises:

def forme(p, xi):
    """Valeurs N_a(xi) et derivees dN_a/dxi de Lagrange, noeuds -1 + 2a/p."""
    xs = [-1 + 2 * a / p for a in range(p + 1)]
    N, dN = [], []
    for a in range(p + 1):
        v, dv = 1.0, 0.0
        for b in range(p

Une seule boucle de Gauss sert ici à la rigidité et à la charge, avec $p + 4$ points: c'est plus que les $p$ points nécessaires pour la rigidité, ce qui ne change rien (une formule exacte au degré $2 p + 7$ l'est a fortiori au degré $2 p - 2$ ), et c'est assez pour que l'erreur de quadrature sur la charge, avec $f = π^{2} sin π x$ , soit invisible: en passant de 8 à 12 points, les valeurs nodales des éléments quadratiques ne bougent pas de plus de . Il reste à mesurer les erreurs. On intègre et élément par élément, avec 8 points de Gauss:

def erreurs(u_h, p, n):
    """Erreurs en energie et L2 pour u = sin(pi x), 8 points de Gauss par element."""
    h, J = 1.0 / n, 0.5 / n
    xg, wg = gauss_legendre(8)
    e_E = e_L = 0.0
    for e in range(n):
        for xi, w in zip(xg, wg):
            N, dN = forme(p, xi)
            x = e * h + (xi

p=1 n=  2 ddl=  1   9.66852e-01   1.5088e-01
p=1 n=  4 ddl=  3   4.98508e-01   3.9284e-02
p=1 n=  8 ddl=  7   2.51182e-01   9.9209e-03
p=1 n= 16 ddl= 15   1.25833e-01   2.4865e-03
p=1 n= 32 ddl= 31   6.29469e-02   6.2202e-04
p=2 n=  2 ddl=  3   1.97190e-01   1.5186e-02
p=2 n=  4 ddl=  7   5.06198e-02   1.9518e-03
p=2 n=  8 ddl= 15   1.27389e-02   2.4568e-04
p=2 n= 16 ddl= 31   3.18999e-03   3.0763e-05
p=2 n= 32 ddl= 63   7.97827e-04   3.8471e-06

Les cinq premières lignes reproduisent la table fixée du cours pour les éléments linéaires — c'est le premier contrôle du programme: un code d'ordre $p$ quelconque doit retrouver le cas $p = 1$ . Les cinq suivantes sont nouvelles.

La table des éléments quadratiques

Disposée avec les rapports entre deux maillages successifs:

$n$	ddl	$∥ u - u_{h} ∥_{E}$	rapport	$∥ u - u_{h} ∥_{L^{2}}$

Les rapports approchent 4 et 8 sans les atteindre: les ordres observés, $lo g_{2}$ des rapports, valent 1,962 puis 1,991, 1,997 et 1,999 en énergie, et 2,960 puis 2,990, 2,997 et 2,999 en norme $L^{2}$ . C'est l'ordre 2 en énergie et l'ordre 3 en $L^{2}$ que prédit l'estimation du chapitre 11 avec $p = 2$ , et c'est un comportement asymptotique: le premier rapport, 3,896, est encore sensiblement au-dessous de 4. Avec des éléments cubiques, le même programme donne des rapports 7,825, 7,956, 7,989, 7,997 en énergie et 15,653, 15,913, 15,978, 15,995 en $L^{2}$ : ordres 3 et 4. Chaque degré supplémentaire ajoute une unité aux deux ordres.

À coût égal

La comparaison honnête n'est pas «un élément quadratique contre un élément linéaire», mais à nombre d'inconnues égal. Or $n$ éléments quadratiques et $2 n$ éléments linéaires ont exactement les mêmes nœuds et les mêmes $2 n - 1$ inconnues. Avec 31 inconnues — 16 éléments quadratiques ou 32 linéaires —, l'erreur en énergie vaut $3, 19 \cdot 1 0^{- 3}$ contre $6, 29 \cdot 1 0^{- 2}$ , soit (19,7 exactement), et l'erreur contre , 20,2 fois moins. Avec 7 inconnues, le gain n'est que d'un facteur 5,0 en énergie et 5,1 en : le gain d'un ordre supérieur croît quand on raffine, puisque les deux erreurs décroissent à des vitesses différentes.

Le coût n'est pas tout à fait le même pour autant. La matrice des éléments linéaires est tridiagonale; celle des éléments quadratiques est pentadiagonale — la ligne d'un nœud d'extrémité couple ce nœud à quatre voisins, à travers les deux éléments qui le partagent —, et sa largeur de bande double. Le travail de la factorisation, proportionnel au carré de la demi-largeur de bande, est donc environ quatre fois plus grand à nombre d'inconnues égal. Pour une précision de $1 0^{- 3}$ en énergie, il faut pourtant 32 éléments quadratiques (63 inconnues), contre quelque 2 000 éléments linéaires si l'on extrapole l'ordre 1 de la table du cours ( $0, 06295 \cdot 32/ n \leq 1 0^{- 3}$ donne $n \geq 2015$ ): aucun facteur de bande ne compense cela. En dimension 2 et 3, le bilan est plus nuancé — la matrice se remplit davantage —, mais la conclusion qualitative est la même pour une solution régulière.

Pour une solution non régulière, en revanche, la hiérarchie s'effondre: près d'un angle rentrant ou d'une charge ponctuelle, $∣ u ∣_{p + 1}$ n'est plus fini, l'ordre chute pour tous les degrés, et monter en degré ne rapporte presque plus rien. Le chapitre 11 traitera ces singularités et le maillage qu'elles demandent.

Les contraintes sont moins précises que les déplacements, et le moins aux nœuds

La table montre une perte d'un ordre entre $u$ (norme $L^{2}$ , ordre 3) et sa dérivée (énergie, ordre 2). Elle ne dit pas où la dérivée est la meilleure. Un calcul ponctuel sur la même solution quadratique le dit: avec 8 éléments, l'erreur maximale sur $u^{'}$ vaut $3, 98 \cdot 1 0^{- 2}$ et elle est atteinte des éléments, tandis qu'aux deux points de Gauss de chaque élément elle ne vaut que — 66 fois moins —, et elle y décroît à l'ordre 3 (rapports 7,62, 7,90 et 7,98 quand double de 4 à 32), contre l'ordre 2 aux nœuds. C'est pourquoi les logiciels calculent les déformations et les contraintes et les extrapolent ensuite aux nœuds: lire une contrainte à un nœud, c'est la lire là où elle est le moins précise.

Question 7.6

Avec 32 éléments quadratiques, l'erreur $L^{2}$ sur le problème modèle vaut $3, 8471 \cdot 1 0^{- 6}$ . En supposant l'ordre 3 exactement atteint, quelle erreur $L^{2}$ prédisez-vous pour 64 éléments quadratiques? Donnez-la en unités de $1 0^{- 7}$ .

Question 7.7

Cette fonction assemble des éléments quadratiques sur le problème modèle, avec la matrice (7.7) et une charge intégrée par Gauss à trois points. Elle numérote mal les nœuds: elle donne à l'élément e les nœuds globaux e, e + 1, e + 2, comme si les éléments se chevauchaient. Corrigez la table de connectivité. Avec f = 1 et deux éléments, le programme doit afficher les cinq valeurs exactes x(1 − x)/2 aux nœuds 0; 0,25; 0,5; 0,75; 1.

assemblage_quadratique.py

from math import 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]))
      M[k], M[p] = M[p], M[k]
      c[k], c[p] = c[p], c[k]
      for i in range(k + 1, n):
          m = M[i][k] / M[k][k]
          for j in range(k, n):
              M[i][j] -= m * M[k][j]
          c[i] -= m * c[k]
  d = [0.0] * n
  for i in range(n - 1, -1, -1):
      s = c[i] - sum(M[i][j] * d[j] for j in range(i + 1, n))
      d[i] = s / M[i][i]
  return d

XI = [-sqrt(3 / 5), 0.0, sqrt(3 / 5)]
W = [5 / 9, 8 / 9, 5 / 9]

def elements_quadratiques(f, n):
  """-u'' = f, u(0) = u(1) = 0, n elements quadratiques; u[i] contient u(i h / 2)."""
  h = 1.0 / n
  nn = 2 * n + 1
  K = [[0.0] * nn for _ in range(nn)]
  F = [0.0] * nn
  ke = [[7 / (3 * h), -8 / (3 * h), 1 / (3 * h)],
        [-8 / (3 * h), 16 / (3 * h), -8 / (3 * h)],
        [1 / (3 * h), -8 / (3 * h), 7 / (3 * h)]]
  for e in range(n):
      noeuds = [e, e + 1, e + 2]
      for xi, w in zip(XI, W):
          x = e * h + (xi + 1) * h / 2
          N = [xi * (xi - 1) / 2, 1 - xi * xi, xi * (xi + 1) / 2]
          for a in range(3):
              F[noeuds[a]] += w * f(x) * N[a] * h / 2
      for a in range(3):
          for b in range(3):
              K[noeuds[a]][noeuds[b]] += ke[a][b]
  Kr = [ligne[1:nn - 1] for ligne in K[1:nn - 1]]
  return [0.0] + resoudre(Kr, F[1:nn - 1]) + [0.0]

u = elements_quadratiques(lambda x: 1.0, 2)
print(["%.6f" % v for v in u])

La quadrature de la charge, encore

L'exemple 7.2 a montré qu'avec un seul élément, intégrer la charge par deux points de Gauss au lieu de la calculer exactement coûtait 24 % au nœud milieu. Sur un maillage, la question est de savoir si l'ordre survit. Le même programme, avec la charge intégrée par un nombre de points imposé, répond:

points pour la charge	$∥ u - u_{h} ∥_{E}$ , $n = 32$	rapports en énergie	rapports en $L^{2}$
8 (référence)	$7, 9783 \cdot 1 0^{- 4}$

Deux points suffisent: l'ordre est intact et l'erreur, au quatrième chiffre près, la même. Un seul point — la règle du point milieu, qui ne charge que le nœud milieu de chaque élément — fait perdre un ordre: les rapports tendent vers 2 et 4, l'erreur en énergie est 40 fois plus grande avec 32 éléments, et l'élément quadratique ne vaut alors pas mieux qu'un élément linéaire. La règle générale, que l'analyse des «crimes variationnels» rend précise (Strang et Fix), est qu'une quadrature exacte au degré $2 p - 2$ préserve l'ordre de convergence en énergie: ici $2 p - 2 = 2$ , que deux points atteignent et qu'un seul n'atteint pas.

La condensation statique du nœud intérieur

Le nœud milieu d'un élément quadratique n'appartient qu'à cet élément. Son équation, dans le système global, ne fait intervenir que les trois nœuds de l'élément. On peut donc la résoudre avant l'assemblage, élément par élément, et n'assembler qu'un système portant sur les nœuds d'extrémité.

Partageons les ddl de l'élément en ddl de bord $d_{b} = (u_{1}, u_{3})$ et ddl intérieur $u_{i} = u_{2}$ :

(k_{bb} k_{ib} k_{bi} k_{ii}) (d_{b} u_{i}) = (f_{b} f_{i}) .

La seconde ligne donne $u_{i} = (f_{i} - k_{ib} d_{b}) / k_{ii}$ ; en la reportant dans la première,

k_{c} (k_{bb} - \frac{k _{bi} k _{ib}}{k _{ii}}) d_{b} = f_{c} f_{b} - \frac{k _{bi}}{k _{ii}} f_{i} . (7.12)

C'est l'élimination de Gauss de l'inconnue intérieure, faite localement. La matrice condensée $k_{c}$ est $2 \times 2$ , symétrique, et on l'assemble comme celle d'un élément linéaire; après la résolution, on récupère $u_{i}$ par la première formule, élément par élément.

Pour l'élément quadratique (7.7) avec $E A$ constant, $k_{bi} = - \frac{8 E A}{3 h} (1, 1)^{T}$ et $k_{ii} = \frac{16 E A}{3 h}$ , d'où

k_{c} = \frac{E A}{3 h} [(7117) - \frac{64}{16} (1111)] = \frac{E A}{h} (1 - 1 - 1 1) .

La matrice condensée de l'élément quadratique est celle de l'élément linéaire. Et pour une charge constante, (7.8) donne $f_{c} = \frac{q h}{6} (1, 1)^{T} + \frac{1}{2} \cdot \frac{4 q h}{6} (1, 1)^{T} = \frac{q h}{2} (1, 1)^{T}$ : la charge de l'élément linéaire. Le système condensé des éléments quadratiques est donc au système des éléments linéaires sur les mêmes nœuds d'extrémité. L'exercice 7.5 montre que cela vaut pour toute charge : c'est l'explication de l'exactitude des nœuds d'extrémité observée dans la question 7.7. Les éléments linéaires y sont exacts (chapitre 4); le système condensé des quadratiques est le même; donc ses solutions le sont aussi. Ce que les quadratiques ajoutent, c'est le nœud milieu et la courbure entre les nœuds — c'est-à-dire tout ce qui fait baisser l'erreur en énergie.

Question 7.8

Remettez dans l'ordre les étapes du calcul de la matrice de rigidité d'un élément isoparamétrique par quadrature de Gauss.

Glissez les éléments pour les mettre dans le bon ordre

Calculer les dérivées

d N_{a} / d x = J^{- 1} d N_{a} / d ξ

et le coefficient

E A

au point

x (ξ_{i})

Calculer le jacobien

J (ξ_{i}) = \sum_{a} x_{a}^{e} d N_{a} / d ξ

et vérifier qu'il est positif

Après le dernier point, assembler

k^{e}

dans

K

par la table de connectivité

En chaque point de Gauss, évaluer les fonctions de forme

N_{a} (ξ_{i})

et leurs dérivées

d N_{a} / d ξ

Ajouter

w_{i} E A (d N_{a} / d x) (d N_{b} / d x) J (ξ_{i})

à chaque coefficient

k_{ab}^{e}

Choisir le nombre de points de Gauss d'après le degré de l'intégrande

La console revisitée

Le chapitre 6 a calculé à la main la matrice de l'élément de poutre et ses charges nodales cohérentes. Avec la quadrature, ce calcul devient mécanique, et l'on peut contrôler le nombre de points.

Un peu d'histoire

Les formules de quadrature qui portent le nom de Gauss datent du début du dix-neuvième siècle; la forme moderne, avec les racines des polynômes de Legendre, est due aux travaux qui ont suivi, notamment de Jacobi. Elles sont restées un outil de calculateur jusqu'à ce que la méthode des éléments finis en fasse un rouage central dans les années 1960, quand les ingénieurs ont voulu des éléments plus riches que le triangle linéaire de Turner, Clough, Martin et Topp.

Les éléments isoparamétriques apparaissent dans cette décennie. L'idée de décrire la géométrie d'un élément avec ses propres fonctions de forme, et d'en calculer les matrices par intégration numérique sur un élément de référence, est associée aux noms de Bruce Irons, à Rolls-Royce puis à Swansea, et d'Olgierd Zienkiewicz; l'article d'Ergatoudis, Irons et Zienkiewicz (1968) sur les quadrilatères isoparamétriques courbes en est une référence classique. Le chapitre 10 en exploitera toute la portée en deux dimensions. Le constat qu'une intégration «trop exacte» rend certains éléments trop raides, et que la sous-intégration peut les assouplir, date du début des années 1970; il a ouvert un long chapitre de recherche sur le verrouillage et les modes parasites, dont les éléments des logiciels actuels sont l'aboutissement.

Synthèse

Tout élément se construit une fois, sur l'élément de référence $[- 1, 1]$ , et se transporte par la transformation isoparamétrique $x (ξ) = \sum_{a} x_{a}^{e} N_{a} (ξ)$ . Le donne les dérivées, , et l'élément de longueur, ; pour une transformation affine, .

Série d'exercices du chapitre 7Exercice 1 sur 5

Question 7.9

Quel est le degré d'exactitude de la formule de Gauss–Legendre à 4 points?

Exercices

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

Exercice 7.1 · Les fonctions de forme quadratiques et leur matrice

On considère l'élément quadratique de référence, de nœuds $ξ_{1} = - 1$ , $ξ_{2} = 0$ , $ξ_{3} = 1$ .

Exercice 7.2 · La formule de Gauss à trois points

On cherche une formule symétrique à trois points, $w g (- a) + w_{0} g (0) + w g (a)$ , exacte pour $1$ , $ξ^{2}$ et . Écrire le système et le résoudre.

Exercice 7.3 · Un élément cubique exact

On résout $- u^{''} = 6 x$ sur $(0, 1)$ , $u (0) = u (1) = 0$ , dont la solution exacte est $u (x) = x - x^{3}$ , avec élément cubique de nœuds , , , .

Exercice 7.4 · Intégration réduite et modes parasites

Calculer la matrice de rigidité de l'élément quadratique ( $E A = 1$ , longueur $h$ , nœud milieu au milieu) intégrée avec un point de Gauss. Quel est son rang? Identifier un vecteur de son noyau qui n'est pas un mode rigide.
On assemble deux tels éléments sur $[0, 1]$ avec $u (0) = u (1) = 0$ . Écrire la matrice réduite des trois inconnues et conclure.

Exercice 7.5 · Condensation statique et exactitude des nœuds d'extrémité

On considère le problème $- u^{''} = f$ sur $(0, 1)$ , $u (0) = u (1) = 0$ , avec $f$ continue quelconque, discrétisé par éléments quadratiques de nœuds milieux au milieu, la charge étant intégrée exactement.

Références

Hughes, T. J. R., The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover, chap. 3 et 4 (élément de référence, éléments isoparamétriques, intégration numérique, intégration réduite et modes parasites).
Zienkiewicz, O. C., Taylor, R. L. et Zhu, J. Z., The Finite Element Method: Its Basis and Fundamentals, Butterworth-Heinemann (fonctions de forme de Lagrange, transformation isoparamétrique, nombre de points d'intégration).
Bathe, K.-J., Finite Element Procedures, Prentice Hall / K.J. Bathe (éléments isoparamétriques, intégration de Gauss, condensation statique).
Cook, R. D., Malkus, D. S., Plesha, M. E. et Witt, R. J., Concepts and Applications of Finite Element Analysis, Wiley (choix du nombre de points, modes à énergie nulle, qualité des éléments, contraintes aux points de Gauss).
Strang, G. et Fix, G., An Analysis of the Finite Element Method, Prentice Hall / Wellesley-Cambridge (effet de la quadrature sur l'ordre de convergence, «crimes variationnels»).
Dhatt, G., Touzot, G. et Lefrançois, E., Méthode des éléments finis, Hermès / Lavoisier (en français; éléments de référence, intégration numérique et programmation).

2	0,96685	—	$1, 509 \cdot 1 0^{- 1}$	—
4	0,49851	1,939	$3, 928 \cdot 1 0^{- 2}$	3,841
8	0,25118	1,985	$9, 921 \cdot 1 0^{- 3}$	3,960
16	0,12583	1,996	$2, 487 \cdot 1 0^{- 3}$	3,990
32	0,06295	1,999	$6, 220 \cdot 1 0^{- 4}$	3,997

éléments	ddl	$u (L)$ [mm]	erreur relative
1 linéaire	1	0,952381	$- 3, 82$ %
2 linéaires	2	0,979592	$- 1, 07$ %
4 linéaires	4	0,987457	$- 0, 278$ %
1 quadratique	2	0,989011	$- 0, 121$ %
2 quadratiques	4	0,990109	$- 0, 0102$ %
1 cubique	3	0,990174	$- 0, 0037$ %

2	3	$1, 9719 \cdot 1 0^{- 1}$	—	$1, 5186 \cdot 1 0^{- 2}$	—
4	7	$5, 0620 \cdot 1 0^{- 2}$	3,896	$1, 9518 \cdot 1 0^{- 3}$	7,780
8	15	$1, 2739 \cdot 1 0^{- 2}$	3,974	$2, 4568 \cdot 1 0^{- 4}$	7,945
16	31	$3, 1900 \cdot 1 0^{- 3}$	3,993	$3, 0763 \cdot 1 0^{- 5}$	7,986
32	63	$7, 9783 \cdot 1 0^{- 4}$	3,998	$3, 8471 \cdot 1 0^{- 6}$	7,997

Objectifs du chapitre

Pourquoi monter en degré

L'élément de référence et la transformation isoparamétrique

Un seul élément pour tous

Dériver et intégrer sur l'élément de référence

Les éléments de Lagrange quadratiques et cubiques

Les fonctions de forme de degré p

La matrice de rigidité de l'élément quadratique

L'intégration numérique de Gauss–Legendre

Pourquoi intégrer numériquement

Les formules de Gauss–Legendre

Combien de points faut-il?

Compter les degrés

Ce que coûte une intégration trop pauvre

Les coefficients variables: une barre conique

Les éléments distordus et le signe du jacobien

Le problème modèle avec des éléments quadratiques

Le programme

La table des éléments quadratiques

À coût égal

La quadrature de la charge, encore

La condensation statique du nœud intérieur

La console revisitée

Un peu d'histoire

Synthèse

Exercices

Références

Le jacobien aux extrémités

La transformation n'est plus affine

Gauss à deux points

La valeur exacte et l'erreur de quadrature