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

HAL
open science

# A dual algorithm for the solution of non linear variational problems via finite element approximation

Daniel Gabay, Bertrand Mercier

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

## To cite this version:

Daniel Gabay, Bertrand Mercier. A dual algorithm for the solution of non linear variational problems via finite element approximation. [Research Report] IRIA-RR-126, IRIA. 1975, pp.43. hal-04716124

HAL Id: hal-04716124

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

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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RIA
laboria

laboratoire de recherche
en informatique
et automatique

A dual algorithm
for the solution
of non linear variational problems
via finite element approximation

D. Gabay, B. Mercier

Institut de Recherche
d'Informatique
et d'Automatique

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

Rapport de Recherche n° 126

Juin 1975
A DUAL ALGORITHM
FOR THE SOLUTION OF NON LINEAR VARIATIONAL PROBLEMS
VIA FINITE ELEMENT APPROXIMATION*

D. Gabay, B. Mercier**

RÉSUMÉ

Pour des problèmes variationnels de la forme :

$$ \text{Inf}_{\text{veV}} \left\{ f(Av) + g(v) \right\} $$

on propose une méthode de dualité qui découple les difficultés relatives aux fonctionnelles f et g d'un éventuel mauvais conditionnement de l'opération linéaire A.

La technique est fondée sur l'emploi d'un lagrangien augmenté qui conduit à introduire un algorithme efficace et simple à implémenter.

On étudie aussi l'approximation par éléments finis de tels problèmes, compatible avec l'emploi de l'algorithme.

La méthode est finalement appliquée à la résolution de plusieurs problèmes de mécanique des milieux continus.

<footnote>* The authors are greatly indebted to Prof. GLOWINSKI, who originates the method, for his friendly encouragements and for very stimulating discussions during the research.</footnote>

<footnote>**Now : Ecole Polytechnique - Centre de Mathématiques Appliquées - 91120 - PALAISEAU.</footnote>
ABSTRACT

For variational problems of the form

$$ \inf_{v \in V} \{f(Av) + g(v)\}, $$

we propose a dual method which decouples the difficulties relative to the functionals f and g from the possible ill-conditioning effects of the linear operator A.

The approach is based on the use of an Augmented Lagrangian functional and leads to an efficient and simply implementable algorithm. We study also the finite element approximation of such problems, compatible with the use of our algorithm. The method is finally applied to solve several problems of continuum mechanics.

Daniel GABAY,
Centre National de la Recherche Scientifique,
Laboratoire d'Analyse Numérique, L.A. 189,
4, Place Jussieu - 75230 PARIS CEDEX 05 - France

Bertrand MERCIER,
IRIA - LABORIA - Domaine de Voluceau -
78150 LE CHESNAY - France.
A DUAL ALGORITHM FOR THE SOLUTION OF NON LINEAR VARIATIONAL
PROBLEMS VIA FINITE ELEMENT APPROXIMATION.

# 1 - INTRODUCTION.

Many problems, in Physics, Mechanics, and Mathematical
Economics, can be formulated as the following variational problem :

$$ (\mathcal{P}) \qquad \inf \{ f(Av) + g(v) \} . \\ v \in V $$

V and Y are two Hilbert spaces, endowed with the norm topology ;
A is a continuous linear operator from V into Y ; f and g
are convex functions defined respectively on Y and V, and
taking their values in ( -∞, +∞] . This formulation, first
used by ROCKAFELLAR [R3] to extend Fenchel's duality theory,
includes, in particular, the ordinary problems of convex programming.
We call v* the solution of (P), when it exists.

$$ \text{Let } \phi = f \circ A + g . $$

A natural approach to solve (P) consists in searching
directly for the minimum of φ(v) over V. If φ is differentiable,
iterative techniques based on the use of its gradient are available
to construct a sequence of points in V, converging to v*.
But the speed of convergence is slow if the operator A is
ill-conditioned [L4].

This difficulty may disappear if we use the following device :
in a first time, we introduce the additional variable y ∈ Y,
linked to the original variable v ∈ V by the constraint
Av - y = 0, and consider the constrained problem (P_c) on
V × Y, obviously equivalent to (P) :

$$ (\mathcal{P}_c) \qquad \inf \{ f(y) + g(v) \mid (v,y) \in V \times Y, \text{Av} - y = 0 \} . $$

We then eliminate the artificial constraint, just introduced,
by a dual approach to solve (P_c).
In the following, (.,..) and ((.,..)) denote the scalar
product on respectively Y and V, while |.| and ||.|| are
the corresponding norms. V' and Y' are the dual spaces of
V and Y. We denote the duality between V and V' by
<.,.>, while we identify Y' to Y.

The classical dual approach to ($\mathcal{P}_c$) consists in intro-
ducing a Lagrange multiplier $\lambda \in Y'$ corresponding to the
constraint and the Lagrangian functional $\mathcal{L}$, defined on
$(V \times Y) \times Y'$ with values in $(-\infty, +\infty]$

$$
(1.1) \qquad \mathcal{L}(v,y; \lambda) = g(v) + f(y) + (\lambda, Av - y).
$$

Duality theory [G1] [R3] tells us that, if there exists a
saddle point (v*,y*; λ*) of L on (V × Y) × Y' (*), then (v*,y*)
is solution of (Pc). In this case, we have also the equality
between

$$
\text{Inf}_{(v,y) \in V \times Y} \{ \text{Sup}_{\lambda \in Y'} \mathcal{L}(v,y; \lambda) \},
$$

where we recognize the primal problem ($\mathcal{P}_c$), and

$$
(\mathcal{D}) \quad \text{Sup}_{\lambda \in Y'} \{ \text{Inf}_{(v,y) \in V \times Y} \mathcal{L}(v,y; \lambda) \},
$$

which constitutes the dual problem.

In order to guarantee the existence of a unique solution to
the inner minimization problem in (D) under the same hypotheses
that insure the existence and uniqueness of a solution of (P),
it is convenient, as originally proposed by GLOWINSKI and
MARROCCO [G3], to define the dual problem with the Augmented
Lagrangian functional L_r. The expression of the Augmented
Lagrangian, introduced by HESTENES [H] and POWELL [P], is
obtained by adding to the standard Lagrangian L(v,y; λ),
expressed in (1.1), a term, depending on a positive parameter r,
penalizing for the violation of the constraint; in our case

$$
(1.2) \qquad \mathcal{L}_r (v,y; \lambda) = f(y) + g(v) + (\lambda, Av - y) + \frac{r}{2} |Av - y|^2 .
$$

$(v^*, y^*; \lambda^*)$ is a saddle point of $\mathcal{L}$ iff

$$
\mathcal{L}(v^*, y^*; \lambda) \le \mathcal{L}(v^*, y^*; \lambda^*) \le \mathcal{L}(v, y; \lambda^*) \quad \forall (v, y), \lambda \in (V \times Y) \times Y'
$$
This formulation presents also the advantage over the ordinary
penalty function (which corresponds to the case where $\lambda$ is
maintained equal to 0) that a minimizing sequence ${v_n}$ for
$(\mathcal{P})$ can be generated without making $r \to +\infty$, thus avoiding
the well-known ill-conditioning of ordinary penalty methods. For
any arbitrary fixed $r > 0$, if $(v^*, y^*; \lambda^*)$ is a saddle point
of $\mathcal{L}_r(v,y; \lambda)$, then $v^*$ is solution of $(\mathcal{P})$.

In practice, we must solve a sequence of unconstrained
minimization problems

$$
(1.3) \qquad \inf_{(v,y) \in V \times Y} \mathcal{L}_r (v,y; \lambda^n) ,
$$

where the multipliers $\lambda^n$ form a maximizing sequence of the
dual functional

$$
(1.4) \qquad d_r(\lambda) = \inf_{(v,y) \in V \times Y} \mathcal{L}_r(v,y;\lambda).
$$

Let $(v^{n+1}, y^{n+1})$ be the solution of (1.3). The
concave functional $d_r$ is differentiable and its gradient is
given by

$$
\nabla d_r(\lambda^n) = A_v^{n+1} - y^{n+1}
$$

We can use the gradient method to maximize $d_r(\lambda)$ by generating
the multipliers $\lambda^n$ according to the iteration

$$
(1.5) \qquad \lambda^{n+1} = \lambda^n + \rho (A_v^{n+1} - y^{n+1}).
$$

A subsequence of ${\lambda^n}$ converges to $\lambda^*$, provided that the
stepsize $\rho$ is chosen $0 < \rho < 2r$ [R3]. This dual algorithm
due to UZAWA [A], is mentioned in the pioneering paper of
Hestenes with $\rho = r$, and is used in [G5] to solve the non-
linear Dirichlet's problem. The convergence results for convex
programming problems, in finite dimensional spaces, have been
extended, by ROCKAFELLAR [R1] and BERTSEKAS [B1][K], to
the case where the minimization in (1.3) is performed
approximately at each cycle and is only asymptotically exact.

An alternative formulation of the dual function (1.4),
provided by FORTIN [F], consists in writing

$$d_r(\lambda) = \operatorname{Inf}_{v \in V} \{ g(v) + \operatorname{Inf}_{p \in Y} [f(Av + p) + \frac{r}{2} |p|^2 - (\lambda, p)] \}.$$

It is then possible to explicit the problem in $p$ with the proximation operator of $Y$ relatively to $f$ [M3]. The resulting problem is still non linear in $v$ and $Av$, and is of the same type as the original problem ($\mathcal{P}$). This approach is particularly interesting when $f$ is a support function of a closed convex set of $Y$, since the nondifferentiable unconstrained problem ($\mathcal{P}$) is replaced by a differentiable sup-inf problem. A classical dual method has been presented earlier for the same type of situation in [L5] [C2]. See also [B2] for an approach similar to [F].

We present, in this paper, a modification of Uzawa's algorithm, hinted in [G3], and experimented in [G5], without proof of convergence.

In this modification, the minimization in (1.3) is neither performed exactly nor subject to a termination criterion like in [B1], [R1]. This proposal is based on the observation that problem (1.3) involves a problem of minimization in $v$ coupled with a problem in $y$ by the term $-r(y, Av)$. An approximate solution $(v^{n+1}, y^{n+1})$ to (1.3) is therefore provided by the following algorithm :

Given $v^n, y^n, \lambda^n$:

$$(1.6) \quad \text{select } v^{n+1}, \text{ solution of } \operatorname{Inf}_{v \in V} \{g(v) + \frac{r}{2} \|Av\|^2 + \langle \lambda^n - ry^n, Av \rangle;$$

$$(1.7) \quad \text{select } y^{n+1}, \text{ solution of } \operatorname{Inf}_{y \in Y} \{f(y) + \frac{r}{2} \|y\|^2 - \langle \lambda^n + rAv^{n+1}, y \rangle;$$

which we complete by iteration (1.5) for the multiplier $\lambda^{n+1}$.

The algorithm is interesting for computation only if problems (1.6) and (1.7) are relatively easy to solve. If $g(\cdot)$ is a linear functional or a quadratic form, (1.6) becomes a quadratic problem. In the finite dimensional case, it can be solved by direct methods of matrix inversion, which are less sensitive to the ill-conditioning effect of $A$ than gradient methods. We observe also that (1.7) is a nonlinear problem in $y$ which is not sensitive at all to the ill-conditioning of $A$. We have
thus achieved a decoupling of f and A which constitutes
a serious improvement upon (P).

We have applied our algorithm to several problems of
continuum mechanics, namely the minimal hypersurfaces problem,
the problem of visco-plastic flow in a cylindrical pipe, and
two problems of elasto-plasticity : torsion and equilibrium.
All these problems are of the special form (P) where g(.)
is a linear functional, and to which we shall restrict our
analysis in the following. In 5.1, for example,
V is taken as the Sobolev space H1(Ω) (see e.g. [L3]),
where Ω is a regular open subset of RN (N = 1, 2 or 3).
Y is taken as [L2(Ω)]N and A is the gradient operator.
Finally, the convex function f can be represented by an
integral over Ω.

Although our algorithm can be applied to problem (P)
in this functional framework, numerical computations must,
in practice, be performed in finite dimensional spaces. We
thus construct internal approximations Vh ⊂ V and Yh ⊂ Y
of finite dimension [G2] via a finite element method [Z].
We introduce a regular triangulation T_h of Ω in a finite
number of "triangles" T of size ≲ h. V_h is chosen, for instance, as
a space of continuous functions which are piecewise polynomial.
Since we are not able, in general, to calculate exactly the
integral defining f, we must use a formula for numerical
integration with a certain number of integration nodes. Naturally,
the approximation A_h of the operator A is defined by :
the values taken by Av_h at each integration node. We thus
choose for Y_h a space of piecewise constant functions of
dimension N= N times the number of integration nodes. We
establish that the solution v_h* of the approximate problem

$$
(\mathcal{P}_h) \quad \text{Inf} \left\{ f (A_h v_h) - \langle b, v_h \rangle \right\} \\
v_h \in V_h
$$

converges towards v* as h → 0, if the numerical integration is
precise enough. Using our dual algorithm to solve (P_h), we
benefit of a powerful decomposition property, [G5], since each iteration
consists in a quadratic problem (1.6) in V_h and h nonlinear
problems (1.7) in $\mathbb{R}^N$ independent of each other. The iteration (1.5) can also be performed component by component.

The article is organized as follows: in Section 2, after a brief review of duality theory, we study the properties of the augmented Lagrangian, and, in particular, the existence and uniqueness of a saddle point. In Section 3, we state a dual algorithm to compute such a saddle point, and we establish its convergence. In Section 4, we turn to the approximation of infinite dimensional problems. A convergence result is given and applied to a finite element approximation compatible with our dual approach. We use this approximation for several problems of continuum mechanics, stated in Section 5, and apply our algorithm. Numerical results are reported in Section 6, and compared with direct methods of solution.
## 2 - CONVEXITY PROPERTIES AND THE AUGMENTED LAGRANGIAN.

### Duality Theory and Saddle-Points :

Duality theory has recently received an elegant formulation [R3][E], based on the consideration of a family of perturbed problems $(\mathcal{P}_p)$ associated to a problem $(\mathcal{P})$: $\text{Inf}_{v \in V} \phi(v)$. We consider a general bifunction $\Phi:V \times Y \to (-\infty, +\infty]$ such that

$$ \Phi(v,0) = \Phi(v) $$

and the problem depending on a perturbation $p \in Y$

$$ (\mathcal{P}_p) \qquad \text{Inf}_{v \in V} \Phi(v,p). $$

In this framework, we define a Lagrangian function $\Lambda: V \times Y' \to [-\infty, +\infty]$ associated to $(\mathcal{P})$ by the relation

[2.1]

$$ \Lambda(v; \lambda) = \text{Inf} \{ \Phi(v,p) - (\lambda,p) \}. \qquad p \in Y $$

We must immediately observe that this is different from the classical Lagrangian function defined for constrained minimization problem like $(\mathcal{P}_c)$, since this function $\Lambda$ is directly associated to the unconstrained problem $(\mathcal{P})$. $\Lambda$ depends, however, on the choice of the perturbation bifunction $\Phi$.

We can verify that

$$ \text{Sup}_{\lambda \in Y'} \Lambda(v; \lambda) = \Phi(v,0); $$

hence the problem

$$ \text{Inf}_{v \in V} \text{Sup}_{\lambda \in Y'} \Lambda(v; \lambda) $$

is nothing else than problem $(\mathcal{P})$ and is independent of $\Phi$.

Parallel to the duality theory for convex programming in terms of mini-max, we define a dual problem to $(\mathcal{P})$

$$ (\mathcal{D}) \qquad \text{Sup} \quad \text{Inf}_{\lambda \in Y'} \Lambda(v; \lambda) \qquad v \in V $$
From now on, let

$$ (\mathcal{P}) \quad \operatorname{Inf}_{v \in V} \{ f(Av) - \langle b, v \rangle \} \quad \text{with} \quad b \in V' $$

We assume that $f$ is the sum of two lower semi-continuous convex functions $f_1$ and $f_2$ from $Y$ into $(-\infty, +\infty]$,

$$ (2.2) \qquad f = f_1 + f_2 $$

$f_1$ is $C^1$. Gateaux-differentiable and its gradient $f_1'$ is weakly continuous on finite dimensional subspaces of $V$ and **strongly monotone**: i.e. there exists $\gamma > 0$ such that

$$ (2.3) \qquad \forall y, z \in Y \quad (f'_1(y) - f'_1(z), y-z) \geq \gamma |y-z|^2. $$

$f_2$ is non empty. The interior in $Y$ of dom $f_2$ is the whole $Y$.

We assume moreover that

the operator $A'A$ is an isomorphism from $V$ onto $V'$; i.e. there exists $\alpha > 0$ such that

$$ (2.5) \qquad |Av|^2 \geq \alpha^2 ||v||^2. $$

The strict convexity of $f^o A$ follows.

We consider the specific perturbation bifunction

$$ \Phi(v,p) = f(Av+p) - \langle b,v \rangle $$

The dual problem is defined by

$$ \begin{array}{lll} \operatorname{Sup} & \operatorname{Inf} & \{ -\langle b, v \rangle + \operatorname{Inf}_{p \in Y} [f(Av+p) - (\lambda, p)] \} = \\ & \lambda \in Y' & v \in V \\ & \varepsilon \in \mathcal{V} & p \in Y \end{array} $$

$$ (2.6) \qquad \begin{array}{l} = \operatorname{Sup} \operatorname{Inf} \{ f(y) + (\lambda, Av-y) - \langle b, v \rangle \}, \\ \lambda \in Y' \quad (v,y) \in \mathcal{V} \times \mathcal{Y} \end{array} $$

after the change of the variable $y = Av + p$. Inside the brackets, we recognize the expression (1.1) of the standard Lagrangian functional $\mathcal{L}(v, y; \lambda)$ corresponding to the constrained problem $(\mathcal{P}_c)$. We have thus established the equivalence between the choice of the perturbation bifunction $\Phi$ and the introduction
of the artificial variable y in problem (P<sub>c</sub>).

We study, now, the existence of solutions to (P) and (D).

**PROPOSITION 2.1** : Under (2.3), (2.5), there exists a unique
solution v* to (P).

Proof : There exists v₀ ∈ V such that f₂(Av₀) = (f₂ ○ A)(v₀) < +∞ otherwise (P) has no meaning.

Let now z = Av₀ and φ(t) = f₁(z+t(y-z)) for y∈Y, using
(2.3) and the equality φ(1) = φ(0) + ∫₁₀ φ'(t)dt , we deduce

$$
(2.7) \qquad f_1(y) \ge f_1(z) + (f'_1(z), y-z) + \frac{\gamma}{2} |y-z|^2
$$

Since $f_2$ is a proper convex function, even if (2.4) is not
satisfied, it has a continuous affine lower bound :
there exists $y_2 \in Y$ and $\beta \in \mathbb{R}$ such that :

$$
(2.8) \qquad f_2(y) \ge (y_2, y) + \beta \quad \forall y \in Y
$$

Then, applying (2.5), we deduce the coercivity of (P):
($\phi(v) \to +\infty$ if $||v|| \to +\infty$), and then, the existence of $v^*$ [E].
Uniqueness follows from strict convexity of $f \circ A$. ■

We now prove the existence of a saddle point of $\mathcal{L}(v,y;\lambda)$.
we recall first that the subgradient of a function $\phi: V \to (-\infty, +\infty]$ at a point $u \in V$ is the set (possibly empty)

$$
(2.9) \quad \partial \phi(u) = \{ \mu \in V' \mid \forall v \in V \quad \phi(v) - \phi(u) \ge \langle \mu, v - u \rangle \} .
$$

**THEOREM 2.1** : *Any saddle point (v*,y*; λ*) of *L*(v,y; λ) over (V × Y) × Y' satisfies :*

(2.10) $v^*$ is solution of $(\mathcal{P})$, $y^* = Av^*$ and $\lambda^* \in \partial f(Av^*)$ with $A'\lambda^* = b^{(*)}$
Conversely, if (2.3), (2.4), (2.5) hold, there exists at least
one such saddle point.

<footnote>(*) A' is the continuous linear operator from Y' → V', dual of A.</footnote>
**Proof :** We first assume that $(v^*, y^*, \lambda^*)$ is a saddle point of

$\mathcal{L}(v, y; \lambda) = f(y) + (\lambda, Av-y) - \langle b, v \rangle$, i.e. which satisfies

$$
(2.11) \qquad \mathcal{L}(v^*, y^*; \lambda) \le \mathcal{L}(v^*, y^*; \lambda^*) \quad \forall \lambda \in Y'
$$

$$
(2.12) \qquad \mathcal{L}(v^*, y^*; \lambda^*) \le \mathcal{L}(v, y; \lambda^*) \quad \forall v, y \in V \times Y.
$$

(2.11) implies immediately $y^* = Av^*$. Let $v = v^*$ in (2.10).
We get

$$
f(y) - f(Av^*) \geq ( \lambda^*, y - Av^* ) \quad \forall y \in Y,
$$

hence the necessary condition

$$
(2.13) \qquad \lambda^* \in \partial f(Av^*).
$$

Let $y = Av$ in (2.12) . We get, this time,

$$
f(Av) - \langle b, v \rangle \geq f(Av^*) - \langle b, v^* \rangle \quad \forall v \in V,
$$

which shows that $v^*$ must be solution of $(\mathcal{J})$. It can also
be written

$$
f(Av) - f(Av^*) \geq \langle b, v - v^* \rangle ,
$$

which means that

$$
(2.14) \qquad b \in \partial (f \circ A)(v^*).
$$

We now prove that such a triple $(v^*, y^*; \lambda^*)$ exists.
Proposition 2.1 concludes to the existence of $v^*$ since
(2.3), (2.5) hold. The qualification hypothesis (2.4) guarantees
that there exists a point of $Y$ where $f$ is finite and continuous.
Then ([E] ch. 1, PR 5.7)

$$
\partial (f \circ A) (v) = A' \partial f (Av) \quad \forall v \in V ;
$$

Therefore (2.14) shows that $\partial f(Av^*)$ is non empty ; we can
find a $\lambda^* \in \partial f(Av^*)$ such that

$A' \lambda^* = h$
since $A'$ is an operator onto $Y'$ by (2.5). We verify easily
that $(v^*, y^*, \lambda^*)$ is actually a saddle point of $\mathcal{L}$. ■

**REMARK 1** : Since $f_1$ is Gateaux-differentiable and because of the qualification hypothesis on $f_2$, we have ($[E]$, ch. 1, PR 5.6)

$$
\partial f = \partial f_1 + \partial f_2
$$

We can therefore replace the condition $\lambda \in \partial f(Av)$ by

$$
(2.15) \qquad \lambda^* - f_1(Av^*) \in \partial f_2(Av^*) . \blacksquare
$$

**REMARK 2** : In the general case where the convex function $g(\cdot)$ is not restricted to be a linear functional, we can prove the existence of a saddle point of the Lagrangian $\Lambda(v,\lambda)$ and hence of $\mathcal{L}(v,y;\lambda)$ under the stronger qualification hypothesis : there exists $v_0 \in \text{dom } g$ such that $Av_0 \in \text{int } \text{dom } f$. ([E], ch. 3) ■

The Augmented Lagrangian :

Although (D) has been shown to have a solution λ*,
the inner minimization problem in (2.6) may not have a bounded
solution for every λ ∈ Y'. For this reason, we switch to
the augmented Lagrangian. Like the standard Lagrangian was
arising from the consideration of a particular perturbation
bifunction, we can derive the _augmented Lagrangian_ from the
class of perturbations defined by [R1][F]

$$
\begin{equation} \tag{2.16}
\begin{aligned}
\Phi_r(v,p) &= \Phi(v,p) + \frac{r}{2} |p|^2 && , \\
&= g(v) + f(Av+p) + \frac{r}{2} |p|^2 .
\end{aligned}
\end{equation}
$$

We define now

$$
(2.17) \qquad \Lambda_r (v ; \lambda) = \inf_{p \in Y} \{ \Phi(v,p) + \frac{r}{2} |p|^2 - (\lambda, p) \} ;
$$

the dual problem

$$
(2.18) \quad \begin{array}{ll} \text{Sup} & \text{Inf} \quad \Lambda_r (v ; \lambda) \\ \mathbb{1} \in \mathcal{Y}' & v \in \mathcal{V} \end{array}
$$
can be written, after the change of variable $y = Av + p$

$$ ( \mathcal{D}_r ) \ \text{Sup} \ \text{Inf} \ \{ f(y) - \langle b, v \rangle + (\lambda, Av-y) + \frac{r}{2} |Av-y|^2 \} , \quad \lambda \in Y' \quad (v,y) \in V \times Y $$

which agrees with $\mathcal{L}_r(v,y;\lambda)$ of (1.2).

The virtues of the augmented Lagrangian result from the following properties.

**THEOREM 2.2** : If $r > 0$, any saddle point of the augmented Lagrangian $\mathcal{L}_r$ is saddle point of the standard Lagrangian $\mathcal{L}$ and conversely.

(Hence, if (2.3), (2.4), (2.5) hold, according to theorem 2.1, there exists a saddle point $(v^*, y^*; \lambda^*)$ of $\mathcal{L}_r$).

*Proof* : By definition (1.2), $\mathcal{L}_r$ differs from $\mathcal{L}$ by a non-negative term ; hence

$$ \mathcal{L}(v,y;\lambda) \le \mathcal{L}_r(v,y;\lambda), $$

with equality if $y = Av$.

If $(v^*, y^*; \lambda^*)$ is saddle point of $\mathcal{L}$, we have $y^* = Av^*$.
Therefore

$$ (2.19) \quad \mathcal{L}_r(v^*, y^*; \lambda) = \mathcal{L}(v^*, y^*; \lambda) \le \mathcal{L}(v^*, y^*; \lambda^*) = \mathcal{L}_r(v^*, y^*; \lambda^*) \le \mathcal{L}(v, y; \lambda^*) \le \mathcal{L}_r(v, y; \lambda^*) $$

and $(v^*, y^*; \lambda^*)$ is also a saddle point of $\mathcal{L}_r$. Conversely, if $(v^*, y^*; \lambda^*)$ is a saddle point of $\mathcal{L}_r$, $y^* = Av^*$, then :

$$ \mathcal{L}(v^*, y^*; \lambda^*) = \mathcal{L}_r(v^*, y^*, \lambda^*) \le \mathcal{L}_r(v, y; \lambda^*) \quad \forall (v, y) \in V \times Y $$

Thus, returning to the definition (1.2) of $\mathcal{L}_r$ :

$$ (2.20) \quad f(y^*) - \langle b, v^* \rangle \le f(y) - \langle b, v \rangle + (\lambda^*, Av-y) + \frac{r}{2} |Av-y|^2 $$

Choosing, first, $v = v^*$ in (2.20), we get, since $y^* = Av^*$:

$$ (2.21) \quad f(y) - f(y^*) + (\lambda^*, y^* - y) + \frac{r}{2} |y-y^*|^2 \ge 0 \quad \forall y \in Y; $$

For any $z \in Y$, we take $y = \theta z + (1-\theta)y^*$ with $\theta \in ]0,1[$ in (2.21).
We get, using the convexity of $f$ :

$$ (2.22) \qquad \theta(f(z) - f(y^*)) - \theta(\lambda^*, z - y^*) + \frac{r}{2}\theta^2|z - y^*|^2 \ge 0; $$

Dividing by $\theta$, and making $\theta \to 0$, we get :

$$ (2.23) \qquad f(y^*) - (\lambda^*, y^*) \ge f(y) - (\lambda^*, y) \quad \forall y \in Y. $$

On the other hand, choosing $y = y^*$ in (2.20), we get, since $y^* = Av^*$:

$$ (2.24) \qquad \frac{r}{2} |A(v-v^*)|^2 + (\lambda^*, A(v-v^*)) - \langle b, v-v^* \rangle \ge 0 \quad \forall v \in V; $$

In the same way, taking $v = \theta w + (1-\theta)v^*$ in (2.24) with $\theta \in ]0,1[$, we get, as $\theta \to 0$ :

$$ (2.25) \qquad (\lambda^*, Av^*) - \langle b, v^* \rangle \le (\lambda^*, Aw) - \langle b, w \rangle \quad \forall w \in V; $$

Finally, adding (2.23) and (2.25) we obtain the desired result :

$$ \mathcal{L}(v^*, y^*; \lambda^*) = f(y^*) - \langle b, v^* \rangle \le f(y) + (\lambda^*, Av - y) - \langle b, v \rangle = \mathcal{L}(v, y; \lambda^*) \quad \forall (v, y) \in V \times Y $$
the other inequality being obvious. \blacksquare

The augmented Lagrangian is useful to evaluate the dual functional :

$$ (2.26) \qquad d_r(\lambda) = \begin{cases} \mathcal{L}(v, y; \lambda) \\ (v, y) \in V \times Y \end{cases} $$

**THEOREM 2.3** : Under (2.3), (2.5), for every $\lambda \in Y'$, there exists a unique solution $(v_\lambda, y_\lambda)$ to the minimization problem in (2.24).

*Proof* : From (2.5), (2.7) and (2.8), we see that $\mathcal{L}_r(v, y; \lambda) \to +\infty$ if $v$ or $y \to +\infty$, so that the problem (2.26) is coercive. Then, it has a solution for any $\lambda \in Y'$.

Uniqueness follows from the strict convexity of $\mathcal{L}_r$, for $r > 0$. ■
### 3 - A DUAL ALGORITHM :

**Description :**

From now on, we assume that the hypotheses (2.2) to (2.5) hold. To solve $(\mathcal{O}_r)$, we can maximize the differentiable function $d_r(\lambda)$ using a gradient algorithm. This forms the basis for Uzawa's algorithm [A]:

$$
\begin{array}{l}
\text{(3.1) UZAWA'S ALGORITHM :} \\
\quad \text{Let } \lambda^0 \in Y. \text{ By induction } \lambda^n \text{ being given} \\
\quad \text{Step 1 - Find } v^{n+1}, y^{n+1} \text{ minimizing on } V \times Y \\
\qquad \mathcal{L}_r(v,y;\lambda^n) = f(y) - \langle b,v \rangle + (\lambda^n, Av-y) + \frac{r}{2} |Av-y|^2 . \\
\quad \text{Step 2 - Make} \\
\qquad \lambda^{n+1} = \lambda^n + \rho (Av^{n+1} - y^{n+1}) . \blacksquare
\end{array}
$$

Because of the convexity of $f$, the problem in step 1 is equivalent to the variational inequality :

$f(y) - f(y^{n+1}) + (\lambda^n + rAv^{n+1}) - ry^{n+1}, Av - (y - y^{n+1}) - \langle b, v \rangle \ge 0 \quad \forall v \in V, \forall y \in Y$;

$$
\text{(3.2)} \quad r(Av^{n+1}, Av) = (ry^{n+1} - \lambda^n, Av) - \langle b, v \rangle \quad \forall v \in V,
$$
and a variational inequality in $Y$, equivalent to

$$
\text{(3.3)} \quad (f_1(y^{n+1}), y - y^{n+1}) + f_2(y) - f_2(y^{n+1}) + (ry^{n+1} - rAv^{n+1} - \lambda^n, y - y^{n+1}) \ge 0 \quad \forall y \in Y,
$$
since $f_1$ is Gateaux-differentiable. These two problems are coupled to each other by the terms $(y^{n+1}, Av)$ in (3.2) and $(Av^{n+1}, y - y^{n+1})$ in (3.3). A simple way to decouple these problems consists in replacing, for instance, $y^{n+1}$ by $y^n$ in (3.2). The resulting system is much simpler to solve and
provides an approximation to the minimization in step 1. We
state the new algorithm :

$$
(3.4) \quad \text{\textbf{A MODIFIED DUAL ALGORITHM}}
\begin{align*}
& \text{Let } (y^0, \lambda^0) \in Y \times Y'. \quad \text{By induction } (y^n, \lambda^n) \text{ being given} \\
& \text{Step 1 : Find } v^{n+1} \text{ such that :}
\end{align*}
$$

$$
(3.5) \quad r(Av^{n+1}, Av) = (ry^n - \lambda^n, Av) - \langle b, v \rangle \quad \forall v \in V.
\quad \text{Step 2 : Find } y^{n+1} \text{ such that :}
$$

$$
(3.6) \quad (f_1'(y^{n+1}), y - y^{n+1}) + f_2(y) - f_2(y^{n+1}) + (ry^{n+1} - \lambda^n - rAv^{n+1}, y - y^{n+1}) \ge 0 \quad \forall y \in Y.
\quad \text{Step 3 : Make}
$$

$$
(3.7) \qquad \lambda^{n+1} = \lambda^n + \rho(Av^{n+1} - y^{n+1}) . \blacksquare
$$

**REMARK 1** : (3.5) is the variational formulation of the linear
system

$$
rA'A v^{n+1} = rA'y^n - A' \lambda^n - b ,
$$

which is invertible by (2.5). At each iteration, Step 1 consists
of the resolution of this system for a different second member. •

**REMARK 2** : (3.6) is a nonlinear variational problem in y
independent of A. An interesting decomposition occurs if
Y is finite dimensional and can be viewed as the product of
k spaces Y_i such that f separates itself in

$$
f(y) = \sum_{i=1}^{k} f_i(y_i) , \text{ with } y = (y_1, \dots, y_k);
$$

(3.6) is then equivalent to k variational problems

$$
(f'_i(y_i^{n+1}), y_i - y_i^{n+1}) + f_{2i}(y_i) - f_{2i}(y_i^{n+1}) + (ry_i^{n+1} - s_i^n, y_i - y_i^{n+1}) \ge 0 \quad \forall y_i \in Y_i
$$
where $s_i^n$ is the component on $Y_i$ of $s^n = \lambda^n + rAv^{n+1}$. ■
### Convergence :

According to theorem 2.2, there exists a saddle point $(v^*, y^*; \lambda^*)$ satisfying (2.8). It verifies also

$$ (3.8) \qquad r(Av^*, Av) = (ry^* - \lambda^*, Av) - \langle b, v \rangle \quad \forall v \in V, $$

$$ (3.9) \qquad (f_1(y^*), y-y^*) + f_2(y) - f_2(y^*) + (ry^* - \lambda^* - rAv^*, y-y^*) \ge 0 \quad \forall y \in Y $$

$$ (3.10) \qquad y^* = Av^*. $$

Subtracting (3.8) from (3.5) yields :

$$ (3.11) \qquad r(A(v^{n+1} - v^*), Av) = (r(y^n - y^*) - (\lambda^n - \lambda^*), Av) \quad \forall v \in V. $$

We introduce the **projection operator** $P$ from $Y$ onto the range of $A$, $R(A)$. Since $R(A)$ is closed, we have $Y = R(A) \oplus R(A)^\perp$. Given any $y \in Y$, $Py$ is the unique element of $R(A)$ such that

$$ (Py, Av) = (y, Av) \quad \forall v \in V. $$

This allows to write (3.11) under the explicit form

$$ (3.12) \qquad A(v^{n+1} - v^*) = P(y^n - y^*) - \frac{1}{r} P(\lambda^n - \lambda^*) . $$

Adding (3.9) with $y = y^{n+1}$ to (3.6) with $y = y^*$, we get
$$ (f_1(y^{n+1}) - f_1(y^*), y^{n+1} - y^*) + (r(y^{n+1} - y^*) - (\lambda^n - \lambda^*) - rA(v^{n+1} - v^*), y^{n+1} - y^*) \le 0 $$

where the terms in $f_2(y^{n+1})$ and $f_2(y^*)$ have cancelled out.
Using the expression (3.12) of $A(v^{n+1} - v^*)$, we obtain

$$ (3.13) \qquad (f_1(y^{n+1}) - f_1(y^*), y^{n+1} - y^*) + r|y^{n+1} - y^*|^2 \le r(y^n - y^*, P(y^{n+1} - y^*)) + (\lambda^n - \lambda^*, (I-P)(y^{n+1} - y^*)) $$

where $I$ is the identity operator on $Y$. A majorant for the first term is given by Schwartz' inequality :

$$ (3.14) \qquad (y^n - y^*, P(y^{n+1} - y^*)) = (P(y^n - y^*), P(y^{n+1} - y^*)) < \frac{1}{2} \{|P(y^n - y^*)|^2 + |P(y^{n+1} - y^*)|^2\}. $$
We deduce from (3.7), (3.10), (3.12)

$$ \lambda^{n+1} - \lambda^* = \lambda^n - \lambda^* - \rho [y^{n+1} - y^* - P(y^n - y^*) + \frac{1}{\rho} P(\lambda^n - \lambda^*)] $$

The projection of the equality of $R(A)$ and $R(A)^\perp$ yields

$$ (3.15) \qquad P(\lambda^{n+1} - \lambda^*) = \left(1 - \frac{\rho}{r}\right) P(\lambda^n - \lambda^*) - \rho P(y^{n+1} - y^n), $$

$$ (3.16) \qquad (I-P)(\lambda^{n+1} - \lambda^*) = (I-P)(\lambda^n - \lambda^*) - \rho (I-P)(y^{n+1} - y^*) \\ \text{and squaring the norm of both members of (3.16)} $$

$$ (3.17) \qquad \left|(I-P)(\lambda^{n+1} - \lambda^*)\right|^2 = \left|(I-P)(\lambda^n - \lambda^*)\right|^2 + \rho^2 \left|(I-P)(y^{n+1} - y^*)\right|^2 - 2\rho (\lambda^n - \lambda^*, (I-P)(y^{n+1} - y^*)) \\ \text{We now use the strong monotonicity of } f; \text{ (2.3) to obtain} $$

$$ (3.18) \qquad (f'_1(y^{n+1}) - f'_1(y^*), y^{n+1} - y^*) \geq \gamma |y^{n+1} - y^*|^2 ; $$

Combining (3.14), (3.17), (3.18) with (3.13) gives

$$ (3.19) \qquad \gamma |y^{n+1}-y^*|^2 + \left(r-\frac{\rho}{2}\right) |(I-P)(y^{n+1}-y^*)|^2 + \frac{r}{2} |P(y^{n+1}-y^*)|^2 + \frac{1}{2\rho} |(I-P)(\lambda^{n+1}-\lambda^*)|^2 \\ \leq \frac{r}{2} |P(y^n-y^*)|^2 + \frac{1}{2\rho} |(I-P)(\lambda^n-\lambda^*)|^2. $$

Adding up (3.19) for $n=0, 1, \dots, N$:

$$ (3.20) \qquad \gamma \sum_{n=0}^{N} |y^{n+1} - y^*|^2 + \left(r - \frac{\rho}{2}\right) \sum_{n=0}^{N} |(I-P)(y^{n+1} - y^*)|^2 + \frac{1}{2\rho} |(I-P)(\lambda^{N+1} - \lambda^*)|^2 \\ \leq \frac{r}{2} |P(y^0 - y^*)|^2 + \frac{1}{2\rho} |(I-P)(\lambda^0 - \lambda^*)|^2. $$

This shows that for any choice of $\rho$ such that $0 < \rho \leq 2r$,
and for any $N$, the series

$$ \sum_{n=0}^{N} |y^{n+1} - y^*|^2 $$

remains bounded and so, convergent; therefore

$$ \lim_{n \to +\infty} |y^n - y^*|^2 = 0 $$

which proves the strong convergence of $\{y^n\}$ to $y^*$ in $Y$.
Squaring the norm of (3.15) and setting $\theta = \frac{\rho}{r}$, we get

$$ (3.21) \qquad |P(\lambda^{n+1}-\lambda^*)| = (1-\theta)^2 |P(\lambda^n-\lambda^*)|^2 + \rho^2 |P(y^{n+1}-y^n)|^2 - 2\rho(1-\theta)(P(\lambda^n-\lambda^*),P(y^{n+1}-y^n)) $$

We use Schwartz' inequality under the form

$$ (P(\lambda^n - \lambda^*), P(y^{n+1} - y^n)) \le \frac{1}{2} \{ \epsilon |P(\lambda^n - \lambda^*)|^2 + \frac{1}{\epsilon} |P(y^{n+1} - y^n)|^2 \} $$

where $\epsilon$ is an arbitrary small positive number; (3.21) gives

$$ (3.22) \qquad |P(\lambda^{n+1}-\lambda^*)|^2 \le ((1-\theta)^2 + \rho\epsilon|1-\theta|)|P(\lambda^n-\lambda^*)|^2 + (\rho^2 + |1-\theta|\frac{\rho}{\epsilon}) |P(y^{n+1}-y^n)|^2 $$

Consider the sequence of positive scalars $u_n$ defined by

$$ (3.23) \qquad u_{n+1} = a u_n + b w_n $$

with $0 \le a < 1$, $b > 0$, $w_n > 0$ and $\sum_{n=0}^{\infty} w_n < +\infty$.

It is easy to verify that for any $N$, $\sum_{n=0}^{N} u_n < +\infty$ and therefore $u_n \to 0$. We observe that the sequence

$$ w_n = |P(y^{n+1} - y^n)|^2 $$

satisfies the previous condition since

$$ \sum_{0}^{\infty} w_n \le 2 \sum_{0}^{\infty} |P(y^n - y^*)|^2 \le 2 \sum_{0}^{\infty} |y^n - y^*|^2, $$

which has been shown to be bounded. We choose $b = (\rho^2 + |1-\theta|\frac{\rho}{\epsilon})$ and $\rho$ such that

$$ (3.24) \qquad a = \left(1 - \frac{\rho}{r}\right)^2 + \rho\epsilon \left|1 - \frac{\rho}{r}\right| < 1 $$

It is easy to check that, if $u_0 = |P(\lambda^0 - \lambda^*)|^2$, $|P(\lambda^n - \lambda^*)|^2 \le u_n$ for every $n$ and $|P(\lambda^n - \lambda^*)| \to 0$. Since $\epsilon$ is arbitrary, the inequality (3.24) is satisfied for all $0 < \rho < 2r$.

(3.12) allows us to conclude that $|A(v^{n+1}-v^*)| \to 0$ and, applying (2.5), we see that $|v^{n+1}-v^*| \to 0$.

We have thus established that the sequence $\{v^n\}$, constructed by algorithm (3.4), converges strongly to $v^*$ in $V$. However, we have only proved that $|P(\lambda^n - \lambda^*)| \to 0$; we can deduce from (3.20) that $|(I-P)(\lambda^n - \lambda^*)|$ remains bounded. In fact $\lambda^*$ may not be unique.
**THEOREM 3.1** : For every stepsize $\rho$ such that $0 < \rho < 2r$, the sequence $\{(v^n, y^n)\}$, constructed by algorithm (3.4), converges strongly in $V \times Y$ to $(v^*, Av^*)$ where $v^*$ is the unique solution of $(\mathcal{P})$. The sequence $\{\lambda^n\}$ remains bounded in $Y$. ■

**REMARK 1** : Assuming only the monotonicity of $f_1'$ instead of (2.3), we obtain a modification of (3.20) :

$$ (r - \frac{\rho}{2})_{n=0}^{N} |(I-P)(y^{n+1} - y^{*})|^2 + \frac{r}{2} |P(y^{N+1} - y^{*})|^2 + \frac{1}{2\rho} |(I-P)(\lambda^{N+1} - \lambda^{*})|^2 \le \frac{r}{2} |P(y^{0} - y^{*})|^2 + \frac{1}{2\rho} |(I-P)(\lambda^{0} - \lambda^{*})|^2 ; $$

thus, for all $0 < \rho < 2r$, the sequence $\{u^n\}$ generated by algorithm (3.4) remains bounded.

To prove theorem 3.1, it is therefore enough to assume the strong monotonicity of $f_1'$ on every bounded set of $Y$. ■

**REMARK 2** : Convergence of the algorithm can be shown for $\rho=r$ under the weaker hypothesis on $f_1'$: suppose that for all C positive, there exists a family of forcing functions $\delta_C$: $\mathbb{R}^+ \to \mathbb{R}^+$ strictly monotone increasing with $\delta_C(0) = 0$ such that

$$ (3.25) \qquad (f_1'(y) - f_1'(z), y-z) \ge \delta_C(|y-z|) \quad \forall y,z \in Y \quad \text{with } |y| \le C \\ \text{and } |z| \le C $$

We can still prove that $y^n \to y^*$; for $\rho = r$. (3.12) and (3.15) give immediately the convergence of $v^n \to v^*$. ■

This assumption is satisfied in the important case where $Y$ is finite dimensional and $f_1'$ strictly convex $[C1]$. ■

**REMARK 3** : For a differentiable problem ($f_2 = 0$), we have the existence and uniqueness of $\lambda^*$. Moreover $\lambda^n \to \lambda^*$. ■
## 4 - APPROXIMATION VIA FINITE ELEMENTS

It is in practice convenient to approximate the infinite dimensional problem $(\mathcal{P})$ by a finite dimensional one.

Let $V_h \subset V$ be a family of internal approximation of $V$:

For $v \in V$, there exists $v_h \in V_h$ such that $v_h \to v$ when $h \to 0$ (In fact, it is sufficient to check it for a dense subset of $V$).

In the same way, we shall approximate $Y$ by $Y_h \subset Y$.

Since $A(v_h)$

may not be included in $Y_h$, we must approximate $A$ by a continuous linear operator $A_h: V_h \to Y_h$. We assume that $A_h$ satisfies:

$$ (4.1) \qquad \left\{ \begin{array}{ll} (i) & \alpha' \|v_h\| \le |A_h v_h| \le M \|v_h\| \quad \forall v_h \in V_h ; \\ (ii) & \text{If } v_h \to v \text{ weakly in } V, \text{ then } A_h v_h \to Av \text{ weakly in } \\ (iii) & \text{If } v_h \to v \text{ strongly in } V, \text{ then } A_h v_h \to Av \text{ strongly in } \\ (iv) & |A_h v_h - Av_h| \to 0 \text{ when } \|v_h\| \le C \text{ and } h \to 0. \end{array} \right. $$

We restrict, from now on, our analysis to the simplified case where

$$ (4.2) \qquad f(y) = \frac{1}{2} |y|^2 + f_2(y) $$

$f_2$ must be compatible with the approximation in the following sense:

$$ (4.3) \qquad \forall v \in V, \exists v_h \in V_h \text{ with } v_h \to v \text{ and } f_2(A_h v_h) \to f_2(Av) $$

we verify that the hypothesis (2.3) is satisfied in this case.

According to proposition (2.1), there exists a unique solution $v^*$ of $(\mathcal{P})$.

We approximate $(\mathcal{P})$ by the following problem:

$$ (\mathcal{P}_h) \qquad \inf_{v, \in V} \{ f(A_h v_h) - \langle b, v_h \rangle \} $$
We now establish the convergence of the approximation
in the following sense :

**THEOREM 4.1** : *Under the hypothesis (4.1), (4.2), (4.5), there is a unique solution $v_h^*$ to $(\mathcal{P}_h)$ which converges to the solution $v^*$ of $(\mathcal{P})$ :*

$v_h^* \to v^*$ in $V_h$ **strongly**.

*Proof* : (4.1)(i) shows the existence and uniqueness of $v_h^*$.
Moreover $v_h^*$ satisfies the variational inequality (equivalent
to $(\hat{f}_h)$)

$$
(4.4) \quad (A_h v_h^*, A_h(v_h-v_h^*)) + f_2(A_h v_h) - f_2(A_h v_h^*) \ge \langle b, v_h - v_h^* \rangle \quad \forall v_h \in V_h.
$$
for $v_h = 0$, we obtain

$$
(4.5) \quad |A_h v_h^*|^2 \le f_2(A_h v_h^*) - f_2(0) + \langle b, v_h^* \rangle .
$$

Applying (2.8) and (4.1)(i) to (4.5), we deduce that

$$
\alpha^2 ||v_h^*||^2 \le (M|y_2| + ||b||_{V^*}) ||v_h^*|| + \beta - f_2(0);
$$

hence $||v_h^*|| \le C$.

Thus, we can extract, from $v_h^*$, a subsequence (which
we still call $v_h^*$ ) converging weakly to $w \in V$, as $h \to 0$.
From (4.1)(ii) and (4.4) we have, choosing, for any $v \in V$,
$v_h \in V_h$ as in (4.3) :

$$
\begin{align*}
& |Aw|^2 + f_2(Aw) \le \liminf_{h \to 0} |A_h v_h^*|^2 + f_2(A_h v_h^*) \\
& \qquad \le f_2(Av) + (Aw, Av) + \langle b, w-v \rangle \quad \forall v \in V
\end{align*}
$$

which shows that $w = v^*$ solution of $(\mathcal{P})$. Thus, $v_h^* \to v^*$ in $V$
weakly (without extraction of any subsequence). In fact, it
converges strongly : by (4.1) (i)

$$
(4.8) \quad \alpha' |v_h^{*} - v^{*}| \le |A(v_h^{*})| \le |Av_h^{*} - A_h v_h^{*}| + |A_h v_h^{*} - Av^{*}|
$$

(4.1)(iv) implies that the first term tends to zero. On the
other hand, (4.4) implies that, for any $v \in V$ with $v_h \to v$
in $V$ strongly :

$$
(4.9) \quad
\begin{aligned}
X_h &= |A_h v_h^* - Av^*|^2 = |A_h v_h^*|^2 - 2(A_h v_h^*, Av^*) + |Av^*|^2 \le \\
    &\le (A_h v_h^*, A_h v_h) + f_2 (A_h v_h) - f_2 (A_h v_h^*) - \langle b, v_h - v_h^* \rangle - 2 (A_h v_h^*, Av^*) + |Av^*|^2
\end{aligned}
$$

Since $f_2$ is convex 1.s.c., we know that $f_2(Av^*) \le \liminf f_2(A_h v_h^*)$
and choosing $v_h$ as in (4.3), we have :

$f_2 (A_h v_h) \rightarrow f_2 (Av) .$

Then, taking the lim sup of both members in (4.9), we get :

$$
(4.10) \quad \limsup X_h \le (Av^*, Av) + f_2(Av) - f_2(Av^*) - \langle b, v - v^* \rangle - |Av^*|^2 .
$$

Choosing $v = v^*$, we get :

$$
(4.11) \qquad X_h \to 0
$$

which, by (4.8) implies the strong convergence of $v_h^*$ toward $v^*$.

**REMARK 1** : If $f_2$ has finite values, it is continuous (see [E]).
Then, applying (4.1)(iv), we see that (4.3) holds. ■

**REMARK 2** : Suppose that $f_2 = \delta_K$, the indicator function of some closed convex set $K \subset Y$ :

$$
(4.12) \quad \delta_K(y) = \begin{cases} 0 & \text{if } y \in K \\ +\infty & \text{otherwise} \end{cases}
$$

Let us define $K_0$ and $K_h$ :

$$
(4.13) \quad K_0 = \{v \in V \mid Av \in K\} ,
$$

$$
(4.14) \quad K_h = \{v \in V_h | A_h v_h \in K\} ;
$$
(4.15)

Hypothesis (4.3) can be formulated as in [G2] :

∀v ∈ K₀, ∃vₕ ∈ Kₕ such that vₕ → v in V strongly.

**REMARK 3** : One could think of approximating problem (P) by

$$ \operatorname{Inf}_{v_h \in V_h} \frac{1}{2} |A_h v_h|^2 + f_2(B_h v_h) - \langle b, v_h \rangle $$

where $A_h$ and $B_h$ are two approximations of A, satisfying (4.1). We see that, in order to apply the penalty-duality algorithm (3.4), it is necessary that $A_h = B_h$. Thus, when $f_2$ is the indicator function of some convex K of Y corresponding to a convex K₀ of V as in (4.13), the choice of $A_h$ automatically determines $K_h$ by (4.14), and we have

$$ \delta_{K_h} (v_h) = \delta_K (A_h v_h) . \quad \blacksquare $$

**Finite element approximation.**

We consider now an open set $\Omega$ in $\mathbb{R}^N$ and specify

$$ V = H_0^1(\Omega), \quad Y = [L^2(\Omega)]^N \quad \text{and} $$

(4.16)

$$ (Av)_i = \sum_{j=1}^{N} a_{ij} \frac{\partial v}{\partial x_j} $$

where the coefficients $a_{ij}$ are smooth functions.

We introduce a triangulation $\mathcal{T}_h$ of $\Omega$ consisting of a finite family of "triangles" $T \subset \Omega$, sufficiently regular, and of size less or equal to the positive parameter $h$ [C3]. Let $P_k$ denote the space of polynomials of degree less or equal to $k$ and define the following spaces : $V_h$ is the space of continuous functions on $\bar{\Omega}$ which are equal to polynomials of $P_k$ on each triangle of $\mathcal{T}_h$, $C_h$ is the space of vectorial functions, continuous on each triangle of $\mathcal{T}_h$, and $Z_h \subset C_h$ is the space of vectorial functions, the restriction of which belongs to $(P_{k-1})^N$ on each triangle of $\mathcal{T}_h$,
$$
(4.17) \qquad
\begin{aligned}
V_h &= \{ v_h \in C^0(\bar{\Omega}) \mid |v_h|_T \in P_k \quad \text{for all } T \in \mathcal{T}_h \}, \\
C_h &= \prod_{T \in \mathcal{T}_h} (C^0(T))^N, \\
Z_h &= \{ z_h \in Y \mid z_h|_T \in (P_{k-1})^N \quad \text{for all } T \in \mathcal{T}_h \}.
\end{aligned}
$$

we notice that $A(V_h) \subset C_h$. But $C_h$ is not finite dimensional.
If the functions $a_{ij}$ are constant, $A(V_h) \subset Z_h$, and the
approximation defined by $(V_h, Z_h, A)$ is convergent, according
to theorem 4.1. If the $a_{ij}$ are not constant, $A(V_h) \not\subset Z_h$;
we now need an approximation formula for the numerical integration
of

$$
\int_{\Omega} |Av_h|^2 dx .
$$

Such a formula is also necessary in any case, if the function $f_2$
is defined by an integral over $\Omega$.

This formula is defined by L integration nodes $\hat{a}_\ell$ on
a reference element $\hat{T}$ of $\mathcal{T}_h$ ; then, for any triangle $T \in \mathcal{T}_h$,
we consider the affine mapping $J_T$ on $\Omega$, transforming
$\hat{T}$ into $T$. To the reference nodes $\hat{a}_\ell$ correspond on $T$
the integration nodes $a_{T,\ell} = J_T(\hat{a}_\ell)$. We thus approximate the
integral $\int_{\Omega} \phi(x) \, dx$ by

$$
\sum_{T \in \mathcal{T}_h} \left| \det J_T \right| \hat{\omega}_\ell \phi(a_{T,\ell}) .
$$

We suppose that the formula is exact for the polynomials of $P_{k'}$, $k' \ge 0$
The use of this integration formula can be interpreted as the
introduction of an approximate scalar product $(\cdot, \cdot)_h$ on $C_h$
defined in the following way : for any $y_1, y_2 \in C_h$

$$
(4.18) \qquad (y_1, y_2)_h = \sum_{T \in \mathcal{T}_h} \left| \det J_T \right| \hat{\omega}_{\ell} \sum_{\ell=1}^{L} y_1(a_{T,\ell}) y_2(a_{T,\ell}) .
$$

Since $k' \ge 0$, the formula is exact for the constant functionals :

$$
\sum_{\ell=1}^{L} \hat{\omega}_{\ell} = \text{meas} (\hat{T}) ;
$$

we can partition $\hat{T}$ into $\hat{T}_{\ell}, \ell = 1, \dots, L$ such that

$$
\bigcup_{\ell=1}^{L} \hat{T}_{\ell} = \hat{T} , \quad \hat{a}_{\ell} \in \hat{T}_{\ell} , \text{ measure } (\hat{T}_{\ell}) = \hat{\omega}_{\ell} .
$$

We extend this partitioning on every triangle $T$ of $\mathcal{T}_h$, defining

$$
T_{\ell} = J_{T}(\hat{T}_{\ell}) . \quad \ell = 1, \dots, L .
$$
We are now able to define

$$
(4.19) \qquad Y_h = \{ y_h \in Y \mid |y_h|_{T_\ell} \in (P_0)^N \text{ for } \ell = 1, \dots, L \text{ and all } T \in \mathcal{T}_h \}.
$$

Even though $Y_h \not\subset C_h$, we can easily extend the scalar
product $(\cdot, \cdot)_h$ on $Y_h$ and moreover

$$
\forall y_1, y_2 \in Y_h \quad (y_1, y_2)_h = (y_1, y_2) .
$$

We now construct the mapping $p_h : C_h \to Y_h$, defined for any $z_h \in C_h$ by the restrictions of $p_h(z_h)$ to each $T_\ell$:

$$
(4.20) \quad \forall z_h \in C_h \quad p_h(z_h)|_{T_\ell} = z_h(a_{T,\ell}) \text{ for } \ell = 1, \dots, L \text{ and all } T \in \mathcal{T}_h .
$$

We obviously have

$$
\forall y_h, z_h \in C_h \quad (y_h, z_h)_h = (p_h(y_h), p_h(z_h)) .
$$

We finally define

$$
(4.21) \qquad A_h = P_h \circ A,
$$

which is an operator from $V_h$ (4.17) into $Y_h$ (4.19), such that

$$
(4.22) \quad \forall v_h, w_h \in V_h \quad (Av_h, Aw_h)_h = (A_h v_h, A_h w_h) .
$$

We now show the convergence of this approximation in the
sense of theorem 4.1 . This results from the following proposition :

**PROPOSITION 4.1** : If $k' > 2k - 3$, the operator $A_h$, defined
**in** (4.21) **satisfies** (4.1) and provides us with a converging
**approximation**.

Proof : Following Nédelec [N] , we show that :

$$
(4.23) \qquad \qquad \qquad \qquad \qquad \qquad \qquad \| |Av_h^*|^2 - |Av_h^*|^2 \| \le C_h h^{k'} - (2k-3) \|v_h^*\|^2 ,
$$

and, for all $z_h \in Z_h$,

$$
(4.24) \qquad | (Av_h^*, z_h)_h - (Av_h^*, z_h) | \le C_2 h^{k'} - (2k-3) \|v_h\| \|z_h\| .
$$
Then, if $k' > 2k - 3$, (4.23) implies (4.1)(i) for $h$ sufficiently small.

To prove (4.1)(ii), suppose that $v_h \to v$, we want to prove that $(A_h v_h, y) \to (Av, y) \quad \forall y \in Y$. But $y$ being given, since $\mathcal{T}_h$ is a regular family of triangulations, we can choose a particular $z_h \in z_h$ such that $z_h$ is constant on each $T \in \mathcal{T}_h$ and $z_h \to y$ in $Y$ strongly. (In fact, $z_h \in Y_h$ too). Since $z_h = p_h(z_h)$, $(A_h v_h, z_h) = (Av_h, z_h)_h$. Then, we can write :

$$|(A_h v_h, y) - (Av, y)| \le |(A_h v_h, y - z_h)| + |(Av_h, z_h)_h - (Av_h, z_h)| + |(Av_h, z_h) - (Av, y)|$$

Since $v_h \to v$, $|A_h v_h|$ is bounded and the first term tends to zero. It is the same for the second term applying (4.24), and for the third one, since $A$ is continuous; (4.1)(ii) is then established.

Taking into account (4.23), (4.1)(iii) and (iv) follow immediately. ■
# 5 - APPLICATION TO SOME PROBLEMS IN CONTINUUM MECHANICS.

In the following, $\Omega$ is a regular open bounded set of $\mathbb{R}^N$ of boundary $\Gamma$. ($N = 1, 2, 3$).

## 5.1 - Bingham fluids.

We consider, following MOSSOLOV-MIASNIKOV [M4] and Duvaut-Lions [D], the flow of a visco-plastic fluid in a cylindrical pipe of cross section $\Omega$. Let $u$ denote the component of the velocity of the fluid parallel to the axis, and $b$ be the pressure drop per unit length of the pipe. The problem is to find the function $u(x) \in H_0^1(\Omega)$ minimizing

$$ (5.1) \qquad \frac{\nu}{2} \int_{\Omega} |\mathrm{grad} \, v|^2 \, dx + g \int_{\Omega} |\mathrm{grad} \, v| \, dx - \int_{\Omega} b v \, dx $$

where $v$ is the viscosity of the material in the fluid regions and $g$ is the plasticity threshold.

This is a nondifferentiable problem of the form $(\mathcal{P})$, where $V = H_0^1(\Omega)$, $Y = (L^2(\Omega))^2$, $A = \mathrm{grad}$, $f_1(y) = \frac{\nu}{2} \int_{\Omega} |y(x)|^2 dx$ and $f_2(y) = g \int_{\Omega} |y(x)| dx$. We check that (2.2) to (2.5) hold.

**The algorithm**: Algorithm (3.4) can be made explicit as follows :

By induction, $y^n, \lambda^n$ being given

Step 1 : Find $v^{n+1}$ solution of $\begin{cases} -r\Delta v^{n+1} + \mathrm{div}(ry^n - \lambda^n) = b$ on $\Omega$ \\ v^{n+1}|_{\Gamma} = 0 \end{cases}$

(5.2)

Step 2 : Find $y^{n+1}$ solution of the variational inequality

$$ (5.3) \qquad (v+r) \int_{\Omega} y^{n+1} (y-y^{n+1}) dx + g \int_{\Omega} (|y|-|y^{n+1}|) dx \geq \int_{\Omega} (\lambda^n + r \operatorname{grad} v^{n+1}) (y-y^{n+1}) dx \quad \forall y \in [L^2(\Omega)] $$

(5.4)

Step 3 : Make

$$ \lambda^{n+1} = \lambda^n + \rho (\mathrm{grad} \, v^{n+1} - y^{n+1}) . \quad \blacksquare $$

(5.4)
In fact, (5.3) is equivalent to

$$
(5.5) \qquad (v+r)y^{n+1} - s^n \in \partial f_2(y^{n+1}) \quad \text{a.e. in } \Omega,
$$

where $s^n = \lambda^n + r \operatorname{grad} v^{n+1}$. Since almost everywhere :

$$
(5.6) \qquad \partial f_2 (y) = \begin{cases} \frac{g_y}{|y|} & \text{if } y \neq 0 \\ gB_2 (0; 1) & \text{if } y = 0 \end{cases}, \quad (5.4) \text{ can be solved explicitly ;}
$$

$$
(5.7) \qquad y^{n+1} = \frac{s^n}{v+r} \quad \text{Max} \ (0, 1 - \frac{g}{|s^n|}) \quad \text{a.e. .}
$$

The algorithm reduces to a sequence of Dirichlet problems (5.2), accompanied by the updating formulae (5.7), (5.4) for $y^{n+1}$ and $\lambda^{n+1}$. Since (2.2) to (2.5) hold, the sequence $\{v_n\}$, constructed by this algorithm, converges to the unique solution $v^*$ of (5.1).

**Approximation**: We approximate ($\mathcal{P}$) by ($\mathcal{P}_h$) as described
in Section 4. Since $Y_h$ is a space of piecewise constant functions
formulae (5.5) to (5.7) are still valid, and the algorithm applied
to problem ($\mathcal{P}_h$) still consists of iterations formed by (5.2),
(5.7), (5.4). Proposition 4.1 shows that, if $k' > 2k - 3$,
$v_h^* \to v^*$ as $h \to 0$

5.2 - Minimal hypersurfaces problem.

We want to minimize the area of an hypersurface supported
by a given contour in $\mathbb{R}^{N+1}$. Given $\Omega \subset \mathbb{R}^N$ and a function $g$
defined on $\Gamma$ (the contour being $\{(x,y) \in \Gamma \times \mathbb{R} | y = g(x)\}$ ),
the problem is to find a solution $u$ to :

$$
(5.8) \qquad \begin{array}{ll} \underset{v \in V}{\operatorname{Inf}} & \int_{\Omega} \sqrt{1 + |\operatorname{grad} v|^2} \, dx \ . \\ v|_{\Gamma} = g & \end{array}
$$

With a slight modification due to the presence of the non-
homogeneous boundary condition, this can be formalized as
a problem $(\mathcal{P})$ with $V = W_6^{1,0}(\Omega)$, $Y = [L^4(\Omega)]^N$, $A = \operatorname{grad}$,
$f_1(y) = \int_{\Omega} \sqrt{1 + |y(x)|^2} dx$, $f_2 \equiv 0$, $b = 0$.

$$
(B_2 (0,1) \equiv \{z \in \mathbb{R}^2 \mid |z| \le 1\})
$$
We see that V and Y are non-reflexive Banach spaces, and
we cannot apply the theory we have developed so far. We refer
to [E] for the study of the existence of solutions or generalized
solutions to (5.8).

**Approximation**: Formulated in a finite dimensional space $V_h$,
these difficulties disappear : $f'$ is strongly monotone, there
exists a unique bounded solution $v_h^*$ to $(\mathcal{P}_h)$ and the
algorithm (3.4) is convergent. We use to this purpose an
approximation of the problem in the functional framework,
by piecewise linear finite elements. We refer to [V] for further
details and for convergence results on the approximation.

**Algorithm**: Algorithm (3.4) can be made explicit as previously
We notice that Step 2 decomposes itself into $\mathcal{N}$ variational
inequalities in $y_i^{n+1}$ where $y_i^{n+1}$ is the component of $y^{n+1}$
on the i.th element $T_\ell$. Each of these problems is equivalent
to

$$ (r + \frac{1}{\sqrt{1 + |y_i^{n+1}|^2}}) y_i^{n+1} = s_i^n \quad \begin{array}{ll} \text{for } i = 1, \dots, \mathcal{N} & \\ \text{for } y_i^{n+1} \in \mathbb{R}^N, & \end{array} $$

where $s^n$ is, as before, given by $s^n = \lambda^n + r \operatorname{grad} v^{n+1}$.
This system of $\mathcal{N}$ equations in $N$ variables can be reduced
to a system of $\mathcal{N}$ equations in one variable $\theta \ge 0$

$$ (r + \frac{1}{\sqrt{1 + \theta^2}}) \theta = |s_i^n| \quad \text{for } i = 1, \dots, \mathcal{N}. $$

This nonlinear equation can be solved efficiently by Newton's
method. Finally, the algorithm is as follows :

**Choose** $y^0, \lambda^0$. **By induction,** $y^n, \lambda^n$ **being given**
**Step 1 :** **Find** $v^{n+1}$ **solution in** $V_h$ **of the non-homogeneous Dirichlet problem**

$$ \left\{ \begin{array}{l} -r \Delta v^{n+1} + \operatorname{div}(ry^n - \lambda^n) = 0 \quad \text{on } \Omega \\ v^{n+1}|_{\Gamma} = g \end{array} \right. \qquad (5.9) $$
**Step 2:** For $\ell = 1, \dots, L$ and for $T \in \mathcal{T}_h$, solve the nonlinear equation

$$ (r + \frac{1}{\sqrt{1 + \theta^2}}) \theta = |\sigma|^{n+1} \\ \text{with } \sigma^n = (\lambda^n + r \operatorname{grad} v^{n+1})|_{T_\ell} \\ \text{and set } y|_{T_\ell}^{n+1} = \theta \frac{\sigma^n}{|\sigma|^n}. $$

**Step 3:** Make

$$ \lambda^{n+1} = \lambda^n + \rho(Av^{n+1} - y^{n+1}) $$

### 5.3 - Elasto-plastic torsion of a cylindrical bar.

We consider a cylindrical pipe of section $\Omega$. We are led to the following minimization problem [D]:

$$ (5.10) \quad \inf_{v \in K_0} \frac{1}{2} \int_{\Omega} |\operatorname{grad} v|^2 - \int_{\Omega} bvdx $$

where $K_0$ is the closed convex set :

$$ (5.11) \quad K_0 = \{v \in H_0^1(\Omega) \mid |\operatorname{grad} v| \le 1 \text{ a.e. } \}, \\ \text{Let } Y = (L^2(\Omega))^N; \quad V = H_0^1(\Omega); \quad A = \operatorname{grad} \text{ and} $$

$$ (5.12) \quad K = \{y \in (L^2(\Omega))^N \mid |y| \le 1 \text{ a.e. } \}. $$

We see that, if $\delta_K$ is the indicator function of $K$ (see (4.12)), (5.10) can be written :

$$ (5.13) \quad \inf_{v \in V} \frac{1}{2} \int_{\Omega} |\operatorname{grad} v|^2 dx + \delta_K(\operatorname{grad} v) - \langle b, v \rangle . $$

It is a particular case of $(\mathcal{P})$ with $f_1(y)_0 = \frac{1}{2}|y|^2$ and $f_2 = \delta_K$. We check that (2.3) and (2.5) hold, but $K$ being empty in $Y$, $f_2$ is not continuous in any point, and (2.4) is not satisfied.
We approximate (5.10) by finite elements, and to prove
convergence, we must check that (4.14) holds. This is performed
in [G2] for the particular case of piecewise linear finite
elements (k = 1).

We then use algorithm (3.4) to solve the approximate
problem ($\mathcal{P}_h$). It converges to the unique solution of (5.10).
Step 1 is again a Dirichlet problem while Step 2 is

$$
(4.12) \qquad (1+r)(y^{n+1}, z - y^{n+1}) - (\lambda^n + rAv^{n+1}, z - y^{n+1}) \ge 0 \quad z \in K
$$

which can be solved explicitly in $Y_h$: in each integration node,
we have ;

$$
(4.13) \qquad y^{n+1} = (\lambda^n + rAv^{n+1}) \min \left( \frac{1}{1+r}, \frac{1}{|\lambda^n + rAv^{n+1}|_{\mathbb{R}^N}} \right)
$$

5.4 - Elasto-plastic displacement problem.

Hencky's law results in displacement problem [D] [M1],
which is, again, a particular case of (P) with

. V = {v ∈ (H¹(Ω))N | v = 0 on a part Γ_u of Γ }

. Y = {e_{ij} ∈ L²(Ω) | i,j = 1,...,N | e_{ij} = e_{ji} }

N
with the scalar product (e,τ) : Σ_{i,j=1}^N ∫_Ω τ_{ij} e_{ij} dx

. A = ε "strain operator", where, by definition,

$$
(5.14) \qquad \epsilon_{ij}(v) \equiv \frac{1}{2} \left( \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i} \right)
$$

. Finally $f_2 = 0$, and $f_1$ is defined in the following
way. (We suppose for the sake of simplicity that the matrix
of the elasticity coefficients is equal to identity ;
otherwise, see [M1])

$$
(5.15) \qquad f_1(e) = \sup_{\tau \in K} (\tau, e) - \frac{1}{2} |\tau|^2
$$

where K is the plasticity convex set

$$
(5.16) \qquad K = \{ \tau \in Y \mid \mathcal{F}(\tau) \le 0 \text{ a.e. } \}
$$
and $\bar{F}$ is a convex function on $\mathbb{R}^{N(N+1)/2}$ called plasticity criterion (it may be Von Misès, or Tresca Criterion).

It can be checked that $f_i$ is differentiable and its gradient is $f_i'(e) = P_K(e)$. Hypothesis (2.3) does not hold. In fact, $f_i$ is not even strictly convex.

We refer to [M2] where approximation by piecewise linear finite elements has been studied, and convergence has been proved using the dual problem (Stress problem) which is coercive. Algorithm (3.4) can again be used to solve the approximate problem. Step 2 of the algorithm can be solved explicitly, element by element, and the algorithm is then reduced to the solution of a sequence of elastic problems (linear systems, the matrix of which is fixed). The method is particularly performing in this case, as we can see in the numerical results, for which we refer to [M2], because the elasticity matrix is ill-conditioned.

## 5.5 - Nonlinear Dirichlet problem. [G3]

This is the following problem : find $u \in W^1, \bar{P}(\Omega)$ such that

$$ (5.17) \qquad - \sum_{i=1}^{N} \frac{\partial}{\partial x_i} \left( \left| \nabla u \right|^{P-2} \frac{\partial u}{\partial x_i} \right) = b $$

Once again, algorithm (3.4) can be applied, with $A = \text{grad}$ and $f(y) = \frac{1}{\rho} \int_{\Omega} |y|^P dx$: this has been performed by GLOWINSKI-MARROCCO [G4], to which we refer for details and numerical results.
# 6 - NUMERICAL RESULTS.

Let $\epsilon = \frac{1}{r}$ in the following

## 6.1 - Mossolov problem.

We treat a 2-dimensional example, the exact solution of which is known [G2] ; $\Omega$ is the disk

$$ \Omega = \{(x,y)\} \in \mathbb{R}^2 \mid x^2 + y^2 < R^2 \} . $$

We used a triangulation with 512 triangles and 225 internal nodes (Fig.1) corresponding to piecewise linear finite elements.

The exact solution is constant on a disk of radius $R/2$ and its value is 0.25. The value of the approximate solution is $0.2506 \pm 5.10^{-5}$. With $\epsilon = 1$ and $\rho = 1$, this value is reached already at the second iteration of the algorithm.

However, the convergence of the algorithm is much slower with respect to the multiplier $\lambda$ (in fact, it is not proved). This can be seen on the quantity $|y^{n+1} - Av^{n+1}|$ which tends to zero quite slowly (Fig.2 and 3).

Little improvements can be obtained by a variation of $\rho$ (Fig.2). The choice of $\epsilon$ is more important, even though it is easy : taking $\epsilon < 1$ accelerates slightly the convergence in $\lambda$, but gives a reverse effect on $u$. (This can be seen on Fig. 3, where $n^*$ is the number of iterations from which $v^n = 0.2506 \pm 5.10^{-5}$). Conversely, when $\epsilon > 1$, both convergence in $\lambda$ and $u$ are slower.

A good compromise seems $\epsilon = 1$ and $\rho = 1$. For this value, the computing time is 2 seconds on IBM 370/168. Comparison with duality methods [G2] or direct non-differentiable methods of Davidon type [L1][L2] shows the efficiency of the algorithm (5.2), (5.3), (5.4).

For quadratic approximations with numerical integrations, it would be interesting to compare with BRISTEAU [B4].

<figure><img src="image_4.png" /><figcaption>Fig. 1 - Triangulation chosen for the unit disk [G2] (512 nodes, 225 internal nodes). (dark regions are plastic).</figcaption></figure>
<table><caption>Fig. 2 - Influence of ρ on the speed of convergence.
(ε=1, 225 independent variables)</caption><thead><tr><th>ρ</th><th>Σ |y<sub>T</sub><sup>n</sup>-Av<sub>T</sub><sup>n</sup>|<br>T∈T<sub>h</sub> (1)</th><th>n<br>(2)</th><th>U<br>(3)</th><th>n<sup>*</sup><br>(4)</th></tr></thead><tbody><tr><td>0.33</td><td>0.7</td><td>5</td><td>0.467</td><td></td></tr><tr><td>0.66</td><td>0.1</td><td>5</td><td>0.2725</td><td></td></tr><tr><td>1</td><td>0.03</td><td>5</td><td>0.2506</td><td>2</td></tr><tr><td>1.1</td><td>0.0022</td><td>10</td><td>0.2506</td><td>7</td></tr><tr><td>1.2</td><td>0.0020</td><td>10</td><td>0.2505</td><td>8</td></tr><tr><td>1.3</td><td>0.0018</td><td>10</td><td>0.2505</td><td>9</td></tr><tr><td>1.33</td><td>0.05</td><td>5</td><td>0.2574</td><td></td></tr><tr><td>1.66</td><td>0.5</td><td>5</td><td>0.4093</td><td></td></tr></tbody></table>

(1) Norm chosen for the stopping test.

(2) n: number of performed iterations.

(3) U: value of the solution at the center of the disk.

(4) n*: number of iterations necessary for the convergence of U.

<table><caption>Fig. 3 - Influence of ε on the speed of convergence.
(ρ = 1/ε)</caption><thead><tr><th>ε</th><th>Σ |y<sub>T</sub><sup>n</sup>-Av<sub>T</sub><sup>n</sup>|<br>T∈T<sub>h</sub></th><th>n</th><th>U</th><th>n<sup>*</sup></th></tr></thead><tbody><tr><td>0.5</td><td>4.10<sup>-4</sup></td><td>20</td><td>0.2505</td><td>19</td></tr><tr><td>1</td><td>7.10<sup>-4</sup></td><td>20</td><td>0.2606</td><td>2</td></tr><tr><td>1.5</td><td>1.2 10<sup>-3</sup></td><td>20</td><td>0.2506</td><td>10</td></tr></tbody></table>

## 6.2 - Minimal surface problem.

We treat, as [V], a 2-dimensional example where Ω is the circular crown bounded by two concentrical circles Γ₁ and Γ₂, of radius 1 and 4 (fig. 4). We take g=0 on Γ₂ and g=C (constant) on Γ₁. As noted by [V], for C ≥ 2.07, the solution has some vertical parts and is no more in W¹,¹(Ω).

We used, first, a regular triangulation with 192 triangles and 72 internal nodes. Then, dividing each triangle in 4 triangles, we got a new triangulation with 768 triangles and 336 variables (Fig. 5). At last, as the solution is irregular near Γ₁, we also used a triangulation, isomorphic to the previous one, but refined around Γ₁ (fig. 6).
The stopping test, chosen for the algorithm, is

$$E_n = \sum_{T \in \mathcal{T}_h} |y_T^n - Av_T^n| \le \alpha ;$$

For $\alpha = 10^{-6}$ and 72 variables, the variation of the number of iterations (n) with respect to $\epsilon$ is small (Fig. 7).

For $\epsilon = 1.8$, little improvements can be obtained by a variation of $\rho$ (Fig. 8).

For 336 variables, the number of iterations increases within reasonable proportions for the regular triangulation, but more for the refined one (Fig. 9).

The computing time is 6 seconds, for 72 variables, and 20 seconds, for 336 variables (IBM 370/168), which is reasonable compared to more classical methods like non linear successive over relaxation [V] or non-linear conjugate gradient [L2].

For another finite element approximation, see JOURON [J].

<figure><img src="image_5.png" /><figcaption>Fig. 5 - Normal triangulation (768 triangles, 336 variables)</figcaption></figure>

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 6 - Refined triangulation
(768 triangles, 336 variables)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7 - Variation of n (number of iterations) with respect to ε (α=10⁻⁶, 72 variables)</figcaption></figure>

<table><caption>Fig. 8 - Variation of n with respect to ρ (α=10⁻⁷, 72 variables for ε =1.8)</caption><thead><tr><td>ρ.ε</td><td>n</td></tr></thead><tbody><tr><td>1</td><td>27</td></tr><tr><td>1.1</td><td>26</td></tr><tr><td>1.2</td><td>27</td></tr><tr><td>1.3</td><td>27</td></tr></tbody></table>
<table><caption>Fig. 8 - Results with 336 variables.</caption><thead><tr><th>type of the<br>triangulation</th><th>ε</th><th>C</th><th>E<sub>n</sub></th><th>n</th></tr></thead><tbody><tr><td rowspan="4">normal</td><td>1.8</td><td>2</td><td>10<sup>-7</sup></td><td>43</td></tr><tr><td>2</td><td>2</td><td>10<sup>-7</sup></td><td>39</td></tr><tr><td>2.2</td><td>2</td><td>10<sup>-7</sup></td><td>36 (12")</td></tr><tr><td>2.2</td><td>3</td><td>2.10<sup>-5</sup></td><td>66</td></tr><tr><td rowspan="3">refined</td><td>1.8</td><td>2</td><td>2.10<sup>-4</sup></td><td>46</td></tr><tr><td>2.2</td><td>2</td><td>3.10<sup>-4</sup></td><td>26</td></tr><tr><td>2.2</td><td>3</td><td>10<sup>-7</sup></td><td>140 (45")</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Fig. 10 - Radial section of the solution (minimal surface)</figcaption></figure>

1. 72 variables

2. 336 variables (normal triangulation)

3. 336 variables (refined triangulation)

4. exact [V].
### 6.3 - Elasto-plastic torsion.

We treat, for the sake of simplicity, a one-dimensional example: $\Omega = ]0,1[$ and a discretization by piecewise linear finite elements on a regular mesh, the step of which is $h$.

When $b$ = constant, we have an exact solution which is parabolic on the elastic zone $]\frac{1}{2} - \frac{1}{b}, \frac{1}{2} + \frac{1}{b}]$ and linear elsewhere on $\Omega$.

Numerical results show that for $b=10$, starting from $y^0=0$ and $\lambda^0=0$, the best choice for $\epsilon$ is the biggest one, limited, however, by the appearance of round-off errors when $\epsilon$ is too big (Fig. 11). For $\rho$, a good choice seems to be $\rho=\frac{1}{5}$ (Fig. 12).

In that case, results are quite independent on $b$, even when $b<2$, which corresponds to a quadratic problem.

A theoretical study of the speed of convergence would show that starting from $\lambda^0=0$, it is better to take $\epsilon$ very big, this resulting from the fact that $\lambda^0 - \lambda^* \approx 0$ is eigenvector of $A(A'A)^{-1}A'$.

On the contrary, if we start from $\lambda^0=(1,1,...,1)$, when $b$ is small (quadratic problem), a good choice is $\epsilon=1$ (Fig. 13).

<table><caption>Fig. 11 - Choice of $\epsilon$ for the elasto-plastic torsion with $h = \frac{1}{20^e}$, $b = 10$ and $\lambda^0 = 0$.</caption><thead><tr><th>ε</th><th>n<br>nb of iterations</th><th>E<sub>n</sub><br>stopping test</th></tr></thead><tbody><tr><td>2</td><td>10</td><td>10<sup>-4</sup></td></tr><tr><td>3</td><td>10</td><td>10<sup>-5</sup></td></tr><tr><td>10</td><td>10</td><td>7.10<sup>-9</sup></td></tr><tr><td>10<sup>2</sup></td><td>7</td><td>10<sup>-10</sup></td></tr><tr><td>10<sup>3</sup></td><td>5</td><td>10<sup>-10</sup></td></tr><tr><td>10<sup>4</sup></td><td>5</td><td>10<sup>-11</sup></td></tr><tr><td>10<sup>6</sup></td><td>4</td><td>10<sup>-13</sup></td></tr><tr><td>10<sup>10</sup></td><td>4</td><td>10<sup>-13</sup></td></tr><tr><td>10<sup>12</sup></td><td>4</td><td>10<sup>-11</sup></td></tr><tr><td>10<sup>20</sup></td><td>10 (round of errors)</td><td>10<sup>3</sup></td></tr></tbody></table>
<table><caption>Fig. 12 - Choice of ρ with h=1/10$^e$, b=2, λ⁰=0, ε=20</caption><thead><tr><th>ρ</th><th>n</th><th>E<sub>n</sub></th><th>U(1/2) (Uex=0.25)</th></tr></thead><tbody><tr><td>0.5</td><td>10</td><td>0.2</td><td>0.26</td></tr><tr><td>0.9</td><td>10</td><td>10<sup>-7</sup></td><td>0.25</td></tr><tr><td>1</td><td>8</td><td>10<sup>-8</sup></td><td>0.25</td></tr><tr><td>1.1</td><td>10</td><td>10<sup>-7</sup></td><td>0.25</td></tr><tr><td>1.8</td><td>10</td><td>12</td><td>-0.35</td></tr><tr><td>2</td><td>10</td><td>90.9</td><td>-4.29</td></tr></tbody></table>

<table><caption>Fig. 13 - Choice of ε when λ⁰=(1,1,...,1)
( h=1/10$^e$, b=2)</caption><thead><tr><th>ε</th><th>n</th><th>E<sub>n</sub></th></tr></thead><tbody><tr><td>10<sup>-3</sup></td><td>10</td><td>Cv. very slow</td></tr><tr><td>0.7</td><td>10</td><td>1.6 10<sup>-2</sup></td></tr><tr><td>1</td><td>10</td><td>9. 10<sup>-3</sup></td></tr><tr><td>10<sup>3</sup></td><td>10</td><td>7.5 ) no Cv.</td></tr><tr><td>10<sup>12</sup></td><td>10</td><td>7.5 ) no Cv.</td></tr></tbody></table>

## CONCLUSION

From the numerical point of view, we can summarize the properties of the penalty-duality algorithm.

It is very easy to implement, the main part of the algorithm being common to every problem and consisting of the solution of a linear system, the matrix of which can be factorized once for all, at the beginning. The only change, from one problem to another, is relative to step 2 of the algorithm, which consists in solving $η$ non linear equations, generally in one variable.

To select $ρ$ and $ε$, the best is to choose $ρ = 1/ε$, and to try $ε=1$ or 2 at the beginning. Then, other values of $ε$ may be tried.

It is better to start from $y⁰=0$ and $λ⁰=0$, at least at the beginning.

The algorithm has revealed to be very efficient, and the comparison with other classical methods [G2],[L2] shows that it is particularly useful when the matrix $A'A$ is ill-conditioned.
REFERENCES :

[A] K.J. ARROW, L.HURWICZ and U. UZAWA : *Studies in nonlinear programming* - Stanford University Press (1968).

[B1] D.P. BERTSEKAS : "On the method of multipliers for convex programming" - *Math. Prog.*, to appear.

[B2] D.P. BERTSEKAS : "Non-differentiable optimization via approximation" - To appear.

[B3] D.P. BERTSEKAS and S.K. MITTER : "A descent numerical method for optimization problem with non-differentiable cost functionals" - SIAM J. on Control, 11 (1973), 637-652.

[B4] M.O. BRISTEAU : Thèse de 3ème cycle - Paris - (1975).

[C1] J. CEA and R. GLOWINSKI : "Sur des méthodes d'optimisation par relaxation" - R.A.I.R.O. - Dec. 1973, R3, 5-32.

[C2] J. CEA, R. GLOWINSKI and J.C. NEDELEC : "Minimisation de fonctionnelles non différentiables" - Conférence on Numerical Analysis, Dundee - ed. by MORRIS : Lectures Notes in Mathematics-Springer Verlag (1971), 19-38.

[C3] P.G. CIARLET and P.A. RAVIART : "The combined effect of curved boundaries and numerical integration in isoparametric finite element method" - In *The Mathematical foundations of the finite element method*, edited by A.K. AZIZ - Academic Press (1972)-409-474.

[D] G. DUVAUT and J.L. LIONS : *Les inéquations en Mécanique et en Physique* - Dunod - Paris - (1972).

[E] I. EKELAND and R. TEMAM : *Analyse convexe et problèmes variationnels* - Dunod - Paris - (1971).

[F] M. FORTIN : "Minimization of some non-differentiable functionals by the Augmented Lagrangian method of HESTENES and POWELL" - To appear.
[G1] D. GEOFFRION : "Duality in nonlinear programming. A simplified applications-oriented development" - SIAM Rev. 13 (1971), 1-37.

[G2] R. GLOWINSKI, J.L. LIONS, R. TREMOLIERES : *Analyse numérique des inéquations variationnelles* - Dunod - Paris - To appear.

[G3] R. GLOWINSKI and A. MARROCCO : "Sur l'approximation, par éléments finis d'ordre 1, et la résolution, par pénalisation-dualité, d'une classe de problèmes de Dirichlet non-linéaires" - C.R. Acad. Sc. Paris, T. 278, série A (1974), 1649-1652.

[G4] R. GLOWINSKI and A. MARROCCO : "On the solution of a class of nonlinear Dirichlet problems by a penalty-duality method and finite element of order one" - To appear in the Proceedings 6th Congress IFIP-TC7 (Optimization) - Springer Verlag.

[G5] R. GLOWINSKI and A. MARROCCO : Rapport N° 74023 - Laboratoire d'Analyse Numérique, L.A. 189 - Université de Paris VI - Paris France (to appear in R.A.I.R.O.).

[H] M. HESTENES : "Multiplier and gradient methods" - J.O.T.A. 4, N° 5 (1969), 303-320.

[J] C. JOURON : "Resolution numérique du problème des surfaces minima" - To appear in Arch. Rat. Mech. Anal.

[K] B.V. KORT and D. BERTSEKAS : "Combined primal-dual and penalty methods for convex programming" - SIAM J. on Control, To appear.

[L1] C. LEMARECHAL : "An extension of Davidon methods to non-differentiable criterions" - To appear in Math. Prog.

[L2] C. LEMARECHAL : Thèse d'Etat - Paris - To appear.

[L3] J.L. LIONS : *Quelques méthodes de résolution des problèmes non-linéaires* - Dunod - Paris (1969).

[L4] D.G. LUENBERGER : *Introduction to linear and nonlinear programming* Addison-Wesley, Reading, Mass. (1973).
[ L5 ] D.G. LUENBERGER : "Control problems with kinks" - I.E.E.E. Trans. on A.C. 15, 570-575.

[ M1 ] B. MERCIER : "Approximation par éléments finis, et résolution, par un algorithme de pénalisation-dualité, d'un problème d'elasto-plasticité" - C.R. Acad. Sc. Paris, T 280, série A (1975), 287-290.

[ M2 ] B. MERCIER : "A finite element approximation of elasto-plastic problem" - Proceedings of International Computing Symposium, Juan-les-Pins - Antibes, June 2-4 1975.

[ M3 ] J.J. MOREAU : "Proximité et dualité dans un espace hilbertien"- Bull. Soc. Math. France, 93 (1965), 273-299.

[ M4 ] P.P. MOSOLOV and V.P. MIASNIKOV : "Variational methods in the theory of the fluidity of a visco-plastic medium" - P.M.M.; 29, (1965), 468-492.

[ N ] J.C. NEDELEC : "Approximation par éléments finis des équations de Riccati" - Université de Rennes - Journées Eléments Finis - May 13-14, 1974.

[ P ] M.J.D. POWELL : "A method for the nonlinear constraints in minimization problems" - Optimization, edited by R. FLETCHER, Academic Press (1972).

[ R1 ] R.T. ROCKAFELLAR : "Augmented Lagrange multipliers functions and duality in non-convex programming" - SIAM J. on Control, 12 (1974), 268-285.

[ R2 ] R.T. ROCKAFELLAR : "A dual approach to solving nonlinear programming problems by unconstrained optimization" - Math. Prog. 5 (1973), 354-373.

[ R3 ] R.T. ROCKAFELLAR : Convex Analysis - Princeton (1970).

[ R4 ] R.T. ROCKAFELLAR : "Solving a nonlinear programming problem by way of a dual problem" - Symposia Mathematica - (to appear)
G. STRANG : "Variational crimes" - from *The mathematical foundations of the finite element method*, edited by A.K. AZIZ, Academic Press, (1972), 689-710.

R. TREMOLLIERES : Thèse d'Etat - Paris (1973).

A. VANENDE-MORREEUW : Thèse de 3ème cycle - Paris (1972).

ZIENKIEWICZ : *The finite element method in engineering science* - Mac Graw Hill (1970).
Special symbols used :

Page 1 and following pages : P capital manuscript letter P

Page 2 and following pages : L capital manuscript letter L

Page 2 and following pages : D capital manuscript letter D

Page 5 and following pages : N capital manuscript letter N

Page 5 and following pages : T capital manuscript letter T