Wavelet-Like Methods in the Design of Efficient Multilevel Preconditioners for Elliptic PDEs

W a v e l e t - L i k e M e t h o d s in t h e D e s i g n of Efficient M u l t i l e v e l P r e c o n d i t i o n e r s for E l l i p t i c P D E s

Panayot S. Vassilevski and Junping Wang

A b s t r a c t . This article introduces a stabilization of the classical hierarchical basis (HB) method by modifying the HB functions using some computationally feasible approximate L~-projections onto finite element spaces of coarser levels. It is proved that the corresponding multilevel additive and multiplicative algorithms give spectrally equivalent preconditioners for model elliptic problems. In practical computation, one action of the proposed preconditioner is of optimal order algebraically. Numerical results are included to illustrate the efficiency and stability of the new preconditioning methods.

w

Introduction

In this paper we are concerned with the construction of efficient numerical methods for matrix problems arising from finite element methods for elliptic partial differential equations. In practical computations, the standard nodal basis for the finite element space is often chosen as the computational basis and the resulting matrices are ill-conditioned. Our objective is to seek a substitution for the standard nodal basis so that the stiffness matrix arising from the new basis is well conditioned. A computationally feasible basis should possess the following properties: (a) the basis functions must be computable and have local support; (b) the resulting stiffness matrix is sparse and well conditioned. This paper will introduce a wavelet-like method proposed by the authors in [41] which can be employed to construct a new basis with the above mentioned features for the finite element application to elliptic problems. Multiscale Wolfgang Copyright

Wavelet Dahmen,

Methods

for P D E s

Andrew

g. K u r d i l a ,

59 and Peter Oswald

(~)1997 by A c a d e m i c P r e s s , Inc.

All rights of reproduction in any form reserved. ISBN 0-12-200675-5

(eds.), pp. 59-105.

60

P. Vassilevski and J. Wang

The method has a very close relation with the multiresolution (or multiscale) decompositions exploited especially in the wavelet literature (e.g., [17, 14]). In our approach, the requirements on the L2-orthogonality and the existence of a single generating function r (more precisely, r generates an orthogonal basis by dilation and translation as in {r k. - i ) , i E zd}) commonly imposed in the wavelet literature (see Mallat [26] for details) are relaxed in order to have a computationally feasible basis. The basic idea lies in approximating the wavelets without deteriorating their stability, yielding a stable Riesz basis for the finite element space under consideration. Attempts in the search of a stable Riesz basis with some restrictions, either on the mesh or on the analysis, have been made by Griebel and Oswald [20], Kotyczka and Oswald [24], and Stevenson [35, 36]. For a recent comparative study on the construction of economical Riesz bases for Sobolev spaces we refer to Lorentz and Oswald [25]. Our method is general and provides a satisfactory answer for most elliptic equations. The method is based on modifying the existing (unstable) hierarchical basis by using operators which are approximations of the L 2projections onto coarse finite element spaces. A similar approach was used by J affard [22] in seeking a wavelet basis for finite element spaces. In [22], the construction starts with the standard nodal basis functions which are transformed to an orthogonal multiresolution basis based on an explicit orthogonalization procedure exploiting wavelets in each subspace Wj of ~ . Here Wj is the L2-orthogonal complement of the finite element space Vj_ 1 in the next finer finite element space k). The latter procedure generally leads to basis functions that are not locally supported but have good decay rates and hence allow for locally-supported approximations. The method proposed by Vassilevski and Wang [41] can be regarded as an approximation of the wavelet basis in [22]; the approximate-wavelet basis functions are locally-supported and, therefore, computationally feasible. There are two alternatives in the implementation of the approximatewavelet basis in practical computation. In the first approach, one first computes an explicit form for each new basis function and then assembles the global stiffness matrix corresponding to the new basis. The major drawback for this approach is the difficulty encountered while assembling the global stiffness matrix; the robustness corresponding to the standard nodal basis is no longer retained for the new basis. In the second approach, one makes use of the standard nodal basis as the computational one and constructs a preconditioner for the stiffness matrix by using the new basis functions. This is possible because most of the iterative algorithms such as the CG method require only the action of the stiffness matrix on vectors. Thus, the technique of basis changing can be employed to compute the action of the matrix corresponding to the new basis. The second approach can be viewed as a "black-box fix" of the conditioning of the standard nodal

Wavelet-Like Multilevel Preconditioning

61

basis stiffness matrices. Details can be found in Section 8. The organization of the present paper is as follows. Section 2 contains some preliminaries and the problem formulation in an abstract Hilbert space setting. The importance of having a stable Riesz basis and its relation to the conditioning of discrete matrix problems, especially those from the finite element discretization, is the topic of Section 3. Specific examples of second order elliptic PDEs are presented in Section 4. Section 5 discusses the basic idea in constructing multilevel direct space decompositions. The hierarchical basis of Yserentant [44] is reviewed in Section 6. Section 7 contains a detailed discussion of the wavelet-like method in constructing a stable basis. The implementation issues, including two basic preconditioning schemes (namely the additive and multiplicative algorithms) and the construction of the approximate L2-projections, are discussed in Section 8. Finally, some numerical results are presented in Section 9 for convectiondiffusion problems, to confirm the theory developed in previous sections. w

P r e l i m i n a r i e s a n d basic p r o b l e m s

This preliminary section introduces the basic problem we address in this paper. Let H be a Hilbert space equipped with the inner product (., .) such that H C V C H', (2.1) where H is dense in another Hilbert space V with inner product (., ")0, and H ~ is the dual of H with respect to the pairing of V defined by the inner product (., ')0. Assume that the imbedding from H to V is continuous. More precisely, there exists a constant ~0> 0 satisfying

Ilullo

,ilull

v

u E/-jr.

(2.2)

Here I1" II- v/U, )indicates the norm in the Hilbert space H and i1" II0 = x / ( , .)0 that of V. Let a(., .) be a bounded bilinear form defined on H x H satisfying the following inf-sup condition of Ladyzhenskaya, Babugka, and Brezzi:

sup

Ilwll

>/31lv[i -

V v E H,

(2.3)

where/3 > 0 is a fixed constant. We are interested in approximate solutions of the following problem: Given f E H ~, find u E H satisfying

a(u, r

f(r

V r E H.

(2.4)

62

P. Vassilevski and J. Wang

It is well known that if a(., .) is symmetric and H-coercive (i.e., there exists a positive c~ satisfying a(v,v) >_ c~[[vl]2 for all v E H), then the problem (2.4) is equivalent to the minimization problem for the quadratic functional if(v) - la(v, v) - f(v). To approximate (2.4), let {gk}~=l be a sequence of finite dimensional subspaces of H satisfying the following approximation property: For any v E H one has lim

(2.5)

inf I I v - r

k---* oa CE Hk

The well-known Galerkin method for (2.4) seeks uk E Hk such that

a(uk, Ok) -- f(Ok) V Ok 6 Hk. (2.6) In order to have a well-posed discrete problem (2.6), we assume the following discrete inf-sup condition" There exists a constant /3k > 0 such that sup a(v, w) > ~k[iv[ [

V v E Hk.

(2.7)

It is known that if (2.7) holds true, then the Galerkin problem (2.6) has a

unique solution in the subspace Hk. Interested readers are referred to [18] and [11] for more details. In practical computation, we usually formulate a matrix problem for (2.6) by choosing a suitable basis for the subspace Hk. More precisely, let {r . i -- 1 , 2 , . . . , n k } be a computational basis of Hk. Expand the approximate solution u(k) in terms of this basis, yielding rtk

F_,

(2.s)

i--1

Let u (k) - (Cl, c 2 , . . . , cnk) T be the coordinates of u (k). It is not hard to see that the vector u (k) is given as the solution of the following linear system" A(k)u(k) _ f(k),

(2.9)

where the right-hand side vector f(k) is defined as follows: f(k) _ (bl, b2,..., b,~) T

with bi - f(r

The matrix is given by A (k) - {a(r k), r

(2.10)

One may view A (k) as

a linear operator on Hk defined, for any v - ~ vi

E Hk, by

i=1 nk

A(k)v- E i=1

nk

cir

where c i - E j=l

a(r k)' r

(2.11)

Wavelet-Like Multilevel Preconditioning

63

Of main interest in this paper, we study techniques for solving the linear system (2.9) by iterative methods. It is known that the conditioning of the matrices A (k) is of great importance in practical computation. Let us see how the condition number of A (k) is related to the choice of the basis (I)(k) - {elk)} i=1 '~k and the discrete inf-sup condition (2.7) For any v E Hk, denote in the bold face the coordinates of v with respect to the basis (I)(k). Namely, the vectors v (Vl, v2, 9 vnk )T and v are related as follows" _

..

_

,

rtk

v-

~

vir k).

(2.12)

i--1

Introduce a new inner-product in Hk by using the basis {r

i=1 as follows:

nk

k,0 - ~

viwi - v . w .

(2.13)

i--1

Denote by I]vl]k,0 - x/V v, the norm induced by the new inner-product. Recall that the condition number of the matrix A (k) is defined by ~ ( A ( k ) ) - IiA(k)ll IlA(k)-~ll ,

(2.14)

where [[A(k)l[ indicates the norm of the matrix A (k) with respect to the new norm I1" []k,0 in the vector space Hk. The following result can be checked easily. T h e o r e m 1. Assume that there are positive constants Ok,1 and Qk,2 satisfying ~ok,l[[v[]k,o _< IIA(k)vlik,o <_ ~ok,2[[v[[k,o. (2.15) ~)k,2 Then, ~(A (k)) _< ~-ZT,~" The best estimate for Qk,1 and Qk,2 is given by

~ok,2 -

sup ""llA(k)vIJk,o, ][vllk'~ ,, ,, []A(k)vl]~,o. inf

v6_Hk,

t~k,1 =

(2.16)

vEHk, [[vilk,o=l

Moreover, one has ~(A (k)) _ ~k'2 for the best estimate of ~k, In fact, the inequality (2.15) implies the following []A(k)[[ _< ~ok,2,

I]A(k)-~]l _< 1/~ok,1,

(2.17)

which verifies the validity of the theorem. Therefore, it suffices to establish an estimate like (2.15) in order to gain some knowledge about the condition number of the matrix A (k).

P. Vassilevski and J. Wang

64 w

M a t r i x c o n d i t i o n i n g a n d s t a b l e Riesz bases

Our objective in this section is to figure out the connection between the condition number of A (k) and the selection of the basis {r for Hk ~=1 To this end, we derive the inequality (2.15) by using the mr-sup condition. First, notice that for any v E Hk,

IIA(k)vllk,0 --

sup

~eH~, ll~llk,o=1

(A(k)v, w)k,o.

(3.1)

The definition of A (k) (see (2.11))implies that

(A(klv, w)k,o - a

vie ~'),

wi

It follows that IIA(k)vllk,0

--

= a(v, w).

a(v, w).

sup wEHk, llwllk,o=l

(3.2)

(3.3)

By letting 0 k , 2 - sup sup

,

(3.4)

one obtains the following (3.5)

IIA(k)vl]k,o < ek,el]vllk,o. Next, by letting Ok 1 '

inf

sup

Ilvllk,011

ll

,0

,

(3.6)

one has the following estimate from below:

IIA(k)vl]k,o >_ 0k,ll]Vl]k,0.

(3.7)

To summarize, we have proved the following result: T h e o r e m 2. If Ok,1 and 0k,2 are given by (3.6) and (3.4), respectively, then the estimate (2.15) holds true. Consequently, the condition number of the matrix A(k) is bounded by ~k,2 ~k,1

"

The question now is about the determination of the constants 0k,1 and 0k,2" We would like to select a basis {r for Hk so that the ratio (or 1 the condition number of A (k)) ek,2 is as small as possible. L0k,1 As our first consideration , we assume that {r ,~k is an orthonormal }i=l basis with respect to the original inner product (., .) of H. It is clear that

65

Wavelet-Like Multilevel Preconditioning

the corresponding discrete norm II. IIk,o is the same as the original norm II. II in Hk. Therefore, the constant 0k,2 is bounded from above by the norm of the bilinear form a(., .). Similarly, the constant 0k,1 is bounded from below by the parameter ~k in the inf-sup condition (2.7). Since in practical problems, the norm of a(., .) and the parameter/3k stay uniformly bounded in terms of k, the condition number of A (k) is uniformly bounded. It is, of course, impractical to assume the existence of an orthonormal basis {r which is also computationally feasible. The next best thing to an orthonormal basis is the Riesz basis in any Hilbert space. We recall that in the Hilbert space H, a Riesz basis is a basis {r }~--1 of H satisfying OO

#lllvll 2 < ~ c,2 < #211vii2

v v e H,

(3.8)

i=1

where v - EiC~=l c i r [l" I] indicates the original norm in H, and 0"1 and 0"2 are two absolute constants. Here we have assumed that the Hilbert space is separable. R e m a r k 1: The strong norm I" II in the Hilbert space H was used in the definition of the Riesz stability (3.8). Notice that H is a subspace of V. Thus, it would be feasible to discuss the Riesz property with respect to the norm !1" II0 of V. In practical applications, V often represents L2(Ft) and I1" II0 denotes the corresponding L2-norm. In the wavelet literature, the Riesz property is commonly studied with respect to the L2-norm, but the ultimate goal in the preconditioning analysis is to establish estimates such as (3.8), since (3.8) implies that the discretization matrix A (k) is well conditioned. We now go back to the finite dimensional subspace Hk of H. As in (3.8), we assume that there exist constants 0"1 (k) and 0"~k) such that

~k)llvll2 < k

c,2 < 0"~k)Ilvll 2

v v e Uk.

(3.9)

i=1

Observe that Condition (3.9) can be rewritten as 0"{k)l]vll2 _< I]vl]~,0 _< 0"~k)l]vll2 Since {r

i=1

g v e H~.

(3.10)

is assumed to be a basis of Hk, the equivalence relation

(3.10) obviously holds for some constants 0"~k) and 0"2(k) with 0"2(k) =

sup I]vll2k , 0 ~eH~, I1~11=1

(3 911)

o'~k) -_

2 inf llvl]k,o. ,,eH~, I1,,11=1

(3 12)

and

66

P. Vassilevski and J. Wang

The important thing is that one should have some control over the ratio

Definition 1. The family { ~(k ) - {r

} iS said to be a uniformly sta-

b•e family of Riesz bases for {Hk} if the quotient a~k)/a~ k) of the Riesz bounds in (3.9) is bounded uniformly with respect to k --~ oo, i.e., if there exists a constant M independent of k such that (k). (k) (3.13) a 2 / a 1 <_ M for any k = 1, 2, ....

The rest of this section explains the importance of having a stable Riesz basis. First we estimate the parameter Qk,2 defined in (3.4). Let [la[[ denote the norm of the bilinear form a(.,-). Then, ~k,2

- sup sup

a(v,~) IIvll~,011wll~,0

< sup sup

a(v.~) IIvllll~ll sup sup II~llll~ll ~ . ~ ~ n ~ II~llk.011~llk.0

vEHk w e H k

-~.~

-I[a[I

~.~

sup

I]v[]l]w[[

(3.14)

< I[all/cr~k)

v,,~eH,, Ilvllk,011wllk,0 --

Next, we estimate the parameter Qk,1 as follows: Qk,1

a(v,~) ,,eHk ,,,eHk Ilvll~,011~ll~,0

-- inf

sup

> inf

sup

>

inf

inf

a(v,w)

,i,;ii;,oii ll

sup

a(v,~)

inf

Ilwil

i1 11 ,o Ilvll

inf

(3.15) Ilwll

-,~eH,~,~eH,. Ilvllll~ll veHk IIvIIk,0 ,,,eHk Ilwllk,0 >_ Z,~/,,'~). a(k) It follows that Ok,2/Qk,1 _< ~l~ko~k).

The result can be summarized as

follows: T h e o r e m 3. Let {r be a computational basis of Hk satisfying(3 10) Then the condition number of the matrix A (k) arising from the basis {r is bounded by 1

~

g(A (k)) _< Ilal] cr~k) /~k ~ k ) '

"

(3.16)

Wavelet-Like Multilevel Preconditioning

67

Consequently, the condition number of A (k) is uniformly bounded for stable Riesz bases {r provided that the discrete Jar-sup condition (2.7) holds true with uniformly bounded constants/~k from below by some/3* > O. The conditioning estimate in Theorem 3 contains two important factors:

a(k) ~k and a~-'~" The first one depends on the norm of the given bilinear form and the stability constant/3k from the discrete Jar-sup condition. In practical computations, the space Hk must be constructed so that it ensures the boundedness of f3k from below by some fixed/3* > 0. This is the case for the model problems and their discretization spaces to be considered in Section 4. The second factor is basis-dependent. More precisely, it is a characterization of the difference between the discrete coefficient norm I1" llk,0 and the continuous norm I1" II. Stable Riesz bases are important because the corresponding discretization matrices are well conditioned. Thus, simple iterative methods such as the conjugate gradient (CG) can be successfully applied to solve the matrix problem from the Galerkin discretization with a geometric rate of convergence. As is well known, the convergence factor 1 where x - n (A(k) ) -< ~#k o(k) is bounded by ~/~-v/-E+I, o~k). For nonsymmetric problems, one could apply the CG method to the normal equation or the GMRES method applied to A (k). For symmetric and indefinite problems, the MINRES (minimum residual) algorithm would be a good choice. The convergence rate is no worse than that of the CG-method applied to A (k)2. There is another practical criterion in the choice of the basis {(~Ik)}n~l . It is of great practical importance to represent the matrix entries of A (k) as sparsely as possible. This is trivially achieved (assuming a(., .) is symmetric and positive definite) if the basis is a(., .)-orthogonal. Thus, the matrix A (k) admits diagonal forms and only nk entries need to be stored. Such a situation is too special and rarely happens in practice. In general, we assume that the basis is computationally feasible in the sense that the basis function are computable and the corresponding matrix A (k) is sparse (the number of nontrivial entries in the matrix is of order O(nk)). This is the case in practice for the standard nodal bases of finite element spaces Hk if the bilinear form arises from partial differential equations. However, this choice will not make a stable Riesz basis for most of the PDEs. The following section contains a detailed discussion of this aspect. w

Model problems

Here we consider some model problems of (2.4) in partial differential equations. Boundary value problems for the second order elliptic equations are

P. Vassilevski and J. Wang

68

of major consideration in this discussion. Finite element methods will be applied to approximate the solution defined on an open bounded domain gt in 11%d with d = 2 or 3. 4.1

S e c o n d - o r d e r elliptic e q u a t i o n s

Consider the homogeneous Dirichlet boundary value problem for the following second-order elliptic equation:

L(u) -- -~7 . (a(x)~Tu) + b(x). ~7u + c(x)u = f(x),

x E f2,

(4.1.1)

where a = a(x) is a symmetric and positive definite matrix with bounded and measurable entries, b = b(x) and c = c(x) are given bounded functions, and f = f(x) is a function in H-l(ft). Note that we do not intend to consider problems which are convectiondominated. Let Hl(f~) be the standard Sobolev space equipped with the norm:

I1 11 - (llullo + II ullo

V u E Hl(ft).

(4.1.2)

Here ll" {i0 stands for the L2-norm. Let H01(ft) be the closed subspace of Hi(f2) consisting of functions with vanishing boundary values. The following relation is well known: H~(ft) C L2(ft)C H-I(f2).

(4.1.3)

A weak form for the Dirichlet problem of (4.1.1) seeks u E H~(~) satisfying b(u, v ) - f(v) V v E H~(~), (4.1.4) where

b(u, v) - / f ~ (a(z)Vu . Vv + b ( z ) . Vu v + c(z)uv) dx

(4.1.5)

and f(v) is the action of the linear functional f on v. Assume that the problem (4.1.4) has a unique solution. Then the inf-sup condition (2.3)is satisfied for the bilinear form b(., .) defined on H01(f2) • H01(f2). Let us approximate (4.1.4) by using the Galerkin method with continuous piecewise polynomials. If Sh denotes the finite element space associated with a prescribed triangulation of ft with mesh size h, then the Galerkin approximation is given as the solution of the following problem: Find Uh E Sh satisfying

b(uh, r

f(r

V r E Sh.

(4.1.6)

69

Wavelet-Like Multilevel Preconditioning

It has been shown that the discrete problem (4.1.6) has a unique solution when the mesh size h is sufficiently small. Thus, the discrete inf-sup condition (2.7) is satisfied for this problem. Details can be found from [32, 33]. Choose the standard nodal basis as the computational basis for the finite element space Sh. Let {r }~=1 be the set of nodal basis functions and Ah be the corresponding discrete matrix (also called the global stiffness matrix). The condition number of Ah can be estimated by using Theorem 3. To this end, let us establish Inequality (3.10) for the standard nodal basis. For any v E Sh, let tl

v - ~

with vi - v(xi),

vir

(4.1.7)

i=1

where xi are the interior nodal points of the finite element partition. Relation (3.10) is equivalent to the following: n

~lllvll21~ ~ v~ ~ llvll ,

v v e Sh

(4.1.8)

i=1

for some positive constants &l and #2. It is not hard to see that n

,lvll02

~ h-

(419)

i=1

It follows that &2 - O(h-d). Also, by the standard inverse inequality one sees that #1 is bounded from below by a constant proportional to h 2-d. Thus, from Theorem 3, the condition number of Ah is bounded from above by Ch -2 for some constant C; the lower bound for its condition number is also bounded from below by some Ch -2.

4.2

Stokes equations

Consider the problem which seeks u e [H~(a)] d and p e L 2 ( a ) s u c h that -Au+Vp V.u u

= = =

f, 0, 0,

inf,, in f~, on Oft,

(4.2.1)

where f E [L2(f~)] d is a given vector-valued function and 0f~ denotes the boundary of f~. A weak form of the problem (4.2.1) involves the following bilinear form: A(u, p; v, w) = (Vu, Vv)0 - ( V . v , p)0 - ( V . u , w)0

(4.2.2)

P. Vassilevski and J. Wang

70

defined on W • W with 142 - [H0i(fl)] d • L20(f~). Here L~(~)is the closed subspace of L2(fl) consisting of functions with vanishing mean value. The weak problem seeks (u, p) E )42 satisfying A(u, p; v, w) - (f, v)

V (v, w) E W.

(4.2.3)

The inf-sup condition is satisfied for the bilinear form defined in (4.2.2). Details can be found from [18, 11]. To apply the finite element method to (4.2.3), we employ the HoodTaylor element [21] which satisfies the discrete inf-sup condition (with a mild restriction on the triangulation). The Hood-Taylor element is a combination of continuous piecewise linear functions for the pressure variable p and continuous piecewise quadratic functions for the velocity variable u. Denote by Wh = Xh x Sh the corresponding finite element space, where Xh contains continuous piecewise quadratic functions and Sh contains continuous piecewise linear functions for the pressure variable. The finite element approximation (Uh, Ph ) E Wh satisfies .A(uh, ph ;',,, w) - (f, v)

V (v, w ) E Wh.

(4.2.4)

If the standard nodal basis is selected in formulating a matrix problem for (4.2.4), then the condition number of the global stiffness matrix can be estimated by using Theorem 3. More precisely, let {r be the that of Sh, as discussed in the standard nodal basis of Xh and {r previous section. For any (v, w) E Xh • Sh, let tl

w - E wir

(4.2.5)

i=1

and

m

v - E

vjCj.

(4.2.6)

j=l

Using the inverse inequality and the Poincar~ inequality one can derive the following relations" n

2

(4.2.7)

v~ _< c2h-d]JvlJ 21,

(4.2.8)

i=1

and

m

ci h2-dl]vi] ~ _< ~ j=l

where ci and c2 are two absolute constants. Thus, we have from Theorem 3 that the condition number of the global stiffness matrix is bounded by Ch -2.

Wavelet-Like Multilevel Preconditioning

71

We emphasize that the poor conditioning for the Stokes problem is caused by Relation (4.2.8), where the Hi-norm of v was approximated by a discrete norm stable in L 2 only. The equivalence (4.2.7) indicates that the standard nodal basis is a good choice for the pressure variable in the Stokes problem. Therefore, attention should be focused on stabilizing the velocity component in the Stokes equation. A direct wavelet approach to the Stokes problem has been developed by Dahmen et al. [16]. One could also use block-diagonal preconditioners for the saddle-point discretization matrices A (k) with one block corresponding to Laplace-like preconditioners for the velocity component and a second block corresponding to massmatrix preconditioners for the pressure unknown in the MINRES method. Details of this approach can be found in Rusten and Winther [31] and Silvester and Wathen [34]. 4.3

Mixed methods

Here we consider the mixed method for the second order elliptic equation (4.1.1). For simplicity, assume that b = 0, c _= 0, and the following Dirichlet boundary condition u =-g on 0f~ (4.3.1) is imposed on the solution. Let (

, V . v C L2(f~)},

H(div; f ~ ) - ~v" v G which is equipped with the following norm: []Vl[H(div; a) -Let c~(x)=

(Iv[ 2 + [V. vl2)dx

a-l(x) be the inverse of the coefficient matrix a = a(x) and .A(q, u; v, w) = (a(x)q, v)o - (V- v, P)o - (V -q, w)o

be a bilinear form defined on W x W, where W = H(div; f~) x L2(fl). Then, a mixed weak form for (4.1.1) with Boundary condition (4.3.1)seeks (q, u) E W satisfying A(q, u; v, w) = (g, v . n)or~ - (f, w)o

V (v, w ) E W,

(4.3.2)

where (., ")0~ denotes the inner product in L2(O~). The inf-sup condition (2.3) can be verified for the bilinear form A(.; .). Furthermore, finite element spaces satisfying the discrete inf-sup condition (2.7) are available for this bilinear form. Details can be found in the book by Brezzi and Fortin [11]. If the standard nodal basis for the finite element

P. Vassilevski and J. Wang

72

space of Raviart and Thomas is employed in practical computation, the condition number of the global stiffness matrix is known to be proportional to h -2. The question is how to stabilize the nodal basis. Observe that in this application, one needs to construct a basis for the mixed finite element space (called Wh) so that the discrete norm is equivalent to the following norm: /l(v; w)ll~v -Ilvll~(div; ~) + IIwli~. One sees from above that the standard nodal basis is a good choice for the pressure unknown. Thus, the difficulty lies in stabilizing the flux component with the norm in H(div; ft). We do not yet have a positive answer yet for this stabilization, but the approach of ttelmholtz decomposition for vectors might provide a partial answer for problems of two-space variables. The method decomposes each vector-valued function v as follows: v - curl r -b Vr

(4.3.3)

where r and r are functions in Hl(f~). A discrete version of (4.3.3) can be studied in order to apply it to the mixed method. Results for the standard conforming and non-conforming finite elements should be investigated first in this direction. We point out that this approach may have some difficulty for problems of three-space variables. For multilevel methods relying on the above Helmholtz decomposition, see Vassilevski and Wang [40] and Arnold, Falk, and Winther [1]. In the following sections we will devote ourselves to modifying the nodal bases in the application to second order elliptic problems. Our goal is to construct some stable Riesz bases that are computationally feasible for elliptic equations. More precisely, the basis should be so constructed that the resulting discretization matrices are both well conditioned and sparse. w

M u l t i l e v e l direct d e c o m p o s i t i o n s

To construct a computationally stable basis for the second-order elliptic and the Stokes equations, we take advantage of the fact that the weak problem (2.4) is discretized on a sequence of finite element subspaces. In particular, a sequence of nested subspaces may be possible in practical computations. Our objective in this section is to exploit ways of constructing stable basis by using the information from each approximating subspace. 5.1

T h e basic idea

The basic idea comes from the fact that an L 2 orthonormal basis of wavelets is also H 1-stable in applications to partial differential equations. There-

Wavelet-Like Multilevel Preconditioning

73

fore, wavelet bases are good candidates in formulating matrix problems for (2.6). Since the conventional wavelet bases have complicated structures which limit their application in the numerical methods, we shall focus our attention on approximations of wavelet bases. Below we present a detailed discussion. Assume that we have a sequence of nested subspaces {l~ }~=0 satisfying

VoCV~C...cvkc....

(5.1.1)

Each vector space V) shall be referred to as a coarse subspace of Vk when j < k. In applications, they are finite element spaces consisting of continuous piecewise polynomials over a sequence of finite element partitions for the domain f~. Upon viewing Vj as a subspace of L2(f~), one has (:x:}

U V~ - L2(a),

(5.1.2)

j=0

where the closure was taken in the strong topology induced by the L2-norm. We assume that V0 is a very coarse subspace of L2(ft) whose dimension is a small number. For every j > 1, define Wj to be the L2-orthogonal complement of Vj-1 in I~. We have Vj -

Vj _ 1 (~ W j

(5.1.3)

and

Wi_I_Wj if i C j ,

(5.1.4)

where we have assumed that W0 - V0. It follows that k

(5.1.5)

j=0

where all these subspaces are orthogonal. By virtue of (5.1.2) and (5.1.5), this implies C~

j=l

which is a decomposition of L2(f~) into mutually orthogonal subspaces. A wavelet basis for L2(f~) can be constructed if one is able to find an orthonormal basis for each subspace Wj. Such a basis would be ideal in preconditioning the discrete problem (2.6), if it is computationally feasible. In practice, it is very hard to find an L2-orthonormal wavelet basis which is also computable. Therefore, the orthogonality requirement in the

74

P. Vassilevski and J. Wang

decomposition (5.1.5) shall be relaxed to allow only a direct decomposition of the following form" - Voe

e Vl

v:,

(5.1.6)

where each Vj1 is a complement of ~ _ 1 in Vj such that the corresponding two-level decomposition is direct. But in order to attain the stability of the wavelet basis, it is crucial to have some approximate orthogonality among the subspaces Vj1. 5.2

A general approach

A general method for deriving the hierarchical complement of each ~ _ 1 in V) is based on the existence of some computationally feasible projections 7rj from C, a dense subspace of the Hilbert space H, onto Vj. In particular, we assume that C D UI>,VI and zrjr = r for any r E Y). Thus, one has 7rjTri - 7ri for j >_ i ifV~ C V). With Vj1 - ( I - T r j _ l ) ~ , one has the following two-level direct decomposition: (5.2.1)

Definition 2. (Multilevel Hierarchical Basis) For j - O, 1 , . . . , k, let {r 1 , . . . , nj } be a c o m p u t a t i o n a l l y feasible basis o f l ~ . A s s u m e that {r

ii --

1 , . . . , n j _ l } U{ r ), i - nj_ 1 "~- 1 , . . . , nj } forms a basis of k~ . A multilevel hierarchical basis for Vk is defined as follows: k

(I) k -- U { ( / - - 7 r j _ l ) r j=O

j)

i--nj

1-~" 1

nj}

(5.2.2)

Here we have assumed that 7r_ 1 - 0 and n-1 - O. We now discuss the stability of the multilevel hierarchical basis. The following result sets a guideline for the selection of the operators 7rj.

Theorem 4. A necessary condition for ~k to be a stable Riesz basis o f Vk is that the projection operators 7rr be uniformly bounded on Vk with respect to r and k for any r <_ k.

75

Wavelet-Like Multilevel Preconditioning

Proofi

Let v E Vk be expanded as follows:

v-

k

cl )Sl i--1

{r

nj

E

Z

j=r+l

i=nj-l+l

'~r~=~

where is the multilevel hierarchical basis of V~ and ~s) - ( I 7rs_l)r ~). Introduce the notation c(~) c~r+l) C --

ci where c~j) -- "( C(j) n j _ l + l , . . . ,c ~) ) T f o r j - r + l , . . . , k ,

c(~) -(c?),...,

and

c~?) ~ .

Since, by assumption, we have a stable Riesz basis, then there exist (rl and cr2 independent of k such that 0" 1

Observe that 7r~v 7r~v, we obtain

(5.2.3)

c T C < Ilvl12 < ~ C r c .

~

. Thus, from (5.2.3), with v replaced by

i=1

Ill, vii ~ < ~c(~)'c(~)

< ~

(

c(~)~c(~) + ~

c~)~c~ ~)

j=r+l

)

< ~-~llvll ~. O'1

Here we have used the lower bound of (5.2.3). This shows the boundedness of the projection operator 7r~ for any r. w

T h e h i e r a r c h i c a l basis

In this section we review the classical hierarchical basis decomposition. First, partition the domain f2 into large elements. Let To denote this initial coarse triangulation and V0 be the corresponding finite element space of continuous piecewise linear functions. The fine finite element space V - Vj corresponds to the triangulation Tj which is obtained by J _> 1 successive refinements of the coarse triangulation To. For problems of two space variables, one can use the triangular element and the refinement of one triangle at level k - 1 will generate four congruent triangles of level k by connecting

P. Vassilevski and J. Wang

76

the midpoints of its edges. Similar techniques are available for problems of three space variables with tetrahedra being used as elements. We refer to Ong [29] for more details of this discussion. It is also possible to use bisection refinement for both two- and three-space problems. Details about this technique can be found from Mitchell [28] and Maubach [27]. Let 7~ denote the finite element partition at level i, and 17/ be the corresponding finite element space of continuous piecewise linear functions. Thus, we obtain a nested sequence of conforming finite element spaces

VoCVxC...cvj, which can be used to discretize the second-order elliptic equation. To describe the classical hierarchical basis, let Afi be the node set of nodal degrees of freedom at level i which consists of vertices of triangles (or tetrahedra in 3-d) in 7~. One has the following natural direct decomposition: = X:

with A/'~ being the set of newly-introduced nodal points. For example, in 2d, the nodal set A/'~ contains the vertices of the triangles at level k that are midpoints of the edges of the triangles from level k - 1. We also introduce the mesh size hi - 2-ih0 for the ith level triangulation Ti. Here, ho stands for the maximum diameter of the elements in To. Recall that the standard nodal basis functions {r k), xi E A/k} are defined to satisfy the condition r k) (xj) = $i,j, where 6i,j is the Kronecker symbol, with xj running over all the nodes in A/'k. One can then define a two-level hierarchical basis of Vk by adding to the nodal basis of Vk-1 the nodal basis functions of Vk corresponding to the complementary nodal set A/'~ - A/'k \ A/'k-1. Consider now the nodal interpolation operator Ik " C(f~) ~ Vk defined by Ikv -v(xi)r k). It is clear that Ik is a projection and Ikr -- r if @ E Vk. x~EA/'k With the choice of Try ~_ Ij, one obtains the classical hierarchical basis by using the general decomposition method discussed in Section 5. We comment briefly on the stability of the classical hierarchical basis in applications to the second order elliptic problems. In these applications, the natural norm for the finite element space is the norm in the Sobolev space H 1(f~). Since the interpolation operator Ik is not bounded in the H 1norm (e.g., [13]), Theorem 4 implies that the classical hierarchical basis is not absolutely stable. Thus, the resulting stiffness matrices computed with respect to the hierarchical basis will not be well conditioned. However, as is well known, the condition number for problems of two-space variables is practically acceptable. In fact, the condition number for the second-order elliptic problems can be verified to be proportional to k 2 at level k. This condition number grows more slowly than that from the standard nodal

Wavelet-Like Multilevel Preconditioning

77

basis, which behaves like O(h~2). For problems of three-space variables, the condition number for the classical hierarchical basis is of order O(h~-1) which is better than the one from the standard nodal basis but worse than in the case d - 2. w

A s t a b l e R i e s z basis b y w a v e l e t m e t h o d

In this section we will construct appropriate projections 7rk which are H Istable and provide a computationally feasible Riesz basis for Vj. The bilinear forms of main interest are those arising from the Hilbert space method for second-order elliptic problems discussed in Section 4. The method presented here was proposed by Vassilevski and Wang in [41]. 7.1

O n t h e basis c o n s t r u c t i o n

Define the L2-projection operators Q k ' L 2 ( ~ ) ~ Vk as follows:

(Qkv, r

- (v,r

vtE

vk.

Also, assume that there are computationally feasible approximations Q~ 9 L2(f~) ~ Vk of Qk such that for some small tolerance r > 0 the following estimate holds:

[I(Q~ - Qf~)vll0 _< ~lIQkvll

V v E L2(f~).

(7.1.1)

The projection operators of major interest are defined as follows: J-1 (7.1.2)

7rk - 1-I (ij + Q;(Ij+x - i i)), j=k

with ~'j - I. It is clear that ~'kr - r i f r E Vk s i n c e / j r 1 6 2 / i t - r for j >_ k based on ( I j + l - / j ) r - 0 and Ijr - r This also implies that 7r~ - 7rk. Note that 7 r k - l ( / k - I k - 1 ) r Q ~ _ l ( X k - Ik-1)r and 7rk - 7rk_l - ( I - Q~_l)(Ik - Ik-1)Trk. Then, the components in the definition (5.2.2) for the wavelet-like multilevel hierarchical basis read as follows" k {r ~ i -

1

no} U { ( / - j=l

Q;_1)r j) i -

nj 1 + 1

hi}

(7.1.3)

The above components { ( I - Q~_l)q~l j), i - nj_ 1 + 1,..., nj) can be seen as a modification of the classical hierarchical basis components based on the interpolation operator Ik, since ( I - Q]_I)r j) - ( I - Q~_I)(Ij -

78

P. Vassilevski and J. Wang

Ij_l)r

j)', the modification of the classical hierarchical basis components

{(Ij

/j_l)r j),

--

+ 1 , . . . , n j } comes from the additional term In other words, the modification was made by sub-

i -- n j _ l

Q ~ _ I ( I j - Ij-1)r

tracting from each nodal hierarchical basis function r its approximate L2-projection Qa_lr onto the coarse level j - 1. Such modifications of the hierarchical basis function r for some particular choices of Q~_I will be shown in Figures 1-3 in Section 8. It can be seen that the modified hierarchical basis functions are close relatives of the known Battle-Lemari~ wavelets [17]. Observe that in the limit case of Q~ - Qk, so that r - 0 in (7.1.1), we get 7rkv - Q k l k + l Q k + l l k + 2 . . . Q j _ l l j v

- QkQk+l...Qj_lv

- Qkv.

Therefore, :rk reduces to the exact L2-projection Qk. As is well known, the L2-projection operators are bounded in both H~(ft) and L2(Ft). This gives us hope that the hierarchical multilevel basis corresponding to the above choice of the operators 7rk may yield a stable Riesz basis if r is sufficiently small.

7.2

Preliminary estimates

For an analysis of the multilevel basis (7.1.3), we need some auxiliary estimates already presented in Vassilevski and Wang [41]. The following result on estimating the error ej - (~rj - Q j ) v will play an important role in our analysis. L e m m a 1. There exits an absolute constant C such that k

j=l

k

h72il~j 11o2_< c~ 2 ~ h~-211(Qj j=l

Q~-~)vll~),

u

(7.2.1)

The estimate (7.2.1) relies on the following recursive relation"

~,_~ - (Q,_~ + R,_~)~, + R,_~(Q, - Q,_~),, where R~-I - ( Q , - 1 e8-1

Q~_I)(/,-1 - / , ) .

It can be seen as follows"

-- 7 r s - I V -- Q s - I v

= (I,_~ + Q L ~ ( I ,

-

~,_~))~,v

-

Q,_~v

= (Q8-1 - Qa_l)Is-17r, v + Qa_ l~r,v - Q , _ l V .

(7.2.2)

79

Wavelet-Like Multilevel Preconditioning Thus, one has e,_~

-- (Qs-1 - Qas_l)Is-l(Trsv - Qsv) + Qa_l(Trsv - Q , v ) +(Q,-1

-

Qa_l)Is-les-1 -b Qa_lIse8

+Qs-I(Is-IQsv

- Qsv) - Q a _ l ( I s _ l Q s v - Qsv)

= (Qs-~ - Q~_~)(18-1 - I8)e8 + (Qs-~ - Q~_~)e8 + Q~_le8

+(Qs_I-Qa_~)(18_I-18)Qsv = [Qs-~ + (Qs-~-

Q~_~)(18-1- 18)] e8

+(Qs_I-Qa_I)(Is_I-I~)Q,v. The latter together with the fact that (Is-1 - I s ) Q s - 1 desired recursive relation (7.2.2).

-

0 implies the

Proofi We now verify L e m m a 1. Let C n be a mesh-independent upper b o u n d of the L2-norm for the operator I8 - I8-1 " V8 ~ Vs. Then, IlR,-xvii0 _< CRr[bv[}0

for all v C Vs.

(7.2.3)

We assume that, C R r _< q 0 - Const < 1.

(7.2.4)

Then, 1

1 + qo = Const < 1 2

(7.2.5)

Next, observe that ek - 0 (since v E Vk). Then a recursive use of (7.2.2) leads to

ILe,-lll0

(1 + CRr)l]e, [io + CRrl[(Q, - Q,-x)v[[o k j=s

P. Vassilevski and J. Wang

80

Therefore, with hj - 2-Jho and h0 being the coarsest mesh size,

II~s-lli0

k

<_ C n r h s _ l ~_, (1 + Cnr)J-sh-21[[(Q i - Qj_l)v[Io j-s k

- Cnrh~_x ~ (1 + Cnr)J-~h-i~hjhy~ll(Q j - Qj-1)v[10 j=s k

1 j-s - CR~-h._~ F, (1 + CR~-)J -~ (~) h-/~ll(Qj

Qj_l)Vll0

-

j=s k

<_ CRvhs-1 ~ qJ-Sh21[l(Q j -Qj-1)v[}o j=s

1 <- CRvhs-11vff-:~-q

=s qJ-Sh-j-2 i](Qj - Qj-1)v[12o

9

The latter inequality shows k E h~---21[[es- 1[[0 s=l

k

<- C~r21-~q E

k

E qJ-Sh;2[l(Qj - Q/-1)v[[02

s=l j=s k

_< c~,2(~_q)~

h;21](Q~ - Qj - ~)vll 2o.

which proves the lemma. The above proof also shows the following corollary. Corollary 1. For any ~r E (0, 1], the following estimate holds: k E

k

-20" ll~,-~ll~ _< C ~ ~ h~_~

1 ~ hj2~ ( 1 - q ) 2 j= 1

Qj 1)vii 2 --

0,

provided that v satisfies the estimate C n r _< q2 ~

1,

(7.2.6)

for a mesh-independent constant q E (0, 1) (actually q > 2-a). R e m a r k 2: Corollary 1 indicates that in order to have the L2-stability of the deviations, one has to assume a level dependence on the tolerance r. More precisely, there exists a 7"o > 0 such that if v < r0J -1, then k

E s--1

[}e~-1[[02 -< C[Iv[[~)

for all v e Vk.

(7.2.7)

Wavelet-Like Multilevel Preconditioning

81

L e m m a 2. Let V~ - ( I - Mk_I)Vt:(1), with V(1) - (Ik - Ik-1)Vk, be the modified hierarchical subspace of level k for any given L2-bounded operator Mk-1. Then, there are positive constants Cl and c2 independent of k such that

c:11r for any r

_ (I-

_< IIr

_< c=11r

r=0,1,

(7.2.8)

Mk_l)r 1 E V1, (/)1 E V(1). Here I1" I1: ~ . , . d ~ tot th~

. o t t o i . th~ Sobo1~v ~ p ~

Hg(~) ~.d

II Iio denotes the L2(f2)-norm.

Proof: The following strengthened Cauchy inequality is valuable: There exists a constant 7 E (0, 1), independent of the mesh size or the level index k such that

(Vr 1, V(~) < "y(Vr 1, Vr

89(V(~, V~) 89 ,

(7.2.9)

for all (~1 E Yk(1) and r E Yk-1. In fact, we shall make use of the following version of (7.2.9). For any r

e V(1) and r e V~-I, one has

(:7(rI + r V(r 1 + 8)) >_ (i - 72)(Vr I,Vr

(7.2.10)

A derivation of (7.2.9) and (7.2.10) can be found from Bank and Dupont [5] or Axelsson and Gustafsson [3]. We first establish (7.2.8) for the case r = 1. With r = - M k _ : r 1 we see from (7.2.10) that

(1

-

~=)11r

< i1r

Thus, the inequality on the left-hand side of (7.2.8) is valid with cl = i - 7 2. To derive the part on the right-hand side, we use the standard inverse inequality to obtain

2 where we have used the L2-boundedness of the linear operator Mk-1. Observe now'that since r E V(1), there exists a constant C such that

IIr ~ ~ Ch~llr

(7.2.11)

It follows that Itr _< ciir for some constant C. This completes the proof of (7.2.8) for r - 1. Similar arguments can be applied to verify the case r = O.

P. Vassilevski and J. Wang

82 Lemma 3. F o r a n y r 1 - ( I - M k _ l ) r w i t h r X1 e V(1) -- (Ik - I k - 1 ) V k ,

1E

V 1 a n d ~o1 - ( I - M k _ l ) X

1 E V~,

del~ne

(A~kl)~l, ~pl) _ a(~l, ~pl) _ / a ( x ) ~ l .

~1.

fl Here, the bilinear form a(., .) is equivalent to the H I - i n n e r product. there are positive constants ri such that

vlh~-2[lr

2 _< (A~kl)r r

r2h~-2llr

Then

2

Since a(.,.) is equivalent to the H~-inner product, there are two

Proof:

positive c o n s t a n t s ~l and ~2 such t h a t

~111~)1[[2 < (A~)$I r m

~

--

h11r

l"

Using the n o r m equivalence (7.2.8), e s t i m a t e (7.2.11), and the inverse inequality, we o b t a i n (with possibly different constants 7"1 and r2),

nh~-2[1r

_< (A~)r

r

(7.2.12)

r2h-~2llr

T h e above inequalities conclude the l e m m a . L e m m a 4. Given v, let v( k)' - (Trk - 7rk-1)v. There exists a sutticiently s m a l l constant 7"0 > 0 such that if the approximate projections Q~ satisfy (7.1.1) with r e (0, 7"0) (see (7.2.4) a n d (7.2.5)), then

J

Ilvll

h;- llv

k=0

II02.

(7.2.13)

Proof." Let v E V. S t a r t i n g with v (J) - v, for s - J down to 1, one defines V(S-- 1) __

( I , _ 1 + Q ~ - I (I~ - I~_ 1)v(~) - 7r._ 1 V. T h e n the d e c o m p o s i t i o n v

in t e r m s of entries in Vs1 - (I -- Q~-I) V(1) _ (i - Q~_ 1)(18 - 18_1 ) v reads as

d v - v (~ + Z

v(/)l'

v(')l - v(S) - v ( s - l ) .

(7.2.14)

j=l F r o m the r e p r e s e n t a t i o n v (~)~ - v (~) - v ( s - l ) - (Q~ - Q , - 1 ) + es - es-1, one arrives at the e s t i m a t e

J k=0

<_ C

d

J

2 E h~21](Q~- Qk-~)vll2o + c E h~2I[ek[Io

k=0

J

_< C(1 + r 2) ~

_< Cllvll

.

k-O

k=0

h ~ 2 [ l ( O k - Qk_l)V[I 2o

Wavelet-Like Multilevel Preconditioning

83

Here we have used the norm equivalence result of Oswald [30] and Lemma 1. Thus, it suffices to establish the upper bound" J

Ilvll~ ~ c ~ h~-211v(k)~1102.

(7.2.15)

k=0

It is possible to give a direct proof for (7.2.15) by using the strengthened Cauchy-Schwarz inequality (see Vassilevski and Wang [41]). Here we would like to adopt an alternative approach by using the following characterization of the H~-norm for finite element functions: J

Ilvll~ ~

inf

,

(7.2.16)

Eh;=llvkll02"

v= ~'~ vk, v k E V k

k=0

k---O

A proof of the above equivalence can be found from [30]. Thus, for the J

particular decomposition v -

~ v (k)l, v (k)~ C Vk, one immediately has k=0 J

ttvllx2 ~ c ~ h~-211v(k)~llo2, k=O

which completes the proof of the lemma. 7.3

Stability analysis

Here we study the Riesz property of the wavelet-like multilevel hierarchical basis defined in (7.1.3). For any v E V, let J

v-

~

co,ir ~ + ~

xi EA/'o

y~

ck,i(I-Q~_l)r

k)

(7.3.1)

k=lx,EX~x)

be its representation with respect to the given wavelet basis. The corresponding coefficient norm of v is given by / 1/2 JJJvrir-

h0

E x, E.Afo

+

9' h d-2 k-i

,

(7.3.2)

x, E./V'(X)

where d - 2 or 3 according to the number of space variables. Our main result in this section is the following norm equivalence"

P. Vassilevski and J. Wang

84

T h e o r e m 5. There exists a small (but fixed) r0 > 0 such that if the approximate projections Q~ satisfy (7.1.1) with 7" E (0, vo), then there are positive constants Cl and c2 satisfying

c~lllvlll 2 _< Ilvll~ _< c2111vlll2

v v E V.

(7.3.3)

In other words, the modified hierarchical basis is a stable Riesz basis for the second order elliptic and Stokes problems. The equivalence relation (7.3.3) shall be abbreviated as IIIvlll2 ~_ Ilvll~, Proofi

We first rewrite (7.3.1) as follows: J

- ~ v(~) ',

(7.3.4)

k=0

where with Q~_ - 0 V(]g)I ----

~/.

c k , i ( I - Q~_ 1)elk) E V~ .

(7.3.5)

xiEAfk(1)

Furthermore, by letting r

E

ck,,r k) e V(1), we see that v(k)' =

xiEAf(1)

( I - Q~_I)r (k). Thus, by using (7.2.8)in Lemma 2 (with r Mk-1 -- Q~-I) we obtain

IIr Since r

~ II~(k)'ll~o .

(7.3.6)

E V(1), then

IIr

~h~

~

c~,,.

x, EA/'( ~)

Combining the above with (7.3.6) yields J

IIIvlll~ ~ ~ h;~ll~(k)'tlo~. k=O

This, together with Lemma 4, completes the proof of the theorem. R e m a r k 3: For any fixed a E (0, 1], define i/2

mvmo -

h ~o- ~~ F~ c~,, + ~ h ~- ~~ F_, c~ ,, xiEAfo

k=l

0 and

x,,E

85

Wavelet-Like Multilevel Preconditioning

with d -

2 or 3. Then, for sufficiently small r, one has

IIIvlll,, ~ Ilvll~,

(7.3.7)

for all finite element functions v E Vj. Here I1" I1~ denotes the H~'(a)norm defined by interpolating H0i(f~) with L2(f~). The constants in the norm equivalence depend on cr as indicated in Corollary 1. For ~r = 0, the equivalence (7.3.7) holds true provided that the tolerance satisfies r < r0J-1 for some small 70. w

I m p l e m e n t a t i o n s for a m o d e l p r o b l e m

In this section we discuss techniques which implement the wavelet-like multilevel hierarchical basis for approximate solutions of PDEs. For simplicity, we consider the self-adjoint second-order elliptic equation discussed in Section 4. The Stokes equation can be covered in a similar manner. The bilinear form under consideration is defined as follows: a(~, r

- f a(x)~z~. X7r

V ~p, r E H0i(Ft).

(S.1)

s

Let d(~, r

f - L V ~ . V r be the Dirichlet form defined

on

Hl(f~)

x

H01(f~).

Since the two bilinear forms a(., .) and d(., .) are equivalent, we have from (7.3.3) that

cl[llvlll 2 _< a(v, v) < c~lllvlll ~ for some positive constants

C1

(8.2)

and c2. Moreover, the following result holds"

L e m m a 5. Let c > 0 be a parameter and ac(~, r

9V r + / ~

r

V ~, r E H~(Ft).

(8.3)

I f the approximation in (7.1.1) is sufficiently accurate such that 7" < voJ -1 for some constant to, then there exist 7"1 and r2 independent of e and the mesh size hk such that J

J

< a,(v, v)< for any finite element function v E Vj. Here ak - eh~ 2 + 1 and ck,i are the coemcients of v in the expansion (7.3.1).

P. Vassilevski and J. Wang

86

The spectral bounds (8.4) can be used for the bilinear form arising from discretizing time-dependent Stokes problems. The appearance of e is due to the time stepping parameter At. Similar bilinear forms can be obtained for the pressure unknown by eliminating the vector unknown u in the steadystate Stokes equation. For more details, we refer to Bramble and Pasciak [7]. We remark that the bilinear form a~(., .) is equivalent to the Dirichlet form d(., .) if the parameter e is bounded away from zero by a fixed positive constant Co (i.e., e >_ Co). In this case, the equivalence (8.4) holds true under the condition of Theorem 5. In practical computations, one has two alternatives in solving the discrete problem (2.6) with the wavelet-like basis presented in the previous sections. The first makes use of the explicit form of the basis functions (I - Q]-1)r j) to assemble the corresponding stiffness matrix. The second uses the stiffness matrix assembled from using the standard nodal basis and then performs a change of basis in the process of iterations. The first approach has a difficulty in that the assembly of the global stiffness matrix is no longer local (element-wise). In general, one is recommended to adopt the second approach because the corresponding stiffness matrix is much easier to assemble. In this case, the wavelet-like hierarchical basis actually provides a preconditioning technique for solving the discrete problem (2.6). The rest of this section will describe some preconditioning procedures for the model problem. A more detailed discussion can be found in [41] and [42]. 8.1

Preconditioners

Here we outline two preconditioners for the elliptic operator A (k) : Vk --* Vk arising from the bilinear form a(.,.). The preconditioners will be constructed by using the following wavelet-like multilevel hierarchical decomposition of Vk:

= Voe vl e v) e . . . e

v:,

where

Vj 1 -- ( I - Q ; _ I ) ( I j - I j _ I ) V . The preconditioners (to be defined below) shall be called AWM-HB (Approximate-Wavelet Modified Hierarchical Basis) preconditioners. The following operators are needed in the construction of the AWM-HB preconditioners: 9 In each coordinate space V~, there exists a discretization operator

A~kl)" Vk1 --* Vk1 as the restriction of A (k) onto the subspace V~ defined by

(A~kl)r ~pl) _ a(r

~1)

~, ~1,

~1 E Vk1.

(8.1.1)

87

Wavelet-Like Multilevel Preconditioning

9 Similarly, we define A ~ )" -

Vk-1 ~

V1

and A ~ )" V~ V~eVk_I

-

---* V k - 1

by

eVk 1.

(8.1.2)

Thus, the operator A (k) naturally admits the following partition" A(k)--[

A~kl) A ~ )

A~kl) A (k-l)

] } Vkl } Vk-1.

(813) " "

Let B~k) 1 be given approximations (symmetric positive definite operators) to A~ ) such that for some positive constant bl the following holds" (A~kl)~pl, ~ 1 ) ~ (B~)~I, ~ 1 ) < (1 + bl)(A~kl)9~1, 9~1), ~/~1 E Vk1. (8.1.4) Let A - A(J) be the discretization operator for which a preconditioner is necessary in practical computation. The following explains how the preconditioners can be constructed based on the block structure of A(~) in (8.1.3). Definition 3. (Multiplicative A W M - H B preconditioners) Let B - B (J) be the multiplicatiye A W M - H B preconditioner of A. It is defined as follows: 9 Set B (~ - A (~ 9 For k -

1 , . . . , J , set

B(k) -

B~ A21)

0 B (k-l)

I 0

Bll

I

.

Definition 4. (Additive A W M - H B preconditioners) Let D - D (J) be the additive A W M - H B preconditioner of A. defined as follows: 9 Set D (~ - A (~ 9 For k -

1 , . . . , J , set 0

0 ]

D (k-l)

"

It is

P. Vassilevski and J. Wang

88 8.2

M a i n r e s u l t s for t h e A W M - H B

preconditioners

A spectral equivalence between A and its preconditioners B and D has been established in [41]. The result can be stated as follows. If the tolerance in (7.1.1) is sufficiently small such that r _< 7"0 for some small r0, then there are two absolute constants cl and c2 satisfying

cl(Sv, v) <_ (Av, v) <_c2(Sv, v)

g v C Vj,

(8.2.1)

where S - B (J) or D (a). The estimate (8.2.1) is based on the following results: (A) There exists a constant (7"N ) 0 such that J

IQovl~ + ~ 22'11(Q,- Q,-~)vll~ ~_ ~NII~II~

V v e V.

(8.2.2)

s--1

(B) There exist constants ~ri > 0 and (5 C (0, 1) (in fact, if h i - 7hi-l,1 then ~ - ~22) such that the following strengthened Cauchy-Schwarz inequality holds for any i < j"

a(~i, ~j)2 _ o.ir

ii,jl12o

V~iCV~, ~j C ~ .

(8.2.3)

Here Aj - O ( h - f 2) is the largest eigenvalue of the operator A(J). The inequalities in (A) and (B) have been respectively verified by Oswald [30] and Yserentant [44, 45]. One important feature of the partition (8.1.3) is that the block A~ ) is well conditioned; this can be seen from Lemma 3. In particular, the block A~ ) is spectrally equivalent to its diagonal part. Thus, the J acobi preconditioner would be a good choice for B ~ ) (see (8.1.4))in approximating A~ ). R e m a r k 4: If one does not assume the strengthened equality (B), then the estimate (8.2.1) for S = B still cl = O(1) and c2 = O(log 2 25). In the case S = D, not required in the equivalence (8.2.1). See Griebel related details. 8.3

Cauchy-Schwarz inholds with constants the condition (B) is and Oswald [19] for

On t h e a p p r o x i m a t e L 2 - p r o j e c t i o n

Denote by Q~-I the approximate L2-projections onto the subspace Vk-1. We begin with describing algorithms for computing the action of Q~-l. For any v E V(1) , let v - [ V0l ] }N'k\A/'k-1 }A/k- 1

beitscoefficientvector

with respect to the standard nodal basis of Vk; the second block-component

89

Wavelet-Like Multilevel Preconditioning

of v is zero since v vanishes on Ark-1. The operator Q~_I can be designed by approximately solving the following equation: (Qk-l V, w) -- (v, w),

V w E

(8.3.1)

Vk-1.

1--.. 0.9---. 0.80.7--. 0.60.5-. 0.40.30.20.12

20 0

0

F i g u r e 1. Plot of a HB function (no modification).

Let Ik_l

--"

(with the abbreviation

}./V'k-1

12-

and /kk-1 -- Ik_rl be the natural coarse-to-fine, and respectively, fine-tocoarse transformation matrices. For example, if the nodal basis coefficient vector of a function v2 E Vk-1 in terms of the nodal basis of Vk-1 is v2, then its coefficient vector with respect to the nodal basis of Vk (note that V2

C Vk-1 C Vk) will be

k

Ik-lV2

--

[

t.

J12v2 V2

}JV'k \ .hi'k-1 }-~k-1

"

Denote by Gk -- {(r r EAf~ the mass (or Gram) matrix at the k-th level. Then (8.3.1) admits the following matrix-vector form:

w2rGk-lV2-

(/k_ 1

w2)TGkv,

V W 2.

Here v2 and w2 are the nodal coefficient vectors of Qk-1 v and w E Va-1 at the ( k - 1)th level respectively. Therefore, one needs to solve the following mass matrix problem" Gk-lV2 - - Ik-lGkv. (8.3.2)

P. Vassilevski and J. Wang

90

In other words, the exact L2-projection Qk-lv is actually given by

G-~l_lI~-lakv. Hence

(6;!11k

IIQ -,vll0 -

l Ck v) T a k _ l ( -1

V)

k

(8.3.3)

-

= IIGS 89I,k- 1G~ vii 2.

1-,

0 . 8 "--

0.6-, 0.4--

0 . 2 ""

0 "-.

-0.2 20

15 ~

10

~ 0

0

10

15

20

Figure 2. Plot of a wavelet-modified HB function; m = 2.

To have a computationally feasible basis, we replace G~-ll by some approximations (~-11 whose action can be computed by simple iterative methods applied to (8.3.2). Such iterative methods lead to the following polynomial approximations of G~-11"

Gk!1 - [ I - ?rm (ak_1)]ak11, where 7rm is a polynomial of degree m >_ 1. The polynomial ~rm also satisfies 7rm(0) - 1 and 0 _< 7r,~(t) < 1 for t E [a, fl], where the latter interval contains the spectrum of the mass matrix Gk-1. Since Gk-1 is

91

Wavelet-Like Multilevel Preconditioning

well conditioned, one can choose the interval [a,/3] independent of k. Thus, the polynomial degree m can be chosen to be mesh-independent so that a given prescribed accuracy r > 0 in (7.1.1) is guaranteed. More precisely, given a tolerance r > 0, one can choose m = re(r) satisfying

i]G~_I (G-kll - O-k:l) Ik-lGkv]l = [[G~_lgrm(Gk-1)Gklllk-lGkv[I

ilQ~_,v-e~-,,ilo

-

< -

=

max

7rm(tlllG-~)lI~-lGkvll

m~x

~',,,(t)ilek-,vllo.

te[o,,~']

tED,#]

Here we have used identity (8.3.2) and the properties of ~rm. The last estimate implies the validity of (7.1.1) with r > max 7rm (t). - tED,~] A simple choice of ~rm(t) is the truncated series m-1

(1 - 7r.~(t))t -1 - p . ~ - l ( t )

- fl-1 E (1 - fit) k,

(8.3.4)

k=0

which yields G~-_I1 - p m - l ( G k - 1 ) . from the following expansion:

We remark that (8.3.4) was obtained

OO

1 - tf1-1E(1

-tfl-1)

k

t E [a fl] ~

9

k=O

With the above choice of the polynomial ~m(t), we have ~m(t) -

1 - tp~_l(t)

- tZ - ~ ~

(1 - Z - ~ t ) k - (1 - Z - ' t ) ~ .

k>m

It follows that max

~m(t)-

1-

In general by careful selection of ~m, we have max 7rm(t) < Cq m for some ' te[~,Z] constants C > 0 and q C (0, 1), both independent of k. Since the restriction on r was that r be sufficiently small, one must have m - O(log r - l ) .

(8.3.5)

The requirement (8.3.5) obviously imposes a very mild restriction on m. In practice, one expects to use reasonably small m (e.g., m - 1, 2). This

P. Vassilevski and J. Wang

92

observation is confirmed by the numerical experiments performed in Vassilevski and Wang [42]. We show in Figure 1 a typical plot of a nodal basis function of V(1), and in Figure 2, we show its approximate-wavelet modification for m - 2. The conjugate gradient method was employed to provide polynomial approximations for the solution of the mass-matrix problem (8.3.2). 8.4

M a t r i x f o r m u l a t i o n s of t h e A W M - H B p r e c o n d i t i o n e r s

We now turn to the description of the multiplicative and additive AWM-HB methods in a matrix-vector form. Let us first derive matrix representations for the operators A~ ), A ~ ), and A(~) introduced in (8.1.1)and (8.1.2). In what follows of this section, capital letters without overhats will denote matrices corresponding to the standard nodal basis of the underlined finite element space. For example, A (k) denotes the standard nodal basis stiffness matrix with entries ra~(k) ~9i , l(k),~ q~j )~,,xj~fk. For any v E Vk and its nodal coefficient vector v, we decompose v as follows:

where w2 E Vk- 1 is uniquely determined as w2 - Ik- 1v + Q~_ 1(Ik - Ik- 1)v. Our goal is to find a vector representation for components of v. Since the above decomposition is direct, it is clear that there are vectors 91 and v2 satisfying

v-(I-I:-lG-~l-lI:-lGk)[ 91 ]0

}Afk-1}A/' \A/'k-1 k + i:_192.

[,1]

(8.4.1)

The vectors 91 and 92 represent the components of our wavelet-modified two-level HB coefficient vector ~ -

"~2

of v.

Now, consider the following problem A v - d,

(8.4.2)

which is in the standard nodal basis matrix-vector form. We transform it into the approximate wavelet modified two-level HB by testing (8.4.2)

(I-I~_IG~llI~-IGk)""[ @1

[ and I~ 1@2 for 0 J arbitrary @1 and @2. By doing so, we get the following two-by-two block system for the approximate wavelet modified two-level HB components of (denoted by Vl and ~r with the two components

k

,1

Wavelet-Like Multilevel Preconditioning

93

where ^

A~I)-[I [I

0]

(

) ( ~" )[I] i - a ~ 2 _ ~ a ; ~ l I 2 -1 A(~) i-i2_~a-i!1~2-~a~ o

O] .(I-

GkI~_IG-~I_II~ -1)

A(k)I~_l,

A~kl) - I~-IA (k) I - I~_~a;~_lI~-lak

o

^ k 1 _ A(k-1) A~k2)I k - 1A(k )Ik_ Note that having computed Vl and v2, the solution v of (8.4.2) can be recovered by using the formula (8.4.1), i.e., V -- g l V l nu Y2"r

where r,

-

I-ILIO;_~II~-*c~

Y~

= I~_1.

,-

~=

o

'

We have, v-Y~,

, Y-[Y1, Y,], Y1-Y1 (~), y~-y~(~)

The transformed right-hand side vectors of (8.4.3) read similarly as follows: dl d2

-

[I O] (I-GkI~_l(l-~llI~ / ~ - l d - yTd.

-1)

d-

yTd,

Therefore, the multiplicative AWM-HB preconditioner B (k) from Definition 3, starting with B(~ = A (~ takes the following block-matrix form:

~ ][ 0 1

844,

The preconditioner B (k) is related to/~(k) in the same way as A (k) to ~(k). More precisely, one has

~(k)_

[y~,y2]rB(k)[yl,y2] ' B(k)-' _ [yl,y2]~(~)-'[yl,y2]r"

We will show below that the inverse actions of B (k) can be computed only via the actions of A (k), Y1, Y2, and yT, yT in addition to the inverse actions o f / ~ ).

'

94

P. Vassilevski and J. Wang

We point out that (8.4.4) has precisely the same form as the algebraic multilevel method studied by Vassilevski [37] (see also Axelsson and Vassilevski [4] and Vassilevski [38]). Observe that, in (8.4.4) , Bii ^(k) is an appropriately scaled approximation of A^~]). We have shown that A~])is well conditioned (see Lemma 3). Thus, it is possible to utilize some simple polynomial approximation B~k1) for A(i]) in the implementation. However, in order to take into account any possible jumps in the coefficient of the differential operator, it would be preferable to compute the diagonal part of .4~]). This is computationally feasible since the basis functions of V~ - ( I - Q~_i)V (i) have reasonably narrow support if m is not too large, which should be the case in practice. Nevertheless, one can employ in actual implementation the CG method to compute reliable approximate actions of ~ ] ) - 1 . A l g o r i t h m 1 (Computing inverse actions of B (k)) The inverse actions of B (k) are computed by solving the system B(k)w = d, with the change of basis w w

-Yi~

dl d2

-Y~d, = Y~d,

^

Y@. Namely, by selling

+Y2~2- [Y~, Y2] ~2

w -- B(k)-ld is computed via the solution of B(k)~r -- ~l as follows: Forward recurrence: 1. compute zi - ~il ~(k)-I ^di, 2. change the basis; z.e., compute z - Yizi, " ^ ~;(k ).. A(k)z), 3. compute d2 "- d 2 - ~2i zi - y T ( d ~. compute ~ r 2 - B(k-i)-l~t2, 5. change the basis, i.e., compute v -

Y2@2.

Backward recurrence: 1. update the fine-grid residual, i.e., compute

di^ "- di^ _ .~(k)~,122 - y T ( d - A(k)Y2@2) - y T ( d - A(k)v),

95

Wavelet-Like Multilevel Preconditioning 2. compute ~r 1 -3. get the solution by w

-

Y l ' ~ r 1 + Y2'~r2 -

Y1'~"1

+ V.

End. Note that the above algorithm requires only the actions of the standard stiffness matrix A (k), the actions of the transformation matrices Y1 and }I2 and their transpositions yT and yT, the inverse action of/3~), and some suitable approximations to the well conditioned matrices .4~). The actions of y - 1 are not required in the algorithm. We now formulate the solution procedure for one preconditioning step using the multiplicative AWM-HB preconditioner B - B(J). A l g o r i t h m 2 (Multiplicative A WM-HB preconditioning) Given the problem Bv- d Initiate d(J) - d. F o r w a r d r e c u r r e n c e . For k -

d down to 1 perform:

1. Compute ~k)_

Ix

O]

(I- GkI~_IG;llI~-1) d (k),

2. Solve 1 __ d~ ^ k), B^(k)~r ll

3. Transform basis

W--(I-- Ik-lGkllik-l~k) [ ~W1

}d~fk_l}'/~fk\'Afk-1

~. Coarse-grid defect restriction d(k-1)

__

= 5. Set k -

k-1.

i~-ld(k) -- ]-i21 7(k)^Wl / ~ - l ( d ( k ) - A(k)w),

If k > O go to (1), else,

6. Solve on the coarsest level A(O)x(O)_ d(O).

P. Vassilevski and J. Wang

96 Backward

recurrence.

I. Interpolate result: Set k := k + 1 and compute x(k)_ i~_lx(k-1),

2. Update fine-grid residual:

._

a~)

Z(~)v(k-~)

-- ~l~k)_[I O](I-Gke2_lS"~llI~-l)A(k)x (k) 3. Solve B1 l^(k)~ 1 _ ~ k ) ,

~. Change the basis 0

5. Finally, set x (k) = x (k) + w.

6. Set k := k + 1. If k < J go to Step (1), else s e t v = x (J).

End. Similarly, one preconditioning solution step for the additive A W M - H B preconditioner D - D (J) takes the following form: 3 (Additive A WM-HB Preconditioning) Given the problem Dv=d

Algorithm

Initiate d (J) -- d. Forward

recurrence.

For k = J down to 1 perform:

1. Compute

~k)_ IX 0] (I--GkX~_lVkll Ik-1) d (k),

97

Wavelet-Like Multilevel Preconditioning 2. Solve

B1 l^(k),~r1 __ ~k)

3. Transform basis

X(k)

--(I--Ikk_lS;!1I~-lGk)[ Wl] 0

\Nk-1 }&-I

~. Coarse-grid defect restriction

d ( k - 1 ) _ i~-ld(k) 5. Set k - k - 1 .

I f k > O go to (1), else

6. Solve on the c o a r s e s t level

A(O)x(O) _ d(O). Backward recurrence. 1. Interpolate result: Set k "- k + 1 and compute

w-

I~_1 x(k-1),

2. Update at level k

x (k) - x (k) + w, 3. Set k "- k + 1. I l k < J go to Step (1), else set

v - x (J). End.

For both the additive and multiplicative preconditioners, it is readily seen that the above implementations require only actions of the stiffness matrices A (k), the mass matrices G (k), and the transformation matrices I~_ 1 and I kk-1 . The approximate inverse actions of A^ ~ ) can be computed via some inner iterative algorithms. Similarly, the action of G~-_I1 can be computed as approximate solutions of the corresponding mass-matrix problem using m steps of some simple iterative methods. Therefore, at each discretization level k, one performs a number of arithmetic operations proportional to the degrees of freedom at that level, denoted by nk. In the case of local mesh refinement, the corresponding operations involve only the stiffness and mass matrices computed for the subdomains where local

P. Vassilevski and J. Wang

98

refinement was made. Hence, even in the case of locally refined meshes, the cost of the AWM-HB methods is proportional to nj. The proportionality constant depends linearly on m = O(log r - l ) , but is independent of J (or h). Some numerical results for the AWM-HB preconditioners can be found from Vassilevski and Wang [42]. A performance comparison with the BPX method [10] and Stevenson's method [36] on more difficult elliptic problems in three dimensions, and in other applications such as interface domain decomposition preconditioning, is yet to be seen. w

Numerical experiments

In this section we present some numerical results to illustrate the efficiency of the method discussed in Section 7. Consider the boundary value problem of seeking u satisfying s

Vu u

-f - g

inFt, on c9~.

(9.1)

Here, / E L2(~), g E H 89 and b - [bib2] are given single-valued or vector-valued functions. We assume that all given functions are sufficiently smooth on their domains. For simplicity, we take ~ to be a square domain and g - 0 on c9~. If u is a solution of (9.1), then it solves the following problem: .A4cu - - S V . (b s

+ s

- - 6 V . (b f) + f

in ~,

(9.2)

subject to the boundary condition u - 0 on 0~. Here 5 > 0 is a parameter. The purpose of considering the problem (9.2) is to get the so-called streamline derivative o~, _ b. Uu in a variational formula for u More precisely, by testing (9.2) against any r E H ] ( ~ ) one obtains

b (u, r

-

r e ) + (r b. b. r e )

+ 6(_b.

b. r e ) (9.3)

-- (f6, r for all r E H~(fl). Here ]'6 - f 9.1

5U. (b f).

Galerkin discretization

The Galerkin method for the approximation of u is based on the variational finite element space of problem (9.3). Let V - Vh be a C~ piecewise polynomials corresponding to a quasiuniform triangulation Th of

Wavelet-Like Multilevel Preconditioning

99

f~. The Galerkin approximation is a function Uh E Vh such that

b~(~h, r

- ~(V~h, r e ) + (r -e5 E

V~) + 6(b. V~, b. r e )

f Auh(b.Vr

TETh T

(9.1.1)

= (f6, r for all r E Vh. For continuous piecewise linear functions, one has AUh -- 0 on each element. It follows that the discrete problem seeks Uh C Vh such that

~(v~h, r e ) + (r b. v~h) + 6(5. v~h, b. r e ) - (f~, r

(9.1.2)

for all r E Vh. The convection term is assumed to satisfy V.b_< 0

in ft.

(9.1.3)

For a convergence analysis of the streamline diffusion finite element approximation Uh, we refer to [23] and [2]. 9.2

Numerical tests

We choose the same test examples as in [2]. Namely, _b- [ ( 1 - x cos a) cos a(1 - y sin a) sin a], for various angles c~. Note that V. _b- - 1 . The right hand side f is chosen so that u - x(1 - x)y(1 - y) is the exact solution. Thus, the right-hand side function f is e-dependent. The stopping criterion is that the relative error of the residual be less than 10 - s in the discrete L2-norm. The objective is to test the number of iterations in the solution procedure by using the wavelet-like hierarchical basis. The matrix form of the discretized problem (9.1.2) reads as follows" Au-f. (9.2.1) Here A is a nonsymmetric matrix. For small e, it is very difficult to find a good preconditioner for A. It was seen in [2] that a block-ILU factorization method turns out to be very robust with respect to arbitrary positive e, though very little is known theoretically about this good performance. For any finite element function v, we use the bold face v to denote the vector with respect to the nodal basis and 9 the vector representation with respect to the modified hierarchical basis. In the implementation, the modified hierarchical basis is employed to provide a preconditioner for the global stiffness matrix A as follows. Let Y be the transformation from -~ to v such that v - Yg. The preconditioner is given by P - ( y y T ) - I with y T being the conjugate of Y. We discuss the computation of Y and y T .

100

P. Vassilevski and J. Wang

Let

and y(k)_ g-l"

(9.2.3)

Here, G (k) stands for the mass matrix at kth discretization level, I2_ 1 is the natural coarse-to-fine interpolation matrix from ( k - 1)th grid to the kth one a n d / ~ - i _ (/~-l)T" Also, (~(8)-1 is an approximate inverse of the mass matrix G(') 9 For example, a good choice for G (8)-1 w would be an approximate solution of G(~)x = w by polynomial iterative methods. In practice, the J acobi and conjugate gradient methods are good candidates. The numerical results in this section are based on the Jacobi iterative method with two iterations. A l g o r i t h m 4 For any given ~ - (9~J) ,''', Y 9 is computed as follows:

~1), ~(0))T,

the action v -

9 Set v(~ - 9(0), 9 Fork-1

toJ

do

9 v - v(J),

Algorithm 5 yTd

For any given d is computed as follows:

(d~ s) , . . - , d~1), d(~ T, the action a

-

. Set d(J) - d (~ 9 Fork-J

down t o l

do

~(k-i) 9 a-(a(J )

a(~ 9

Once the preconditioner P - ( y y T ) - I is known, one can solve the discrete problem (9.2.1) by using the generalized conjugate gradient method employed in [2]. We comment that this simple preconditioner may not work well for convection-dominated diffusion problems. This fact can be seen from the numerical results illustrated in Tables 1-4.

Wavelet-Like Multilevel Preconditioning

101

T a b l e 1. Iteration counts for h -1 = 64, a = 75 ~

iter

c-1 6 - 0.001 30

e-O.1 6 - 0.01

28

c-

c - 10 -2 6-0.1 42

10 -3 5 - 1.0 59

...........

c - - 10 -4 5--10

74

T a b l e 2. Iteration counts for h -1 --64, a = 105 ~

iter

c-1 5 - 0.001 31

c-O1 6 - 0.01 27

c-lO

2

c-lO -a 6 - 1.0

6 - 0.1 70

e - - 10 . 4 6--10

........

*

.

....

T a b l e 3. Iteration counts for h -1 = 128, c~ = 75 ~

iter

e-1 5 - 0.0025

c -0.1 5 - 0.025

32

29

....

c - 10 -2 5 -- 0.25 48

T a b l e 4. Iteration counts for h -1 = 128, a = 105 ~

iter

c-1 6 - 0.0025 32

c-O.1 5 - 0.025 29

c-

10 -2

6 - 0.25 71

In Tables 1-4, " , " is used to indicate a non-convergence in 100 iterations. It is clear t h a t the use of the wavelet-like hierarchical basis gives an efficient p r e c o n d i t i o n e r if the p r o b l e m is not c o n v e c t i o n - d o m i n a t e d . T h e present i m p l e m e n t a t i o n of the modified hierarchical basis is c o m p a rable with the additive p r e c o n d i t i o n i n g m e t h o d discussed in [42]. T h e inform a t i o n is also c o n t a i n e d in A l g o r i t h m 3 of this p a p e r with B"(k) l l _ Ch.~2ik , where Ik s t a n d s for the identity m a t r i x . Acknowledgments. T h e research of Vassilevski is s u p p o r t e d in p a r t by the U.S. NSF g r a n t INT-95-06184 and also in p a r t by B u l g a r i a n M i n i s t r y for E d u c a t i o n g r a n t M M - 4 1 5 , 1994. T h e research of W a n g is s u p p o r t e d in p a r t by the NSF g r a n t INT-93-09286. T h e a u t h o r s are also grateful to Dr.

P. Vassilevski and J. Wang

102

Oswald for his careful comments and remarks on the paper. References [1] Arnold, D. N., R. S. Falk, and R. Winther, Preconditioning in H(div) and applications, Math. Comp., to appear. [2] Axelsson, O., V. Eijkhout, B. Polman, and P. S. Vassilevski, Incomplete block-matrix factorization iterative methods for convectiondiffusion problems, BIT 29 (1989), 867-889.

[3]

Axelsson, O. and I. Gustafsson, Preconditioning and two-level multigrid methods of arbitrary degree of approximations, Math. Comp. 40 (1983), 219-242.

[4]

Axelsson, O. and P. S. Vassilevski, Algebraic multilevel preconditioning methods, II, SIAM J. Numer. Anal. 27 (1990), 1569-1590.

[5]

Bank, R. E. and T. Dupont, Analysis of a two-level scheme for solving finite element equations, TechnicM Report CNA-159, Center for Numerical Analysis, The University of Texas at Austin, 1980.

[6] Bramble, J. H., Multigrid Methods, Pitman Res. Notes Math. Ser., vol. 294, Longman, Harlow, U.K., 1995 (2nd Ed.). [7] Bramble, J. H. and J. E. Pasciak, Iterative techniques for time dependent Stokes problems, ISC Report, Texas A&M University, 1994, preprint.

[8]

Bramble, J. H., J. E. Pasciak, J. Wang and J. Xu, Convergence estimates for product iterative methods with applications to domain decomposition, Math. Comp. 57 (1991), 1-21.

[9]

Bramble, J. H., J. E. Pasciak, J. Wang, and J. Xu, Convergence estimates for multigrid algorithms without regularity assumptions, Math. Comp. 5 7 (1991), 23-45.

[10]

Bramble, J. H., J. E. Pasciak, and J. Xu, Parallel multilevel preconditioners, Math. Comp. 55 (1990), 1-22.

[11] Brezzi, F. and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, New York, 1991. [12] Carnicer, J. M., W. Dahmen, and J. M. Pefia, Local decompositions of refinable spaces and wavelets, Appl. Comput. Harmon. Anal. 3 (1996), 127-153.

Wavelet-Like Multilevel Preconditioning

103

[13] Ciarlet, P., The Finite Element Method for Elliptic Problems, North Holland, Amsterdam, 1977. [14] Chui, C. K., An Introduction to Wavelets, Academic Press, Boston, 1992. [15] Dahmen, W. and A. Kunoth, Multilevel preconditioning, Numer. Math. 63 (1992), 315- 344. [16] Dahmen, W., A. Kunoth, and K. Urban, A wavelet Galerkin method for the Stokes equations, Computing 56 (1996), 259-301. [17] Daubechies, I., Ten Lectures on Wavelets, SIAM, Philadelphia, 1992. [18] Girault, V. and P.-A. Raviart, Finite Element Methods for NavierStokes Equations, Springer-Verlag, Berlin, 1986. [19] Griebel, M. and P. Oswald, On the abstract theory of additive and multiplicative Schwarz algorithms, Numer. Math. 70 (1995), 163-180. [20] Griebel M. and P. Oswald, Tensor product type subspace splittings and multilevel iterative methods for anisotropic problems, Adv. Comput. Math. 4 (1994), 171-206. [21] Hood, P. and C. Taylor, A numerical solution of the Navier-Stokes equations using the finite element technique, Comput. Fluids 1 (1973), 73-100. [22] Jaffard, S., Wavelet methods for fast resolution of elliptic problems, SIAM J. Numer. Anal. 29 (1992), 965-986. [23] Johnson, C., Numerical Solution of Partial Differential Equations by the Finite Element Method, Cambridge Univ. Press, Cambridge, 1994. [24] Kotyczka, U. and P. Oswald, Piecewise linear prewavelets of small support, in Approximation Theory VIII, vol. 2, C. K. Chui and L. L. Schumaker (eds.), World Scientific, Singapore, 1995, pp. 235-242. [25] Lorentz, R. and P. Oswald, Constructing economical Riesz bases for Sobolev spaces, presented at the 9th Domain Decomposition Conference held in Bergen, Norway, June 3-8, 1996. [26] Mallat, S., Multiresolution approximations and wavelet orthonormal bases of L2(IR), Trans. Amer. Math. Soc. 315 (1989), 69-88. [27] Maubach, J. M., Local bisection refinement for n-simplicial grids generated by reflection, SIAM J. Sci. Comput. 16 (1995), 210-227.

104

P. Vassilevski and J. Wang

[28] Mitchell, W. F., Optimal multilevel iterative methods for adaptive grids, SIAM J. Sci. Star. Comput. 13 (1992), 146-167. [29] Ong, M.-E. G., Hierarchical basis preconditioners in three dimensions, SIAM J. Sci. Comput., to appear. [30] Oswald, P., Multilevel Finite Element Approximation: Theory and Applications, Teubner Skr. Numer., Teubner, Stuttgart, 1994. [31] Rusten, T. and R. Winther, A preconditioned iterative method for saddle-point problems, SIAM J. Matrix Anal. Appl. 13 (1992), 887904. [32] Schatz, A. H., An observation concerning Ritz-Galerkin methods with indefinite bilinear forms, Math. Comp. 28 (1974), 959-962. [33] Schatz, A. H. and J. Wang, Some new error estimates for RitzGalerkin methods with minimal regularity assumptions, Math. Comp. 65 (1996), 19-27. [34] Silvester, D. S. and A. J. Wathen, Fast iterative solution of stabilized Stokes systems, part II: using general block preconditioners, SIAM J. Numer. Anal. 31 (1994), 1-16. [35] Stevenson, R., Robustness of the additive and multiplicative frequency decomposition multilevel method, Computing 54 (1995), 331-346. [36] Stevenson, R., A robust hierarchical basis preconditioner on general meshes, Technical Report # 9533, Department of Mathematics, University of Nijmegen, Nijmegen, The Netherlands, 1995. [37] Vassilevski, P. S., Nearly optimal iterative methods for solving finite element elliptic equations based on the multilevel splitting of the matrix, Technical Report # 1989-09, Institute for Scientific Computation, University of Wyoming, Laramie, WY, USA, 1989. [38] Vassilevski, P. S., Hybrid V-cycle algebraic multilevel preconditioners, Math. Comp. 58 (1992), 489-512. [39] Vassilevski, P. S., On two ways of stabilizing the hierarchical basis multilevel method, SIAM Rev., to appear. [40] Vassilevski, P. S. and J. Wang, Multilevel iterative methods for mixed finite element discretizations of elliptic problems, Numer. Math. 63 (1992), 235-248.

Wavelet-Like Multilevel Preconditioning

105

[41] Vassilevski, P. S. and J. Wang, Stabilizing the hierarchical basis by approximate wavelets, I: Theory, Numer. Linear Algebra Appl. (1996), to appear. [42] Vassilevski, P. S. and J. Wang, Stabilizing the hierarchical basis by approximate wavelets, II: Implementation and numerical experiments, SIAM J. Sci. Comput. (1996), submitted. [43] Xu, J., Iterative methods by space decomposition and subspace correction, SIAM Reg. 34 (1992), 581-613. [44] Yserentant, H., On the multilevel splitting of finite element spaces, Numer. Math. 49 (1986), 379-412. [45] Yserentant, H., Old and new convergence proofs for multigrid algorithms, in Acta Numerica, Cambridge Univ. Press, New York, 1993, pp. 285--326.

Panayot S. Vassilevski Center of Informatics and Computing Technology Bulgarian Academy of Sciences "Acad. G. Bontchev" street, Block 25 A 1113 Sofia, Bulgaria p an ayot @iscbg.ac ad. bg Junping Wang Department of Mathematics University of Wyoming Laramie, Wyoming 82071 [email protected]