IRIA
laboria

laboratoire de recherche
en informatique
et automatique

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

Numerical applications
of a dual iterative method
for solving
a finite element approximation
of the biharmonic equation

Jean-François
J.F. Bourgat

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
78150 - Le Chesnay
Tél. 954.90.20

Rapport de Recherche n°156

Janvier 1976
NUMERICAL APPLICATIONS OF A DUAL ITERATIVE METHOD
FOR SOLVING A FINITE ELEMENT APPROXIMATION OF THE
BIHARMONIC EQUATION

J.F. Bourgat

IRIA/LABORIA

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

**Abstract :**

The purpose of this paper is a presentation of numerical applications of an iterative method due to Ciarlet and Glowinski for solving a finite element approximation. The main feature of this method is that it reduces the biharmonic problem to a sequence of Dirichlet problems for the operator -Δ. Therefore, in numerical applications, finite element programs for solving second order problems can be used, and this is an interesting feature of the method.

**Résumé :**

On décrit dans cet article quelques applications numériques d'une méthode itérative due à Ciarlet et Glowinski pour résoudre une approximation par éléments finis du problème de Dirichlet pour l'opérateur biharmonique. La caractéristique de cette méthode est de ramener le problème biharmonique à une suite de problèmes de Dirichlet pour l'opérateur -Δ ce qui est intéressant dans les applications numériques puisqu'il suffit alors d'utiliser un programme d'éléments finis résolvant les problèmes du second ordre.
# I. INTRODUCTION.

We consider the homogeneous Dirichlet problem for the biharmonic operator :

$$ (1) \qquad \Delta^2 u = f \quad \text{in } \Omega, $$

$$ (2) \qquad u = \frac{\partial u}{\partial n} = 0 \quad \text{on } \Gamma, $$

where $\Omega$ denotes an open bounded set of $\mathbb{R}^2$ with boundary $\Gamma$, and $f$ is in the space $L^2(\Omega)$.

It is shown in CIARLET-RAVIART [1] that the pair $(u, -\Delta u)$ can be found in the product space $H_0^1(\Omega) \times L^2(\Omega)$. Then, generalizing a method due to GLOWINSKI [2], it is shown in CIARLET-GLOWINSKI [3] that one can use Uzawa's algorithm for solving problem (1)-(2). In this algorithm, one approximates a saddle-point of the Lagrangian of this problem, each iteration consisting in the solution of two Dirichlet problems for the operator $-\Delta$.

The pair $(u, -\Delta u)$ belonging to $H_0^1(\Omega) \times L^2(\Omega)$, we use in the discretization finite elements of class $C^0$ instead of finite elements of class $C^1$ used in the conforming methods, such as the 21-degree-of-freedom triangle of ARGYRIS [4]. The continuous problem and the discrete problem are studied in CIARLET-RAVIART [1] and the duality method for solving these problems is studied in CIARLET-GLOWINSKI [3].

In numerical applications, it suffices to use a finite element program for solving second order problems. We shall use here linear and quadratic elements and also second order isoparametric elements.

This method is well-suited for problems in fluid dynamics (Stokes problem) since it gives a continuous approximation not only of the stream function $u$, but also of the vorticity $\Delta u$, whereas standard
approximations using finite elements of class C¹ give a discontinuous approximation of the vorticity.

The Dirichlet problem for the operator $-\Delta$ are solved by Cholesky's method, which is appropriate in this case because it suffices to factorize the matrix once and for all and then to solve a sequence of linear systems whose matrices are triangular.

We have also used the over-relaxation method but only for a small number of iterations for each problem. The computing time is twice longer than for Cholesky's method, but the memory space needed in the computer is smaller. This is an important improvement for such fourth-order problems where a large number of variables is usually needed in order to have a good accuracy.

The method is tested on problems whose analytical solutions are known, so that it is possible to compute the error on the solution and its vorticity.

Computations were performed on a CII-IRIS 80.
## 2. THE CONTINUOUS PROBLEM.

Assuming $f \in L^2(\Omega)$, the standard formulation of problem (1)-(2) as a minimization problem is:

Find $u \in H_0^2(\Omega)$ such that

$$ (3) \qquad J(u) = \inf_{v \in H_0^2(\Omega)} J(v) $$

where

$$ (4) \qquad J(v) = \frac{1}{2} \int_{\Omega} |\Delta v|^2 dx - \int_{\Omega} fv dx $$

This problem has a unique solution $u$.

The pair $(u, -\Delta u)$ can be searched in the space $H_0^1(\Omega) \times L^2(\Omega)$ by minimizing the functional

$$ (5) \qquad J(v,q) = \frac{1}{2} \int_{\Omega} |q|^2 dx - \int_{\Omega} fv dx $$

over the set of functions $v \in H_0^2(\Omega)$ and $q \in L^2(\Omega)$ satisfying $-\Delta v = q$. This set can be described in a variational form as the space

$$ (6) \qquad v = \{ (v,q) \in H_0^1(\Omega) \times L^2(\Omega) \; ; \; \beta((v,q),\mu) = 0 \quad \forall \mu \in H^1(\Omega) \} $$

with

$$ \beta((v,q),\mu) = \int_{\Omega} \operatorname{grad} v \cdot \operatorname{grad} \mu \, dx - \int_{\Omega} q \mu \, dx. $$

We have the following result (CIARLET-RAVIART [1]).

**Theorem 1.** Let $u$ denote the solution of the problem (3)-(4). Then we also have:

$$ (7) \qquad J(u, -\Delta u) = \min_{(v,q) \in V} J(v,q). $$
### 3. THE DISCRETE PROBLEM.

In this section, we recall the main results about the finite element approximation of the problem (7) obtained in CIARLET-RAVIART [1] and in CIARLET-GLOWINSKI [3].

#### 3.1. Formulation of the discrete problem.

Let us suppose, for simplicity, that $\Omega$ is a polygon in the plane. We denote by $\mathcal{L}_h$ a triangulation of $\Omega$ with finite elements $K$ such that:

$$ (8) \qquad \Omega = \bigcup_{K \in \mathcal{L}_h} K. $$

The spaces $L^2(\Omega)$ and $H^1(\Omega)$ are approximated by the same space $V_h$ defined by:

$$ (9) \qquad V_h = \left\{ v_h \in C^0(\Omega) ; v_h/K \in P_k, \forall K \in \mathcal{L}_h! \right\}. $$

where $P_k$ denotes the set of all polynomials of degree $\le k$ in two variables.

The spaces $H^1(\Omega)$ and $V$ are respectively approximated by the spaces

$$ (10) \qquad V_{oh} = \left\{ v_h / v_h \in V_h , \ v_h = 0 \text{ on } \Gamma \right\} , $$

$$ (11) \qquad V_h = \left\{ (v_h, q_h) \in V_{oh} \times V_h ; \beta((v_h, q_h), \mu_h) = 0, \forall \mu_h \in V_h \right\}. $$

Then the discrete problem is defined by: Find a pair $(u_n, p_h) \in V_h$ such that:

$$ (12) \qquad J(u_h, p_h) = \min_{(v_h, q_h) \in V_h} J(v_h, q_h). $$
It was proved in [1] that this problem has a unique solution and that we have the following result of convergence :

Suppose that ($\mathcal{L}_h$) is a regular family of triangulations, i.e. for some constants $\alpha$ and $\beta$ independent of $h$, we have :

$$ (13) \qquad \max_{K \in \mathcal{L}_h} \frac{h(K)}{\rho(K)} \le \alpha, $$

$$ (14) \qquad \beta \max_{K \in \mathcal{L}_h} h(K) \le \min_{K \in \mathcal{L}_h} h(K), $$

where $h(K) = \text{diameter of } K$, $\rho(K) = \sup\{\text{diameters of inscribed spheres in } K\}$. Let $h = \max_{K \in \mathcal{L}_h} h(K)$, then

$$ (15) \qquad \|u-u_h\|_{H^1(\Omega)} + \|(-\Delta u)-p_h\|_{L^2(\Omega)} \le C\|u\|_{H^{k+2}(\Omega)} h^{k+1}, $$

where C is a constant independent of h and u.

**Remark 1.** If $k=1$, $(u_h, p_h)$ does not converge to $(u, -\Delta u)$ in general. However if the triangulation $h$ is generated by three families of parallel straight lines, it is proved in GLOWINSKI [3] that the convergence holds.

**Remark 2.** As a consequence of (15), if $u \in H^4(\Omega)$, an $O(h)$ convergence is obtained with polynomials of degree 2 only.

## 2.2. Solution of the discrete problem by a duality method.

For solving the problem (12) we use a classical technique for problems with linear constraints which amounts to approximating the saddle-point of a Lagrangian.
Let $\mathcal{M}_h$ be a given subspace of the space $V_h$ such that we may
write the direct sum.

$$
(16) \qquad V_h = V_{oh} \oplus \mathcal{M}_h.
$$

We define the space :

$$
(17) \quad W_h = \{ (v_h, q_h) \in V_{oh} \times V_h \ ; \ \beta((v_h, q_h), \mu_h) = 0 \ \forall \mu_h \in V_{oh} \} .
$$

and the Lagrangian

$$
(18) \qquad \mathcal{L}((v_h, q_h), \mu_h) = J(v_h, q_h) + \beta((v_h, q_h), \mu_h),
$$

then we have the following result (CIARLET-GLOWINSKI [3]) :

**Theorem 2.** Let $(u_h, p_h)$ denote the solution of problem (12) and let $p_{oh}$ be the unique function of $V_{oh}$ such that $p_h - p_{oh} \in \mathcal{M}_h$. Then the triple $(u_h, p_h, p_h - p_{oh})$ is a saddle-point of the Lagrangian $\mathcal{L}$ over the space $W_h \times \mathcal{M}_h$.

We next describe <u>Zawa's algorithm</u> for computing the saddle-point
((u<sub>h</sub>, p<sub>h</sub>), p<sub>h</sub> - p<sub>oh</sub>).

Let there be given $\lambda_h^1$ in $\mathcal{M}_h$.

Let us suppose $\lambda_h^n$ is known, then we compute successively $p_h^n, u_h^n$,
$\lambda_h^{n+1}$ as follows :

(i) Find the function $p_h^n \in V_h$ which satisfies

$$
p_h^n - \lambda_h^n \in V_{oh},
$$

$$
(19) \quad \int_{\Omega} \operatorname{grad} p_h^n \operatorname{grad} v_h \,dx = \int_{\Omega} f v_h \,dx, \quad \forall v_h \in V_{oh}.
$$
(ii) Find the function $u_h \in V_{oh}$ which satisfies

$$ (20) \qquad \int_{\Omega} \operatorname{grad} u_h \operatorname{grad} v_h \,dx = \int_{\Omega} p_h^n v_h \,dx, \quad \forall v_h \in V_{oh}. $$

(iii) Find the function $\lambda_h^{n+1} \in \mathcal{M}_h$ which satisfies

$$ (21) \qquad (\lambda_h^{n+1} \cdot \mu_h) \mathcal{M}_h = (\lambda_n^n \cdot \mu_n) \mathcal{M}_h + \rho \beta ((u_h^n, p_h^n), \mu_h), \quad \forall \mu_h \in \mathcal{M}_h, $$

where $(\cdot, \cdot) \in \mathcal{M}_h$ is any inner-product in the space $\mathcal{M}_h$.

Let us define the mapping $B_h : V_h \to \mathcal{M}_h$ as follows :

for any function $q_h \in V_h$, the function $B_h q_h$ satisfies

$$ (B_h q_h \cdot \mu_h) \mathcal{M}_h = \beta ((A_h q_h, q_h), \mu_h), \quad \forall \mu_h \in \mathcal{M}_h, $$

where $A_h q_h$ is the solution in $V_{oh}$ of the variational equations :

$$ \int_{\Omega} \operatorname{grad} v_h \operatorname{grad} \mu_h \,dx = \int_{\Omega} q_h \mu_h \,dx, \quad \forall \mu_h \in V_{ch}. $$

Then we have the following result (CIARLET-GLOWINSKI [3]) :

**Theorem 3** : Uzawa's algorithm converges, i.e. $u_h^n \to u_h$ in $V_{oh}$ and $p_h^n \to p_h$ in $V_h$, provided that

$$ (22) \qquad 0 < \rho < 2\sigma_h^2, $$

where

$$ (23) \qquad \sigma_h = \sup_{q_h \in V_n} \frac{\|q_h\|_{L^2(\Omega)}}{\|B_h q_h\|_{L^2(\Gamma)}}. $$

**Remark 3.** Uzawa's method converges for any $\mathcal{M}_h$ such that $V_h = V_{oh} \oplus \mathcal{M}_h$ and for any inner-product in the space $\mathcal{M}_h$. In the sequel we consider the following "canonical" choice for the space $\mathcal{M}_h$ : If $V_h$ is a space made up of piecewise polynomials defined by their values at nodes, $\mathcal{M}_h$ is defined as the space of functions in the space $V_h$ whose values are zero at the interior nodes.
Estimation of σ<sub>h</sub>.

It is difficult to obtain an estimation of σ<sub>h</sub> by using the formula (23), however we have the following "limit" result (CIARLET-GLOWINSKI, [3]).

**Theorem 4 :** With the canonical choice for the space $\eta_h$, assume that the inner-product (., .) $\eta_h$ is the inner-product of the space $L^2(\Gamma)$. Then we have :

$$
(24) \qquad \lim_{h \to 0} \sigma_h = \sigma,
$$

where

$$
(25) \qquad \sigma = \inf_{v \in H^2(\Omega) \cap H_o^1(\Omega)} \frac{\| \Delta v \|_{L^2(\Omega)}}{\| \frac{\partial v}{\partial n} \|_{L^2(\Gamma)}} .
$$

Thus the quantity σh may be estimated in practise, at least for
simple domains Ω . In particular we can find in SMITH [5] the
following results : for a rectangle with sides a and b ,
we have :

$$
(26) \qquad \sigma^2 \le \frac{\pi^2 (a^2 + b^2)^2}{4ab(a^3 + b^3)},
$$

and for a disk of radius r we have

$$
(27) \qquad \sigma^2 = \frac{2}{r} .
$$
## 4. NUMERICAL APPLICATIONS.

In the numerical applications, we shall take $k=1,2$, i.e. linear and quadratic elements.

### 4.1. Choice of the inner-products for the space $\mathcal{M}_h$. Numerical integration.

To use formula (24) we need the inner-product of $L^2(\Gamma)$, i.e.

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

however, we can use other equivalent inner-products for which formula (24) is still true (see GLOWINSKI-LIONS-TREMOLIERES, ch. 4).

Those inner-products are obtained by using numerical integration formulas. For example, if $S$ denotes a side of the finite element $K$, $S_1$ and $S_2$ the vertices (see fig.1), $\lambda_i$ the value of function $\lambda_h$ at the node $S_i$, we can take for $k=1$:

$$ (29) \qquad (\lambda_h, \mu_h) = \frac{1}{2} \sum_{S \subset \Gamma} \text{mes}(S) (\lambda_1 \mu_1 + \lambda_2 \mu_2) $$

and for $k=2$, if $S_3$ denotes the middle of the side $S$:

$$ (30) \qquad (\lambda_h, \mu_h) = \frac{1}{6} \sum_{S \subset I} \text{mes}(S) (\lambda_1 \mu_1 + 4\lambda_3 \mu_3 + \lambda_2 \mu_2). $$

<figure><img src="image_3.png" /><figcaption>Formulas (29), (30)</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Formula (31)</figcaption></figure>

Fig. 1. Nodes of numerical integration formulas.
With the inner-products (29) and (30), the linear system associated
with the problem (21) has a diagonal matrix, so that it is easy and
fast to solve. In the sequel we shall use the inner-products (29) and
(30).

$$
\text{For computing the terms } \int_{\Omega} f v_h \, dx, \int_{\Omega} p_h^n v_h \, dx, \int_{\Omega} p_h^n \mu_h \, dx
$$

in the case of quadratic elements, we use a numerical integration
formula which is exact for products of polynomials of degree $\le 2$,
whose nodes (see fig.1) are the three vertices and the three middles
of sides (see GLOWINSKI [2]):

$$
\int_T f v \, dx = \frac{\text{aire}(T)}{3} \left[ \begin{aligned}
& \frac{1}{10} (f_1 v_1 + f_2 v_2 + f_3 v_3) + \frac{8}{15} (f_4 v_4 + f_5 v_5 + f_6 v_6) - \\
& \quad - \frac{1}{60} (f_1 v_2 + f_1 v_3 + f_2 v_1 + f_2 v_3 + f_3 v_1 + f_3 v_2) + \\
& \quad + \frac{4}{15} (f_4 v_5 + f_4 v_6 + f_5 v_4 + f_5 v_6 + f_6 v_4 + f_6 v_5) - \\
& \quad - \frac{1}{15} (f_1 v_5 + f_2 v_6 + f_3 v_4 + f_4 v_3 + f_5 v_1 + f_6 v_2)
\end{aligned} \right]
$$

4.2. Methods for solving the second order problems.

We remark that solving the problems (19) and (20) for each iteration
of the Uzawa's algorithm is equivalent to solving a sequence of linear
systems whose matrices are identical.

Let us note A this matrix, we use the Cholesky's method to
factorize the matrix A, once and for all, in a product of two triangular
matrices :

$$
A = LL^t.
$$

To solve each subproblem (19) or (20) is then equivalent to solve
two triangular linear systems :

$$
Ly = b,
$$

$$
L^t x = y,
$$
where b denotes the right hand side of equation (19) or (20).

After a suitable renumbering of the nodes (see MERCIER [7]) the sparse matrix A will have a band structure with the smaller bandwidth possible. Then using programs for solving linear systems by Cholesky's method, adapted to band matrices, we reduce both the computing time and the memory space.

Though Cholesky's method is well-suited for solving problems (19) and (20), we have also used the over-relaxation method. We compute the optimal over-relaxation parameter for the matrix A and we start the algorithm for each subproblem with the solution of the corresponding subproblem at the preceding iteration of the Uzawa's algorithm. It is not necessary to solve each subproblem with a good accuracy, so we use a progressive test (see Section 3.3) to stop the iterations of the over-relaxation algorithm. In this fashion, we need about ten iterations for solving each sub-problem.

In the most favorable case, i.e. renumbering of the nodes and simple precision for the computation, the memory space for Cholesky's method is 1.5 larger than for the over-relaxation method, but its computing time is twice shorter.

### 4.3. Initialization and stopping rules of the algorithms.

Uzawa's algorithm is initialized with $\lambda^1 = 0$.

We denote by N the number of variables of each subproblem and we write

$$u^n = (u_1^n, \dots, u_N^n).$$

We define the mean relative difference between two iterates at the n-th iteration by

$$(32) \qquad \mathrm{DRM}(u^n) = \frac{\sum_{i=1}^{N} |u_i^n - u_{i-1}^{n-1}|}{\sum_{i=1}^{N} |u_i^n|}.$$

Uzawa's algorithm is stopped when

$$(33) \qquad \mathrm{DRM}(u^n) \le 10^{-4}.$$
For the overrelaxation algorithm, we define a stopping rule which is progressively finer. If $n$ denotes an iteration of the Uzawa's algorithm, and $p$ denotes an iteration of the overrelaxation algorithm, we define :

$$ (34) \qquad \mathrm{DRM}(u^n, p) = \frac{\sum_{i=1}^{N} |u_i^n, p - u_{i-1}^{n, p-1}|}{\sum_{i=1}^{N} |u_i^n, p|} , $$

and we stop the overrelaxation algorithm when

$$ (35) \qquad \mathrm{DRM}(u^{n+1}, p) \le \min(10^{-3}, \mathrm{DRM}(u^{n-1}), 10). $$

The overrelaxation algorithms are initialized with $u^{1,1} = 0$ and $u^{n+1} = u^{n-1}$.

### 4.4. Test problems.

We choose for numerical applications two problems, the analytical solutions of which are known.

**Test problem on the square. Problem (S).**

$$ \begin{align*} \Omega &= ]0,1[ \times ]0,1[ , \\ u &= 256x^2(1-x^2)y^2(1-y^2) , \\ \Delta u &= 512 [ x^2(1-x^2)(1-6y+6y^2) + y^2(1-y^2)(1-6x+6x^2) ] , \\ f &= 6144 [ x^2(1-x^2)+y^2(1-y^2) ] + 512(2-12x+12x^2)(2-12y+12y^2) . \end{align*} $$

**Test problem on the disk. Problem (D).**

$$ \begin{align*} \Omega &= \text{disk of radius } 0.5 \text{ and centered at the point } (0.5,0.5), \\ u &= 10(r^2 - 0.25)^2, \\ \Delta u &= 160 r^2 - 20 , \\ f &= 640. \end{align*} $$
### 4.5. Triangulations.

The triangulations chosen for the square satisfy the assumptions of Remark 1. They are generated by three families of parallels straight lines. We use the following triangulations (see table 1).

<table><caption>Table 1. Triangulations of the square.</caption><thead><tr><th rowspan="2">Triangulation</th><th rowspan="2">Number of triangles</th><th colspan="2">Linear elements</th><th colspan="2">Quadratic elements</th></tr><tr><th>interior nodes</th><th>boundary nodes</th><th>interior nodes</th><th>boundary nodes</th></tr></thead><tbody><tr><td>ℓ<sub>h1</sub></td><td>32</td><td>--</td><td>--</td><td>49</td><td>32</td></tr><tr><td>ℓ<sub>h2</sub></td><td>128</td><td>49</td><td>32</td><td>225</td><td>64</td></tr><tr><td>ℓ<sub>h3</sub></td><td>512</td><td>225</td><td>64</td><td>961</td><td>128</td></tr><tr><td>ℓ<sub>h4</sub></td><td>2048</td><td>961</td><td>128</td><td>--</td><td>--</td></tr></tbody></table>

The triangulation $ℓ_{h(i+1)}$ is obtained by joining the middles of the sides of all triangles of the triangulation $ℓ_{hi}$. Nodes for linear elements on triangulation $ℓ_{h(i+1)}$ and nodes for quadratic elements on triangulation $ℓ_{hi}$ are identical. The triangulation $ℓ_{h2}$ used with quadratic elements is shown on figure 2.

For the disk, we use the following triangulations (see Table 2).

<table><caption>Table 2. Triangulations of the disk.</caption><thead><tr><th rowspan="2">Triangulation</th><th rowspan="2">Number of triangles</th><th colspan="2">Linear elements</th><th colspan="2">Quadratic elements</th></tr><tr><th>interior nodes</th><th>boundary nodes</th><th>interior nodes</th><th>boundary nodes</th></tr></thead><tbody><tr><td>ℓ<sub>h1</sub></td><td>96</td><td>37</td><td>24</td><td>169</td><td>48</td></tr><tr><td>ℓ<sub>h2</sub></td><td>384</td><td>169</td><td>48</td><td>721</td><td>96</td></tr><tr><td>ℓ<sub>h3</sub></td><td>1536</td><td>721</td><td>96</td><td>--</td><td>--</td></tr></tbody></table>

The triangulation $ℓ_{h2}$ used with quadratic elements is shown on Fig.3.
<figure><img src="image_5.png" /><figcaption>Fig. 2. Triangulation $\mathcal{T}_{h2}$ of the square</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3 Triangulation $\mathcal{T}_{h2}$ of the disk.</figcaption></figure>
### 4.6. Error.

Let $(u, -\Delta u)$ denote the exact solution of a test problem and let $(u_h, p_h)$ denote the discrete solution obtained with the triangulation $\mathcal{T}_h$. Let us write:

$$u = (u_1, \dots, u_N) ; u_n = (u_{h1}, \dots, u_{hN}) ; -\Delta u = ((-\Delta u)_1, \dots, (-\Delta u)_N).$$

We define the mean relative error on $u$ by:

$$(36) \qquad \mathrm{ERM}(u) = \frac{\sum_{i=1}^{N} |u_i - u_{hi}|}{\sum_{i=1}^{N} |u_i|},$$

and on its Laplacian by:

$$(37) \qquad \mathrm{ERM}(-\Delta u) = \frac{\sum_{i=1}^{N} |(-\Delta u)_i - r_{hi}|}{\sum_{i=1}^{N} |(-\Delta u)_i|}.$$

We define the relative error in $L^2(\Omega)$ on $u$ by:

$$(38) \qquad \mathrm{ERL2}(u) = \frac{\|u - u_h\|_{L^2(\Omega)}}{\|u\|_{L^2(\Omega)}},$$

and on its Laplacian by:

$$(39) \qquad \mathrm{ERL2}(-\Delta u) = \frac{\| -\Delta u - p_h \|_{L^2(\Omega)}}{\| \Delta u \|_{L^2(\Omega)}}.$$

### 4.7. Results for the test problems.

In Tables 3 and 4 we have indicated the values of the exact solution and of the discrete solutions at the points shown on Fig. 4 along with the errors on $u$ and $\Delta u$.

<figure><img src="image_7.png" /><figcaption>Fig. 4: Points where exact and discrete solutions are compared.</figcaption></figure>
<table><caption>Table 4. Results of the test problem on the disk.</caption><thead><tr><th rowspan="2">Point</th><th rowspan="2">u</th><th colspan="3">Linear elements</th><th colspan="4">Quadratic elements</th></tr><tr><th>ε<sub>h1</sub></th><th>ε<sub>h2</sub></th><th>ε<sub>h3</sub></th><th colspan="2">Straights</th><th colspan="2">Isoparametric</th></tr><tr><th></th><th></th><th></th><th></th><th></th><th>ε<sub>h1</sub></th><th>ε<sub>h2</sub></th><th>ε<sub>h1</sub></th><th>ε<sub>h2</sub></th></tr></thead><tbody><tr><td>1</td><td>0.625</td><td>0.610</td><td>0.621</td><td>0.625</td><td>0.608</td><td>0.622</td><td>0.595</td><td>0.625</td></tr><tr><td>2</td><td>0.549</td><td>0.535</td><td>0.546</td><td>0.549</td><td>0.533</td><td>0.546</td><td>0.521</td><td>0.549</td></tr><tr><td>3</td><td>0.352</td><td>0.339</td><td>0.349</td><td>0.351</td><td>0.337</td><td>0.349</td><td>0.326</td><td>0.351</td></tr><tr><td>4</td><td>0.120</td><td>0.106</td><td>0.116</td><td>0.119</td><td>0.112</td><td>0.118</td><td>0.1</td><td>0.119</td></tr><tr><td>5</td><td>0.</td><td>0.</td><td>0.</td><td>0.</td><td>0.</td><td>0.</td><td>0.</td><td>0.</td></tr><tr><td colspan="9">- Δu</td></tr><tr><td>1</td><td>20.</td><td>20.46</td><td>20.12</td><td>20.04</td><td>19.86</td><td>19.95</td><td>19.72</td><td>20.</td></tr><tr><td>2</td><td>17.5</td><td>17.96</td><td>17.64</td><td>17.53</td><td>17.36</td><td>17.45</td><td>17.22</td><td>17.5</td></tr><tr><td>3</td><td>10.</td><td>10.36</td><td>10.13</td><td>10.04</td><td>9.85</td><td>9.95</td><td>9.71</td><td>10.</td></tr><tr><td>4</td><td>-2.5</td><td>-3.02</td><td>-2.54</td><td>-2.49</td><td>-2.32</td><td>-2.55</td><td>-2.53</td><td>-2.49</td></tr><tr><td>5</td><td>-20.</td><td>-25.38</td><td>-23.63</td><td>-21.38</td><td>-25.8</td><td>-19.29</td><td>-25.17</td><td>-21.07</td></tr><tr><td>ERM(u)</td><td></td><td>0.031</td><td>0.0086</td><td>0.0013</td><td>0.045</td><td>0.0075</td><td>0.087</td><td>0.0018</td></tr><tr><td>ERL2(u)</td><td></td><td>0.028</td><td>0.0077</td><td>0.0012</td><td>0.042</td><td>0.0067</td><td>0.075</td><td>0.0015</td></tr><tr><td>ERM(-Δu)</td><td></td><td>0.07</td><td>0.026</td><td>0.0056</td><td>0.14</td><td>0.025</td><td>0.17</td><td>0.013</td></tr><tr><td>ERL2(-Δu)</td><td></td><td>0.06</td><td>0.024</td><td>0.0065</td><td>0.11</td><td>0.024</td><td>0.13</td><td>0.016</td></tr></tbody></table>

### 4.8. Analysis of the results.

On the whole, the approximation of u is good over all the domain and the approximation of -Δu is good except near the boundary where a great number of nodes is necessary to get an acceptable accuracy.

Although we have reduced a fourth order problem to a sequence of second order problems, more variables are needed than in the case of a "reasonable" second order problem in order to get a comparable accuracy.

The numerical integration error has a great influence on the approximation of -Δu near the boundary. To estimate this influence, we consider
the problem (S) discretized with quadratic elements on the triangulation $\mathcal{T}_h$. First we integrate the terms $\int_{\Omega} f v_h dx$, $\int_{\Omega} p_h^n v_h dx$, $\int_{\Omega} p_h^n \mu_h dx$ with formula (31) and next we use the following second-order formula :

$$ (41) \qquad \int_{\Omega} f v dx = \frac{\text{aire}(T)}{3} \sum_{i=1}^{3} f_i v_i, $$

where $f_i$ and $v_i$ are the values of $f$ and $v$ at the mid-points of the sides. We obtain for the point 5 of Figure 4 the results shown on Table 5.

<table><caption>**Table 5.** Comparison of numerical integration formulas</caption><thead><tr><th rowspan="2">Exact value of - Δu</th><th colspan="2">Computed value of p<sub>h</sub> with formula (31)</th></tr><tr><th>formula (31)</th><th>formula (41)</th></tr></thead><tbody><tr><td>- 32</td><td>- 31</td><td>- 35.45</td></tr></tbody></table>

From the results of Tables 4 and 5, one cannot conclude that one type of finite element is better than another for our problem.

In the case the boundary is not approximated (problem (S)), quadratic finite elements give better results than linear finite elements except for the approximation of -Δu on the boundary. If the approximation of the boundary is not exact (problem (D)), linear and isoparametric finite elements give better results than straight quadratic elements.

We remark that discrete solutions converge to the exact solution even in the case of the disk discretized with linear elements, though the triangulation is not generated by three families of parallel straight lines. We can explain heuristically the convergence since it appears (see Fig. 3) that this triangulation could be split into two subtriangulations which have this property.

### 4.9. Convergence of Uzawa's algorithm. Choice of the parameter $\rho$.

From formulas (22), (24), (26), (27), we deduce that we have to take :

$$ (42) \qquad 0 < \rho < 2\pi, $$

for the problem on the square,
and,

$$ (43) \qquad 0 < \rho < 8, $$

for the problem on the disk.

We indicate on Fig. 5 the number of iterations of Uzawa's algorithm for the problem (D) with linear elements. We have taken the stopping rules

$$ \mathrm{DRM}(u^n) \le 10^{-4} \quad \text{and} \quad \mathrm{DRM}(u^n) \le 3.10^{-5}. $$

<figure><img src="image_8.png" /><figcaption>Fig. 5. Number of iterations of the Uzawa's algorithm in function of $\rho$.</figcaption></figure>
For values of $\rho$ greater than 5 there are oscillations in the convergence which are progressively larger as $\rho$ approaches 8.

So we have chosen $\rho = 5$ for the disk and $\rho = 4$ for the square, i.e. just above the mid-value.

Uzawa's algorithm converges rapidly for a proper choice of the parameter $\rho$. In this case, about twenty iterations are enough in practice as we infer from Table 6.

In this table, we study the behavior of the solution at point 1 of the disk in Figure 4 as well as the behavior of the error for the problem (D) with linear elements and triangulation $\ell_{h2}$.

<table><caption>Table 6. Convergence of Uzawa's algorithm.</caption><thead><tr><th>Iteration</th><th>2</th><th>4</th><th>6</th><th>10</th><th>20</th><th>50</th></tr></thead><tbody><tr><td>Computed value of u at the center</td><td>0.316</td><td>0.604</td><td>0.622</td><td>0.622</td><td>0.621</td><td>0.620</td></tr><tr><td>Computed value of p at the center</td><td>15.25</td><td>19.85</td><td>20.13</td><td>20.13</td><td>20.12</td><td>20.11</td></tr><tr><td>DRM(u<sup>n</sup>)</td><td>9.45</td><td>0.22</td><td>0.011</td><td>0.0004</td><td>0.0001</td><td>0.00001</td></tr><tr><td>ERM(u<sup>n</sup>)</td><td>2.84</td><td>0.16</td><td>0.03</td><td>0.01</td><td>0.009</td><td>0.01</td></tr></tbody></table>

It is important to use an appropriate stopping rule so as to take into account the fact that DRM(u^n) decreases very slowly after a small number of iterations. For example if we stop the iterations with $DRM(u^n) \le 10^{-5}$ instead of $DRM(u^n) \le 10^{-4}$ then we need thirty more iterations as it is shown in table 6. For the different examples treated on the square the numbers of iterations obtained are indicated in Table 7. The stopping rule is $DRM(u^n) \le 10^{-4}$.

<table><caption>Table 7. Number of iterations of Uzawa's algorithm for problems on the square.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Linear elements</th><th colspan="3">Quadratic elements</th></tr><tr><th>ℓ<sub>h2</sub></th><th>ℓ<sub>h3</sub></th><th>ℓ<sub>h4</sub></th><th>ℓ<sub>h1</sub></th><th>ℓ<sub>h2</sub></th><th>ℓ<sub>h3</sub></th></tr></thead><tbody><tr><td>Number of iterations</td><td>31</td><td>31</td><td>31</td><td>22</td><td>19</td><td>18</td></tr></tbody></table>
### 4.10. The Overrelaxation Algorithm.

For each iteration of Uzawa's algorithm, it is not necessary to have many iterations of the overrelaxation algorithm. Then, using stopping rule (35) we obtain the results shown in Table 8 for the problem (D) with quadratic elements and triangulation $e_{h1}$.

<table><caption>Table 8. Number of iterations for the overrelaxation method.</caption><thead><tr><th>Iteration of Uzawa's algorithm</th><th></th><th>2</th><th>6</th><th>10</th><th>20</th><th>30</th></tr></thead><tbody><tr><td rowspan="2">Number of iterations of overrelaxation</td><td>for problem (19)</td><td>34</td><td>14</td><td>9</td><td>10</td><td>11</td></tr><tr><td>for problem (20)</td><td>38</td><td>16</td><td>7</td><td>5</td><td>5</td></tr></tbody></table>

We can see in Table 8 that about ten iterations of the overrelaxation method are sufficient for each sub-problem, assuming that we use the optimal overrelaxation parameter.

### 4.11. Computing time.

For the same number of nodes, the computing times obtained for each kind of finite element are approximately of the same order. We indicate in Table 9 the computing times obtained with a CII-IRIS 80 computer for the problem on the square with linear elements.

<table><caption>Table 9. Mean computing time by iteration of Uzawa's algorithm.</caption><thead><tr><th></th><th>e<sub>h3</sub><br>512 triangles<br>289 nodes</th><th>e<sub>h4</sub><br>2048 triangles<br>1089 nodes</th></tr></thead><tbody><tr><td>Cholesky's method</td><td>1.7 s.</td><td>9 s.</td></tr><tr><td>Overrelaxation</td><td>3 s.</td><td>16 s.</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>Fig. 6. Isovalue lines of $u_h$ for the problem (44).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. Isovalue lines of $p_h$ for the problem (44).</figcaption></figure>
### 4.12. Some other numerical examples.

In problems (S) and (D) the right hand side $f$ belongs to the space $L^2(\Omega)$. However it is still possible to use Uzawa's algorithm described in Section 2.2. if the right hand sides are Dirac measures (see [6] ch. 4). For example, we have solved a problem corresponding to the following data:

$\Omega$ = disk of radius 0.5 centered at point $C = (0.5,0.5)$

$f = 100 \delta(c)$. $\delta(c)$ is the Dirac measure at point $c$.

$$u = \frac{100}{\pi} \left( \frac{1}{8} - r^2 \left( \frac{1}{2} + \log \frac{1}{2r} \right) \right).$$

$$\Delta u = - \frac{100}{2\pi} \log \frac{1}{2r} + \frac{100}{4\pi}.$$

The quality of the approximation is similar to that obtained for problems (S) and (D) except of course at the center $C$, where $-\Delta u$ is unbounded.

On Figures 6 and 7, the isovalue lines of the function $u_h$ and $p_h$ are shown respectively, for a problem in the square $[0,1] \times [0,1]$ [with

$$f = 200 \delta(x_o) - 100 \delta(x_1) - 100 \delta(x_2) - 100 \delta(x_3) - 100 \delta(x_4)$$

where

$$x_o = (1/2, 1/2); x_1 = (1/4, 1/4); x_2 = (3/4, 1/4); x_3 = (3/4, 3/4); x_4 = (1/4, 3/4).$$

## CONCLUSION.

It follows from the previous results that the method under consideration, which reduces the biharmonic equation to a sequence of Dirichlet problems for the operator $-\Delta$, is a convenient and effective method for solving the biharmonic problem. It suffices to use linear or quadratic elements, which yield a good approximation of $u$ and $\Delta u$ (except along the boundary for $\Delta u$, when the number of nodes is too small). In case of domains with arbitrary shapes, the value of the parameter $\rho$ will have to be estimated through preliminary trials with coarse meshes.

## ACKNOWLEDGEMENTS.

I would like to thank Professors P.G. Ciarlet and R. Glowinski for their helpful suggestions.
REFERENCES

[1] P.G. CIARLET, P.A. RAVIART. A mixed finite element method for the biharmonic equation. Proceedings of the Symposium on Mathematical Aspects of Finite Elements in Partial Differential Equations, Mathematics Research Center, University of Wisconsin, Madison, April 1-3, 1974.

[2] R. GLOWINSKI. Approximations externes, par éléments finis de Lagrange d'ordre un et deux, du problème de Dirichlet pour l'opérateur biharmonique. Méthodes itératives de résolution des problèmes approchés. Topics in Numerical Analysis, (J.H.H.Miller Ed.) Academic Press, London (1973),pp.123-171.

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

[4] J.H. ARGYRIS, I. FRIED, D.W. SCHARPF. The TUBA family of plate elements for the matrix displacement method, The Aeronautical J.R. Ae. S. 72 (1968), pp. 701-709.

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

[6] R. GLOWINSKI, J.L. LIONS, R. TREMOLIERES. Analyse Numérique des Inéquations Variationnelles, Paris, Dunod, (1976).

[7] B. MERCIER. On the frontal method for finite elements and associated questions of re-ordering. Rapport LABORIA, à paraître.

[8] R.S. VARGA. Matrix Iterative Analysis. Prentice Hall, Englewood Cliffs, New Jersey, (1962).