# Sur la résolution, via une approximation par éléments finis mixtes, du problème de dirichlet pour l'opérateur biharmonique, par une méthode quasi-directe et diverses méthodes itératives

Roland Glowinski, Olivier Pironneau

► To cite this version:

Roland Glowinski, Olivier Pironneau. Sur la résolution, via une approximation par éléments finis mixtes, du problème de dirichlet pour l'opérateur biharmonique, par une méthode quasi-directe et diverses méthodes itératives. [Rapport de recherche] IRIA-RR-197, IRIA. 1976, pp.31. hal-04716341

HAL Id: hal-04716341

https://inria.hal.science/hal-04716341v1

Submitted on 1 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
RIA
laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B. P. 5 78150 - Le Chesnay
France
Tél.: 954 90 20

laboratoire de recherche
en informatique
et automatique

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

SUR LA RÉSOLUTION,
VIA UNE APPROXIMATION
PAR ÉLÉMENTS FINIS MIXTES,
DU PROBLÈME DE DIRICHLET
POUR L'OPÉRATEUR
BIHARMONIQUE,
PAR UNE MÉTHODE
"QUASI - DIRECTE"
ET DIVERSES
MÉTHODES ITÉRATIVES

R. GLOWINSKI - O. PIRONNEAU

Rapport de Recherche N° 197
Novembre 1976
SUR LA RÉSOLUTION,
VIA UNE APPROXIMATION PAR ÉLÉMENTS FINIS MIXTES,
DU PROBLÈME DE DIRICHLET POUR L'OPÉRATEUR BIHARMONIQUE,
PAR UNE MÉTHODE "QUASI-DIRECTE"
ET DIVERSES MÉTHODES ITÉRATIVES

R. Glowinski *) - O. Pironneau **)

Résumé :

Dans ce travail on montre que le problème de Dirichlet pour l'opérateur biharmonique peut être résolu par une méthode de décomposition équivalente à une formulation mixte. Le problème est alors équivalent à la résolution de $2n + 4$ problèmes de Dirichlet pour l'opérateur $\Delta$, où $n$ est le nombre de nœuds frontières, plus un petit système linéaire symétrique très bien conditionné. Ce petit système peut être résolu de manière directe ou itérative. La méthode permet donc d'obtenir d'une manière très simple un code éléments finis pour $\Delta^2$ à partir d'un code éléments finis pour $\Delta$.

Abstract:

In this report it is shown that the Dirichlet problem for the biharmonic operator can be decomposed into $2n + 4$ Dirichlet problems for the Laplace operator $\Delta$ plus a small linear system symmetric positive definite and well conditioned ($n$ is the number of boundary nodes). It is also shown that this method is equivalent to a mixed variational formulation of the problem. The Dirichlet problems for $\Delta$ are discretized with finite elements and the small linear system is solved directly or by an iterative process. Therefore this method leads in a straightforward fashion from a finite element code for $\Delta$ to the construction of a finite element code for $\Delta^2$.

<footnote>*) Université de Paris VI, Analyse Numérique, L. A. 189, Tour 55.65, 5ᵉ étage, 4 Place Jussieu, 75230 PARIS CEDEX 05, FRANCE.</footnote>

<footnote>**) IRIA-LABORIA, Domaine de Voluceau, B. P. 105, 78150 LE CHESNAY, FRANCE.</footnote>
# 1, INTRODUCTION. ORIENTATION.

Soit $\Omega$ un domaine borné de $\mathbb{R}^2$, de frontière $\Gamma$ régulière. On présente dans ce travail diverses méthodes permettant de résoudre numériquement le Problème de Dirichlet pour l'Opérateur Biharmonique, soit

$$ (1.1) \qquad (P_o) \quad \begin{cases} \Delta^2 \psi = f \text{ dans } \Omega, \\ \psi|_{\Gamma} = g_1, \\ \frac{\partial \psi}{\partial n}|_{\Gamma} = g_2. \end{cases} $$

Une partie des résultats exposés ci-après (sans démonstration) ont été annoncés dans GLOWINSKI-PIRONNEAU [1], [2] ; pour des démonstrations complètes on renvoie à GLOWINSKI-PIRONNEAU [3], [4]. Pour l'extension des méthodes décrites dans ce travail aux problèmes d'évolution et aux équations de Navier-Stokes dans la formulation {$\psi, \omega$}^{(1)}, on renvoie à des articles en préparation de GLOWINSKI-PIRONNEAU et BOURGAT-GLOWINSKI-PIRONNEAU. Si $\Omega$ est p-connexe ($p \ge 1$) la formulation en {$\psi, \omega$} des équations de Stokes et Navier-Stokes est plus compliquée, mais les méthodes ci-après peuvent être généralisées à ces problèmes (cf. GLOWINSKI [2] et les articles en préparation ci-dessus).

Décrivons brièvement le contenu de ce travail :

On donne au N° 2 un certain nombre de résultats relatifs au problème continu. Aux N° 3, 4 il est montré que par utilisation d'une approximation par éléments finis mixte de $(P_o)$ adéquate (cf. CIARLET-RAVIART [1]), résoudre $(P_o)$ équivaut à résoudre un nombre fini de problèmes de Dirichlet approchés pour $-\Delta$, plus un "petit" système linéaire de matrice symétrique et définie positive. Ce dernier système provient de l'approximation variationnelle d'une "équation intégrale" sur $\Gamma$. Au N° 5 on décrit diverses méthodes itératives de résolution de $(P_o)$, dont celle du gradient conjugué. L'approche qui a été retenue dans ce travail semble fournir un cadre naturel, pour résoudre sur des domaines $\Omega$ assez généraux, le problème $(P_o)$ comme un système couplé de problèmes de Dirichlet pour $-\Delta$. Elle contient en particulier plusieurs méthodes connues de résolution de $(P_o)$ sur des rectangles par des méthodes de différences finies ; en effet le classique schéma à 13 points de discrétisation de $(P_o)$ est un cas particulier (en fait le plus

<footnote>(1) $\psi$: fonction de courant, $\omega$: fonction de tourbillon.</footnote>
simple de tous) de l'approximation par éléments finis mixte utilisée ci-après
(Cf. N° 4.7, Remarque 4.7. et aussi GLOWINSKI [1], GLOWINSKI-LIONS-TREMOLIERES
[1, Ch. 4], GLOWINSKI-PIRONNEAU [3], etc...).

La plupart des applications des méthodes décrites dans ce travail sont relatives
à l'analyse numérique d'écoulement de fluides, visqueux, incompressibles, modélisés
par les équations de Navier-Stokes ; on a, en particulier, en vue l'application
à des problèmes de circulation dans des lacs ou des océans.

Le lecteur essentiellement intéressé par les applications peut, en première lecture, aller directement au N° 3.

## 2. LE PROBLEME CONTINU.

### 2.1. Cadre fonctionnel et notations.

Les espaces fonctionnels ci-dessous sont essentiels dans l'étude de (P₀) :

$$H^2(\Omega) = \{v | v \in L^2(\Omega), \frac{\partial v}{\partial x_i} \in L^2(\Omega), \frac{\partial^2 v}{\partial x_i \partial x_j} \in L^2(\Omega), 1 \le i, j \le 2\},$$

$$V = H^2(\Omega) \cap H_0^1(\Omega) = \{v | v \in H^2(\Omega), v=0 \text{ sur } \Gamma\},$$

$$H_0^2(\Omega) = \overline{D(\Omega)} H^2(\Omega) = \{v | v \in H^2(\Omega), v = \frac{\partial v}{\partial n} = 0 \text{ sur } \Gamma\}.$$

L'espace $H^2(\Omega)$ est un espace de Hilbert pour le produit scalaire

$$(u,v)_{H^2(\Omega)} = (u,v)_{L^2(\Omega)} + \sum_{1 \le i \le 2} \left( \frac{\partial u}{\partial x_i}, \frac{\partial v}{\partial x_i} \right)_{L^2(\Omega)} + \sum_{1 \le i,j \le 2} \left( \frac{\partial^2 u}{\partial x_i \partial x_j}, \frac{\partial^2 v}{\partial x_i \partial x_j} \right)_{L^2(\Omega)}.$$

Par ailleurs, $\Omega$ étant borné et de frontière régulière on a la

**Proposition 2.1** : L'application $v \to \|\Delta v\|_{L^2(\Omega)}$ définit sur $V$ une norme équivalente à celle induite par $H^2(\Omega)$. ■

Les espaces ci-dessous sont également essentiels dans l'étude de (P₀) :

$$H(\Omega; \Delta) = \{v | v \in L^2(\Omega), \Delta v \in L^2(\Omega)\},$$

$$\mathcal{H} = \{v | v \in H(\Omega; \Delta), \Delta v = 0\}.$$

L'espace $H(\Omega;\Delta)$ est un espace de Hilbert pour le produit scalaire
$$ (u,v)_{H(\Omega;\Delta)} = (u,v)_{L^2(\Omega)} + (\Delta u, \Delta v)_{L^2(\Omega)} $$

et pour la norme associée

$$ (2.1) \qquad \|v\|_{H(\Omega;\Delta)} = (\|v\|_{L^2(\Omega)}^2 + \|\Delta v\|_{L^2(\Omega)}^2)^{1/2}. $$

On déduit de (2.1) la

**Proposition 2.2** : *Sur H les topologies induites par H(Ω;Δ) et L²(Ω) sont identiques.*

## 2.2. Propriétés des traces.

On note $\gamma_0$ et $\gamma_1$ les deux applications traces définies par

$$ \gamma_0 v = v|_{\Gamma}, \quad \gamma_1 v = \frac{\partial v}{\partial n}|_{\Gamma}. $$

Il résulte alors de LIONS-MAGENES [1] les

**Proposition 2.3** : *L'application $\{\gamma_0, \gamma_1\}$ est linéaire, continue, surjective de H(Ω;Δ) sur $H^{-1/2}(\Gamma) \times H^{-3/2}(\Gamma)$.*

**Proposition 2.4** : *L'application $\{\gamma_0, \gamma_1\}$ est linéaire, continue, surjective de $H^2(\Omega)$ sur $H^{3/2}(\Gamma) \times H^{1/2}(\Gamma)$.*

**Proposition 2.5** : *L'application $\gamma_1$ est surjective de $V$ sur $H^{1/2}(\Gamma)$.*

**Proposition 2.6** : *L'application $\gamma_0$ restreinte à H est un isomorphisme (algébrique et topologique) de H sur $H^{-1/2}(\Gamma)$.*

## 2.3. Formule de Green.

On note $\langle \cdot, \cdot \rangle$ (resp. $\ll \cdot, \cdot \gg$) la forme bilinéaire de la dualité entre $H^{1/2}(\Gamma)$ et $H^{-1/2}(\Gamma)$ (resp. $H^{3/2}(\Gamma)$ et $H^{-3/2}(\Gamma)$) qui prolonge $(\cdot, \cdot)_{L^2(\Omega)}$ i.e. $\langle v, w \rangle = \int_{\Gamma} v w d\Gamma$ $\forall v \in H^{1/2}(\Gamma)$, $\forall w \in L^2(\Gamma)$ (resp. $\ll v, w \gg = \int_{\Gamma} v w d\Gamma$ $\forall v \in H^{3/2}(\Gamma)$, $\forall w \in L^2(\Gamma)$). On a alors (cf. LIONS-MAGENES, loc. cit.) la formule de Green suivante
$$ (2.2) \quad \left\{ \begin{array}{l} \displaystyle \int_{\Omega} \Delta u v dx - \int_{\Omega} u \Delta v dx = \langle \gamma_1 u, \gamma_0 v \rangle - \langle \gamma_0 u, \gamma_1 v \rangle \\[1em] \forall u \in H^2(\Omega), \forall v \in H(\Omega; \Delta). \end{array} \right. $$

## 2.4. Résultats d'existence, d'unicité et de décomposition pour $(P_0)$.

Si on suppose que dans (1.1) on a

$$ (2.3) \quad f \in L^2(\Omega), g_1 \in H^{3/2}(\Gamma), g_2 \in H^{1/2}(\Gamma), $$

il résulte alors de LIONS-MAGENES [1] le

**Théorème 2.1** : Le problème $(P_0)$ admet une solution unique dans $H^2(\Omega)$. ■
On démontre aisément la

**Proposition 2.7** : $(P_0)$ équivaut à

$$ (2.4) \quad \left\{ \begin{array}{l} - \Delta \omega = f, \\ - \Delta \psi = \omega, \\ \gamma_0 \psi = g_1, \quad \gamma_1 \psi = g_2. \end{array} \right. \quad \blacksquare $$

**Remarque 2.1** : La décomposition (2.4) est bien connue en Mécanique des Fluides où

$\omega$ est la fonction de **tourbillon**,

$\psi$ est la fonction de **courant**. ■

Dans la suite la **trace** de $\omega$ sur $\Gamma$ jouera un rôle fondamental, théoriquement et numériquement.

L'existence de $\gamma_0 \omega$ résulte, via le N° 2.2, de la

**Proposition 2.8** : Si $f,g_1,g_2$ vérifient (2.3), $\omega$ admet sur $\Gamma$ une **trace** $\gamma_0 \omega \in H^{-1/2}(\Gamma)$. ■

## 2.5. Relation entre $\gamma_0 \omega$ et $\gamma_1 \psi$.

Un certain nombre de méthode itératives de résolution de $(P_0)$ (cf. CIARLET-GLOWINSKI [1], BOSSAVIT [1], GLOWINSKI [1], GLOWINSKI-PIRONNEAU [3], le N° 5 de cet article, etc...), ainsi que la méthode quasi-directe du N° 4, sont en fait basées sur les résultats de ce N° 2.5.
Le lemme qui suit est essentiel et se démontrerait à partir des résultats des N° 2.1, 2.2, 2.3, 2.4 :

**Lemme 2.1** : Soit $\lambda \in H^{-1/2}(\Gamma)$ ; alors :

(i) Le problème

$$ (2.5) \qquad \left\{ \begin{array}{l} \Delta^2 \psi = 0, \\ \psi|_{\Gamma} = 0, \\ - \Delta \psi|_{\Gamma} = \lambda \end{array} \right. $$

admet une solution unique dans $V = H^2(\Omega) \cap H_0^1(\Omega)$.

(ii) Soit $\psi$ la solution dans $V$ de (2.5) ; alors l'opérateur linéaire $A$ défini par

$$ (2.6) \qquad A\lambda = - \gamma_1 \psi $$

est un isomorphisme de $H^{-1/2}(\Gamma)$ dans $H^{1/2}(\Gamma)$.

(iii) La forme bilinéaire $a(\cdot, \cdot) : H^{-1/2}(\Gamma) \times H^{-1/2}(\Gamma) \to \mathbb{R}$ définie par

$$ (2.7) \qquad a(\lambda, \mu) = \langle A\lambda, \mu \rangle $$

est continue, symétrique et $H^{-1/2}(\Gamma)$-elliptique. ■

Considérons de nouveau $(P_o)$ ; on a vu précédemment (cf. Théorème 2.1) que si $f \in L^2(\Omega)$, $g_1 \in H^{3/2}(\Gamma)$, $g_2 \in H^{1/2}(\Gamma)$ alors $(P_o)$ admet une solution unique dans $H^2(\Omega)$. De plus $\omega = - \Delta\psi$ admet une trace $\lambda = \gamma_0\omega$, dans $H^{-1/2}(\Gamma)$. On va montrer que $\lambda$ est solution d'une équation variationnelle linéaire dans $H^{-1/2}(\Gamma)$.

Soit $\bar{\psi}$ l'unique solution dans $V$ de

$$ (2.8) \qquad \left\{ \begin{array}{l} \Delta^2 \bar{\psi} = 0, \\ \bar{\psi}|_{\Gamma} = 0, \\ - \Delta \bar{\psi}|_{\Gamma} = \lambda; \end{array} \right. $$
le problème (2.8) équivaut à

$$
(2.9) \quad \left\{
\begin{array}{l}
- \Delta \bar{\omega} = 0, \\
\bar{\omega}|_{\Gamma} = \lambda, \\
- \Delta \bar{\psi} = \bar{\omega}, \\
\bar{\psi}|_{\Gamma} = 0.
\end{array}
\right.
$$

De même, soit $\psi_o$ l'unique solution dans $H^2(\Omega)$ de

$$
(2.10) \quad \left\{
\begin{array}{l}
\Delta^2 \psi_o = f, \\
\psi_o|_{\Gamma} = g_1, \\
- \Delta \psi_o|_{\Gamma} = 0;
\end{array}
\right.
$$

le problème (2.10) équivaut à

$$
(2.11) \quad \left\{
\begin{array}{l}
- \Delta \omega_o = f, \\
\omega_o|_{\Gamma} = 0, \\
- \Delta \psi_o = \omega_o \\
\psi_o|_{\Gamma} = g_1.
\end{array}
\right.
$$

Si $\psi$ est la solution de $(P_o)$, et si $\omega = -\Delta\psi$, on a évidemment

$$
\begin{align*}
\psi &= \psi_o + \bar{\psi}, \\
\omega &= \omega_o + \bar{\omega},
\end{align*}
$$

et on remarque que _le calcul de ψ_o demande la résolution de deux problèmes de_
_Dirichlet pour -Δ_._ On a un résultat analogue pour ψ̄, dès que λ est connu.

Le *résultat fondamental* de ce N° est énoncé dans le

**Théorème 2.2** : *Soit $\psi$ la solution de $(P_o)$ ; la trace $\lambda$ de $-\Delta\psi$ sur $\Gamma$ est l'unique solution du problème variationnel linéaire*

$$
(2.12) \quad \left\{
\begin{array}{ll}
\langle A\lambda, \mu \rangle = \left\langle \frac{\partial \psi_o}{\partial n} - g_2, \mu \right\rangle & \quad \forall \mu \in H^{-1/2}(\Gamma), \\
\\
\lambda \in H^{-1/2}(\Gamma). & \quad \blacksquare
\end{array}
\right.
$$
**Remarque 2.2** : La forme $a(\cdot, \cdot)$ étant symétrique, (2.12) équivaut au problème de minimisation

$$
(2.13) \quad \left\{
\begin{array}{l}
J(\lambda) \le J(\mu) \quad \forall \mu \in H^{-1/2}(\Gamma), \\
\lambda \in H^{-1/2}(\Gamma),
\end{array}
\right.
$$

où, dans (2.13)

$$
J(\mu) = \frac{1}{2} \langle A\mu, \mu \rangle - \left\langle \frac{\partial \psi_o}{\partial n}, g_2, \mu \right\rangle . \blacksquare
$$

**Remarque 2.3** : Considérant la condition aux limites $\frac{\partial \psi}{\partial n}|_{\Gamma} = g_2$ comme une contrainte linéaire on peut associer à $(P_o)$ le Lagrangien $\mathcal{L}: H^2(\Omega) \times H^{-1/2}(\Gamma) \to \mathbb{R}$ défini par

$$
\mathcal{L}(v, \mu) = \frac{1}{2} \int_{\Omega} |\Delta v|^2 dx - \int_{\Omega} fvdx + \langle \frac{\partial v}{\partial n} - g_2, \mu \rangle .
$$

Si $\tilde{V} = \{v \in H^2(\Omega), v = g_1 \text{ sur } \Gamma\}$ on peut montrer que $\{\psi, -\gamma_0 \Delta \psi\}^{(2)}$ est l'unique
point-selle de $\mathcal{L}$ sur $V \times H^{-1/2}(\Gamma)$ et que

$$
J(\mu) = - \min_{v \in \tilde{V}} \mathcal{L}(v, \mu).
$$

Il en résulte que (2.13) est le **problème dual de** $(P_o)$ associé au Lagrangien $\mathcal{L}$.
On renvoie à CIARLET-GLOWINSKI [1], GLOWINSKI [1], GLOWINSKI-LIONS-TREMOLIERES
[1, Ch. 4] pour une étude plus détaillée de $(P_o)$, par des méthodes de dualité
associées à des Lagrangiens de même type que $\mathcal{L}$.
■

**Remarque 2.4** : Les données $f$ et $g_1$ interviennent dans (2.12) via $\psi_o$ (cf. (2.10)).

**Remarque 2.5** : Soit $\tilde{\mu}$ un relèvement de $\mu$ dans $\Omega$. Il résulte alors de la **formule de Green** que

$$
(2.14) \quad \left\{
\begin{aligned}
a(\lambda, \mu) &= \langle A\lambda, \mu \rangle = - \left\langle \frac{\partial \psi}{\partial n}, \mu \right\rangle = \\
&= - \int_{\Omega} \Delta\psi \tilde{\mu} dx - \int_{\Omega} \nabla\psi \cdot \nabla \tilde{\mu} dx = \int_{\Omega} \omega \tilde{\mu} dx - \int_{\Omega} \nabla\psi \cdot \nabla \tilde{\mu} dx,
\end{aligned}
\right.
$$

<footnote>(2) $\psi$ est bien entendu la solution de $(P_o)$.</footnote>
où, dans (2.14), $\psi$ est la solution de (2.5) et où $\omega = -\Delta\psi$. On a de même

$$ (2.15) \quad \left\langle \frac{\partial \psi_o}{\partial n}, \mu \right\rangle = \int_{\Omega} \nabla \psi_o \cdot \nabla \tilde{\mu} dx + \int_{\Omega} \Delta \psi_o \tilde{\mu} dx = \int_{\Omega} \nabla \psi_o \cdot \nabla \tilde{\mu} dx - \int_{\Omega} \omega_o \tilde{\mu} dx. $$

où $\{\omega_o, \psi_o\}$ est la solution de (2.11).

Il faut remarquer que (2.14), (2.15) ne sont que formelles, en général, $\tilde{\mu}$ n'étant pas assez régulière. Cependant si $\mu$ est suffisamment régulière ($\mu \in H^{1/2}(\Gamma)$ par exemple) on peut prendre $\tilde{\mu} \in H^1(\Omega)$ et (2.14), (2.15) sont vraies.

L'intérêt de (2.14), (2.15) est de fournir des expressions des deux membres de (2.12) dans lesquelles $\frac{\partial \psi}{\partial n}$ et $\frac{\partial \psi_o}{\partial n}$ n'apparaissent pas explicitement. On tirera parti de cette remarque aux N° 3 et 4, lorsque $(P_o)$ et (2.12) seront approchés à l'aide d'une approximation par éléments finis mixte.

## 2.6. Récapitulation sur la résolution de $(P_o)$.

Soit $\psi$ la solution de

$$ (P_o) \quad \left\{ \begin{array}{l} \Delta^2 \psi = f \text{ dans } \Omega, \\ \psi|_{\Gamma} = g_1, \\ \frac{\partial \psi}{\partial n}|_{\Gamma} = g_2, \end{array} \right. $$

et $\omega = -\Delta\psi$, $\lambda = \omega|_{\Gamma}$.

On a vu, au N° 2.5, que résoudre $(P_o)$ équivaut à résoudre séquentiellement les problèmes suivants :

$$ (2.16) \quad \left\{ \begin{array}{l} -\Delta\omega_o = f \text{ dans } \Omega, \\ \omega_o|_{\Gamma} = 0, \end{array} \right. $$

$$ (2.17) \quad \left\{ \begin{array}{l} -\Delta\psi_o = \omega_o \text{ dans } \Omega, \\ \psi_o|_{\Gamma} = g_1, \end{array} \right. $$

$$ (2.18) \quad A\lambda = \frac{\partial \psi_o}{\partial n} - g_2, $$
$$
(2.19) \quad \left\{
\begin{array}{l}
- \Delta \omega = f \text{ dans } \Omega, \\
\omega|_{\Gamma} = \lambda,
\end{array}
\right.
$$

$$
(2.20) \quad \left\{
\begin{array}{l}
-\Delta\psi = \omega \text{ dans } \Omega, \\
\psi|_{\Gamma} = g_1.
\end{array}
\right.
$$

Il faut donc résoudre **quatre problèmes de Dirichlet pour $-\Delta$ et une "équation intégrale" sur $\Gamma$** dont la formulation variationnelle est donnée par (2.12). Dans les N° qui suivent on s'intéressera plus particulièrement à l'approximation de (2.18).

**Remarque 2.6** : L'opérateur A est en fait un **opérateur pseudo-différentiel qui n'est pas explicitement connu, en général.**

2.7. **Forme explicite de l'isomorphisme A dans le cas d'un domaine circulaire.**
Les résultats de ce N° ne sont pas essentiels pour la suite et sont surtout donnés
à titre de curiosité. On suppose dans ce N° que

$$
(2.21) \quad \Omega = \{x \in \mathbb{R}^2, x_1^2 + x_2^2 < R^2\},
$$

on a alors le

**Théorème 2.3** : *Soit A l'isomorphisme de $H^{-1/2}(\Gamma)$ sur $H^{1/2}(\Gamma)$ défini au N° 2.5. Alors si $\Omega$ est donné par (2.21) et si $\lambda$ est une fonction définie sur $\Gamma$ et suffisamment régulière ($\lambda \in L^2(\Gamma)$ par exemple), on a*

$$
(A\lambda)(x) = \int_{\Gamma} A(x,y)\lambda(y)d\Gamma(y), \forall x \in \Gamma,
$$

le noyau $A(x,y)$ étant donné par

$$
\left\{
\begin{array}{l}
A(x,y) = \frac{1}{2\pi} \left[ \left(1 - \frac{|y-x|}{2R^2}\right) \text{Log} \frac{R}{|y-x|} + \frac{|y-x|}{R} \sqrt{1 - \frac{|y-x|}{4R^2}^2} \dots \right. \\
\qquad \left. \dots \text{Arc cos} \frac{|y-x|}{2R} - \frac{1}{2} \right],
\end{array}
\right.
$$

avec $|y-x| = \text{distance}(x,y).$ ■
### 3. APPROXIMATION DE $(P_o)$ PAR UNE METHODE D'EQUIMENTS FINIS MIXTES.

On suppose dans ce N° que $\Omega$ est un domaine polygonal de $\mathbb{R}^2$ mais les résultats ci-après peuvent être généralisés au cas d'une frontière $\Gamma$ courbe, par utilisation d'éléments finis de type isoparamétrique (cf. CIARLET-RAVIART [2]).

#### 3.1. Triangulation de $\Omega$. Espaces fondamentaux.

Soit $\mathcal{T}_h$ une triangulation de $\Omega$ qui vérifie

$$ (3.1) \qquad \mathcal{T}_h \text{ est finie, } T \subset \bar{\Omega}, \forall T \in \mathcal{T}_h, \bigcup_{T \in \mathcal{T}_h} T = \bar{\Omega}, $$

$$ (3.2) \qquad \left\{ \begin{array}{l} T, T' \in \mathcal{T}_h, T \neq T' \Rightarrow \mathring{T} \cap \mathring{T}' = \emptyset \text{ et } T \cap T' = \emptyset \text{ ou } T \text{ et } T' \\ \text{n'ont qu'un sommet commun ou tout un côté commun.} \end{array} \right. $$

On désignera par $h$ la longueur du plus grand côté des $T \in \mathcal{T}_h$. Soit $P_k$ l'espace des polynômes de degré $\le k$ ; on considère alors les espaces de dimension finie suivants

$$ (3.3) \qquad V_h = \{v_h \in C^0(\bar{\Omega}), v_h|_T \in P_k, \forall T \in \mathcal{T}_h\} $$

$$ (3.4) \qquad V_{oh} = V_h \cap H_o^1(\Omega) = \{v_h \in V_h, v_h = 0 \text{ sur } \Gamma\} $$

$$ (3.5) \qquad M_h \subset V_h \text{ tel que } V_h = V_{oh} \oplus M_h, $$

$$ (3.6) \qquad \left\{ \begin{aligned} W_{gh} &= \{(v_h, q_h) \in V_h \times V_h, v_h|_{\Gamma} = g_{1h}, \int_{\Omega} \nabla v_h \cdot \nabla u_h dx = \\ &= \int_{\Omega} q_h u_h dx + \int_{\Gamma} g_{2h} u_h d\Gamma \quad \forall u_h \in V_h\}. \end{aligned} \right. $$

Dans (3.5), $M_h$ n'est pas défini de façon unique, mais un choix naturel pour $M_h$ sera donné ultérieurement. Dans (3.6), $g_{1h}$ est une approximation de $g_1$ appartenant à $\gamma_o V_h$ et $g_{2h}$ est une approximation de $g_2$ telle que $\int_{\Gamma} g_{2h} u_h d\Gamma$ soit "facile" à calculer.
### 3.2. Approximation de ($P_o$).

On approche ($P_o$) par

$$ (P_h) \quad \left\{ \begin{array}{ll} j_h(\psi_h, \omega_h) \le j_h(v_h, q_h) & \forall (v_h, q_h) \in W_{gh}, \\ (\psi_h, \omega_h) \in W_{gh}, & \end{array} \right. $$

où

$$ (3.7) \quad j_h(v_h, q_h) = \frac{1}{2} \int_{\Omega} |q_h|^2 dx - \int_{\Omega} f_h v_h dx. $$

Dans (3.7), $f_h$ est une approximation de f telle que $\int_{\Omega} f_h v_h dx$ soit facile à calculer.

Une telle approximation par éléments finis est dite **mixte** (cf. CIARLET-RAVIART [1], CIARLET-GLOWINSKI [1]). On a la

**Proposition 3.1** : ($P_h$) admet une solution unique.

### 3.3. Résultats de convergence.

On suppose que les angles de $\mathcal{C}_h$ sont **bornés inférieurement**, **uniformément en h**, par $\theta_o > 0$, et que

$$ (3.8) \quad \left\{ \begin{array}{ll} \displaystyle \frac{\max_{T \in \mathcal{T}_h} h(T)}{\min_{T \in \mathcal{T}_h} h(T)} & \le \tau, \forall \mathcal{T}_h, \tau \text{ indépendant de } h, \\ \\ \displaystyle \frac{\max_{T \in \mathcal{T}_h} h(T)}{\min_{T \in \mathcal{T}_h} h(T)} & \ge \tau, \forall \mathcal{T}_h, \tau \text{ indépendant de } h, \end{array} \right. $$

où, dans (3.8), $h(T)$ est la longueur du plus grand côté de T. Alors si $k \ge 2$ il résulte de CIARLET-RAVIART [1] que

$$ (3.9) \quad \|\psi_h^{-\psi}\|_{H^1(\Omega)}^2 + \|\omega_h^{-(\Delta\psi)}\|_{L^2(\Omega)}^2 \le C \|\psi\|_{H^{k+2}(\Omega)}^{k-1} $$

avec C indépendant de h et $\psi$.$^{(3)}$

<footnote>(3) Cette estimation suppose que $f,g_1,g_2$ ont été correctement approchés.</footnote>
Pour le cas $k = 1$ on renvoie à GLOWINSKI [1], GLOWINSKI-LIONS-TREMOLIERES [1,ch. 4]. On renvoie également à SCHOLZ [1] où, sous les hypothèses ci-dessus sur $\mathcal{C}_h$, il est démontré que pour $k \ge 3$ on a

$$ \|\psi_h - \psi\|_{L^2(\Omega)}^2 + h^2 \|\omega_h - (-\Delta\psi)\|_{L^2(\Omega)}^2 \le C h^{k+1} \|\psi\|_{H^{k+1}(\Omega)}^2 $$

avec $C$ indépendant de $h$ et $\psi$.

**Remarque 3.1** : Les résultats ci-dessus valent encore si au lieu d'éléments finis *triangulaires* on utilise des éléments finis *quadrilatéraux*.

### 3.4. Décomposition de $(P_h)$.

On a $V_h = V_{oh} \oplus M_h$ ; soit $\{\psi_h, \omega_h\}$ la solution de $(P_h)$ et $\lambda_h$ la composante de $\omega_h$ dans $M_h$. On a alors

$$ (3.10) \quad \left\{ \begin{array}{l} \omega_h = (\omega_h - \lambda_h) + \lambda_h, \quad \text{avec} \\ \\ \omega_h - \lambda_h \in V_{oh}, \quad \lambda_h \in M_h. \end{array} \right. $$

Dans CIARLET-GLOWINSKI [1], GLOWINSKI-PIRONNEAU [3] est démontré le

**Théorème 3.1** : Soit $\{\psi_h, \omega_h\}$ la solution de $(P_h)$ et $\lambda_h$ la composante de $\omega_h$ dans $M_h$ ; on peut alors caractériser $\{\psi_h, \omega_h, \lambda_h\}$ comme étant l'unique élément de $V_h \times V_h \times M_h$ tel que

$$ (3.11) \quad \left\{ \begin{array}{l} \displaystyle \int_{\Omega} \nabla \omega_h \cdot \nabla v_h dx = \int_{\Omega} f_h v_h dx \quad \forall v_h \in V_{oh}, \\ \\ \omega_h - \lambda_h \in V_{oh}, \end{array} \right. $$

$$ (3.12) \quad \left\{ \begin{array}{l} \displaystyle \int_{\Omega} \nabla \psi_h \cdot \nabla v_h dx = \int_{\Omega} \omega_h v_h dx \quad \forall v_h \in V_{oh}, \\ \\ \psi_h \in V_h, \quad \psi_h|_{\Gamma} = g_{1h}, \end{array} \right. $$

$$ (3.13) \quad \int_{\Omega} \nabla \psi_h \cdot \nabla \mu_h dx = \int_{\Omega} \omega_h \mu_h dx + \int_{\Gamma} g_{2h} \mu_h d\Gamma \quad \forall \mu_h \in M_h. \blacksquare $$
**Remarque 3.2** : Les relations (3.11)-(3.13) sont les analogues en dimension finie de (2.19), (2.20) et de

$$ (3.14) \qquad \left\{ \begin{array}{l} \displaystyle \int_{\Omega} \nabla \psi \cdot \nabla \mu dx = \int_{\Omega} \omega \mu dx + \int_{\Gamma} g_2 \mu d\Gamma \quad \forall \mu \in M, \text{ avec} \\[2ex] \qquad M \subset H^1(\Omega), \quad H^1(\Omega) = H_o^1(\Omega) \oplus M. \quad \blacksquare \end{array} \right. $$

### 3.5. Analogue discret du Lemme 2.1.

Soit $\lambda_h \in M_h$ ; on définit alors $\omega_h$ et $\psi_h$ comme étant les solutions des deux problèmes de Dirichlet approchés ci-dessous

$$ (3.15) \qquad \left\{ \begin{array}{l} \displaystyle \int_{\Omega} \nabla \omega_h \cdot \nabla v_h dx = 0 \quad \forall v_h \in V_{oh}, \\[2ex] \omega_h \in V_h, \quad \omega_h - \lambda_h \in V_{oh}. \end{array} \right. $$ 
$$ (3.16) \qquad \left\{ \begin{array}{l} \displaystyle \int_{\Omega} \nabla \psi_h \cdot \nabla v_h dx = \int_{\Omega} \omega_h v_h dx \quad \forall v_h \in V_{oh}, \\[2ex] \psi_h \in V_{oh}. \end{array} \right. $$

On définit alors la forme bilinéaire $a_h : M_h \times M_h \to \mathbb{R}$ par

$$ (3.17) \qquad a_h(\lambda_h, \mu_h) = \int_{\Omega} \omega_h \mu_h dx - \int_{\Omega} \nabla \psi_h \cdot \nabla \mu_h dx \quad \forall \mu_h \in M_h. $$

On remarque que la définition de $a_h(\cdot, \cdot)$ suit le point de vue de la Remarque 2.5. On démontre aisément le

**Lemme 3.1.** : La forme bilinéaire $a_h(\cdot, \cdot)$ est symétrique et définie positive ; on a en outre

$$ (3.18) \qquad a_h(\lambda_{1h}, \lambda_{2h}) = \int_{\Omega} \omega_{1h} \omega_{2h} dx \quad \forall \lambda_{1h}, \lambda_{2h} \in M_h $$

où $\omega_{jh}$ ($j=1,2$) est la solution de (3.15) correspondant à $\lambda_{jh}$.

### 3.6. Application du Lemme 3.1 à la résolution de $(P_h)$.

Soit $\{\psi_h, \omega_h\}$ la solution de $(P_h)$ et $\lambda_h$ la composante de $\omega_h$ dans $M_h$.

Soit $\bar{\omega}_h$ et $\bar{\psi}_h$ les solutions respectives de
$$
(3.19) \quad \left\{
\begin{array}{l}
\displaystyle \int_{\Omega} \nabla \bar{\omega}_h \cdot \nabla v_h \, dx = 0 \quad \forall v_h \in V_{oh}, \\
\bar{\omega}_h \in V_h, \quad \bar{\omega}_h - \lambda_h \in V_{oh},
\end{array}
\right.
$$

$$
(3.20) \quad \left\{
\begin{array}{l}
\displaystyle \int_{\Omega} \nabla \bar{\psi}_h \cdot \nabla v_h \, dx = \int_{\Omega} \bar{\omega}_h v_h \, dx \quad \forall v_h \in V_{oh}, \\
\bar{\psi}_h \in V_{oh}.
\end{array}
\right.
$$

Soit $\omega_{oh}$ et $\psi_{oh}$ les solutions respectives de

$$
(3.21) \quad \left\{
\begin{array}{l}
\displaystyle \int_{\Omega} \nabla \omega_{oh} \cdot \nabla v_h \, dx = \int_{\Omega} f_h v_h \, dx \quad \forall v_h \in V_{oh}, \\
\omega_{oh} \in V_{oh},
\end{array}
\right.
$$

$$
(3.22) \quad \left\{
\begin{array}{l}
\displaystyle \int_{\Omega} \nabla \psi_{oh} \cdot \nabla v_h \, dx = \int_{\Omega} \omega_{oh} v_h \, dx \quad \forall v_h \in V_{oh}, \\
\psi_{oh} \in V_h, \quad \psi_{oh} = g_{1h} \quad \text{sur } \Gamma.
\end{array}
\right.
$$

On a $\psi_h = \bar{\psi}_h + \psi_{oh}$, $\omega_h = \bar{\omega}_h + \omega_{oh}$ et les relations (3.19)-(3.22) sont les analogues en dimension finie de (2.8)-(2.11).

L'analogue discret du Théorème 2.2. est le

**Théorème 3.2** : Soit {$\psi_h, \omega_h$} la solution de $(P_h)$ et $\lambda_h$ la composante de $\omega_h$ dans $M_h$. Alors $\lambda_h$ est l'unique solution du problème variationnel linéaire

$$
(3.23) \quad \left\{
\begin{array}{l}
a_h(\lambda_h, \mu_h) = \int_{\Omega} \nabla \psi_{oh} \cdot \nabla \mu_h \, dx - \int_{\Omega} \omega_{oh} \mu_h \, dx - \int_{\Gamma} g_{2h} \mu_h \, d\Gamma \quad \forall \mu_h \in M_h, \\
\lambda_h \in M_h
\end{array}
\right.
$$

qui équivaut à un système linéaire de matrice symétrique, définie positive.

**Remarque 3.3** : Pour calculer le second membre de (3.23) il faut résoudre les
deux problèmes de Dirichlet approchés (3.21), (3.22). De même, une fois que $\lambda_h$
est connu il faut résoudre les deux problèmes de Dirichlet approchés (3.11), (3.12)
pour obtenir $\omega_h$ et $\psi_h$.
### 3.7. Sur le conditionnement de $a_h(\cdot, \cdot)$.

Dans les applications pratiques le système linéaire équivalent à (3.23) sera résolu par des méthodes directes ou itératives. Pour cette raison il est important de connaître le conditionnement de la matrice de ce système (cette matrice sera explicitée au N° 4). Le Théorème 3.3 ci-après fournira une estimation du nombre de conditionnement de cette matrice (cf. N° 4.6).

**Théorème 3.3** : On suppose que $\Omega$ est convexe ; sous les hypothèses sur $\mathcal{C}_h$ du N° 3.3 on a alors pour $k \ge 1$ et $h$ suffisamment petit

$$ (3.24) \qquad \alpha h \|\gamma_0 \lambda_h\|_{L^2(\Gamma)}^2 \le a_h(\lambda_h, \lambda_h) \le \beta \|\gamma_0 \lambda_h\|_{L^2(\Gamma)}^2 \quad \forall \lambda_h \in M_h, \\ \text{avec } \alpha, \beta > 0 \text{ et indépendants de } h \text{ et } \lambda_h. \blacksquare $$

**Remarque 3.4** : On peut déduire du Théorème ci-dessus

$$ (3.25) \qquad \lim_{h \to 0} \sup_{\lambda_h \in M_h - \{0\}} \frac{a_h(\lambda_h, \lambda_h)}{\|\gamma_0 \lambda_h\|_{L^2(\Gamma)}^2} = \sup_{\lambda \in L^2(\Gamma) - \{0\}} \frac{a(\lambda, \lambda)}{\|\lambda\|_{L^2(\Gamma)}^2} = |A| $$

où $|A|$ est la plus grande valeur propre de $A$. ■

### 3.8. Récapitulation sur la résolution du problème approché.

Soit $\{\psi_h, \omega_h\}$ la solution de $(P_h)$ et $\lambda_h$ la composante de $\omega_h$ dans $M_h$. Le vecteur $\lambda_h$ est solution d'un **système linéaire de matrice symétrique et définie positive**. Une formulation variationnelle équivalente est donnée par (3.23) mais la forme bilinéaire $a_h(\cdot, \cdot)$ n'est pas explicitement connue. Il en est malheureusement de même pour la matrice du système linéaire ci-dessus. La construction de cette matrice et les méthodes de résolution du système linéaire correspondant seront décrites au N° 4 ci-après. Des méthodes itératives de résolution évitant la construction de la matrice seront décrites au N° 5.
## 4. CONSTRUCTION ET RESOLUTION DU SYSTEME LINEAIRE EQUIVALENT A (3.23).

### 4.1. Généralités.

Soit $N_h = \dim(M_h)$ et $\mathcal{B}_h = \{w_i\}_{i=1}^N$ une base de $M_h$. Si $\lambda_h \in M_h$ on a

$$ (4.1) \qquad \lambda_h = \sum_{j=1}^{N_h} \lambda_j w_j. $$

On démontrerait alors la

**Proposition 4.1.** : Le problème (3.23) est équivalent au système linéaire en $\{\lambda_1, \lambda_2, \dots, \lambda_{N_h}\}$

$$ (E_h) \qquad \left\{ \begin{array}{l} \displaystyle \sum_{j=1}^{N_h} a_h(w_j w_i)^\lambda_j = \int_{\Omega} \nabla \psi_{oh} \cdot \nabla w_i \,dx - \int_{\Omega} \omega_{oh} w_i \,dx - \int_{\Omega} g_{2h} w_i \,d\Gamma, \\ 1 \le i \le N_h. \quad \blacksquare \end{array} \right. $$

Notons $a_{ij} = a_h(w_j w_i)$ et $A_h = (a_{ij})_{1 \le i,j \le N_h}$ ; on démontrerait alors aisément la

**Proposition 4.2.** : La matrice $A_h$ est une matrice $N_h \times N_h$ symétrique et définie positive. ■

On va étudier dans les N° suivants la construction de $A_h$ et du second membre de $(E_h)$, à partir d'un choix convenable de $M_h$ et $\mathcal{B}_h$.

### 4.2. Choix de $M_h$

Il convient de choisir $M_h$ de telle façon que, toutes choses égales par ailleurs, le calcul des $a_{ij}$ et des seconds membres soit peu coûteux. Ce sera le cas si les fonctions $w_i \in \mathcal{B}_h$ ont un "petit" support. Il résulte de CIARLET-GLOWINSKI [1], GLOWINSKI-LIONS-TREMOLIERES [1, Ch. 4] qu'un choix optimal pour $M_h$ semble correspondre à

$$ (4.2.) \qquad \left\{ \begin{array}{l} M_h : \text{complémentaire de } V_{oh} \text{ dans } V_h, \\ v_h \in M_h \Rightarrow v_{h|T} = 0 \quad \forall T \in \mathcal{T}_h, T \cap \Gamma = \emptyset. \end{array} \right. $$

Si, en particulier, $V_h$ est défini à partir d'éléments finis de Lagrange (cf. Figure 4.1., ci-dessous pour $k=2$), $M_h$ est l'espace des fonctions de $V_h$ nulles
sur les noeuds de $\mathcal{T}_h$ n'appartenant pas à $\Gamma$. On aura dans ce cas

$$N_h = \dim(M_h) = \mathrm{Card} (\Sigma_h),$$

où

$$\Sigma_h = \{P \mid P \in \Gamma, P \text{ noeud de } \mathcal{T}_h\}$$

et le choix canonique pour $\mathcal{B}_h$ correspondra à

$$\mathcal{B}_h = \{w_i\}_{i=1}^{N_h}$$

où la fonction $w_i$ est définie par

$$(4.3) \qquad \left\{ \begin{array}{l} w_i \in V_h, \\ w_i(P_i) = 1, P_i \in \Sigma_h, w_i(Q) = 0 \ \forall Q \text{ noeud de } \mathcal{T}_h, Q \neq P_i. \end{array} \right.$$

Dans (4.3) on suppose que $\Sigma_h$ a été ordonné de 1 à $N_h$. Avec ce choix de $M_h$ et $\mathcal{B}_h$ les coefficients $\lambda_j$ dans le développement (4.1) de $\lambda_h$ sont précisément les valeurs prises par $\lambda_h$ aux noeuds frontières $P_j$, $1 \le j \le N_h$ d'où

$$\lambda_j = \lambda_h(P_j) \quad \forall P_j \in \Sigma_h, \ 1 \le j \le N_h.$$

<figure><img src="image_3.png" /><figcaption>Figure 4.1.</figcaption></figure>

(k=2 ; on a entouré les noeuds frontières et hachuré les supports de $w_4$ et $w_j$).
### 4.3. Calcul des seconds membres de $(E_h)$.

Soit $b_h = \{b_1, b_2, \dots, b_{N_h}\}$ le vecteur second membre de $(E_h)$. On a donc

$$ (4.5) \qquad b_i = \int_{\Omega} \nabla \psi_{oh} \cdot \nabla w_i \, dx - \int_{\Omega} \omega_{oh} w_i \, dx - \int_{\Gamma} g_{2h} w_i \, dx, \quad 1 \le i \le N_h. $$

Le calcul de $b_h$ nécessite donc la détermination de $\omega_{oh}$ et $\psi_{oh}$ par résolution de (3.37), (3.38).

**Remarque 4.1** : Le calcul du vecteur $b_h$ sera d'autant plus rapide que le support de $w_i$ est petit (cf. N° 4.2). Par ailleurs si $M_h$ vérifie (4.2), il suffit pour calculer $b_h$ de connaître $\omega_{oh}$ et $\psi_{oh}$ sur la réunion des $T \in \mathcal{T}_h$ tels que $T \cap \Gamma \neq \emptyset$, ce qui permet éventuellement, de réduire l'occupation mémoire de l'ordinateur utilisé. ■

### 4.4. Calcul de la matrice $A_h$.

Soit $w_j \in \mathcal{B}_h$ ; omettant l'indice h on note $w_j$ et $\psi_j$ les solutions de (3.15), (3.16) correspondant à $w_j$, soit

$$ (4.6) \qquad \left\{ \begin{array}{l} \displaystyle \int_{\Omega} \nabla w_j \cdot \nabla \mu_h \, dx = 0 \quad \forall \mu_h \in V_{oh}, \\ \omega_j \in V_h, \quad \omega_j - w_j \in V_{oh}, \end{array} \right. $$

$$ (4.7) \qquad \left\{ \begin{array}{l} \displaystyle \int_{\Omega} \nabla \psi_j \cdot \nabla \mu_h \, dx = \int_{\Omega} \omega_j \mu_h \, dx \quad \forall \mu_h \in V_{oh}, \\ \psi_j \in V_{oh}. \end{array} \right. $$

Compte tenu de (3.17) on en déduit

$$ (4.8) \qquad \left\{ \begin{array}{l} a_{ij} = a_h(w_j, w_i) = \int_{\Omega} \omega_j w_i \, dx - \int_{\Omega} \nabla \psi_j \cdot \nabla w_i \, dx, \\ 1 \le i, j \le N_h. \end{array} \right. $$

Il résulte de (4.8) que pour calculer la $j$-colonne de $A_h$ il faut résoudre les deux problèmes de Dirichlet approchés (4.6), (4.7), puis $w_i$ décrivant $(\mathcal{B}_h)_h$, évaluer les intégrales figurant dans (4.8). La Remarque 4.1. vaut encore à propos du calcul des $a_{ij}$. Il convient également de noter que compte tenu de la symétrie de $A_h$, il suffit, lorsque l'on construit la $j$-colonne, de se limiter au calcul des éléments $a_{ij}$ tels que $1 \le j \le i$. ■
On verra au N° 4.5.2. l'intérêt que l'on peut tirer des propriétés de $A_h$ mentionnées ci-dessus lorsque l'on résoud $(E_h)$ par la Méthode de Cholesky.

**Remarque 4.2.** : Il résulte de (3.18) que

$$ (4.9) \qquad a_{ij} = \int_{\Omega} \omega_i \omega_j dx \quad \forall 1 \le i,j \le N_h. $$

Il semble donc que l'on puisse calculer $A_h$ en résolvant seulement $N_h$ problèmes de Dirichlet approchés au lieu de $2N_h$ comme c'est le cas lorsque l'on emploie la formule (4.8). En fait cette simplification n'est qu'apparente. En effet, l'utilisation de (4.9) nécessite d'une part le stockage en mémoire de $\omega_1, \omega_2, \dots, \omega_{N_h}$, ce qui est toujours possible par utilisation de fichiers-bandes ou disques, mais cela aura pour effet d'augmenter assez considérablement le temps de calcul de $A_h$. D'autre part les intégrales intervenant dans (4.9) sont à effectuer sur $\Omega$ tout entier et non pas au voisinage de $\Gamma$ comme c'est le cas lorsque l'on utilise (4.8) (cf. Remarque 4.1.).

## 4.5. Résolution de $(E_h)$.

### 4.5.1. Généralités

Soit $\lambda_h \in \mathbb{R}^{N_h}$ le vecteur $\{\lambda_1, \lambda_2, \dots, \lambda_{N_h}\}$ ; $(E_h)$ s'écrit alors

$$ (4.10) \qquad A_h^{\lambda_h} = b_h. $$

La matrice $A_h$ étant **symétrique** et **définie positive** on peut utiliser pour résoudre (4.10) une **méthode directe** de type **Cholesky**. On peut également utiliser une **méthode itérative** du type S.O.R., S.S.O.R. (cf. VARGA [1], D.M. YOUNG [1]) ou du type **gradient** et **gradient conjugé** (cf. J.W. DANIEL [1], CEA [1], E. POLAK [1], CONCUS-GOLUB [1]). On détaillera au N° 4.5.2. la méthode de Cholesky qui semble particulièrement bien adaptée à la résolution de $(E_h)$. Quant aux méthodes du type **gradient** et **gradient conjugué** leur mise en oeuvre ne nécessite pas en fait, la connaissance de $A_h$. On reviendra là-dessus au N° 5.

### 4.5.2. Résolution de $(E_h)$ par la méthode de Cholesky.

La matrice $A_h$ étant **symétrique** et **définie positive**, il est classique qu'il existe une **matrice triangulaire inférieure** $L_h$ et une seule telle que
$$
(4.11) \quad \left\{
\begin{array}{l}
A_h = L_h L_h^t, \\
\ell_{ii} > 0, \ 1 \le i \le N_h,
\end{array}
\right.
$$

où les $\ell_{ii}$ sont les éléments de la **diagonale principale** de $L_h$. On notera donc
$\ell_{ij}$ les éléments de $L_h$, on a

$$
(4.12) \qquad \ell_{ij} = 0 \quad \forall 1 \le i < j \le N_h.
$$

Rappelons les formules de base de la méthode de Cholesky :

Pour $j=1$, on a :

$$
(4.13) \quad \left\{
\begin{array}{l}
l_{11} = \sqrt{a_{11}}, \\
l_{i1} = \frac{a_{i1}}{l_{11}} \quad \forall 2 \le i \le N_h.
\end{array}
\right.
$$

Pour $2 \le j \le N_h$, on a :

$$
(4.14) \quad \left\{
\begin{array}{@{}l@{\qquad}l@{}}
\ell_{jj} & \sqrt{a_{jj} - \sum_{k=1}^{j-1} \ell_{jk}^2}, \\
\ell_{ij} & \displaystyle\frac{a_{ij} - \sum_{k=1}^{j-1} \ell_{ik} \ell_{jk}}{\ell_{jj}}, \quad \forall j+1 \le i \le N_h.
\end{array}
\right.
$$

Compte-tenu des formules (4.13), (4.14) il apparaît que pour construire $L_h$ il est inutile de stocker $A_h$ (ou du moins les $a_{ij}$ avec $1 \le j \le i \le N_h$). En effet supposons connues les (j-1) premières colonnes de $L_h$ ; pour déterminer la $j^{\text{ê}}$ colonne de $L_h$, on détermine la solution $\{\omega_j, \psi_j\}$ de (4.6), (4.7), puis successivement, $a_{jj}$ par (4.8), $\ell_{ij}$ par (4.14), et pour $j+1 \le i \le N_h$, $a_{ij}$ par (4.8), $\ell_{ij}$ par (4.14). Les mêmes principes s'appliquent également à la construction de la première colonne de $L_h$ à partir de la première colonne de $A_h$. On notera également que pour déterminer la $j^{\text{ê}}$ colonne de $L_h$ il est inutile de stocker les $a_{ij}$ tels que $j \le i \le N_h$, puisque si les (j-1) premières colonnes de $L_h$ sont connues le calcul de $\ell_{ij}$ peut se faire dès que $a_{ij}$ est connu, à condition de commencer par calculer $a_{jj}, \ell_{jj}$ (cf. ci-dessus).■

Dès que $L_h$ est connu il est standard que la détermination de $\lambda_h$ se ramène à la résolution des deux systèmes triangulaires
$$ (4.15) \quad \left\{ \begin{array}{l} L_h y_h = b_h , \\ L_h^\lambda y_h = y_h . \end{array} \right. $$

La détermination de $\lambda_h$ à partir de $\lambda_h^\vee$ étant triviale on calcule ensuite $\omega_h$ et $\psi_h$ par la résolution des deux problèmes de Dirichlet approxérés (3.11), (3.12).■

**Remarque 4.3.** : Lorsque $L_h$ a été déterminée on peut facilement résoudre des problèmes $(E_h)$, correspondant à diverses valeurs de $f,g_1,g_2$, les autres paramètres restant inchangés. En effet dans la résolution d'un système linéaire par la méthode de Cholesky, la phase la plus coûteuse est la factorisation, i.e. le calcul de $L_h$ ; ceci est encore plus vrai dans le cas présent puisque pour calculer $A_h$, donc $L_h$, il faut résoudre $2N_h$ problèmes de Dirichlet approxérés. ■

### 4.5.3. Comptabilité des sous-problèmes linéaires dans la méthode de Cholesky.

La résolution de $(P_h)$, via $(E_h)$, par la méthode de Cholesky exige donc la résolution

* des deux problèmes de Dirichlet approxérés (3.21), (3.22) pour déterminer $b_h$,

* des $2N_h$ problèmes de Dirichlet approxérés (4.6), (4.7), $1 \le j \le N_h$ pour déterminer $L_h$,

* des deux systèmes linéaires (4.15) pour déterminer $\lambda_h$,

* des deux problèmes de Dirichlet approxérés (3.11), (3.12) pour déterminer $\omega_h$ et $\psi_h$.

Il faudra donc, en particulier, résoudre $2N_h+4$ problèmes de Dirichlet approxérés.

### 4.6. Conditionnement de $A_h$.

On rappelle que le nombre de conditionnement $v(M)$ d'une matrice carrée $N \times N$, inversible, $M$ est donné par

$$ (4.16) \quad v(M) = ||M|| \ ||M^{-1}|| , $$

où, dans (4.16), la norme matricielle considérée est celle associée à la norme vectorielle euclidienne canonique de $\mathbb{R}^N$. On rappelle également que si $M$ est symétrique et définie positive alors

$$ (4.17) \quad v(M) = \frac{\mu_{\max}}{\mu_{\min}} $$
où $\mu_{max}$ (resp. $\mu_{min}$) est la plus grande (resp. plus petite) valeur propre de M.

Le système $(E_h)$ sera d'autant plus facile à résoudre-que ce soit par des méthodes directes ou itératives - que le nombre de conditionnement de la matrice $A_h$ sera petit. Sous les hypothèses sur $\tau_h$ du N° 3.3. on déduit assez aisément du Théorème 3.3 et de (4.17) le

**Théorème 4.1.** : Dans le cas d'une approximation par éléments finis de Lagrange, sous les conditions de validité du Théorème 3.3 et pour $k \ge 1$ on a

$$ (4.18) \qquad v(A_h) = O\left(\frac{1}{h}\right) . \blacksquare $$

**Remarque 4.4** : On notera que les matrices d'approximation classique, par différences ou éléments finis, de $\Delta(\text{resp. } \Delta^2)$ ont un nombre de conditionnement en $O(\frac{1}{h^2})$ (resp. $O(\frac{1}{h^4})$) et sont donc nettement moins bien conditionnées que $A_h$. Cela étant ces matrices sont creuses alors que $A_h$ est une matrice pleine. ■

### 4.7. - Remarques diverses

**Remarque 4.5.** : Les $2N_h+4$ problèmes de Dirichlet approchés mis en évidence au N° 4.5. sont de la forme

$$ (4.19) \qquad (-\Delta)_h U_h = C_h $$

où $(-\Delta)_h$ est une matrice $N'_h \times N'_h$ **symétrique, définie positive** (approchant $-\Delta$), avec $N'_h = \text{dim}(V_{oh})$.

Les $2N_h+4$ systèmes (4.19) ne diffèrent que par les seconds membres et la matrice $(-\Delta)_h$ étant symétrique et définie positive, on peut effectuer une factorisation de Cholesky tenant compte de la structure **creuse** de $(-\Delta)_h$, soit

$$ (-\Delta)_h = \Lambda_h \Lambda_h^t . $$

La matrice $\Lambda_h$ est **triangulaire inférieure, inversible**.

La matrice $\Lambda_h$ étant calculée une fois pour toutes, la résolution des $2N_h+4$ systèmes linéaires (4.19) se réduit à celle (peu coûteuse) de systèmes linéaires dont les matrices sont $\Lambda_h$ et $\Lambda_h^t$. ■
**Remarque 4.6.** : Si on résoud (4.19) par une **méthode itérative**, il n'est pas déraisonnable, dans le calcul de $\{\omega_j, \psi_j\}$ à partir de la fonction de base $w_j$ (cf. N° 4.4), d'initialiser l'algorithme utilisé par $\omega_{j-1}$ (resp. $\psi_{j-1}$) si les supports de $w_{j-1}$ et $w_j$ sont "voisins" (on peut toujours stocker $\omega_{j-1}$ et $\psi_{j-1}$ dans les tableaux destinés aux solutions finales $\omega_h$ et $\psi_h$ de $(P_h)$).■

**Remarque 4.7.** : Tout ce qui a été dit ci-dessus reste valable si aux N° 3.1,3.2 on utilise des méthodes d'intégration numérique pour définir $W_{gh}$ et $(P_h)$. En particulier si $k=1$ et pour des triangulations particulières, l'approximation de $\int_{\Omega} q_h \mu_h dx$ par

$$ (4.20) \qquad \left\{ \begin{array}{ll} \frac{1}{3} \sum_{T \in \mathcal{T}_h} \text{Mesure}(T) \sum_{i=1}^{3} q_h(M_{iT}) \mu_h(M_{iT}) , & M_{iT}, 1 \le i \le 3, \text{ sommets de T,} \\ \\ \int_{\Omega} q_h \mu_h dx & \end{array} \right. $$

redonne le classique **schéma à 13 points** d'approximation par **différences finies** du problèmes de Dirichlet pour $\Delta^2$ (cf. GLOWINSKI [1], GLOWINSKI-LIONS-TREMOLIERES [1, ch. 4]).

## 5. REMARQUES SUR L'UTILISATION DE METHODES ITERATIVES. METHODES DU GRADIENT CONJUGUE.

L'utilisation de méthodes de gradient ou gradient conjugué permet de résoudre $(E_h)$ sans connaître explicitement $A_h$. Il faudra, en fait, résoudre "seulement" deux problèmes de Dirichlet approchés à chaque itération. Etant donné le nombre de pages limité dont nous disposerons, nous ne donnerons qu'une brève description de ces méthodes, renvoyant à GLOWINSKI-PIRONNEAU [3], BOURGAT-GLOWINSKI-PIRONNEAU [1], BOURGAT [1] pour une analyse plus détaillée. En vue de décrire ces algorithmes il est commode d'introduire $r_h : M_h \to \mathbb{R}^{N_h}$ défini comme suit :

Soit, comme au N° 4.1., $\mathcal{B}_h = \{w_i\}_{i=1}^{N_h}$ une base de $M_h$. Si $\mu_h \in M_h$ on a

$$ (5.1) \qquad \mu_h = \sum_{i=1}^{N_h} \mu_i w_i ; $$

$r_h$ est alors l'**isomorphisme** défini par

$$ (5.2) \qquad r_h \mu_h = \{\mu_1, \mu_2, \dots, \mu_{N_h}\} \quad \forall \mu_h \in M_h. $$
Si on note $(\cdot, \cdot)_h$ le produit scalaire euclidien canonique de $\mathbb{R}^{N_h}$ et par $\|\cdot\|_h$ la norme correspondante, on a

$$
(5.3) \qquad a_h(\lambda_h, \mu_h) = (A_h r_h \lambda_h, r_h \mu_h)_h \quad \forall \lambda_h, \mu_h \in M_h,
$$

$$
(5.4) \quad \int_{\Omega} \nabla \psi_{oh} \cdot \nabla \mu_h \, dx - \int_{\Omega} \omega_{oh} \mu_h \, dx - \int_{\Gamma} g_{2h} \mu_h \, d\Gamma = (b_h, r_h \mu_h)_h \quad \forall \mu_h \in M_h,
$$

où $A_h$ et $b_h$ ont été respectivement définis aux N° 4.1 et 4.3.

5.2. Méthodes de gradient à pas constant.

5.2.1. Description de l'algorithme.

Soit $s_h : M_h \times M_h \rightarrow \mathbb{R}$ une forme bilinéaire, symétrique et définie positive et soit $\rho > 0$. La méthode du gradient à pas constant est alors définie par

$$
(5.5) \qquad \lambda_h^0 \in M_h, \text{ arbitrairement choisi,}
$$

$\lambda_h^n$ connu, on calcule $\lambda_h^{n+1}$ par

$$
(5.6) \quad \left\{
\begin{array}{l}
\displaystyle \int_{\Omega} \nabla \omega_h^n \cdot \nabla v_h \, dx = \int_{\Omega} f_h v_h \, dx \quad \forall v_h \in V_{oh} \\[1em]
\omega_h^n - \lambda_h^n \in V_{oh},
\end{array}
\right.
$$

$$
(5.7) \quad \left\{
\begin{array}{l}
\displaystyle \int_{\Omega} \nabla \psi_h^n \cdot \nabla v_h \, dx = \int_{\Omega} \omega_h^n v_h \, dx \quad \forall v_h \in V_{oh}, \\[1em]
\psi_h^n \in V_h, \quad \psi_h^n = g_{1h} \quad \text{sur } \Gamma,
\end{array}
\right.
$$

$$
(5.8) \quad \left\{
\begin{array}{l}
s_h(\lambda_h^{n+1}, \mu_h) = s_h(\lambda_h^n, \mu_h) + \rho \left( \int_{\Omega} \nabla \psi_h^n \cdot \nabla \mu_h \, dx - \int_{\Omega} \omega_h^n \mu_h \, dx - \int_{\Gamma} g_{2h} \mu_h \, d\Gamma \right) \\
\forall \mu_h \in M_h, \quad \lambda_h^{n+1} \in M_h.
\end{array}
\right.
$$

La forme $s_h(\cdot, \cdot)$ étant symétrique et définie positive, il existe une matrice $S_h$
symétrique, définie positive, unique telle que

$$
(5.9) \qquad s_h(\lambda_h, \mu_h) = (S_h r_h \lambda_h, r_h \mu_h)_h .
$$

On en déduit que (5.8) peut encore s'écrire

$$
(5.10) \quad r_h \lambda_h^{n+1} = r_h \lambda_h^n + \rho S_h^{-1} (A_h r_h \lambda_h^n - b_h).
$$
Il résulte de (5.10) que le calcul de $\lambda_h^{n+1}$ exigera, à chaque iteration, la
solution d'un système linéaire de matrice $S_h$.

5.2.2. Choix de $s_h$.

Le choix de $s_h(\cdot, \cdot)$ est un point important et on renvoie à GLOWINSKI-PIRONNEAU [3],
BOURGAT-GLOWINSKI-PIRONNEAU [1] pour une discussion détaillée de ce problème.

Si on suppose que $M_h$ est défini par (4.2) et si on utilise des **éléments finis de Lagrange**, on peut prendre pour $s_h(\cdot, \cdot)$ (cf. CIARLET-GLOWINSKI [1], GLOWINSKI-LIONS-TREMOLIERES [1, ch. 4])

$$
(5.11) \qquad s_h(\lambda_h, \mu_h) = \int_{\Gamma} \lambda_h \mu_h d\Gamma,
$$

$$
(5.12) \qquad s_h(\lambda_h, \mu_h) = \int_{\Omega} \lambda_h \mu_h \, dx,
$$

$$
(5.13) \qquad s_h(\lambda_h, \mu_h) = \int_{\Omega} \nabla \lambda_h \cdot \nabla \mu_h \, dx.
$$

Les matrices associées à (5.11), (5.12), (5.13) sont creuses mais non diagonales.
Par utilisation d'un procédé d'intégration numérique il est facile d'approcher
les formes (5.11), (5.12) par des formes bilinéaires pour lesquelles $S_h$ est
diagonale (cf. GLOWINSKI-PIRONNEAU [3], LIONS-GLOWINSKI-TREMOLIERES, loc. cit.).

5.2.3. <u>Convergence de l'algorithme de gradient</u> (5.5)-(5.8).

On peut démontrer le

Théorème 5.1. : Soit {$\lambda_h^n$}$_{n}$ la suite définie par (5.5)-(5.8) et $\lambda_h$ la solution de ($E_h$). Alors $\forall \lambda_h^o \in M_h$

$$
\lim_{n \to +\infty} \lambda_h^n = \lambda_h
$$

<u>si et seulement si</u>

$$
0 < \rho < \frac{2}{\Lambda_{N_h}},
$$

où $\Lambda_{N_h}$ est la plus grande valeur propre de $S_h^{-1} \Lambda_{N_h} \cdot \blacksquare$
### 5.3. Méthodes de gradient à pas variable.

Lors de l'utilisation de l'algorithme (5.5)-(5.8) le choix optimal de $\rho$ est un pro-
blème important et non trivial à résoudre. Dans le but d'éviter ce type de diffi-
cultés on peut utiliser, au prix de certains calculs supplémentaires, les deux
variantes ci-dessous de l'algorithme (5.5)-(5.8) :

* Méthode de descente optimale :

$$ (5.14) \qquad \lambda_h^O \in M_h, \text{ arbitrairement choisi,} $$

$$ (5.15) \qquad g_h^O = A_h r_h \lambda_h^O - b_h, $$

puis pour $n \ge 0$

$$ (5.16) \qquad \rho_n = \frac{(S_h^{-1} g_h^n, g_h^n)}{(A_h S_h^{-1} g_h^n, S_h^{-1} g_h^n)}, $$

$$ (5.17) \qquad r_h \lambda_h^{n+1} = r_h \lambda_h^n - \rho_n S_h^{-1} g_h^n, $$

$$ (5.18) \qquad g_h^{n+1} = g_h^n - \rho_n A_h S_h^{-1} g_h^n. \blacksquare $$

* Méthode du résidu minimal :

C'est le même algorithme que ci-dessus, à l'exception de (5.16) qu'il faut
remplacer par

$$ (5.19) \qquad \rho_n = \frac{(A_h S_h^{-1} g_h^n, S_h^{-1} g_h^n)}{(S_h^{-1} A_h S_h^{-1} g_h^n, A_h S_h^{-1} g_h^n)}, \blacksquare $$

Dans ces deux algorithmes, l'étape (5.18) exigera, si $A_h$ n'est pas connue, la
résolution de deux problèmes de Dirichlet approchés.

**Remarque 5.1.** : Par utilisation des méthodes de (par exemple) MARCHOUK-KUZNETSOV [1], il serait intéressant d'estimer le taux de convergence des deux algorithmes ci-dessus. ■
### 5.4. Résolution de $(E_h)$ par une méthode de gradient conjugué.

La matrice $A_h$ étant symétrique et définie positive il est naturel de songer à la méthode du gradient conjugué pour résoudre $(E_h)$. On rappelle que cette méthode converge surlinéairement (i.e. plus rapidement que toute suite géométrique) et qu'en l'absence d'erreurs d'arrondis la convergence a lieu en un nombre fini d'itérations.

Dans GLOWINSKI-PIRONNEAU [3] on montre que l'algorithme du gradient conjugué appliqué à la résolution de $(E_h)$ prend la forme ci-dessous :

$$ (5.20) \qquad \lambda_h^0 \in M_h, \text{ choisi arbitrairement,} $$

$$ (5.21) \qquad g_h^0 = A_h r_h \lambda_h^0 - b_h, $$

$$ (5.22) \qquad z_h^0 = g_h^0, $$

et pour $n \ge 0$

$$ (5.23) \qquad \rho_n = \frac{(z_h^n, g_h^n)_h}{(A_h z_h^n, z_h^n)_h}, $$

$$ (5.24) \qquad r_h \lambda_h^{n+1} = r_h \lambda_h^n - \rho_n z_h^n, $$

$$ (5.25) \qquad g_h^{n+1} = g_h^n - \rho_n A_h z_h^n, $$

$$ (5.26) \qquad \gamma_n = \frac{(g_h^{n+1}, g_h^{n+1})_h}{(g_h^n, g_h^n)_h}, $$

$$ (5.27) \qquad z_h^{n+1} = g_h^{n+1} + \gamma_n z_h^n. \blacksquare $$

On remarque que pour $n \ge 0$ et si on suppose que $\lambda_h^n, g_h^n, z_h^n$ sont connus, il faut résoudre deux problèmes de Dirichlet approchés pour calculer $A_h z_h^n$. Une fois ce vecteur connu on calcule $\rho_n, \lambda_h^{n+1}, g_h^{n+1}$, puis $\gamma_n$ et $z_h^{n+1}$. ■
Remarque 5.2 : Pour $N_h$ suffisamment grand la convergence surlinéaire est plus intéressante que la convergence théorique en un nombre fini d'itérations. ■

## 6. Commentaires.

Il résulte de nos expériences numériques que les deux méthodes de résolution les plus efficaces sont

(i) La **méthode du gradient conjugué** décrite au N° 5, si $(P_h)$ ne doit être résolu qu'une seule fois.

(ii) La méthode "quasi-directe" du N° 4, si l'on désire un "solver" biharmonique devant être utilisé un grand nombre de fois. Ce sera le cas lorsque l'on résoudra par une méthode itérative les équations de Navier-Stokes dans la formulation {$\psi, \omega$}, ou encore un problème de Stokes (ou Navier Stokes) dépendant du temps, puisque dans ce dernier cas les méthodes décrites au N° 4.5 peuvent être facilement généralisées (cf. BOURGAT-GLOWINSKI-PIRONNEAU [3], GLOWINSKI-PIRONNEAU [4]). ■

Pour conclure nous voudrions signaler qu'il existe dans la littérature un certain nombre d'articles relatifs à la résolution du problème de Dirichlet pour $\Delta^2$, via une suite de couples de problèmes de Dirichlet pour $-\Delta$. Citons entre autres, J. SMITH [1], [2], [3], BOSSAVIT [1], ERLICH [1], [2], [3], Mc LAURIN [1], ERLICH-GUPTA [1]. Ces travaux sont relatifs à des approximations, par **différences finies sur des rectangles** et n'utilisent pas le fait que le problème approché équivaut à résoudre un système linéaire de matrice **symétrique**, définie positive dont la solution est la **trace de la fonction de tourbillon approchée**. Il résulte de ceci, qu'à notre connaissance, la méthode "quasi-directe" du N° 4, les méthodes de gradient à pas variable et de gradient conjugué du N° 5 semblent être nouvelles.

Il est bien évident que l'outil de base pour obtenir ces algorithmes était l'approximation par éléments finis mixtes du N° 3.

Des applications de la méthode de gradient à pas constant du N° 5.2. sont donnés dans BOURGAT [1], à partir de l'article de CIARLET-GLOWINSKI [1].
BIBLIOGRAPHY.

BOSSAVIT, A. [1] Une méthode de décomposition de l'opérateur biharmonique. Note HI 585/2, Electricité de France, (1971).

BOURGAT, J.F. [1] - Numerical study of a dual iterative method for solving a finite element approximation of the biharmonic equation, Rapport LABORIA N° 156, et à paraître dans Comp. Meth. Applied Mech. Eng.

BOURGAT J.F., GLOWINSKI R., PIRONNEAU O. [1] - Numerical methods for the Dirichlet problem for the biharmonic equation and applications (à paraître)

CEA J. [1] Optimisation. Théorie et Algorithmes, Dunod, 1971.

CIARLET P.G., GLOWINSKI R. [1] - Dual iterative techniques for solving a finite element approximation of the biharmonic equation, Comp. Meth. Applied Mech. Eng. 5, (1975), pp. 277-295.

CIARLET P.G., RAVIART P.A. [1] - A mixed finite element method for the biharmonic equation, in Mathematical aspects of finite elements in partial differential equations. C. de Boor, Ed. Acad. Press, (1974), pp. 125-145.

[2] - Interpolation theory over curved element with application to finite element methods. Comp. Meth. Applied Mech. Eng. 1, (1972), pp. 217-249.

CONCUS P., GOLUB G.H. [1] - Monography on conjugate gradient (à paraître).

DANIEL J.W. [1] - The approximate minimization of functionals. Prentice Hall (1970).

EHRLICH L.W. [1] - Solving the biharmonic equation as coupled difference equations Siam J. Num. Anal. 8 (1971), pp. 278-287.

[2] - Coupled harmonic equations, SOR and Chebyshev acceleration, Math. Comp. 26, (1972), pp. 335-343.

[3] - Solving the biharmonic equation in a square. Comm. ACM, 16 (1973), pp. 711-714.

EHRLICH L.W., GUPTA M.M. [1] - Some difference schemes for the biharmonic equation Siam J. Num. Anal. 12, (1975), pp. 773-790.

GLOWINSKI R. [1] - Approximations externes par éléments finis d'ordre un et deux du problème de Dirichlet pour $\Delta^2$. In Topics in Numerical Analysis J.J.H. Miller Ed., Academic Press, (1973), pp. 123-171.

[2] - Sur l'équation biharmonique dans un domaine multi-connexé C.R.A.S. Paris. (à paraître).

GLOWINSKI R., LIONS J.L., TREMOLIERES R. [1] Analyse Numérique des Inéquations Variationnelles (Tome 2), Dunod-Bordas, (1976).

GLOWINSKI R., PIRONNEAU O. [1] Sur la résolution numérique du problème de Dirichlet pour l'opérateur biharmonique par une méthode "quasi-directe". C.R.A.S. Paris, t. 282 A, pp. 223-226, (1976).
[2] Sur la résolution numérique du problème de Dirichlet pour $\Delta^2$ par la méthode du gradient conjugué. Applications. C.R.A.S. Paris. t. 282 A, p. 1315-1318 (1976).

[3] Sur la résolution par une méthode "quasi-directe", et par diverses méthodes itératives, d'une approximation par éléments finis mixtes du problème de Dirichlet pour $\Delta^2$. Rapport 76010, Laboratoire d'Analyse Numérique, Université Paris 6, (1976).

[4] Rapport de Stanford University (à paraître).

LIONS J.L., MAGENES E. [1] - *Problèmes aux limites non homogènes*, (T.1), Dunod, 1968.

MARCHOUK G.I., KUZNETSOV J.A., [1] - Méthodes itératives et fonctionnelles quadratiques, dans : *Sur les Méthodes Numériques en Sciences Physiques et Economiques*, LIONS J.L., MARCHOUK G.I., Ed., Dunod, 1974, p. 1-132.

Mc LAURIN J.W. [1] - A general coupled equation approach for solving the biharmonic boundary value problem, *Siam J. Num. Anal.* 11, (1974) pp. 14-33.

POLAK E. [1] *Computational Methods in Optimization*, Acad. Press, 1971.

SHOLZ R. [1] Approximation Von sattelpunkten mit finiten elementen (à paraître).

SMITH J. [1] The coupled equation approach to the numerical solution of the biharmonic equation by finite differences, I. *SIAM J. Num. Anal.* 5 (1968), pp. 323-339.

[2] The coupled equation approach to the numerical solution of the biharmonic equation by finite differences, II. *SIAM J. Num. Anal.* 7 (1970), pp. 104-111.

[3] On the approximate solution of the first boundary value problem for $\nabla^4 u = f$. *SIAM J. Num. Anal.*, 10 (1973), pp. 967-982.

VARGA R.S. [1] *Matrix iterative Analysis*, Prentice-Hall, 1962.

YOUNG D.M. [1] *Iterative solution of large linear systems*, Acad. Press, 1971.