• Nous introduisons alors le tableau XY des coordonnées des sommets du maillage :

$$
(37) \quad X = XY \ (II, 1) \quad Y = XY \ (II, 2)
$$

expriment que les coordonnées X,Y du sommet II ont pour adresses (II, 1) et (II, 2) dans le
tableau XY. Notons que pour éviter de dupliquer l'information, on fait appel aux numéros
globaux des sommets, mais on peut aussi introduire les coordonnées X(K,I) et Y(K,I) du
noeud de numéro local I dans l'élément K, à travers la double dépendance suivante :

$$
(38) \quad X(K,I) = XY \ (\text{NELT}(K,I), 1)
$$

$$
(39) \quad Y(K,I) = XY \ (\text{NELT}(K,I), 2).
$$

L'explicitation, coordonnées par coordonnées, de la relation (36) prend alors la
forme suivante : un point M = (X,Y) de K est paramétré par un point M = (x, y) de K, grâce
aux coordonnées baricentriques (31) (32) (33) et aux relations :

$$
(40) \quad X = \sum_{j=1}^{3} \hat{\lambda}_{j}(\hat{M}) X(K, j)
$$

$$
(41) \quad Y = \sum_{j=1}^{3} \hat{\lambda}_j(\hat{M}) Y(K, j).
$$

On le voit, tous les calculs avec des nombres réels se font "élément par élément"
dans la méthode des éléments finis. C'est cette remarque qui motive la notion de **fonction de base locale** puis de **matrice élémentaire**.

• On se donne un élément K du maillage τ ainsi qu’un numéro local de sommet (ou de degré de liberté puisqu’on s’intéresse à l’élément P1 où les noeuds sont biunivoquement associés aux sommets du maillage) I. On appelle **fonction de base locale** et on note $\tilde{\phi}_{I,K}$ la restriction à l’élément K de la fonction de base globale $\phi_{II}$, où $II = \text{NELT}(K,I)$ et le numéro global du degré de liberté associé. Les relations de définition (9) et (10) sont inchangées ; nous ne les réécrirons pas.

L'expression algébrique de $\tilde{\phi}_{I,K}$, (fonction de base locale) est très simple ; nous avons la :
**Proposition**

La fonction de base locale $\tilde{\phi}_{I,K}$, est identiquement égale à la $j^0$ coordonnée barycentrique $\lambda_j$ relativement à l'élément K :

$$
(42) \quad \tilde{\phi}_{j,K}(M) = \lambda_j(M) \qquad \forall M \in K.
$$

Démonstration

Un point $M \in K$ admet un unique triplet de coordonnées barycentriques $\lambda_1, \lambda_2, \lambda_3$
qui par définition sont les poids qu'il faut donner aux sommets $A_1, A_2, A_3$ pour que le
barycentre pondéré soit exactement égal à M. Par définition, la somme de ces poids vaut 1.
Nous avons donc :

$$
(43) \quad \sum_{j=1}^{3} X(K,j) \lambda_j(M) = X
$$

$$
(44) \quad \sum_{j=1}^{3} Y(K,j) \lambda_j(M) = Y
$$

$$
(45) \quad \sum_{j=1}^{3} \lambda_j (M) = 1
$$

La fonction $M \mapsto \lambda_j (M)$ est clairement fonction linéaire du triplet $(X,Y,1)$,
c'est-à-dire fonction affine de $(X,Y)$. De plus, prenant dans le second membre la valeur
particulière $X = X(K,i)$, $Y = Y(K,i)$ (coordonnées au sommet $A_i$), il est évident qu'alors $\lambda_j$ (j non coordonnée barycentrique) vaut 1 pour $j = i$ et 0 si $j \neq i$ (le vérifier sur les relations (43)-(45) si vous n'êtes pas convaincus !):

$$
(46) \quad \lambda_j (A_i) = \delta_{ij} \dots
$$

La fonction $\lambda_j$ est donc affine, vaut 1 pour le sommet de numéro local $j$ et 0 pour les deux autres sommets ; elle est donc identiquement égale à la restriction dans l'élément $K$, de la fonction de base associée à ce sommet. C'est exactement ce qu'exprime la relation (42) !

L'expression algébrique de $\tilde{\phi}_{j,K}$ s'obtient alors facilement en fonction de $X,Y$, par
résolution du système (43) (45) de trois équations à trois inconnues. Nous explicitons ce
calcul élémentaire, en adoptant une notation simplifiée le temps du calcul algébrique :

$$
(47) \quad x_j \equiv X(K,j), \quad x \equiv X, \quad y \equiv Y.
$$
On a :

$$
(48) \quad x_1 \lambda_1 + x_2 \lambda_2 + x_3 \lambda_3 = x
$$

$$
(49) \quad y_1 \lambda_1 + y_2 \lambda_2 + y_3 \lambda_3 = y
$$

$$
(50) \quad \lambda_1 + \lambda_2 + \lambda_3 = 1
$$

on retranche $x_3$ de (48), en l'écrivant $x_3 (\lambda_1 + \lambda_2 + \lambda_3)$ et on procède de façon analogue pour $y_3$ avec la relation (49). Il vient :

$$
(51) \quad (x_1-x_3)\lambda_1 + (x_2-x_3)\lambda_2 = x-x_3
$$

$$
(52) \quad (y_1 - y_3) \lambda_1 + (y_2 - y_3) \lambda_2 = y - y_3
$$

et, en remarquant que :

$$
(53) \quad \left\{
\begin{aligned}
&(y_2-y_3)(x_1-x_3) - (x_2-x_3)(y_1-y_3) = \\
&= (y_2-y_1)(x_1-x_3) - (x_2-x_1)(y_1-y_3) \\
&= (x_2-x_1)(y_3-y_1) - (y_2-y_1)(x_3-x_1)
\end{aligned}
\right.
$$

on trouve l'expression suivante pour $\lambda_1$ :

$$
(54) \quad \lambda_1 (x,y) = \frac{(y_2-y_3)(x-x_3) - (x_2-x_3)(y-y_3)}{(x_2-x_1)(y_3-y_1) - (y_2-y_1)(x_3-x_1)}
$$

Les autres coordonnées barycentriques se déduisent de (54) par permutation
circulaire des indices 1, 2, 3. On remarque que le dénominateur est le double de la surface du
triangle K, puisque égal au produit mixte suivant :

$$
(55) \quad S_K = \frac{1}{2} \left( \overrightarrow{A_1A_2} \times \overrightarrow{A_1A_3}, k \right)
$$

ceci dans l'hypothèse où l'orientation des sommets $A_1, A_2, A_3$ est "directe", ce qui revient
à dire que le déterminant de $F_K$ est positif. Cette remarque constitue en fait une règle de
numérotation locale. Il n'y a donc pas, même avec cette convention, une seule numérotation
locale possible ; toute permutation paire de 1,2,3 permet de retrouver une nouvelle
numérotation locale cohérente des noeuds de l'élément K.
• La notion de matrice élémentaire est fondée, comme dans le cas monodimensionnel, sur la propriété suivante :

**Proposition :**

L'intégrale $\int_K \nabla\varphi_{II} \cdot \nabla\varphi_{JJ} dx$ ne peut être non nulle que si les sommets II et JJ sont tous deux des sommets appartenant au bord de l'élément K.

La preuve donnée pour le cas monodimensionnel est inchangée. La définition des matrices élémentaires AELT (K,I,J) est encore obtenue par la relation (11), mais I et J varient cette fois entre 1 et 3. Le calcul de cette matrice élémentaire est un exercice laissé au lecteur (compte tenu de l'explicitation de la relation (54) donnant les fonctions de base $\tilde{\phi}_{j,K}$). Nous pouvons de plus en donner un résultat sous forme géométrique.

**Proposition :**

Avec les notations introduites à la figure ci-contre,
on a :

$$ (56) \quad \int_K \nabla \tilde{\phi}_{1,K} \cdot \nabla \tilde{\phi}_{2,K} \, dx = - \frac{1}{2 \operatorname{tg} \alpha_3} $$

$$ (57) \quad \int_K |\nabla \tilde{\phi}_{1,K}|^2 dx = \frac{\sin \alpha_1}{2 \sin \alpha_2 \sin \alpha_3} $$

<figure><img src="image_55.png" /></figure>

**Proposition**

Les matrices élémentaires sont assemblées pour former la matrice globale A(II,JJ).
L'algorithme (16) qui consiste à "éclater" la matrice élémentaire AELT (K,•,•) au sein de la
matrice assemblée est inchangé. De même, on peut définir un second membre élémentaire
BELT (K,I) à l'aide de la relation (19). L'assemblage du second membre s'obtient encore à
l'aide de la relation (20).
### 4) Condition de Dirichlet non homogène

* Nous avons jusqu'ici supposé que la condition de Dirichlet à prendre en compte est homogène : pour les relations (2) et (22), la donnée de $u(\cdot)$ sur le bord de $\Omega$ est la fonction nulle. Nous supposons dans ce paragraphe que ce n'est plus le cas, c'est-à-dire que $u$ est donnée non nulle sur une partie $\Gamma_1$ du bord.

$$ (58) \qquad u = \bar{u}_0 \quad \text{sur } \Gamma_1 $$

Sur la partie $\Gamma_2$ complémentaire de $\Gamma_1$, nous adoptons une condition de Neumann.

$$ (59) \qquad \frac{\partial u}{\partial n} = g \quad \text{sur } \Gamma_2 $$

et nous continuons à résoudre un problème de Poisson dans le domaine d'étude.

$$ (60) \qquad -\Delta u = f \quad \text{dans } \Omega $$

Les degrés de liberté du problème discrétisé sont les sommets $A_j$ du maillage $\tau$. On numérote en premier les sommets qui n'appartiennent pas à $\Gamma_1$.

$$ (61) \qquad \text{sommets n'appartenant pas à } \Gamma_1 : A_{1}, A_{2}, \dots, A_{N} $$

<figure><img src="image_56.png" /></figure>

Choix d'une numérotation pour les noeuds ; on numérote en dernier les sommets de $\Gamma_1$

puis ensuite ceux qui appartiennent à $\Gamma_1$.

$$ (62) \qquad \text{sommets de } \tau \text{ sur } \Gamma_1 : A_{N+1}, \dots, A_{N+M} $$

On dispose de ce fait de deux espaces de dimension finie $V_h$ et $W_h$ : les fonctions de ces deux espaces sont continues sur $\bar{\Omega}$, affines dans chaque triangle $K$ ; les fonctions de $V_h$ sont
supposées nulles sur la partie $\Gamma_1$ du bord, alors que nous ne faisons pas cette hypothèse pour les fonctions de $W_h$. Avec les notations usuelles pour les fonctions de base, nous avons donc :

$$
(63) \quad V_h = \operatorname{Vect} \langle \phi_1, \dots, \phi_N \rangle
$$

$$
(64) \quad W_h = \operatorname{Vect} \langle \phi_1, \dots, \phi_N, \phi_{N+1}, \dots, \phi_{N+M} \rangle
$$

et bien entendu, les dimensions de ces deux espaces s'en déduisent :

$$
(65) \quad \dim V_h = N, \quad \dim W_h = N+M.
$$

• Ces préliminaires étant posés, nous allons développer deux approches pour la prise en compte de la condition de Dirichlet non homogène (58) : l'une théorique qui consiste à se ramener au cas où $\bar{u}_o$ est nulle grâce à un relèvement de cette condition non homogène, l'autre pratique où l'on "élimine les noeuds bloqués" qui est mise en oeuvre au sein des algorithmes de calcul.'

Le relèvement dans l'espace W_h de la condition de Dirichlet non homogène u̅_o
consiste à chercher une fonction ũ_o telle que

$$
(66) \quad \tilde{u}_o \in W_h
$$

$$
(67) \quad \tilde{u}_o (A_k) = \bar{u}_o (A_k) \quad \forall k = N+1, \dots, N+M:
$$

la fonction Ű₀ est définie par ses valeurs aux sommets du maillage, est donnée égale à Ű₀
pour les sommets de Γ₁ (de numéros N+1 à N+M) et n'est pas déterminée a priori pour les
autres sommets. Nous pouvons par exemple choisir :

$$
(68) \quad \tilde{u}_o (A_j) = 0 \quad \forall j = 1, \dots, N
$$

mais tout autre choix est possible aussi. Nous faisons le choix (68) car il simplifie un peu les calculs. La fonction $\tilde{u}_o$ peut donc se calculer sans difficulté compte tenu des conditions (66) (67) (68) et la condition de Dirichlet (58) peut se réécrire de façon équivalente :

$$
(69) \quad u = \tilde{u}_o + z, \quad z \in V_h
$$

ce qui ne fait qu'exprimer que u s'obtient à partir de Ű₀ grâce à une fonction z nulle sur Γ₁. La formulation variationnelle s'obtient en multipliant la relation (60) par une fonction test v appartenant à V_h et en intégrant par parties, ce qui permet de prendre en compte la condition de Neumann (59) Il vient :
$$
(70) \quad \int_{\Omega} \nabla u \cdot \nabla v \, dx = \int_{\Omega} fv \, dx + \int_{\Gamma_2} gv \, d\gamma, \quad \forall v \in V_h
$$

En prenant en compte la relation (69), on obtient une formulation variationnelle
discrète pour la fonction z.

$$
(71) \quad z \in V_h
$$

$$
(72) \quad \int_{\Omega} \nabla z \cdot \nabla v \, dx = \int_{\Omega} fv \, dx + \int_{\Gamma_2} gv \, d\gamma - \int_{\Omega} \nabla \tilde{u}_o \cdot \nabla v \, dx, \quad \forall v \in V_h
$$

Le système linéaire décrit par (71) (72) prend la forme classique.

$$
(73) \quad A_1 Z = F + G - B \bar{u}_0
$$

avec des notations naturelles rappelées ci-dessous :

$$
(74) \quad (A_1)_{ij} = \int_{\Omega} \nabla \varphi_i \cdot \nabla \varphi_j \, dx \qquad 1 \le i, j \le N
$$

$$
(75) \quad F_i = \int_{\Omega} f \varphi_i \, dx
$$

$$
(76) \quad G_i = \int_{\Gamma_2} g \varphi_i \, d\gamma
$$

$$
(77) \quad B_{ik} = \int_{\Omega} \nabla \varphi_i \cdot \nabla \varphi_k \, dx, \quad \begin{cases} 1 \le i \le N \\ N+1 \le k \le N+M \end{cases}
$$

• Tout est bien établi du point de vue théorique ; du point de vue pratique, le calcul explicite de la fonction ũ₀ n'est jamais réalisé. On adopte une autre démarche, fondée sur l'écriture de la relation (70) non pas pour v ∈ V_h, mais pour w ∈ W_h. L'écriture des termes de bord est un peu différente et nous avons en toute rigueur

$$
(78) \quad \int_{\Omega} \nabla u \cdot \nabla w \, dx = \int_{\Omega} fw \, dx + \int_{\Gamma_2} gw \, d\gamma + \int_{\Gamma_1} \frac{\partial u}{\partial n} w \, d\gamma, \quad \forall w \in W_h
$$

expression qui est identique à (70) si nous particularisons le choix (70). La matrice de ce
système (78), qui suppose u ∈ W_h mais oublie (pour un moment) la relation (69), est donnée
par :
$$
(79) \quad A_{ij} = \int_{\Omega} \nabla \varphi_i \cdot \nabla \varphi_j \, dx \quad 1 \le i, j \le N+M
$$

et l'on notera que la seule différence entre (79) et (74) est la borne supérieure de variation du
double indice (i, j). De façon plus précise, la matrice A se décompose en blocs sous la forme :

$$
(80) \quad A = \begin{pmatrix} A_1 & B_1 \\ t_{B_1} & A_2 \end{pmatrix}.
$$

Quand on explicite la relation (78), en recherchant l'inconnue u sous la forme :

$$
(81) \quad u = z + \sum_{k=N+1}^{N+M} \lambda_k \varphi_k
$$

(sans savoir pour l'instant que $\lambda_k = \bar{u}_o$ ($A_k$) !), on obtient un système qui s'écrit par blocs :

$$
(82) \quad \begin{pmatrix} A_1 & B_1 \\ t_{B_1} & A_2 \end{pmatrix} \begin{pmatrix} Z \\ \Lambda \end{pmatrix} = \begin{pmatrix} F+G \\ * \end{pmatrix}
$$

On reconnait la relation (70) pour les N premières lignes de (82) alors qu'on n'a pas explicité
le second membre des M dernières.

• **L'élimination des noeuds bloqués** consiste à transformer le système (82) en écrivant brutalement la relation (69) à la place de la seconde ligne, c'est-à-dire ici

$$
(83) \quad \Lambda = \bar{u}_0,
$$

en éliminant d'une part les matrices blocs $^tB_1$ et $A_2$ et en les remplaçant par 0 et la matrice $\mathrm{Id}$ respectivement, en modifiant le second membre pour "forcer" la relation (83) d'autre part. Il vient :

$$
(84) \quad \begin{pmatrix} A_1 & B_1 \\ 0 & \mathrm{Id} \end{pmatrix} \begin{pmatrix} Z \\ \Lambda \end{pmatrix} = \begin{pmatrix} F+G \\ \bar{u}_0 \end{pmatrix}
$$

qui a l'inconvénient de faire apparaître une matrice non symétrique. Ce défaut peut être
corrigé en éliminant $\Lambda$ de la première équation et nous déduisons :
$$ (85) \quad \begin{pmatrix} A_1 & 0 \\ 0 & \mathrm{Id} \end{pmatrix} \begin{pmatrix} Z \\ \Lambda \end{pmatrix} = \begin{pmatrix} F+G - B_1 \bar{u}_o \\ \bar{u}_o \end{pmatrix} $$

qui est parfaitement équivalent à la relation (73).

L'intérêt de cette dernière approche est qu'elle est généralisable facilement au cas où l'on n'a pas adopté de numérotation particulière pour les degrés de liberté sur $\Gamma_1$ (cf. relations (61) et (62)). On assemble la matrice A et le second membre F+G grâce aux algorithmes (16) et (20). Puis on élimine les noeuds bloqués, opération qui consiste à passer de (82) à (85). En pratique, on a donc l'algorithme suivant :

boucle sur les noeuds : $1 \le i \le N+M$

i est-il un noeud sur $\Gamma_1$ ?

non → passer à i+1

oui

boucle sur les noeuds : $1 \le j \le N+M$

j est-il un noeud sur $\Gamma_1$ ?

non :

$B_j \leftarrow B_j - A_{ji} \bar{u}_o (A_i)$

$A_{ji} \leftarrow 0$

oui → passer à j+1

boucle sur les noeuds : $1 \le j \le N+M$

j = i ?

non :

$A_{ij} \leftarrow 0$

oui :

$A_{ii} \leftarrow 1$
Notons enfinque l'adressage matriciel est rarement à double indice comme écrit à l'algorithme (86), ce qui en pratique en rend la lecture plus délicate ...

## 5) Étude d'un exemple, dit "Hadhri"

On considère le domaine $\Omega = ]0,3[ \times ]0,3[$, sa frontière $\partial\Omega$ est décomposée en $\Gamma_1 = \{0\} \times ]0,3[ \cup ]0,3[ \times \{3\}$, et $\Gamma_2 = ]0,3[ \times \{0\} \cup \{3\} \times ]0,3[$ (voir figure ci-dessous). On veut résoudre le problème.

$$ (87) \quad \left\{ \begin{array}{ll} -\Delta u = f & \Omega \\ u = 0 & \Gamma_1 \\ \displaystyle\frac{\partial u}{\partial n} = g & \Gamma_2 \end{array} \right. $$

de façon approchée, à l'aide de la méthode des éléments finis, avec le maillage (de 18 éléments triangulaires) proposé à la figure ci-dessous. Les sommets de numéros 1 à 9 engendrent l'espace $V_h$ alors que les sommets 10 à 16 sont situés sur $\Gamma_1$.

<figure><img src="image_57.png" /></figure>

maillage en triangles de :

$$ \Omega = ]0,3[ \times ]0,3[ $$
* Compte tenu de la simplicité du maillage, le calcul des matrices élémentaires se réduit aux deux cas de figure suivants :

$$
\text{1) y} \quad \phi_1 (x,y) = 1-y \\
\text{K} \quad \phi_2 (x,y) = x \\
\text{x} \quad \phi_3 (x,y) = y-x
$$
$$
\text{matrice élémentaire } \text{AELT} = \begin{pmatrix} \frac{1}{2} & 0 & -\frac{1}{2} \\ 0 & \frac{1}{2} & -\frac{1}{2} \\ -\frac{1}{2} & -\frac{1}{2} & 1 \end{pmatrix}
$$

$$
\text{2) y} \quad \phi_1 (x,y) = 1-x \\
\text{K} \quad \phi_2 (x,y) = x-y \\
\text{x} \quad \phi_3 (x,y) = y
$$
$$
\text{matrice élémentaire } \text{AELT} = \begin{pmatrix} \frac{1}{2} & -\frac{1}{2} & 0 \\ -\frac{1}{2} & 1 & -\frac{1}{2} \\ 0 & -\frac{1}{2} & \frac{1}{2} \end{pmatrix}
$$

Il est utile aussi de calculer la matrice élémentaire de masse, ie :

$$
(88) \quad \text{MELT} (K, I, J) = \int_K \bar{\phi}_{I,K} \bar{\phi}_{J,K} \, dx.
$$

Dans les deux cas précédents, on obtient :

$$
(89) \quad \text{MELT} = \frac{1}{24} \begin{pmatrix} 2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{pmatrix} .
$$
* Si on interpole les fonctions f et g dans l'espace $W_h$ (avec les notations précédentes), on doit calculer les seconds membres :

$$ (90) \quad (\mathrm{MF})_i = \sum_{j=1}^{16} \left( \int_{\Omega} \varphi_i \varphi_j \, dx \right) f_j, \quad 1 \le i \le 9 $$

et

$$ (91) \quad (\Gamma G)_i = \sum_{j=1}^{16} \left( \int_{\Gamma_2} \varphi_i \varphi_j d\gamma \right) g_j, \quad 1 \le i \le 9 $$

Le calcul (89) de la matrice élémentaire et la relation (90) conduisent à une matrice de masse M d'ordre 9 x 16 (9 lignes et 16 colonnes) dont nous donnons ci-dessous les éléments non nuls.

$$ (92) \quad M = \frac{1}{24} \begin{bmatrix} 6 & 1 & 0 & 2 & 2 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 6 & 1 & 0 & 2 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 2 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 2 & 0 & 0 & 12 & 2 & 0 & 2 & 2 & 0 & 2 & 2 & 0 & 0 & 0 & 0 & 0 \\ 2 & 2 & 0 & 2 & 12 & 2 & 0 & 2 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 2 & 1 & 0 & 2 & 6 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 2 & 0 & 0 & 24 & 2 & 0 & 0 & 2 & 2 & 0 & 2 & 2 & 0 \\ 0 & 0 & 0 & 2 & 2 & 0 & 2 & 24 & 2 & 0 & 0 & 0 & 0 & 0 & 2 & 2 \\ 0 & 0 & 0 & 0 & 2 & 1 & 0 & 2 & 6 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix} $$

La relation (91) est non triviale pour $i = 1,2,3,6,9$ (ce sont les numéros des sommets sur le bord $\Gamma_2$) et nous avons finalement :

$$ (93) \quad \left\{ \begin{aligned} (\Gamma G)_1 &= \frac{2}{3} g_1 + \frac{1}{6} g_2 + \frac{1}{6} g_{10} \\ (\Gamma G)_2 &= \frac{1}{6} g_1 + \frac{2}{3} g_2 + \frac{1}{6} g_3 \\ (\Gamma G)_3 &= \frac{1}{6} g_2 + \frac{2}{3} g_3 + \frac{1}{6} g_6 \\ (\Gamma G)_6 &= \frac{1}{6} g_3 + \frac{2}{3} g_6 + \frac{1}{6} g_9 \\ (\Gamma G)_9 &= \frac{1}{6} g_6 + \frac{2}{3} g_9 + \frac{1}{6} g_{16} \end{aligned} \right. $$
puisque nous avons :

$$
(94) \quad \int_{\text{côté}} \varphi_1^2 d\gamma = \frac{1}{3} \quad \text{et} \quad \int_{\text{côté}} \varphi_1 \varphi_2 d\gamma = \frac{1}{6}.
$$

• Le système d'équations à résoudre s'écrit :

$$
(95) \quad AU = MF + \Gamma G
$$

avec A obtenu également par assemblage des matrices élémentaires de rigidité, ce qui conduit
à :

$$
A =
\begin{bmatrix}
2 & -\frac{1}{2} & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\
-\frac{1}{2} & 2 & -\frac{1}{2} & 0 & -1 & 0 & 0 & 0 & 0 \\
0 & -\frac{1}{2} & 1 & 0 & 0 & -\frac{1}{2} & 0 & 0 & 0 \\
-1 & 0 & 0 & 4 & -1 & 0 & -1 & 0 & 0 \\
0 & -1 & 0 & -1 & 4 & -1 & 0 & -1 & 0 \\
0 & 0 & -\frac{1}{2} & 0 & -1 & 2 & 0 & 0 & -\frac{1}{2} \\
0 & 0 & 0 & -1 & 0 & 0 & 4 & -1 & 0 \\
0 & 0 & 0 & 0 & -1 & 0 & -1 & 4 & -1 \\
0 & 0 & 0 & 0 & 0 & -\frac{1}{2} & 0 & -1 & 2
\end{bmatrix}
$$

Nous laissons au lecteur le soin d'expliciter la matrice B₁ (cf. relation (73)) pour
transformer le système (95) lorsqu'une condition de Dirichlet non homogène est à prendre en
compte sur Γ₁.
# IX - COMPLÉMENTS

## i) Formule d'intégration par parties

Nous établissons la formule d'intégration par parties :

$$ (1) \quad \int_{\Omega} u \frac{\partial v}{\partial x_j} dx = - \int_{\Omega} \frac{\partial u}{\partial x_j} v dx + \int_{\partial \Omega} uv n_j d\gamma $$

dans le cas où $\Omega$ est un ouvert de $\mathbb{R}^2$ de frontière polygonale.

* Nous commençons par écrire la relation (1) sous la forme équivalente suivante :

$$ (2) \quad \int_{\Omega} \frac{\partial w}{\partial x_j} dx = \int_{\partial \Omega} w n_j d\gamma $$

* En effet, prenant $u \equiv 1$, $v \equiv w$ dans la relation (1), nous établissons (2) et réciproquement, si (2) est vrai nous posons $w = uv$, ce qui donne (1) par un simple calcul de la dérivée d'un produit de deux fonctions.

Un ouvert de $\mathbb{R}^2$ à frontière polygonale est réunion de triangles donc nous commençons par le cas où $\Omega$ est le **triangle de référence** : $\hat{T}$

$$ (3) \quad \hat{T} \{ (x,y) \mid x \ge 0, y \ge 0, x+y \le 1 \} $$

<figure><img src="image_58.png" /><figcaption>Triangle de référence</figcaption></figure>
La normale extérieure sur le bord de $\hat{T}$ est donnée par :

(4) $$ n = (0, -1) \quad \text{sur } [0,1] \times \{0\} $$

(5) $$ n = \left(\frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}}\right) \quad \text{sur} \ \{(x,y), x \ge 0, y \ge 0, x + y = 1\} $$

(6) $$ n = (-1, 0) \quad \text{sur} \ \{0\} \times [0,1] $$

Nous établissons la relation (2) pour $j = 1$, la propriété s'en déduisant pour $j = 2$ en échangeant les rôles de $x$ et $y$. Nous montrons donc :

(7) $$ \iint_{\hat{T}} \frac{\partial w}{\partial x} \, dx \, dy = \int_{\partial \hat{T}} w n_x \, d\gamma $$

avec :

(8) $$ \int_{\partial \hat{T}} w n_x \, d\gamma = \int_{0}^{1} w (1-t, t) \frac{1}{\sqrt{2}} (\sqrt{2}dt) - \int_{0}^{1} w (0,t) \, dt $$

Nous calculons l'intégrale double du membre de gauche de (7) par la formule de Fubini, en commençant par l'intégrale en $y$ :

(9) $$ \iint_{\hat{T}} \frac{\partial w}{\partial x} \, dx \, dy = \int_{0}^{1} dy \int_{0}^{1-y} dx \frac{\partial w}{\partial x} (x,y) $$

(10) $$ \iint_{\hat{T}} \frac{\partial w}{\partial x} \, dx \, dy = \int_{0}^{1} dy [w(1-y, y) - w(0, y)] $$

qui est, à la notation près, identique au second membre de (8), ce qui établit la propriété dans ce cas.

• Nous établissons la relation (2) dans le cas où $\Omega$ est un **triangle quelconque** $T = (A,B,C)$. Un tel triangle est image de $\hat{T}$ par la transformation suivante :

(11) $$ M = (1 - \hat{x} - \hat{y}) \, A + \hat{x} \, B + \hat{y} \, C, \quad \hat{M} = (\hat{x}, \hat{y}) \in \hat{T} $$

dont la jacolienne $J = \frac{\partial M}{\partial \hat{M}}$ est une constante égale au double de la surface du triangle $T$ :

(12) $$ J = \frac{\partial M}{\partial \hat{M}} = 2 |T| $$
Nous calculons une intégrale sur T par changement de variable en la ramenant sur $\hat{T}$ :

$$ (13) \quad \int_{T} u(M) \, dx = J \int_{\hat{T}} u(M(\hat{M})) \, d\hat{x} $$

donc :

$$ (14) \quad \int_{T} \frac{\partial W}{\partial x_j} \, dx = J \int_{\hat{T}} \sum_{k=1}^{2} \frac{\partial w}{\partial \hat{x}_k} \frac{\partial \hat{x}_k}{\partial x_j} \, d\hat{x} $$

$$ (15) \quad \int_{T} \frac{\partial w}{\partial x_j} dx = J \sum_{k=1}^{2} \frac{\partial \hat{x}_k}{\partial x_j} \int_{\partial \hat{T}} w \hat{n}_k \, d\hat{\gamma} $$

compte tenu du fait que la formule de Green est vraie dans le triangle $\hat{T}$. Nous devons maintenant transformer l'expression de $\hat{n}_k$ et $d\hat{\gamma}$ afin de faire apparaître une intégrale sur le bord du triangle T. Les directions tangentielles de $\partial\hat{T}$ sont images de celles de $\partial T$ dans la transformation $T \to \hat{T}$ inverse de (11). Nous avons donc :

$$ (16) \quad \hat{\tau}_k = \alpha \sum_{k=1}^{2} \frac{\partial \hat{x}_k}{\partial x_m} \tau_m $$

à un coefficient de proportionalité près qui assure que $\hat{\tau}$ et $\tau$ sont unitaires. Cette relation montre que l'élément de longueur $d\hat{\gamma} = \sqrt{d\hat{M}^2}$ peut s'écrire en fonction de $d\gamma = \sqrt{dM^2}$ le long de $\partial T$ selon le changement de variables

$$ (17) \quad d\hat{\gamma} = \left\{ \sum_{k,l,m} \frac{\partial \hat{x}_k}{\partial x_l} \frac{\partial \hat{x}_l}{\partial x_m} \tau_l \tau_m \right\}^{\frac{1}{2}} d\gamma $$

<figure><img src="image_59.png" /><figcaption>Vecteur tangent $\tau$ le long de $\partial T$ et vecteur normal $n$ vers l'extérieur</figcaption></figure>
ou en échangeant les rôles de T et $\hat{T}$ :

$$ (18) \quad d\hat{\gamma} = \left\{ \sum_{k,l,m} \frac{\partial x_k}{\partial \hat{x}_l} \frac{\partial x_k}{\partial \hat{x}_m} \hat{\tau}_k \hat{\tau}_m \right\}^{-\frac{1}{2}} d\gamma $$

De plus, puisque le produit scalaire $\hat{n} \cdot \hat{\tau}$ est nul, nous avons :

$$ (19) \quad \sum_{m=1}^{2} \left\{ \sum_{k=1}^{2} \frac{\partial \hat{x}_k}{\partial x_m} \hat{n}_k \right\} \tau_m = 0 $$

ce qui montre que la mième composante $n_m$ du vecteur normal n est proportionnelle à l'expression entre parenthèses dans la relation (19), puisque $n \cdot \tau$ est nul :

$$ (20) \quad n_m = \beta \sum_k \frac{\partial \hat{x}_k}{\partial x_m} \hat{n}_k $$

avec

$$ (21) \quad \beta^2 \left\{ \sum_m \left( \sum_k \frac{\partial \hat{x}_k}{\partial x_m} \hat{n}_k \right)^2 \right\} = 1 $$

puisque n est un vecteur unitaire. Nous supposons pour simplifier les calculs (et ne pas écrire de signe qui in fine arriveront en nombre pair !) que la transformation $\hat{M} \rightarrow M$ est directe, ie conserve l'orientation. Une direction extérieure à $\hat{T}$ se transforme donc en une direction extérieure à T puisque $\hat{\tau}$ se transforme en $\tau$ et $(\hat{n}, \hat{\tau})$ d'une part, $(n, \tau)$ d'autre part, sont deux angles droits de $+\frac{\pi}{2}$. Donc $\beta$ est positif. Nous obtenons donc, en utilisant la relation (20) au sein du second membre de la relation (15) :

$$ (22) \quad \int_T \frac{\partial w}{\partial x_j} dx = \int_{\partial T} J \frac{1}{\beta} n_j w \frac{d\gamma}{\left\{ \sum_{k,l,m} \frac{\partial x_k}{\partial \hat{x}_l} \frac{\partial x_k}{\partial \hat{x}_m} \hat{\tau}_k \hat{\tau}_m \right\}^{\frac{1}{2}}} $$

compte tenu de la relation (18). La fin de la preuve résulte d'un calcul élémentaire sur l'inverse d'une matrice 2x2. Nous notons provisoirement la jacobienne de la transformation (11) sous la forme :

$$ (23) \quad \frac{\partial M}{\partial \hat{M}} = \begin{pmatrix} a & c \\ b & d \end{pmatrix} $$
et détaillons le calcul de $\frac{1}{\beta^2}$ à l'aide de (21), en tenant compte du fait que :

$$ (24) \qquad \hat{n}_1 = -\hat{\tau}_2, \quad \hat{n}_2 = \hat{\tau}_1 $$

et de la propriété classique :

$$ (25) \qquad \frac{\partial \hat{M}}{\partial \hat{M}} = \frac{1}{J} \begin{pmatrix} d & -c \\ -b & a \end{pmatrix} $$

Nous avons :

$$ \begin{aligned} \frac{1}{\beta^2} &= \frac{1}{J^2} \left[ (d(-\hat{\tau}_2) - c\hat{\tau}_1)^2 + ((-b)(-\hat{\tau}_2) + a\hat{\tau}_1) \right] \\ &= \frac{1}{J^2} \left( (a\hat{\tau}_1 + b\hat{\tau}_2)^2 + (c\hat{\tau}_1 + d\hat{\tau}_2)^2 \right) \end{aligned} $$

c'est-à-dire :

$$ (26) \qquad \frac{1}{\beta^2} = \frac{1}{J^2} \sum_{k \ell m} \frac{\partial x_k}{\partial \hat{x}_\ell} \frac{\partial x_k}{\partial \hat{x}_m} \hat{\tau}_k \hat{\tau}_m $$

La relation (22) se simplifie alors en :

$$ (27) \qquad \int_T \frac{\partial w}{\partial x_j} dx = \int_{\partial T} w n_j d\gamma $$

ce qui prouve la propriété dans le cas où $\Omega$ est un triangle quelconque du plan.

• Nous abordons maintenant le cas général. L'ouvert $\Omega$ étant polygonal, on peut le recouvrir par un maillage de type éléments finis, formé de triangles T :

$$ (28) \qquad \overline{\Omega} = \bigcup_{I \in \tau} \overline{T} $$

où l'intersection $\overline{T} \cap \overline{T}'$ de deux triangles de la triangulation $\tau$ est formée de l'ensemble vide, d'un sommet de T et $T'$, d'une arête commune à T et $T'$ ou est égale à T et $T'$. On découpe alors l'intégrale (2) en autant d'intégrales sur les triangles du maillage.

$$ (29) \qquad \int_{\Omega} \frac{\partial w}{\partial x_j} dx = \sum_{T \in \tau} \int_{\Gamma} \frac{\partial w}{\partial x_j} dx $$
pour lesquelles la relation (2) est applicable :

$$ (30) \quad \int_{\Omega} \frac{\partial w}{\partial x_j} dx = \sum_{I \in \tau} \int_{\partial I} w n_j d\gamma $$

<figure><img src="image_60.png" /><figcaption>Triangulaire d'un ouvert Ω polygonal du plan</figcaption></figure>

Le membre de droite de la relation (30) est une somme d'intégrales sur les arêtes a du maillage. Si l'arête a est intérieure au domaine, elle appartient au bord de deux triangles T et T' donc contribue par les deux intégrales suivantes :

$$ (31) \quad \int_{a \cap \partial T} w n_j d\gamma + \int_{a \cap \partial T'} w n_j d\gamma = 0 $$

dont la somme est nulle car les normales extérieures à T et T' sont opposées.

<figure><img src="image_61.png" /><figcaption>Arête a intérieure au domaine Ω</figcaption></figure>

Si l'arête a est au bord du domaine, la normale extérieure à T coincide avec la normale extérieure à ∂Ω et l'arête n'apparaît qu'une seule fois dans la somme (30) puisqu'elle n'appartient qu'à un seul triangle du maillage. Nous avons donc :
$$
(32) \quad \int_{\Omega} \frac{\partial w}{\partial x_j} dx = \sum_{\substack{a \in \tau \\ a \subset \partial\Omega}} \int_a w n_j d\gamma
$$

qui peut s'écrire plus simplement (puisque la réunion des arêtes du bord de Ω constitue le bord
lui-même) :

$$
(33) \quad \int_{\Omega} \frac{\partial w}{\partial x_j} dx = \int_{\partial \Omega} w n_j d\gamma,
$$

ce qui montre la propriété.
## ii) Algorithme du gradient conjugué

A désigne une matrice carrée d'ordre n symétrique définie positive et b un vecteur de $\mathbb{R}^n$.

### 1) Montrer que la résolution du système linéaire

(1) $A x = b$

équivaut à rechercher le minimum, pour $x$ appartenant à $\mathbb{R}^n$, de la fonctionnelle.

(2) $J(x) = \frac{1}{2} (x, A x) - (b, x)$

### 2) Description de l'algorithme du gradient conjugué

• Initialisation : $x^0 \in \mathbb{R}^n$, $w^0 = g^0 = A x^0 - b$.

• Itération de l'algorithme.

On suppose connus l'état $x^{k-1}$ après la kième itération ainsi que la direction de descente $w^{k-1}$. L'état suivant $x^k$ de l'algorithme est issu de $x^{k-1}$ via un incrément proportionnel à la direction de descente :

(3) $x^k = x^{k-1} + \rho^k w^{k-1}$

de façon à minimiser la fonctionnelle $J$ :

(4) $J(x^k) \le J(x^{k-1} + \rho^k w^{k-1}) \quad \forall \rho \in \mathbb{R}$

On introduit également le gradient de $J$ au point $x^k$.

(5) $g^k = A x^k - b$.

2.1) Calculer le coefficient $\rho^{k-1}$ en fonction de $g^{k-1}, w^{k-1}$ et $A$.

2.2) Montrer que l'on a :

(6) $(g^k, w^{k-1}) = 0$.
La nouvelle direction de descentre $w^k$ est recherchée sous la forme :

$$ (7) \quad w^k = g^k + \alpha^k w^{k-1} $$

de sorte que $w^k$ soit orthogonal à la direction $w^{k-1}$ pour le produit scalaire associé à la matrice A, c'est-à-dire :

$$ (8) \quad (w^k, A w^{k-1}) = 0 $$

### 2.3) Calculer la valeur de $\alpha^k$ en fonction de $g^k, w^{k-1}$ et A.

On notera que le choix plus simple $\alpha^k = 0$ conduit à la méthode du gradient simple, où la direction de descente de l'algorithme $w^k$ correspond à la ligne de plus grande pente $g^k$.

## 3) L'algorithme converge en au plus n étapes

3.1) Remarquer que si les gradients successifs $g^0, g^1, ..., g^{k-1}$ sont non nuls mais que $g^k$ est nul, alors l'état $x^k$ est la solution du système linéaire (1).

3.2) On suppose dans cette question que tous les gradients $g^i$ sont non nuls jusqu'à l'étape m incluse. Montrer qu'on a alors les relations d'orthogonalité suivantes :

$$ (9) \quad (g^k, w^j) = 0, \quad 0 \le j < k \le m $$

$$ (10) \quad \rho^j \neq 0 \text{ et } w^j \neq 0, \quad 0 \le j < k \le m $$

$$ (11) \quad (g^k, g^j) = 0, \quad 0 \le j < k \le m $$

$$ (12) \quad (w^k, A w^j) = 0, \quad 0 \le j < k \le m $$

On pourra raisonner par récurrence sur $k$ pour l'ensemble des quatre propriétés et démontrer les relations (9) à (12) dans l'ordre où elles apparaissent.

3.3) Montrer que l'algorithme converge en au plus n itérations.
### 4) Propriétés auxiliaires

4.1) Montre que l'on a :

$$ (13) \quad \alpha^k = \frac{\|g^k\|^2}{\|g^{k-1}\|^2} $$

En déduire que le calcul très précis des produits scalaires est crucial pour assurer le succès pratique (informatique) de la convergence de la méthode.

4.2) Montre que l'état $x^k$ réalise le minimum de la fonctionnelle J sur le sous-espace affine $x^0 + \langle w^0, \dots, w^{k-1} \rangle$. Quel commentaire pouvez-vous faire ?

#### Solution

1) La fonctionnelle a bien un minimum car la matrice A est symétrique définie positive. L'équation d'Euler au minimum s'écrit ici :

$$ J'(x) \equiv Ax - b = 0, $$

ce qui montre la propriété.

$$ 2.1) \quad (14) \quad \rho^{k-1} = - \frac{\begin{pmatrix} g^{k-1} & w^{k-1} \\ w^{k-1} & A w^{k-1} \end{pmatrix}}{\begin{pmatrix} g^{k-1} & w^{k-1} \\ w^{k-1} & A w^{k-1} \end{pmatrix}} $$

2.2) La condition précédente a été obtenue en écrivant qu'au point $x^k$, la fonctionnelle J est minimale si on la restreint à la droite affine passant par $x^{k-1}$ et dirigée par le vecteur $w^{k-1}$, ce qui s'exprime en annulant le gradient de J au point $x^k$ dans la direction $w^{k-1}$, ce qu'exprime exactement la condition (6).

$$ 2.3) \quad (15) \quad \alpha^k = - \frac{\begin{pmatrix} g^k & A w^{k-1} \\ w^{k-1} & A w^{k-1} \end{pmatrix}}{\begin{pmatrix} g^k & w^{k-1} \\ w^{k-1} & A w^{k-1} \end{pmatrix}} $$

3.1) Dès que l'un des gradients $g^m$ est nul, on est en un point $x^m$ qui est solution du système (1) et l'algorithme ne présente plus alors aucun intérêt puisque l'on a résolu le problème posé.
### 3.2) • La propriété est vraie pour $k=1$

On remarque que $w^0 = g^0$ est non nul en vertu de la question précédente. On a ensuite $(g^1, w^0) = 0$ compte tenu de la relation (6). Si on considère la relation (7) pour $k=1$ et qu'on la multiplie scalairement par $g^1$, on obtient $(g^1, w^1) = \|g^1\|^2$ ce qui montre que $w^1$ est non nul puisque $g^1$ est non nul et le numérateur comme le dénominateur de l'expression (14) qui permet de calculer $\rho^0$ sont non nuls, donc $\rho^0$ est non nul et la relation (10) est vraie à l'ordre $j=0$. La relation (11) est une conséquence simple du choix de la direction de descente initiale : $(g^1, g^0) = (g^1, w^0) = 0$ compte tenu de la relation (9). Enfin, la relation (12) exprime simplement la relation (8) pour $k=1$.

* On suppose les relations (9) à (12) vérifiées jusqu'à l'ordre $k$ inclus et on les étend pour l'indice $k+1$, en supposant $g^{k+1}$ non nul.

On remarque d'abord que la relation (3) entraîne clairement :

$$ (16) \quad g^{k+1} = g^k + \rho^k A w^k $$

On a d'une part $(g^{k+1}, w^k) = 0$ compte tenu de la relation (6) et pour $j < k$, on a d'autre part :

$$ \begin{aligned} (g^{k+1}, w^j) &= (g^{k+1} - g^k, w^j) + (g^k, w^j) \\ &= \rho^k (A w^k, w^j) \end{aligned} $$

compte tenu de (16) et de l'hypothèse de récurrence (9). Or cette dernière expression est nulle en vertu de l'hypothèse de récurrence (12), donc la relation (9) est établie.

Si on multiplie scalairement l'identité (7) écrite au rang $k+1$ par $g^{k+1}$, la relation (9) que nous venons de démontrer entraîne que l'on a :

$$ (17) \quad (g^{k+1}, w^{k+1}) = \|g^{k+1}\|^2 $$

et la relation (10) est alors une conséquence simple de la relation qui permet de calculer $\rho^k$.

Exprimons $g^j$ grâce à la relation (7) (considérée avec $k=j$). Il vient :

$$ (g^{k+1}, g^j) = (g^{k+1}, w^j) - \alpha^j (g^{k+1}, w^{j-1}) $$
et cette expression est nulle compte tenu de la relation (9) prise à l'ordre k+1. Ceci montre la relation (11).

La relation (12) est la conséquence de la relation (8) pour $j=k$ et du calcul suivant $j<k$ :

$$ 
\begin{align*} 
(w^{k+1}, A w^j) &= (g^{k+1}, A w^j) \text{ compte tenu de (7) et de (12),} \\ 
&= \left(g^{k+1}, \frac{1}{\rho^j} (g^{j+1} - g^j)\right) \text{en vertu de (16),} \\ 
&= 0 \text{ compte tenu de (11).} 
\end{align*} 
$$

La propriété est donc démontrée par récurrence.

3.3) Compte tenu de la relation (11), la suite de gradients $g^0, g^1, ..., g^k$ est composée de vecteurs orthogonaux deux à deux jusqu'à un ordre m et ceci a lieu dans l'espace $\mathbb{R}^n$. La conclusion est alors claire.

4.1) Dans la relation (15), on exprime $A w^{k-1}$ au numérateur et au dénominateur à l'aide de la relation (16) et on utilise la relation (11) pour le numérateur. Il vient :

$$ \alpha^k = \frac{\|g^k\|^2}{(g^{k-1}, w^{k-1})} $$

et la relation (13) est alors conséquence directe de la relation (17).

4.2) Si on écrit l'inéquation d'Euler qui caractérise le minimum de la fonction convexe dérivable J sur le convexe fermé $x^0 + \langle w^0, \dots, w^{k-1} \rangle$, atteint au point y, il vient $J'(y) \cdot w^j = 0$ pour tout $j = 1, 2, \dots, k-1$. Mais le choix $y = x^k$ vérifie ces relations puisque d'une part $x^k$ appartient au convexe $x^0 + \langle w^0, \dots, w^{k-1} \rangle$ compte tenu de l'initialisation de l'algorithme et de la relation constitutive (3), et d'autre part $J'(x^k) \cdot w^j = (g^k, w^j)$ lequel est nul en vertu de la relation (9).

Cette propriété montre que le point $x^k$, conçu initialement à la relation (4) pour minimiser la fonctionnelle J le long de la droite affine passant par $x^{k-1}$ et dirigée selon le vecteur $w^{k-1}$ la minimise en fait sur tout un sous-espace de dimension k ! On touche là au génie de Hestenes et Stiefel qui ont inventé la méthode du gradient conjugué en 1952.
### iii) Étude du schéma de Newmark

On se propose de définir une famille à deux paramètres pour approcher numériquement la solution du système différentiel.

$$ (D) \quad M \frac{d^2 q}{dt^2} + C \frac{dq}{dt} + K q = p(t) $$

où $q(t)$ est un vecteur d'état appartenant à $\mathbb{R}^m$ (modèle de structure à m degrés de liberté). On suppose connu à l'instant $t_n = n\Delta t$ le vecteur d'état $U = \begin{pmatrix} \dot{q} \\ q \end{pmatrix}$ et on cherche à calculer ce même vecteur à l'instant $t_{n+1} = (n+1)\Delta t$.

1) Écrire la formule de Taylor avec reste intégral pour $\dot{q}_{n+1}$ et $q_{n+1}$, en faisant apparaître la dérivée seconde $\ddot{q}(\theta)$ ($t_n < \theta < t_{n+1}$).

2) Quel schéma obtient-on si on suppose que dans la relation précédente, on remplace $\ddot{q}(\theta)$ par la moyenne, ie $\ddot{q}(\theta) \approx \frac{1}{2}(\ddot{q}_n + \ddot{q}_{n+1})$ ?

3) Même question, mais avec $\ddot{q}(\theta)$ remplacé par l'interpolé affine entre les valeurs $\ddot{q}_n$ et $\ddot{q}_{n+1}$ :

$$ \ddot{q}(\theta) = \ddot{q}_n + \frac{1}{\Delta t} (\theta - t_n) (\ddot{q}_{n+1} + \ddot{q}_n) $$

Le schéma de Newmark consiste à écrire :

$$ (N) \quad \begin{cases} \dot{q}_{n+1} = \dot{q}_n + (1-\gamma) \Delta t \ddot{q}_n + \gamma \Delta t \ddot{q}_{n+1} \\ q_{n+1} = q_n + \Delta t \dot{q}_n + \Delta t^2 (\frac{1}{2} - \beta) \ddot{q}_n + \Delta t^2 \beta \ddot{q}_{n+1} \end{cases} $$

4) A quelles valeurs de $(\beta, \gamma)$ correspondent les schémas calculés aux questions 2 et 3 ?
5) Montrer que si $\gamma \neq \frac{1}{2}$, le schéma est d'ordre 1 quelque soit $\beta$ et pour $\gamma = \frac{1}{2}$, le schéma est d'ordre 2 si $\beta \neq \frac{1}{6}$ et d'ordre 3 si $\beta = \frac{1}{6}$.

6) Proposer un algorithme de résolution du schéma (N) qui prenne en compte la loi de la dynamique (D). (On pourra résoudre un système linéaire du type $S \dot{\ddot{q}}_{n+1} = f$ avec $S = M + \gamma \Delta t C + \beta \Delta t^2 K$ et f un second membre qu'on précisera).

7) Pour étudier la stabilité du schéma, on se restreint dans la suite au cas où $M = 1$, $C = 0$, $K = \omega^2$, $p = 0$ dans le modèle (D). Montrer qu'alors le schéma de Newmark (N) s'écrit :

$$U_{n+1} = A U_n$$

où A est une matrice 2x2 qu'on précisera, on pourra introduire $\xi$ tel que :

$$\xi^2 = \frac{\omega^2 \Delta t^2}{1 + \beta \omega^2 \Delta t^2}$$

8) Montrer, en étudiant le module des racines complexes conjuguées de l'équation caractéristique relative à la matrice A que :

* le schéma est instable pour $\gamma < \frac{1}{2}$,

* le schéma est stable inconditionnellement pour $\gamma \ge \frac{1}{2}$ et $\beta \ge \frac{1}{4} (\gamma + \frac{1}{2})^2$,

* le schéma est stable sous la condition

$$(5) \quad \left(\gamma + \frac{1}{2}\right)^2 - 4\beta \le \frac{4}{\omega^2 \Delta t^2}$$

$$\text{lorsque } \gamma \ge \frac{1}{2} \text{ et } \beta < \frac{1}{4} \left(\gamma + \frac{1}{2}\right)^2$$

On tracera dans le plan $(\gamma,\beta)$ ces domaines de stabilité.

9) Quel schéma "optimum" proposez-vous de choisir ?
10) Montrer que l'erreur relative d'amplitude des ondes d'une part calculées par le schéma de Newmark et d'autre part exacte est donnée par le développement.

$$ \rho - 1 = -\frac{1}{2} \left(\gamma - \frac{1}{2}\right) \omega^2 \Delta t^2 + O(\Delta t^4) $$

et que l'erreur de relative de périodicité, dans le cas $\gamma = \frac{1}{2}$, est donnée par :

$$ \frac{\Delta T}{T} = \frac{\omega \Delta t}{\varphi} - 1 = \frac{1}{2} \left(\beta - \frac{1}{2}\right) \omega^2 \Delta t^2 + O(\Delta t^3) $$

où $\lambda = \rho \exp (\pm i\varphi)$ sont les racines (complexes) de l'équation caractéristique de la matrice A.

{plus difficile}.

Corrigé

$$ \begin{array}{l} \text{1) } \dot{q}_{n+1} = \dot{q}_n + \int_{t_n}^{t_{n+1}} \ddot{q}(\theta) \, d\theta \\ q_{n+1} = q_n + \Delta t \dot{q}_n + \int_{t_n}^{t_{n+1}} (t_{n+1} - \theta) \ddot{q}(\theta) \, d\theta . \end{array} $$

2) On remplace $\ddot{q}(\theta)$ par $\frac{1}{2}(\ddot{q}_n + \ddot{q}_{n+1})$ dans les relations précédentes. On obtient :

$$ \begin{aligned} \dot{q}_{n+1} &= \dot{q}_n + \frac{\Delta t}{2} (\ddot{q}_n + \ddot{q}_{n+1}) \\ q_{n+1} &= q_n + \Delta t \dot{q}_n + \Delta t^2 (\frac{1}{4} \ddot{q}_n + \frac{1}{4} \ddot{q}_{n+1}) \end{aligned} $$

3) Même méthode mais pour $\ddot{q}(\theta) = \ddot{q}_n + \frac{1}{\Delta t}(\theta - t_n)(\ddot{q}_{n+1} + \ddot{q}_n)$. On trouve :

$$ \begin{aligned} \dot{q}_{n+1} &= \dot{q}_n + \frac{\Delta t}{2} (\ddot{q}_n + \ddot{q}_{n+1}) \\ q_{n+1} &= q_n + \Delta t \dot{q}_n + \Delta t^2 (\frac{1}{3} \ddot{q}_n + \frac{1}{6} \ddot{q}_{n+1}) \end{aligned} $$
4) Question 2 : $\gamma = \frac{1}{2}$, $\beta = \frac{1}{4}$ ; question 3 : $\gamma = \frac{1}{2}$, $\beta = \frac{1}{6}$

5) On suppose que $q(t)$ vérifie l'équation différentielle (D), c'est-à-dire le système différentiel suivant pour $U(t)$ :

$$ \frac{dU}{dt} = \begin{bmatrix} -M^{-1}C & -M^{-1}K \\ 1 & 0 \end{bmatrix} U + \begin{pmatrix} M^{-1}p(t) \\ 0 \end{pmatrix} $$

On commence donc par mettre le schéma (N) sous une forme analogue, en remplaçant la dérivée $\frac{dU}{dt}$ par le quotient des différences finies.

$$ \frac{1}{\Delta t} (U_{n+1} - U_n) = \begin{pmatrix} (1-\gamma) \ddot{q}_n + \gamma \ddot{q}_{n+1} \\ \dot{q}_n + \Delta t (\frac{1}{2} - \beta) \ddot{q}_n + \Delta t \beta \ddot{q}_{n+1} \end{pmatrix} $$

Puis l'erreur de troncature est ici vectorielle :

$$ \tau = \frac{1}{\Delta t} (U_{n+1} - U_n) - \begin{pmatrix} (1-\gamma) \ddot{q}_n + \gamma \ddot{q}_{n+1} \\ \dot{q}_n + \Delta t (\frac{1}{2} - \beta) \ddot{q}_n + \Delta t \beta \ddot{q}_{n+1} \end{pmatrix} $$

Compte tenu des formules de Taylor établies à la question 1, il suffit d'évaluer les restes dans les formules de quadrature.

$$ (1) \quad \frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} \ddot{q}(\theta) \, d\theta \approx (1-\gamma) \ddot{q}_n + \gamma \ddot{q}_{n+1} $$

$$ (2) \quad \frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} (t_{n+1} - \theta) \ddot{q}(\theta) \, d\theta \approx \Delta t \left[ \left(\frac{1}{2} - \beta\right) \ddot{q}_n + \beta \ddot{q}_{n+1} \right] $$

pour trouver l'ordre demandé.

* La formule (1) est a priori du premier ordre, sauf si $\gamma = \frac{1}{2}$ et elle est alors du 2ème ordre (formule des trapèzes). L'ordre global du schéma est donc 1 si $\gamma \neq \frac{1}{2}$.
Si $\gamma = \frac{1}{2}$, la formule des trapèzes est d'ordre 2, et un surcroît de précision peut être donné par la formule de quadrature (2). Celle-ci est d'ordre deux si $\ddot{q}(\theta)$ est approchée par n'importe quelle constante compte tenu de l'infiniment petit ($t_{n+1}-\theta$) à intégrer, et elle est d'ordre trois si $\dddot{q}(\theta)$ est approchée par l'interpolation affine proposée à la question 3, ie pour $\beta = \frac{1}{6}$ (cf question 4).

## 6) Algorithme

Dans les relations (N), les dérivées secondes sont calculées à l'aide de l'équation (D). Mais si $\beta \neq 0$, la présence de $\dddot{q}_{n+1}$ au second membre de la seconde équation rend le schéma implicite. Il convient donc de calculer d'abord ce terme, en tenant compte du schéma (N) et de l'équation (D) écrite à l'instant $t_{n+1}$.

$$M \ddot{q}_{n+1} + C \dot{q}_{n+1} + K q_{n+1} = p_{n+1}$$

On introduit dans cette relation $\dot{q}_{n+1}$ et $q_{n+1}$ calculés par le schéma :

$$M \ddot{q}_{n+1} + C \{ \dot{q}_n + (1-\gamma) \Delta t \ddot{q}_n + \gamma \Delta t \ddot{q}_{n+1} \} + K \{ q_n + \Delta t \dot{q}_n + \Delta t^2 (\frac{1}{2} - \beta) \ddot{q}_n + \Delta t^2 \beta \ddot{q}_{n+1} \} = p_{n+1}$$

D'où le système à résoudre :

$$S \ddot{q}_{n+1} = f$$

$$où \quad S = M + \gamma \Delta t C + K \beta \Delta t^2$$

$$f = p_{n+1} - C \left( \dot{q}_n + (1-\gamma) \Delta t \ddot{q}_n \right) - K \left( q_n + \Delta t \dot{q}_n + \left(\frac{1}{2} - \beta\right) \Delta t^2 \ddot{q}_n \right)$$

* Algorithme

* $q_o, \dot{q}_o$ connus.

Calculer $\ddot{q}_o$ par résolution du système (D).
• $q_n, \dot{q}_n, \ddot{q}_n$ connus

$$
(i) \quad \begin{aligned}
\dot{q}^* &= \dot{q}_n + (1-\gamma) \Delta t \ddot{q}_n \\
q^* &= q_n + \Delta t \dot{q}_n + \left(\frac{1}{2} - \beta\right) \Delta t^2 \ddot{q}_n
\end{aligned}
\quad \text{(prédiction)}
$$

$$
(ii) \quad S \ddot{q}_{n+1} = p_{n+1} - C \dot{q}^* - K q^* \quad \text{(résolution)}
$$

$$
(iii) \quad
\begin{aligned}
\dot{q}_{n+1} &= \dot{q}^* + \Delta t \gamma \ddot{q}_{n+1} \\
q_{n+1} &= q^* + \Delta t^2 \beta \ddot{q}_{n+1}
\end{aligned}
\qquad \text{( incrémentation)}
$$

7) Dans le cas simplifié d'un seul oscillateur non amorti, l'algorithme précédent s'explicite simplement :

$$
\begin{cases}
\dot{q}_{n+1} = \dot{q}_n + (1-\gamma) \Delta t (-\omega^2 q_n) + \gamma \Delta t (-\omega^2 q_{n+1}) \\
q_{n+1} = q_n + \Delta t \dot{q}_n + \Delta t^2 (\frac{1}{2} - \beta) (-\omega^2 q_n) + \Delta t^2 \beta (-\omega^2 q_{n+1})
\end{cases}
$$

c'est-à-dire :

$$
\begin{pmatrix} 1 & \gamma \Delta t \omega^2 \\ 0 & 1+\beta \Delta t^2 \omega^2 \end{pmatrix} \begin{pmatrix} \dot{q} \\ q \end{pmatrix}_{n+1} = \begin{pmatrix} 1 & -(1-\gamma) \Delta t \omega^2 \\ \Delta t & 1-(\frac{1}{2}-\beta) \Delta t^2 \omega^2 \end{pmatrix} \begin{pmatrix} \dot{q} \\ q \end{pmatrix}_n
$$

D'où le résultat, avec :

$$
A = \begin{pmatrix} 1 - \gamma \xi^2 - \omega^2 \Delta t (1 - \frac{\gamma}{2} \xi^2) \\ \frac{\xi^2}{\omega \Delta t} - \frac{1}{2} \xi^2 \end{pmatrix}
$$

et
$$
\xi^2 = \frac{\Delta t^2 \omega^2}{1 + \beta \omega^2 \Delta t^2}
$$
$$8) \quad \text{On a det } (A - \lambda I) = \lambda^2 - \lambda \left(2 - \left(\gamma + \frac{1}{2}\right) \xi^2\right) + 1 - \left(\gamma - \frac{1}{2}\right) \xi^2$$

Les racines sont complexes si le discriminant réduit $\Delta'$ est négatif. Or :

$$\begin{align*} \Delta' &= \left(1 - \frac{1}{2}\left(\gamma + \frac{1}{2}\right) \xi^2\right)^2 - \left(1 - \left(\gamma - \frac{1}{2}\right) \xi^2\right) \\ &= \xi^2 \left(\frac{1}{4} \left(\gamma + \frac{1}{2}\right)^2 \xi^2 - 1\right) \end{align*}$$

$\Delta' \le 0$ équivaut à :

$$(\gamma + \frac{1}{2})^2 \omega^2 \Delta t^2 \le 4 \left(1 + \beta \omega^2 \Delta t^2\right)$$

c'est-à-dire :

$$(\gamma + \frac{1}{2})^2 - 4\beta \le \frac{4}{\omega^2 \Delta t^2} \quad (S)$$

Cette condition limite le pas de temps, sauf si $(\gamma + \frac{1}{2})^2 - 4\beta \le 0$, auquel cas elle ne donne aucune condition. Par ailleurs, le produit des deux racines est égal au carré de leur module commun, c'est-à-dire $1 - (\gamma - \frac{1}{2}) \xi^2$. La condition $|\lambda|^2 \le 1$ équivaut à $(\gamma - \frac{1}{2}) \xi^2 \ge 0$ ie $\gamma \ge \frac{1}{2}$.

Les deux racines sont réelles si $\Delta' \ge 0$. La condition :

$$(\gamma + \frac{1}{2})^2 - 4\beta > \frac{4}{\omega^2 \Delta t^2}$$

donne une valeur minimale du pas de temps, ce qui indique que pour les valeurs toutes petites de $\Delta t$, le schéma est instable. Ce type de comportement est en général rejeté par les utilisateurs, et nous ne considérons pas ce cas dans la discussion qui suit.

* Si $\gamma \le \frac{1}{2}$, le schéma est instable si $\gamma \ge \frac{1}{2}$ et $(\gamma + \frac{1}{2})^2 < 4\beta$, le schéma est stable lorsque le pas de temps $\Delta t$ satisfait à la condition (S).
<figure><img src="image_62.png" /></figure>

9) Le schéma optimal correspond à la limite de stabilité inconditionnelle, c'est-à-dire
$\gamma = \frac{1}{2}$, $\beta = \frac{1}{4}$, qui donne un schéma d'ordre **deux**.

10) Dans le cas de deux valeurs propres complexes pour l'équation caractéristique, on a :

$$
\lambda = \rho \exp (i\phi)
$$

$$
\text{avec : } \rho^2 = 1 - \left(\gamma - \frac{1}{2}\right) \xi^2
$$

$$
\operatorname{tg} \varphi=\frac{\xi \sqrt{1-\frac{1}{4}\left(\gamma+\frac{1}{2}\right)^{2} \xi^{2}}}{1-\frac{1}{2}\left(\gamma+\frac{1}{2}\right) \xi^{2}}
$$

L'erreur relative de phase et d'amplitude s'obtient en comparant la solution exacte
de l'équation q̈ + ω² q = 0 à la solution numérique calculée par le schéma, c'est-à-dire λⁿ
à une constante près. On doit donc comparer (e^(iω/n))^n à λⁿ, avec λ donné plus haut.

L'erreur d'amplitude s'obtient en comparant les modules, or :

$$
\begin{align*}
\rho-1 &= \frac{1}{\rho+1} (\rho^2-1) = \frac{1}{2} \left( -\left(\gamma-\frac{1}{2}\right) \xi^2 \right) + 0 \left(\xi^2\right) \\
&= -\frac{1}{2} \left(\gamma-\frac{1}{2}\right) \Delta t^2 \omega^2 + 0 \left(\Delta t^4\right)
\end{align*}
$$
ce qui montre que le choix $\gamma = \frac{1}{2}$ conduit à une erreur d'ordre 4 en $\Delta t$. Pour l'erreur de phase, on compare l'écart de périodes entre solution exacte et numérique, ramenée à la période exacte.

$$ \frac{\Delta T}{T} = \frac{\omega}{2\pi} \left( \frac{2\pi}{\varphi} - \frac{2\pi}{\omega} \right) $$

$$ \text{ie} \quad \frac{\Delta T}{T} = \frac{\omega h}{\varphi} - 1 $$

$$
\begin{align*}
\text{or } \varphi &= \operatorname{Arctg}(\operatorname{tg} \varphi) \\
&= (\operatorname{tg} \varphi) - \frac{1}{3} (\operatorname{tg} \varphi)^3 + 0 (\operatorname{tg} \varphi)^5 \\
&= \xi \left(1 + \frac{1}{2} \left(\gamma + \frac{1}{2}\right) \xi^2\right) \left(1 - \frac{1}{8} \left(\gamma + \frac{1}{2}\right) \xi^2\right) - \frac{1}{3} \xi^2 + 0 (\xi^5) \\
&= \xi \left(1 + \left[\frac{1}{2} \left(\gamma + \frac{1}{2}\right) - \frac{1}{8} \left(\gamma + \frac{1}{2}\right)^2 - \frac{1}{3}\right] \xi^2\right) + 0 (\xi^5) \\
&= \omega \Delta t \left(1 - \frac{\beta}{2} \omega^2 \Delta t^2\right) \left(1 + \frac{1}{24} \omega^2 \Delta t^2\right) + 0 (\Delta t^4) \quad \text{si } \gamma = \frac{1}{2} \\
&= \omega \Delta t \left(1 - \frac{1}{2} \left(\beta - \frac{1}{12}\right) \omega^2 \Delta t^2\right) + 0 (\Delta t^4) \\
\frac{1}{\varphi} &= \frac{1}{\omega \Delta t} \left(1 + \frac{1}{2} \left(\beta - \frac{1}{12}\right) \omega^2 \Delta t^2\right) + 0 (\Delta t^3) \\
\frac{\omega \Delta t}{\varphi} &= 1 + \frac{1}{2} \left(\beta - \frac{1}{12}\right) \omega^2 \Delta t^2 + 0 (\Delta t^3)
\end{align*}
$$

ce qui montre le résultat proposé.

<figure><img src="image_63.png" /></figure>