Analysis of the Inexact Uzawa Algorithm for Saddle Point Problems.
In this paper, we consider the so-called "inexact Uzawa" algorithm for iteratively solving linear block saddle point problems.  Such saddle point problems arise, for example, in finite element and finite difference discretizations of Stokes equations, the equations of elasticity, and mixed finite element discretization of second-order problems.  We consider both the linear and nonlinear variants of the inexact Uzawa iteration. We show that the linear method always converges as long as the preconditioners defining the algorithm are properly scaled.  Bounds for the rate of convergence are provided in terms of the rate of convergence for the preconditioned Uzawa algorithm and the reduction factor corresponding to the preconditioner for the upper left-hand block.  In the case of nonlinear iteration, the inexact Uzawa algorithm is shown to converge provided that the nonlinear process approximating the inverse of the upper left-hand block is of sufficient accuracy.  Bounds for the nonlinear iteration are given in terms of this accuracy parameter and the rate of convergence of the preconditioned linear Uzawa algorithm.  Applications to the Stokes equations and mixed finite element discretization of second-order elliptic problems are discussed and, finally, the results of numerical experiments involving the algorithms are presented.
Introduction
. This paper provides a new analysis for the inexact Uzawa
method applied to the solution of saddle point systems which arise in the discretization
of various systems of partial differential equations. Such systems typically are obtained
when "multiplier" or mixed discretization techniques are employed. Examples of these
include the discrete equations which result from approximation of elasticity problems,
Stokes equations and sometimes linearizations of Navier-Stokes equations [4], [14], [15],
[16]. In addition, these systems result from Lagrange multiplier [2], [3], [24] and mixed
formulations of second order elliptic problems [8], [21], [24].
We shall consider iterative solution of an abstract saddle point problem. Let H 1
and H 2 be finite dimensional Hilbert spaces with inner products which we shall denote
by (\Delta; \Delta). There is no ambiguity even though we use the same notation for the inner
products on both of these spaces since the particular inner product will be identified by
This manuscript has been authored under contract number DE-AC02-76CH00016 with the U.S.
Department of Energy. Accordingly, the U.S. Government retains a non-exclusive, royalty-free license
to publish or reproduce the published form of this contribution, or allow others to do so, for U.S.
Government purposes. This work was also supported in part under the National Science Foundation
Grant No. DMS-9007185 and by the U.S. Army Research Office through the Mathematical Sciences
Institute, Cornell University.
y Department of Mathematics, Cornell University, Ithaca, New York 14853 and Department of Math-
ematics, Texas A&M University, College Station,
z Brookhaven National Laboratory, Upton, NY 11973.
x Department of Mathematics, Texas A&M University, College Station,
the type of functions appearing. We consider the abstract saddle point problem:
!/
Y
F
G
are given and are the unknowns. Here
assumed to be a linear, symmetric, and positive definite operator. In
addition, the linear is the adjoint of B: H 1 7! H 2 . Applying block
elimination to (1.1) yeilds
G:
is symmetric and nonnegative and a straightforward computation
shows that
(AU;U)
Consequently, a necessary and sufficient condition for the unique solvability of (1.1) is
that the Ladyzhenskaya-Babu-ska-Brezzi condition hold, i.e.
sup
(AU;U)
for some positive number c 0 . Here k \Delta k denotes the norm in the space H 2 (or H 1
corresponding to the inner product (\Delta; \Delta).
One could iteratively solve (1.2) for Y by conjugate gradient (or preconditioned
conjugate gradient) iteration [12]. Then X is obtained by
Uzawa method [1] is a particular implementation of a linear iterative method for solving
(1.2). One common problem with the methods just described is that they require
the evaluation of the action of the operator A \Gamma1 in each step of the iteration. For
many applications, this operation is expensive and is also implemented as an iteration.
The inexact Uzawa methods replace the exact inverse in the Uzawa algorithm by an
"incomplete" or "approximate" evaluation of A \Gamma1 . These algorithms are defined in
Section 2 and 4. They were also studied in [11].
There are other general iterative techniques for solving saddle point problems of the
form of (1.1), e.g., [5], [26]. In [5], a preconditioner for A is introduced and system (1.1)
is reformulated as a well conditioned symmetric and positive definite algebraic system
which may be solved efficiently by applying the conjugate gradient algorithm. In [26],
the authors consider the convergence properties when the minimal residual algorithm
is applied to a more direct preconditioned reformulation of (1.1). Both of the above
mentioned techniques incorporate preconditioning and avoid the inversion of A.
There is also a variety of application specific techniques that depend strongly on
the particular approximation spaces, geometry of the domain etc. In the case of the
mixed approximation of second order problems, those include domain decomposition
techniques [17], a reduction technique involving the use of additional Lagrange multipliers
[9], as well as an indefinite preconditioner [13].
The inexact Uzawa algorithms are of interest because they are simple and have
minimal computer memory requirements. This could be important in large scale scientific
applications implemented for todays computing architectures. In addition, an
Uzawa algorithm implemented as a double iteration can be transformed trivially into
an inexact Uzawa algorithm. It is not surprising that the inexact Uzawa methods are
widely used in the engineering community.
In this paper we present new estimates for the inexact Uzawa algorithm both in
the linear and nonlinear case. In the former case, the evaluation of A \Gamma1 is replaced by
the inverse of a linear preconditioner. Theorem 3.1 shows that the resulting algorithm
always converges and gives bounds on the rate of convergence provided that the preconditioner
is properly scaled. The inexact Uzawa algorithm in the nonlinear case replaces
the evaluation of A \Gamma1 by some approximate nonlinear process. Theorem 4.1 shows that
the resulting algorithm converges provided that the nonlinear approximation to A \Gamma1 is
suitably accurate. More restrictive results for variants of the inexact Uzawa algorithms
have already appeared in the literature [11], [23].
The outline of the remainder of the paper is as follows. In Section 2, we define
and motivate the linear version of the inexact Uzawa algorithm. Section 3 provides an
analysis of this algorithm. In Section 4, the nonlinear version of the inexact Uzawa
algorithm is defined and analyzed. Section 5 discusses a model application to the
Stokes problem while Section 6 considers a model application to a mixed finite element
discretization of a second order problem boundary value problem. Finally, the results of
numerical experiments involving the inexact Uzawa algorithms are given in Section 7.
A comparison with some other methods is presented as well.
2. The abstract inexact Uzawa algorithm. In this section, we define
the inexact Uzawa method when linear preconditioners are used. This algorithm is
motivated by first considering the Uzawa iteration [1] which can be defined as follows.
Algorithm 2.1 (Uzawa). For the sequence
with - a given real number.
be the iteration error generated by the above method. It is easy
to show that
Let c 1 denote the largest eigenvalue of BA converges to Y if - is chosen
such that
In this case, X i and Y i converge respectively to X and Y with a rate of convergence
bounded by ae i .
One problem with the above method is that it may converge slowly if BA is not
well conditioned. Thus, it is natural to introduce a preconditioner QB We
assume that QB is linear, symmetric and positive definite and define the preconditioned
Uzawa algorithm as follows.
Algorithm 2.2 (Preconditioned Uzawa). For
the sequence f(X
For convenience of notation, we have absorbed the parameter - into the preconditioner
QB . Accordingly, we assume that QB is scaled so that
Note that since QB is positive definite, it follows that
holds for some fl in the interval [0; 1). In practice, effective preconditioners satisfy (2.4)
with bounded away from one.
generated by (2.2) then
is symmetric with respect to the inner product
Let k\Deltak QB denote the corresponding norm
Then by (2.3) and (2.4),
Here and in the sequel, for a symmetric and positive definite linear operator L on H j ,
will denote the norm (L\Delta; \Delta) 1=2 .
One problem with the above algorithms is that they require the computation of the
action of the operator A \Gamma1 at each step of the iteration. For many of the applications,
this is an expensive operation which is also done iteratively. This leads to a two level
iteration, an inner iteration for computing the action of A \Gamma1 coupled with the outer
Uzawa iteration (2.1) or (2.2). The inexact Uzawa method replaces the action of A \Gamma1
by a preconditioner. A preconditioner QA is a linear operator which is
symmetric and positive definite. In practice, good preconditioners are relatively cheap
to invert. For example, the computational cost for one evaluation of Q
A should be
comparable with the cost of evaluating the action of A (not A \Gamma1 ). The inexact Uzawa
algorithm is then given as follows (this algorithm was also studied in [11]).
Algorithm 2.3 (Inexact Uzawa). For the
sequence
A
One step of the inexact Uzawa algorithm involves an evaluation of each of the op-
erators,
A and Q
B . In contrast to Krylov space minimization algorithms
such as conjugate residual, there are no discrete inner products involved in the itera-
tion. This makes this algorithm very well suited for implementation on contemporary
massively parallel computer architectures.
3. Analysis of the inexact Uzawa algorithm. In this section, we investigate
the stability and convergence rate of the inexact Uzawa algorithm defined above. The
main theorem will show that the inexact Uzawa algorithm will always converge provided
that the preconditioners are properly scaled. By this we mean that (2.3) holds and that
for all W 2 H 1 with W 6= 0. The strict inequality above will be replaced by
in a subsequent corollary.
Bounds for the rates of iterative convergence will be provided in terms two natural
parameters. The first parameter has already been defined and is the convergence factor
(see (2.4)) for the preconditioned Uzawa algorithm. The second parameter is the rate
ffi at which the preconditioned iteration
converges to the solution of
A A)E A
A A is a symmetric operator with respect to the inner product (QA \Delta; \Delta) and
hence the convergence rate ffi is the largest eigenvalue of I
A A. Alternatively, ffi is
the smallest number for which the inequality
is satisfied. It will sometimes be convenient to rewrite (3.3) as
The first convergence estimate will be provided in terms of a norm on H 1 \Theta H 2
which we shall now define. Consider the bilinear form on H 1 \Theta H 2 given by
"/
U
R
By (3.1), [\Delta; \Delta] generates a norm on H 1 \Theta H 2 which we shall denote by
We can now state the main result of this section.
Theorem 1. ssume that (2.3) and (3.1) hold and that fl and ffi satisfy (2.4) and
respectively. Let fX; Y g be the solution pair for (1.1), fX be defined by the
inexact Uzawa algorithm and set
Then, for
where
Remark 3.1. It is elementary to see that
Thus the inexact Uzawa method converges if (2.3) and (3.1) hold. As expected, the
convergence rate deteriorates as either fl or ffi approach one. In addition, if ffi tends to
zero (and thus, QA tends to (defined by (3.7)) tends to fl, the convergence
rate of the preconditioned Uzawa algorithm.
Proof.[Theorem 1] We first derive a relationship between the errors e i+1 and e i . The
components of the corresponding errors are denoted by E X
From (1.1) and (2.5) we see that the errors satisfy the recurrence
A A
Replacing
i+1 in the second equation with its expression from the first givesB @
A
A A
iC A
This can be rewritten as
The proof of the theorem will be complete if we can show that the operator norm
is bounded by ae given by (3.7).
The operator M can be written in the form
A
A A
It is straightforward to check that both E and M 1 are symmetric in the [\Delta; \Delta]-inner
product. Let M   denote the adjoint of M with respect to [\Delta; \Delta]. Then we have
and
2:
Consequently,
[M   Mx;x]
Therefore, to estimate the norm of M, it suffices to bound the spectrum oe(M 1 ) of M 1 .
Since M 1 is symmetric with respect to the [\Delta; \Delta] inner product, its eigenvalues are real.
We shall bound the positive and negative eigenvalues of M 1 separately.
We first provide a bound for the positive eigenvalues of M 1 . The operator I \Gamma Q
A A
is symmetric with respect to the inner product ((QA \Gamma A)\Delta; \Delta). Moreover, it follows from
(3.1) that it is positive definite and its positive square root is well defined. Let
A
It follows from (3.1) that D is invertible and from (3.3) that
\GammaffiI
A
A
The largest eigenvalue -m of M 1 satisfies
We used (3.11) for the last inequality above. Since both D and M 1 are symmetric with
respect to [\Delta; \Delta], it follows that N is also. Consequently, -m is bounded by the largest
eigenvalue of N .
Let - be a nonnegative eigenvalue of N with corresponding eigenvector f/
i.e.,
Eliminating / 1 in the above equations gives
\Gamma-L   L/
and hence
By (3.3) and (2.4), it follows that
is nonnegative, we see from the first equation in (3.13) that if /
then Consequently, / 2 is not equal to zero. Thus, from (3.14) and (3.15), we
get
from which it follows that - ae where ae is given by (3.7). This provides the desired
bound for the positive eigenvalues of M 1 .
We next estimate the negative eigenvalues of M 1 . Let - be a negative eigenvalue
of M 1 with corresponding eigenvector (/
A A
A A
The first equation in (3.16) together with (2.4) imply that if /
Consequently, any eigenvector must have a nonzero component / 1 .
Multiplying the first equation of (3.16) by Q
B B from the left and adding it to the
second one yields
Substituting (3.17) into the first equation of (3.16) and taking an inner product with
which we rewrite as
For any V 2 H 1 ,
For the last inequality above we used (2.3). Applying (3.19) to the left hand side of
and (3.4) on the right hand side of (3.18) gives
or
This implies that - \Gamma
since / 1 is nonzero. It is elementary to check that
where ae is defined by (3.7). This completes the proof of the theorem.
The proof of Theorem 3.1 depended on (3.1) so that the inner product [\Delta; \Delta] induced
a norm. The next result shows that the inexact Uzawa method converges even when
only (3.2) is assumed. It also provides an estimate for the error E X
in a more
natural norm.
Corollary 3.1. Assume that (2.3) and (3.2) hold and that fl and ffi satisfy (2.4)
and (3.3) respectively. Let fX; Y g be the solution pair for (1.1), let fX be defined
by the inexact Uzawa algorithm and set E X
where ae is given by (3.7). In addition,
The above inequalities hold for
Proof. Taking the (\Delta; \Delta)-inner product of the first equation of (3.8) with QA e X
applying the Schwarz inequality, and (2.3) gives
QB
Thus, applying (3.2) gives
and by (3.3),
for consider the sequence of iterates fX ffl;i ; X ffl;i g
generated by the inexact Uzawa algorithm which replaces QA in the first equation of
(2.5) by Q A;ffl . Applying Theorem 3.1 gives that the error
satisfies
[je ffl;i j] ffl - ae i
"/
U
R
and
Clearly,
ffl;i
Inequality (3.20) results from combining (3.24) and (3.25) and taking the limit as ffl
tends to zero.
In a similar manner we prove (3.21). Taking the limit in (3.24) as ffl tends to zero
gives
[je
Combining (3.22) and (3.26) gives (3.21) and completes the proof of the corollary.
Remark 3.2. More restrictive convergence results (in these norms) were obtained
by Queck [23]. He proved a convergence result which required stronger conditions with
respect to the scaling of QA and QB . In particular, there are cases which fail to satisfy
the hypothesis of the theory of [23] yet convergence is guaranteed by the corollary
above. In addition, there are many cases when the convergence estimates given above
are substantially better than those of [23].
4. Analysis of the nonlinear inexact Uzawa algorithm. As was pointed
out in Section 2, the Uzawa algorithm is often implemented as a two level iterative
process, an inner iteration for computing A \Gamma1 coupled with the outer Uzawa iteration
or (2.2). In this section we investigate the stability and convergence rate of
an abstract inexact Uzawa algorithm where the computation of the action of A \Gamma1 is
replaced with that of an approximation to A \Gamma1 which results from applying a nonlinear
iterative process for inverting A. Two examples of such approximations come from
defining the approximate inverse by a preconditioned conjugate gradient iteration or
the operator which results from the application of a multigrid cycling algorithm with a
nonlinear smoother.
The nonlinear approximate inverse is described as a map
\Psi(OE) is an "approximation" to the solution - of
We shall assume that our approximation satisfies
for some ffi ! 1. As will be seen below, (4.2) is a reasonable assumption which is
satisfied by the approximate inverse associated with the preconditioned conjugate gradient
algorithm. It also can be shown that (4.2) holds under reasonable assumptions
for approximate inverses defined by one sweep of a multigrid algorithm with conjugate
gradient smoothing.
Perhaps the most natural example of a nonlinear approximate inverse is defined in
terms of the preconditioned conjugate gradient procedure [22]. Let QA be a symmetric
and positive definite operator on H 1 and consider applying n steps of the conjugate
gradient algorithm preconditioned by QA to solve the problem (4.1) with a zero starting
iterate. We define is the resulting approximation to -. The
preconditioned conjugate gradient algorithm (PCG) provides the best approximation
(with respect to the norm corresponding to the (A\Delta; \Delta)-inner product) to the solution -
in the space
A
It is well known that this implies [6]
where
A A) is the condition number of
A A. Note that ffi n is a decreasing function of n and ffi 1 is less than one. Thus, (4.2)
holds in the PCG example. In fact,
A
A
tends to zero as n tends to infinity, it is possible to make ffi n as small as we
want by taking a suitably large number PCG iterations.
The variant of the inexact Uzawa algorithm we investigate in this section is defined
as follows.
Algorithm 4.1 (Nonlinear Uzawa). For the
sequence
Clearly, (4.4) reduces to the preconditioned Uzawa algorithm (2.2) if \Psi(f
reduces to the inexact Uzawa algorithm if \Psi is a linear operator
A .
We will provide bounds for the rate of convergence for the above algorithm in terms
of two parameters, the convergence factor fl for the preconditioned Uzawa algorithm
defined in (2.4) and the parameter ffi of (4.2). The main result of this section provides a
sufficient condition on ffi for convergence of the nonlinear Uzawa algorithm and bounds
for the resulting rate of convergence.
Theorem 2. ssume that (2.3) and (4.2) hold and that fl satisfies (2.4). Let fX; Y g
be the solution pair for (1.1) and fX be defined by the nonlinear Uzawa algorithm
(4.4). Then X i and Y i converge to X and Y respectively if
In this case the following inequalities hold:
and
where
Remark 4.1. The result of Theorem 2 is somewhat weaker than the results obtained
in Section 3 for the linear case due to the threshold condition (4.5) on ffi. In the case of
PCG, it is possible to take sufficiently many iterations n so that (4.5) holds for any fixed
A A). In applications involving partial differential equations, fl and -(Q \Gamma1
A
may depend on the discretization parameter h. If, however, -(Q
A can be bounded
and fl can be bounded away from one independently of h then by Theorem 2, a fixed
number (independent of h) of iterations of PCG are sufficient to guarantee convergence
of the nonlinear Uzawa algorithm.
Remark 4.2. An analysis of (4.4) is given in [10] and [11] in the case of applications
to Stokes problems. The sufficient condition for convergence derived there is that the
iterate
A
where - is independent of the mesh size. The above norms are not natural for procedures
such as PCG and multigrid with nonlinear smoothing. PCG does not give rise to
monotone error behavior in the norm k\Deltak even though convergence is guaranteed by the
canonical bound (4.3),
and equivalence of norms in finite dimensional spaces. Such norm equivalences depend
on the mesh parameter h. A second problem with the requirement (4.9) is that the
norm on the right hand side converges to zero as X i converges to the solution X. This
means that even thought - is fixed independent of h, considerably more iterations of
PCG may be required to satisfy (4.9) as the approximate solution converges.
Proof.[Theorem 2] We start by deriving norm inequalities involving the errors E X
. As in (3.8),
The first equation above can be rewritten
It follows from the triangle inequality, (4.2) and (2.3) that
Using (4.11) in the second equation of (4.10), we obtain
is a symmetric operator in the ! \Delta; \Delta ?-inner product, it follows
from (2.4) that
Thus, by the triangle inequality, (2.3), (3.19) and (4.2),
Let us adopt the notation
for vectors of nonnegative numbers x Repeated
application of (4.12) and (4.13) gives@
where M is given by
We consider two dimensional Euclidean space with the inner product
$/
!%
A trivial computation shows that M is symmetric with respect to the b\Delta; \Deltac-inner prod-
uct. It follows from (4.14) that
where ae is the norm of the matrix M with respect to the b\Delta; \Deltac-inner product. Since
M is symmetric in this inner product, its norm is bounded by its spectral radius. The
eigenvalues of M are the roots of
It is elementary to see that the spectral radius of M is equal to its positive eigenvalue
which is given by (4.8).
Examining the expression for ae given by (4.8) we see that ae is an increasing function
of ffi for any fixed fl 2 [0; 1]. Moreover,
This completes the proof of the (4.6).
To prove (4.7) we apply the arithmetic-geometric mean inequality to (4.12) and get
for any positive j,
A
A
QB
Inequality (4.7) follows taking applying (4.6). This completes the proof
of the theorem.
5. Application to a Stokes problem . In this section we consider an application
of the theory developed in the previous sections to solving indefinite systems
of linear equations arising from finite element approximations of the Stokes equations.
For simplicity we restrict our discussion to the following model problem: Find u and p
such that
Z\Omega
where\Omega is the unit cube in R d , d=2, 3, \Delta is the componentwise Laplace operator, u is
a vector valued function representing the velocity, and the pressure p is a scalar func-
tion. Generalizations to domains with more complex geometry and variable coefficients
equations are possible.
Let L 2(\Omega\Gamma be the set of functions in L
2(\Omega\Gamma with zero mean value
on\Omega and H
denote the Sobolev space of order one
on\Omega (cf., [18], [20]). The space H 1
consists
of those functions
in\Omega whose traces vanish on @
\Omega\Gamma the boundary of \Omega\Gamma Also, (H 1
d
will denote the product space consisting of vector valued functions with each vector
component in H 1
In order to derive the weak formulation of (5.1) we multiply the first two equations
of (5.1) by functions in (H 1
d and L 2
respectively and integrate
over\Omega to get
Here (\Delta; \Delta) is the L 2
inner product and D(\Delta; \Delta) denotes the vector Dirichlet form for
vector functions
on\Omega defined by
d
Z\Omega
We next identify approximation subspaces of (H 1
d and L 2
0(\Omega\Gamma4 In order to
avoid unnecessary complexity of the presentation only a two dimensional example will
be considered. The discussion here is very closely related to the examples given in [4]
and [5] where additional comments and other applications can be found. We
Fig. 5.1. The square mesh used for ~
the support
(shaded) and values for a typical OE ij
into 2n \Theta 2n square shaped elements, where n is a positive integer and define
Each of the square elements is
further partitioned into two triangles by connecting the lower right corner to the upper
left corner. Let S h be the space of functions that vanish on
@\Omega and are continuous and
piecewise linear with respect to the triangulation thus defined. We set H 1
. The choice of H 2 is motivated by the observation [19] that the space ~
H 2 of
functions that are piecewise constant with respect to the square elements together with
defined above form an unstable pair of approximation spaces. This means that
the functions from H 1 \Theta ~
H 2 do not satisfy (1.4) with a constant c 0 independent of the
discretization parameter h. To overcome this problem, one may consider a smaller space
defined as follows. Let j kl for k; 2n be the function that is one on the square
element [x
by
(see

Figure

5.1). The space H 2 is then defined by
The now satisfies (1.4) with a constant c 0 independent of h [19]. Moreover,
the exclusion of the functions OE i;j does not change the order of approximation for the
space since the H 2 still contains the piecewise constant functions of size 2h.
The approximation to the solution of (5.2) is defined as the unique pair (X; Y
(r
Obviously, (5.3) is a system of linear equations whose unique solvability is guaranteed
by (1.4).
The system (5.3) can be reformulated in terms of operators as follows. Let
It follows that the solution (X; Y ) of (5.3) satisfies (1.1) with F equal to the L
projection of f into H 2 and G equal to the (L
projection of g into H 1 .
It is straightforward to check that (2.3) holds for A, B, and B T as above. Moreover,
it follows from (1.4) that (2.4) holds with fl independent of the mesh size h.
Remark 5.1. It appears from the definition of the above operators that one has
to invert Gram matrices in order to evaluate the action of A, B T and B on vectors
from the corresponding spaces. In practice, the H 1 Gram matrix inversion is avoided
by suitable definition of the preconditioner QA . For the purpose of computation, the
evaluation of Q
A f for f 2 H 1 is defined as a process which acts on the inner product
data (f; is the basis for H 1 . Moreover, from the definition of the Uzawa-
like algorithms in the previous sections, it is clear that every occurrence of A or B T is
followed by an evaluation of Q
A . Thus the inversion of the Gram matrix is avoided
since the data for the computation of Q
A , ((B T Q; can be computed
by applying simple sparse matrices. In the case of this special choice of H 2 , it is possible
to compute the operator B in an economical way (see Remark 5 of [5]) and we can take
QB to be the identity. For more general spaces H 2 , the inversion of Gram matrices can
be avoided by introducing a preconditioner QB whose inverse is implemented acting on
inner product data as in the H 1 case above.
We still need to provide preconditioners for A. However, A consists of two copies
of the operator which results from a standard finite element discretization of Dirichlet's
problem. There has been an intensive effort focused on the development and analysis
of preconditioners for such problems. In our examples Section 7, we will use a preconditioning
operator which results from a V-cycle variational multigrid algorithm. Such
a preconditioner is known to be scaled so that both (3.2) holds and (2.4) holds with fl
bounded away from one independently of the mesh parameter h.
6. Applications to mixed finite element discretizations of elliptic prob-
lems. In this section we discuss applications of the algorithms analyzed in Section 3
to solving indefinite systems arising from mixed finite element discretizations of second
order partial differential equations. For this application, it will be relatively easy to
construct preconditioners QA while the development of a suitable operator QB is more
difficult.
The basic problem we consider here is
\Gammar
i;j=1 is a symmetric positive definite matrix whose entries are bounded
functions of the spatial
variable,\Omega is a bounded domain with polygonal or polyhedral
boundary in d-dimensional Euclidean space for or 3. This is a classical model
problem in continuum mechanics or fluid flow in porous media.
Introducing a new variable u, (6.1) can be written as a first order system as follows:
In the typical applications K is the elasticity/permeability tensor, u usually represents
the stress/velocity, p is the displacement/pressure. The mixed method naturally takes
into account constraints that appear in the variational formulation of a given differential
provides direct approximations to the two variables of
interest: u and p. Often these features are more attractive then those corresponding to
the standard finite element method.
Then the weak formulation of (6.2) is
The space H div is the set of vector functions in (L
2(\Omega\Gamma8 d whose divergences are
also in L 2
(\Omega\Gamma4 Here, as in the previous section, (\Delta; \Delta) denotes the L
or (L
inner product. The mixed discretizations involve the introduction of two approximation
subspaces, H 1 ae H div for the velocities and H 2 ae L
2(\Omega\Gamma for the pressures. To
this type of application, we will only discuss the simplest mixed finite element
discretization of (6.1), namely the lowest order Raviart-Thomas spaces [24]. We assume
some familiarity with the mixed approximation approach and only give limited detail.
Detailed development can be found in [8], [21], and [24].
Let T h be a partitioning
of\Omega into simplices of quasi-uniform size h. The space H 1
is defined to be the vector valued functions which are linear on the simplices and have
a continuous constant normal component on each of the face of the mesh. The space
consists of the set functions which are constant (discontinuous across the faces) on
each of the simplices. The mixed finite element approximation is defined to be the pair
(r
The operators B and B T are defined as in the previous section. However, for this
application, the operator defined by
In terms of these operators, we get a discrete system of linear equations of type (1.1)
with is the L
orthogonal projection of f into H 2 .
The operator A is well conditioned and hence a simple multiple of the identity
provides an effective QA . On the other hand, the operator BA is not uniformly
well conditioned. In fact, it exhibits a condition number growth like h \Gamma2 and should
be preconditioned in order to get an efficient algorithm of type (2.2) or (2.5). It is well
known that BA behaves like a discretization of a second order operator. In some
applications, it can be preconditioned by cell centered techniques [25], multigrid [7], or
incomplete Choleski factorization of BB T [26].
7. Numerical examples. In this section we present the results from numerical
experiments that illustrate the theory developed in the earlier sections. We also report
similar results obtained from applying the conjugate gradient algorithm for saddle point
problems introduced in [5].
Even though the most effective algorithms result from the use of good precondi-
tioners, we shall initially present results using one of the worst possible preconditioners,
the identity operator. This is important since in some engineering applications, good
preconditioners may not be readily available. We also report results when effective
preconditioners are employed.
The test problem was (5.1)
its exact
solution was zero for both pressure and velocity. We started the iterations with an
arbitrary but fixed initial iterate. All of the iterative methods considered are functions
of the error and thus, iterating for a problem with a zero solution and a nonzero starting
guess is equivalent to solving a related problem with a nonzero solution and a zero initial
guess. We used the discretization described in Section 5.
Our objectives in conducting the numerical experiments were to establish experimentally
the conclusions from the theoretical analysis of the algorithms tested and to
assess their effectiveness in terms of error reduction after fixed number of iterations.
The same nonzero initial iterate was used for all algorithms. As discussed in Section 5,
we used QB j I. The experimental results are organized in four tables.
In

Table

7.1 we give results for three algorithms using QA equal to an appropriate
multiple of the identity. The algorithms are described as follows.
: The algorithm (2.5) with I and
- max is an upper
bound for the largest eigenvalue of A.
: The algorithm (4.4) with I and \Psi defined by one step of the
steepest descent method (SDM) applied to approximate the action of A \Gamma1 .
preconditioned conjugate gradient algorithm for saddle point problems
given in [5] with
- min is a lower bound for the smallest
eigenvalue of A and Notice that the scaling required by Theorem 1 of
[5] is in the opposite direction of (3.1).
The reported error values in Table 7.1 represent the relative error norm after i iterations
computed by
Clearly, this is not the norm which appears in the theory and one cannot expect the
errors to behave in a monotone way. This explains the increase in the reported error
for UID when 1=64. That the USTD method appears convergent
for h - 32 is surprising since (4.5) is not satisfied for these applications. The BPID
method converges considerably faster in these examples since the saddle point method
of [5] is known to give a rate of convergence which exhibits square root acceleration in
cases when poor preconditioners are employed. As expected, all methods deteriorate
due to lack of preconditioning as the mesh size is decreased.

Table
Errors in UID, USTD and BPID by (7.1)
h 200 iterations
UID USTD y BPID
1/8 4.2\Theta10 \Gamma3 5.1\Theta10 \Gamma6 z 6.5\Theta10 \Gamma12
y one SDM step per inexact Uzawa iteration.
z for 109 BPID iterations.
In order to establish experimentally the convergence of UID and USTD, we ran
these two algorithms for 2000 iterations. The results are shown in Table 7.2. Even
though improved convergence is observed in all cases when compared to Table 7.1, the
UID algorithm still appears unstable for 1=64. We ran UID for 10000 iterations
and observed an error of :0048. Although convergent, the inexact Uzawa method with
such a poor preconditioner converges too slowly to be of practical use.
The above results may at first appear to contradict the validity of the theory of
Section 4. The reason that the methods appear divergent at a relatively low numbers
of iterations is that the theorems guarantee monotonicity of the errors in norms which
are different from those used in (7.1). Our next experiment was designed to illustrate
the monotone convergence of UID and BPID predicted by Theorem 1 and Theorem 1
in [5]. Accordingly, we measured the errors in the norms appearing in the theorems. In
the case of UID, we use
(7.
Table
Errors in UID and USTD by (7.1)
h 2000 iterations
UID USTD y
1/32 2.5\Theta10 \Gamma2 2.1\Theta10 \Gamma4
y one SDM step per inexact
Uzawa iteration.
In the case of BPID, we used
The convergence results in these norms are reported in Table 7.3. Note that all of
the reported errors are less than one. We made additional runs at lower number of
iterations. All runs reflected the monotone error behavior in these norms as guaranteed
by the theory.

Table
Errors in UID and BPID by (7.1) and (7.3)
h 200 iterations
UID BPID
1/8 4.29\Theta10 \Gamma3 z 2.1\Theta10 \Gamma12
z for 109 BPID iterations.
The last experiment given in this section is intended to illustrate the performance
of the algorithms when effective preconditioners are used. In this case, we define Q
A to
be the operator which corresponds to one V-cycle sweep of variational multigrid with
point Gauss-Seidel smoothing. The order of points in the Gauss-Seidel iteration was
reversed in pre- and post-smoothing. Note that QA automatically satisfies (3.2) and
satisfies (3.3) with ffi independent of h. We consider the following two algorithms:
: The algorithm (2.5) with I and Q
A being the action of multigrid.
The algorithm from [5] with the A block preconditioned by :5Q
A and

Table

7.4 contains the error reductions for this example. The effect of applying a better
preconditoner QA is clearly seen when we compare the results from UID (Tables 7.1
and 7.2) with those from UMG. Notice that the UMG data in Table 7.4 show little,
if any, deterioration as the mesh size becomes small.

Table
Errors in UMG and BPMG by (7.1)
iterations
UMG y BPMG y
y one multigrid V-cycle per
iteration.
In all of the reported results, the reformulation method of [5] shows faster con-
vergence. Nevertheless, the inexact Uzawa methods are still of interest since they are
robust, simple to implement, have minimal memory requirements and avoid the necessity
of computing inner products. These properties may make the inexact Uzawa
methods attractive in certain applications.



--R

Studies in Nonlinear Programming

The Lagrange multiplier method for Dirichlet's problem
Iterative techniques for time dependent Stokes problems.

Pitman Research Notes in Mathematics Series
The analysis of multigrid algorithms with non-nested spaces or non-inherited quadratic forms
Mixed and Hybrid Finite Element Methods
Balancing domain decomposition for mixed finite elements
Multigrid and Krylov subspace methods for the discrete Stokes equations.
Inexact and preconditioned Uzawa algorithms for saddle point prob- lems
Computational aspects of mixed finite element methods
Preconditioning indefinite systems arising from mixed finite element discretization of second-order elliptic problems
An analysis of the finite element method using Lagrange multipliers for the stationary Stokes equations
estimates for mixed methods
Finite Element Approximation of the Navier-Stokes Equations
Domain decomposition and mixed finite element methods for elliptic problems
Elliptic Problems in Nonsmooth Domains

Probl'emes aux Limites non Homog'enes et Applications
Elements finis mixtes incompressibles pour l'equation de Stokes dans R 3

The convergence factor of preconditioned algorithms of the arrow-hurwicz type
A mixed finite element method for 2-nd order elliptic problems
Finite element and finite difference methods for continuous flows in porous media
A preconditioned iterative method for saddle point problems
--TR

--CTR
Mingrong Cui, A sufficient condition for the convergence of the inexact Uzawa algorithm for saddle point problems, Journal of Computational and Applied Mathematics, v.139 n.2, p.189-196, February 2002
Zhi-Hao Cao, Fast uzawa algorithm for generalized saddle point problems, Applied Numerical Mathematics, v.46 n.2, p.157-171, August
Howard C. Elman, Preconditioners for saddle point problems arising in computational fluid dynamics, Applied Numerical Mathematics, v.43 n.1-2, p.75-89, October 2002
Walter Zulehner, Analysis of iterative methods for saddle point problems: a unified approach, Mathematics of Computation, v.71 n.238, p.479-505, April 2002
Iterative Stokes solvers in the harmonic Velte subspace, Computing, v.67 n.1, p.13-33, July 2001
Qiya Hu , Jun Zou, Substructuring preconditioners for saddle-point problems arising from Maxwell's equations in three dimensions, Mathematics of Computation, v.73 n.245, p.35-61, January 2004
M.-R. Cui, Analysis of iterative algorithms of Uzawa type for saddle point problems, Applied Numerical Mathematics, v.50 n.2, p.133-146, August 2004
Z. Chen , R. E. Ewing, Degenerate Two-Phase Incompressible Flow IV: Local Refinement and Domain Decomposition, Journal of Scientific Computing, v.18 n.3, p.329-360, June
Yanqiu Wang, Overlapping Schwarz preconditioner for the mixed formulation of plane elasticity, Applied Numerical Mathematics, v.54 n.2, p.292-309, July 2005
Angela Kunoth, Fast Iterative Solution of Saddle Point Problems in Optimal Control Based on Wavelets, Computational Optimization and Applications, v.22 n.2, p.225-259, July 2002
Franois Bertrand , Philippe A. Tanguy, Krylov-based Uzawa algorithms for the solution of the Stokes equations using discontinuous-pressure tetrahedral finite elements, Journal of Computational Physics, v.181 n.2, p.617-638, 20 September 2002
uzawa algorithm;indefinite systems;stokes equations;saddle point problems;iterative methods;preconditioners