Sparse-Grid Finite-Volume Multigrid for 3D-Problems 1 ... - CiteSeerX
Sparse-Grid Finite-Volume Multigrid for 3D-Problems P. W. Hemker Abstract We introduce a multigrid algorithm for the solution of a second order elliptic equation in three dimensions. For the approximation of the solution we use a partially ordered hierarchy of nite-volume discretisations. We show that there is a relation with semicoarsening and approximation by more-dimensional Haar wavelets. By taking a proper subset of all possible meshes in the hierarchy, a sparse grid nite-volume discretisation can be constructed. The multigrid algorithm consists of a simple damped point-Jacobi relaxation as the smoothing procedure, while the coarse grid correction is made by interpolation from several coarser grids levels. The combination of sparse grids and multigrid with semi-coarsening leads to a relatively small number of degrees of freedom, N , to obtain an accurate approximation, together with an O(N ) method for the solution. The algorithm is symmetric with respect to the three coordinate directions and it is t for combination with adaptive techniques. To analyse the convergence of the multigrid algorithm we develop the necessary Fourier analysis tools. All techniques, designed for 3D-problems, can also be applied for the 2D case, and {for simplicity{ we apply the tools to study the convergence behaviour for the anisotropic Poisson equation for this 2D case.
1 Introduction In this paper we describe the approximation of a function on a nite-volume sparse grid, and a multigrid algorithm for the solution of partial di erential equations in three dimensions. The algorithm is intended for the solution of ow problems described by conservation laws, and therefore nite volumes are a natural choice for the discretisation. But to introduce the main principles, we will restrict the treatment here to second order elliptic equations, and in particular to the anisotropic Poisson equation. In contrast to the usual multigrid approach, we do not use a sequentially ordered set of discretisations on di erent meshes, but we use a partially ordered hierarchy of `semicoarsened' grids as proposed e.g. by Mulder [6, 7] and Naik-VanRosendale [8] or Zenger et al. [3, 9]. As indicated in [9], adaptive `sparse grid' discretisations can be constructed by taking a suitable subset of all possible discretisations in such a hierarchy. However, in contrast to the sparse grid approximation proposed in [3, 9], we base our approximation on nite volumes rather than on nite elements. The multigrid algorithm consists of damped Jacobi relaxation as a smoothing procedure and a coarse grid correction constructed by extrapolation from simultaneous corrections on several coarser grids levels. The algorithm is completely symmetric with respect to the three coordinate directions and it is suitable for combination with adaptive techniques. A description of a data structure to implement such adaptive three-dimensional algorithms is given elsewhere [5]. 1
2 Finite volume sparse grids In this section we introduce nite-volume sparse grids. We show the relation between the approximation by Haar wavelets (when this notion is extended to more dimensions) and the sparse-grid approximation. For the theory of wavelets, multiresolution analysis (MRA) etc. we refer to Daubechies [2].
2.1 The more-dimensional MRA A multidimensional multiresolution analysis of L2 ( ), = IR3 , is a partially ordered set of closed linear subspaces fVnjVn L2( )gn2ZZ3 with the properties:
T
S
(1) n Vn = f0g; n Vn dense L2 ( );3 m (2) f (x)2Vn , f (2 x)2Vn+m 8n2ZZ ; m2E; (3) f (x)2Vn , f (x ? 2?n k)2Vn 8k2ZZ 3 ; n2E; (4) 9 2V0 : f (x ? k)gk2ZZ 3 is a Riesz basis for V0 :
(1)
Here n = (n1 ; n2 ; n3)2ZZ 3 , and we denote jnj = n1 +n2 +n3 ; 2n = (2n1 ; 2n2 ; 2n3 ). We also use the notation 0 = (0; 0; 0)2IN 3 ; x = (x1 ; x2 ; x3 )2IR3 ; 2nx = (2n1 x1 ; 2n2 x2 ; 2n3 x3 ). Further we introduce in IN 3 the unit vectors ek , k = 1; 2; 3; as follows: e1 = (1; 0; 0); e2 = (0; 1; 0); e3 = (0; 0; 1), and we use e = (1; 1; 1). Finally we de ne E = fe1 ; e2; e3g. Although we are particularly interested in the three-dimensional case, generalisation to a di erent number of space dimensions is straightforward. The function (x) in (1.4) is called the the scaling function of the multiresolution analysis.
2.2 Piecewise constant function spaces Let either = IR3 be the three-dimensional Euclidean space, or let = (0; 1)3 IR3 be the open unit cube. For any n2ZZ 3 we introduce the function space Vn , i.e. the space of piecewise constant functions on a uniform grid with meshsize h = (2?n1 ; 2?n2 ; 2?n3 ). These grids are uniformly spaced in each of the three coordinate directions, but possibly with a di erent meshsize in the di erent directions. The volume of these cells is denoted by h3 = 2?jnj . The functions in Vn are all constant in each cell
n;k = [k12?n1 ; (k1 + 1)2?n1 ] [k22?n2 ; (k2 + 1)2?n2 ] [k3 2?n3 ; (k3 + 1)2?n3 ] ; and this family of cells forms the grid n = f n;k j n;k ; k2ZZ 3 g. The family of cell centres or cell nodes is denoted by n = fzn;k j zn;k = (k+ e=2)2?n ; k2ZZ 3 ; zn;k 2 g. The number of these nodes is equal to the dimension of Vn. Apparently, all grids are identi ed by a triple n; the number jnj is called the level of the grid n. Notice that {di erent from classical multigrid{ here and later in our multigrid algorithm, there is a clear distinction between the grid-identi cation index n and the level number jnj. Because for j = 1; 2; 3; (2) Vn Vn+e ; j
we see, by construction, that nesting relations exist between spaces Vn and that the 2
nesting provides a partial ordering1 : Vn Vm , n m: (3) Spaces Vn and Vm or grids n and m that satisfy this nesting relation n < m are called related. The construction of the spaces Vn shows that even a stronger relation holds than (2), namely Vn?e \ Vn?e = Vn?e ?e ; j; k = 1; 2; 3; j 6= k: (4) We see also that for = IR3 the spaces fVngn2ZZ 3 form a multiresolution analysis and that in this case the characteristic function on the unit cube, 2V0 , (x) = 01 ifif xx622
0;0 ;; (5) 0;0 serves as the scaling function. The set f nk j nk (x) = (2n x ? k); k2ZZ 3 g forms a basis in Vn , which corresponds with the usual Haar-basis for the one-dimensional case. In the case = (0; 1)3 we restrict ourselves to Vn with n1 ; n2 ; n3 0 and we see dim(Vn ) = 2jnj . Formally, for = (0; 1)3 and n1 , n2 or n3 negative we de ne Vn = Vn1;n2;n3 by Vn = f0g. For all spaces Vn we introduce the restriction operator Rn : L2 ( ) ! Vn, the L2 projection such that for u2L2 ( ) we have Rn u2Vn and j
j
k
k
(Rn u)(zn;k ) = 2jnj
Z
n k
u(z) dz :
(6)
;
2.3 More-dimensional wavelets
We introduce the wavelet space Wn Vn which consists of all functions in Vn that are not represented in any of the related function spaces on the next coarser level, i.e. they are in Vn but not in Span(Vn?e1 ; Vn?e2 ; Vn?e3 ), or (7) Vn = Wn Span(Vn?e1 ; Vn?e2 ; Vn?e3 ); (8) f0g = Wn \ Span(Vn?e1 ; Vn?e2 ; Vn?e3 ): This means that Wn contains the `di erence information' that is available in the ne grid Vn but not in the the span of the coarser grids Vn?e1 , Vn?e2 and Vn?e3 . In our case, where Vn contains the piecewise constant functions, it is simple to construct the spaces Wn such that (9) Wn ? Span(Vn?e1 ; Vn?e2 ; Vn?e3 ): This makes Wn the orthogonal complement of Span(Vn?(1;0;0) ; Vn?(0;1;0); Vn?(0;0;1)) in Vn , and (8) follows immediately from (9). For IR3 the relation (9) allows a straightforward decomposition of Vn. In the case = IR3 we have2 n1 n2 n3 M M M
Vn (IR3 ) =
j1 =?1 j2 =?1 j3 =?1 1 We de ne the inequalities between triples by
and
n<m,n nm,n
Wj ;
1
< m1 ; n2 < m2 ; n3 < m3 ;
1
m1 ; n2 m2 ; n3 m3 :
(10)
2 Notice that here, in the more-dimensional case, it is convenient to choose an indexing that is di erent from the usual indexing in the well-known one-dimensional case.
3
where all Wj are orthogonal to each other. To handle the bounded domain = (0; 1)3 we introduce the functions Vn0 ( ) Vn( ) which have a zero mean value on , i.e. Vn0 ( ) = fu2Vn( ) j R0 (u) = 0g, and we have3
Vn0 ( ) = and hence
n1 M n2 M n3 M
j1 =0 j2 =0 j3 =0
Vn ( ) = V0 Vn0 = V0
Wj ;
M
(11)
Wj :
(12)
0j n The `di erence information' between the approximations of a function f 2L2 ( ) on two successive levels, Rn f 2Vn on the one hand and Rn?e f 2Vn?e , j = 1; 2; 3; on the other hand, is given by the orthogonal projection Qnf of f onto the orthogonal complement Wn of Vn?e in Vn . This is described in the following theorem. Theorem 2.1 Let the operator Qn : L2 ( ) ! Wn be the orthogonal projection onto Wn , then it follows that Qnu = Rn u ? Rn?e1 u ? Rn?e2 u ? Rn?e3 u + Rn?e2 ?e3 u + Rn?e1 ?e3 u + Rn?e1 ?e2 u ? Rn?e u ; or, equivalently, X ? Qnu = Rnu ? Rn?eu + (13) Rn?e+e u ? Rn?e u : j
j
j
j
j
j =1;2;3
Proof: From (10) or (12) it follows that (possibly neglecting functions in V0 if = (0; 1)3 )
Vn =
so that
Rn u =
and
j n
X
j n
Wj ; Qj u ;
(14)
Rn ? Rn?e1 ? Rn?e2 ? Rn?e3 + Rn?e2 ?e3 + Rn?e1 ?e3 + Rn?e1 ?e2 ? Rn?e =
2
M
X
Qj ?
X
Qj ?
X
Qj ?
X
Qj
j n j n?e3 j n?e2 X j n?e1 X X X + Qj + Qj + Qj ? Qj j n?e2 ?e3 j n?e1 ?e3 j n?e1?e2 X X X Xj n?e = Qj + Qj + Qj + Qj n?e<j X n n?e1 <j n X n?e2<j n X n?e3<j n ? Qj ? Qj ? Qj n?e2?e3<j n n?e1?e3<j n n?e1?e2 <j n = Qn :
3 In the case of a bounded domain, ! = (0; 1)3 , W (0;0;0) ( ) is the zero-function and the functions W(j1 ;j2 ;j3 ) with some zero-index (e.g. W(j1 ;j2 ;0) or W(0;0;j3 ) ) have a shape di erent from those with j > 0.
4
Remarks: In the right-hand-side of equation (13) we recognise the information that can be represented on the levels jnj, jnj ? 3, jnj ? 2, jnj ? 1, respectively. In (13) the information on the level jnj ? 2 and jnj ? 3 can directly be derived from the information on level jnj ? 1, e.g. by Rn?e2 ?e3 u = Rn?e2 ?e3 (Rn?e3 u). Thus, expression (13) describes the di erence information between Rn u and its approximation on the related next coarser grids. Notice that in the two-dimensional case the relation (13) reduces to Qnu = Rnu ? Rn?e1 u ? Rn?e2 u + Rn?e u ; where e = (1; 1), and in the one-dimensional case we have
Qnu = Rn u ? Rn?e u :
(15) (16)
First, in the remainder of this section we restrict ourselves to the case of the unbounded domain = IR3 . The four relations (1.1) to (1.4) imply that also the spaces Wn are scaled versions of one space W0 , f (x)2Wn , f (2?n x)2W0 ; 8n2ZZ 3 ; (17) and, moreover, that they are translation invariant for the discrete translations 2?n ZZ 3 ,
f (x)2W0 , f (x ? k)2W0 ; 8n2ZZ 3 : (18) The relations (7) and (9) make that they are mutually orthogonal spaces, generating all functions of L2 (IR3 ), LWn ?3 WWm for nL6=2 (IRm3;): (19)
n2ZZ n
dense
Summarising, we obtain a nesting between the spaces fVn g and fWn g, 9 9 Vn?e1 > >> >
>
>> >> > = Vn > >> Vn?e3 > >> > >; = Wn Vn+e ; > > Vn+e ?e ?1 > > > > Vn+e ?e +1 > > > > Wn+e ; > Vn?e2 > =
j
j
j
j
j
j
that is essentially more complex than in the classical one-dimensional case, where there is a sequential ordering of the spaces fVn gand fWn g. As soon as we nd a function (x) with the property that (x ? k), k2ZZ 3 , is a basis of We , then by a simple rescaling, we see that (2n x ? k), yields a basis of Wn+e . Such a function is the more-dimensional generalisation Since L2 (IR3 ) is the directo n n+e of a wavelet. n + e sum of these Wn+e , the full collection k (x) j k (x) = (2nx ? k); n; k2ZZ 3 is a basis of L2 (IR3 ). 5
It is easy to check that the more-dimensional wavelet (x)2We , corresponding with the scaling function (x)2V0 , from the previous section, is the three-dimensional checkerboard basis function4 given by (5): 8 = 0 if x62 0;0 ; < (x) = : = +1 if x2 0;0 ; x2 e;k ; jkj even; = ?1 if x2 0;0 ; x2 e;k ; jkj odd: This function is the three-dimensional generalisation of the Haar-wavelet. In wavelet theory the spaces Wn are labelled channels, and the distinct channels are linearly independent. The rst decomposition of an arbitrary function from L2 ( ) P 3 with = IR consists in writing u(x) = n un (x), where each un belongs to the corresponding channel Wn with n2ZZ 3 . Similarly, we can write for functions de ned on = (0; 1)3 the relation (12) and make a similar decomposition in channels. Each subspace Wn+e , n 0, has its natural basis, the standard basis5 n n+e n+e(x) = (2nx ? k)o (x) j k k of functions with a minimal support. The basis function kn+e is a scaled, elementary checkerboard function, that may be characterised either by its support which is a single cell in n or by the centerpoint of this cell, zn;k . For = (0; 1)3 , the exceptions related with the boundary are found in the spaces Wn with a zero index (i.e. n1 n2 n3 = 0). These Wn have basis functions with di erent shapes. They are derived from the corresponding functions for the unbounded case, but T their support is restricted to n?e . Their corresponding nodal points zn?e;k are found on the boundary @ = n , n e, n 6= 0. Taking this into account, both for
= (0; 1)3 and for = IR we may write for each u2L2 ( ) a wavelet expansion X u(x) = an;k (2nx ? k): (20) n;k
2.4 Approximation results
The decompositions (10) or (12) clearly allow the approximation of a suciently smooth function in L2 ( ) by a series with elements in Wj . To obtain an impression of the quality of these expansions we derive some error estimates. As the case where boundaries are present is the more general one, we take = (0; 1)3 . To quantify the error of approximation on , we introduce for u2C 3 ( ) the seminorm6 @ 3 u(x) juj = sup @x @x @x + (21) x2 1 2 3 @ p @ q @ r max sup p;q;r=0;1 x2@ @x1 @x2 @x3 u(x) : Now we derive the following Theorem 2.2 If we consider an expansion of a C 3 ( )-function, u, in piecewise constant functions on the grid n, for an arbitrary n2ZZ 3 , n > 0, and if we write X Rn u = v0 + uj ; 0j n
Notice 2We Ve is a function piecewise constant on e . Notice that in more dimensions we use the indexing kn+e , whereas in the one-dimensional case one usually writes kn . 6 The necessity of the boundary terms in this seminorm is seen immediately if we want to approximate in L2 ( ) smooth functions u2C 3 ( ) that do not satisfy homogeneous boundary conditions. 4
5
6
with v0 2V0 and uj 2Wj , 0 j n, then
kuj kL2( ) 2?3jj j=2 juj ; and we get an estimate for the approximation error
ku ? Rn ukL2( ) 8 7?3=2 (2?3n1 + 2?3n2 + 2?3n3 )1=2 juj 8 7?3=2 (h31=2 + h32=2 + h33=2 ) juj :
(22)
Proof: We take the normalised f ~kj g = f2jj ?ej=2 kj g as a basis in Wj , 0 j n, j= 6 0. Clearly, all these functions are orthogonal to all functions v0 2V0 and mutually they form an orthonormal set (an orthonormal Haar basis) in Wj L2 ( ). We see this as follows
2We ; e kj 2We ; k 2 We ;
) = 0;0 ; e kj ) = 0;k ; k ) = j ?e;k ; or, in other words, kj 2Vj , but kj scales like a basis function in Vj ?e . Hence Z j d = 0 for k 6= m ; 2jj ?ej=2 kj 2jj ?ej=2 m and
Z
2jj ?ej=2 kj 2jj ?ej=2 kj d = 2jj ?ej
Thus, we nd with Now
support( support( support(
Z Z j j jj ?ej d
= 2 k k
Rn u = v0 + uj =
X
X
d = 1:
;
uj ;
0j n
X ajk ~kj = (u; ~kj ) ~kj : k k
ajk = (u; ~kj ) =
Z
u ~kj d =
By Taylor expansion around zj ?e;k , we have
Z
j ?e k
(23)
Z
j ?e k ;
u ~kj d :
u ~kj d 2?2jj jjuj :
j ?e k For j e the point zj ?e;k lies in the interior of and the estimate holds with ;
(24)
@ 3 u(x) juj = sup @x @x @x : x2 1 2 3
For j 6 e, i.e. for kj with a j-component equal to zero, the point zj ?e:k lies on the boundary and the function kj is constant in one direction over the whole domain , and 7
it is of Haar-wavelet type for the non-zero indices (or index). In this situation the same estimate (24) holds with, e.g. if j1 = 0, with
@ 2 u(x) juj = sup @x @x : 2 3
For j = 0 the relation (24) is trivially satis ed. Hence, the estimate (24) holds for j 0 if we use the seminorm (21), and we nd jaj ;k j 2?2jj j juj ; X X kuj k2 = jajk j2 2?4jj jjuj2 = 2?3jj jjuj2 ; k k so that kuj k 2?3jj j=2juj ; and
0 BB BB X X X X 2 2 2 + + ku ? Rn uk = kuj k juj B BB BB j1 > n1 j1 0 j 0 j1 > n1 @ j2 0 j21 > n2 j2 0 or j2 > n2 j3 0
or j3 > n3
And it follows that
j3 0
3 ? 78 2?3n1 + 2?3n2 + 2?3n3 juj2 :
j3 > n 3
1 CC CC CC 2?3jj j CC CA
3=2 2?3n1 =2 + 2?3n2 =2 + 2?3n3=2 juj ku ? Rn uk 78
2
3=2 h31=2 + h32=2 + h32=2 juj = 78
If we have no further a-priori knowledge about u, the most ecient approximation will be one with h1 = h2 = h3 because this equalises the three main terms in the error bound. We see that X Rn = Qj ; j n and the truncation error for u ? Rnu is neither particularly promising nor surprising: the major part of the error is produced by the largest meshwidth: (max(h1 ; h2 ; h3 ))3=2 , whereas the total number of degrees of freedom for an element in Vn is 2jnj . Following the idea of sparse grids, as introduced for nite elements in [3, 9], a better accuracy per degree of freedom is obtained for the approximation operator X R^m = Qj (25) jj jm with m2ZZ . 8
Theorem 2.3 For the approximation operator (25) we have the truncation error estimate ku ? R^muk < M (")2?(3?")m=2 juj = M (") (h1 h2 h3 )(3?")=2 juj : (26)
for some arbitrary small constant " and a constant M ("), depending on ". Proof: Following the same lines as in Theorem 2.2, and because kuj k2 2?3jj j juj, we get ku ? R^muk2 = Pjj j>m kuj k2 juj2 Pjj j>m 2?3jj j = juj2 Pl>m Pjj j=l 2?3jj j P (l+1)(l+2) 2?"l2?(3?")l (27) = juj2 l>m P 2 ?(3?")l 2 2 juj M (") l>m 2 = juj2 M (")2 2?m(3?") Hence ku ? R^m uk M (")2?m(3?")=2juj (28) ? m where 2 = h1 h2 h3 is the volume of the smallest cells in the sparse grid used for the approximation of u. 2 In (22) all hj need to be small and in (26) only their product. This means that for convergence in (22) all meshsizes should tend to zero, whereas in (26) only the area should vanish. Further, the estimate (26) is of the same order of accuracy as (22), except for a logarithmic small factor. However, the number of degrees of freedom for the approximation (26) is signi cantly less. Namely, in the unit cube, for Rn u the number of degrees of freedom is 2jnj , whereas for R^m u it is (m2 + m +2)2m ? 1. Because signi cantly less degrees of freedom are involved in the approximation R^m u than in the approximation of R(m;m;m) u, i.e. less coecients aj ;k and less gridpoints zj ;k , following [3], we call the approximation R^ m u the sparse grid approximation and
n
o
?m = zj ;k j zj ;k 2 ?j ; jjj m is the sparse grid or sparse box grid for this approximation on level m. In this paper we are interested in the approximate solution of PDEs discretised on a sparse grid, i.e. we are looking for an approximation of the solution of these equations in the space M M Sm (IR3) = Vn = Wn ; or, for = (0; 1)3 , in the space
Sm ( ) =
jnjm
M 0jnjm
jnjm
Vn = V0
M 0jnjm
Wn :
We call Sm ( ) the m-th level sparse-grid space.
3 The multigrid algorithm The basis principle of multigrid for the solution of the discrete equation Lh uh = fh is that the high frequencies in the error are reduced by relaxation on a ne grid, whereas the low frequencies are taken care of by coarse-grid discrete equations. The classical coarse grid correction (CGC) step is described by u(new) = u(old) + PhH L?H1 RHh (fh ? Lhu(old) ) ; 9
where LH is the coarse-grid discrete operator and PhH and RHh are the grid-transfer operators from the coarse-to-the- ne and ne-to-the-coarse grid respectively. Usually the coarse-grid mesh size is twice the mesh size on the next ner grid. The coarse grid problem is approximately solved by means of the recursive application of the same algorithm on the coarser level. In this classical procedure a linearly ordered sequence of ne and coarse discretisations is required. In the case of our sparse-grid nite-volume approximation, a discretisation should exist for all grids n , jnj m, ne and coarse. On each of these grids we can obtain approximations to un 2Vn, the solution of the discrete problem
Ln un = fn on n : (29) These discretisations, however, don't o er an ordered sequence. Nevertheless, the multidimensional wavelet decomposition of un2Vn, un = v0 +
X
wj ; with wj 2Wj ;
j n allows us to distinguish a high-frequency component, wn , that cannot be represented on coarser grids, and all other components, v0 and wj , j n, j 6= n, which can be present in coarser grid representations. Therefore we may consider the grid n to be solely responsible for the accurate (and ecient) representation of wn . This component is clearly a high-frequency function (in fact a checkerboard-type function), of which an error can be eciently reduced by a simple relaxation procedure as e.g. damped Jacobi. The decomposition (13) in Theorem 2.1 shows us how a CGC can be obtained from these coarser grids in n?e , j = 1; 2; 3, j
P
u(new) n + Pj=1;2;3 Pn;n?e L?n1??e1Rn?e ;n rn n = u(old) ? j=1;2;3 Pn;n?e+e Ln?e+e Rn?e+e ;n rn + Pn;n?e L?n1?e Rn?e;n rn ; j
j
j
j
with
j
j
(30)
(31) rn = fn ? Lnu(old) n : Here we denote by Rm;n : Vn ! Vm, m n, the restriction operator de ned by Rm;nun = Rmun for all un2Vn L2( ). The prolongation operator Pn;m : Vm ! Vn can be de ned e.g. as the adjoint of Rm;n. The third remark following Theorem 2.1 shows how the two- or one-dimensional case can be treated similarly and we see that {for the one-dimensional case{ our approach reduces to the classical scheme. The approach (30) would imply three coarser levels to be active for a CGC, and {as was shown in the remark after Theorem 2.1{ we can do with only one coarser level by deriving the information on the levels jnj ? 3 and jnj ? 2 from the information on level jnj ? 1. If we consider the corrections from level jnj ? 1, (32) cn?e = L?n1?e Rn?e ;n rn; j = 1; 2; 3; j
j
j
as approximating a single (but unknown) correction function cn 2Vn, the corrections from the levels jnj ? 2 and jnj ? 3 can be computed as the mean values ? cn?e+e = 21 Rn?e+e ;n?e +1 cn?e +1 + Rn?e+e ;n?e ?1 cn?e ?1 ; (33) j = 1; 2; 3, and X (34) Rn?e;n?e+e cn?e+e : cn?e = 31 j
j
j
j
j
j
j =1;2;3
10
j
j
j
This is justi ed by the fact that the restrictions are commutative, i.e.
m n l ) Rm;nRn;l = Rm;l and the following (trivial) lemma.
Lemma 3.1 If all discrete operators Ln are stable and relatively consistent, i.e. kRn;n+e Ln+e ? Ln Rn;n+e k O(2?jnjp) then
j
j
j
kcn?e ? Rn?e ;n cnk O(2?jnjp) : j
j
The consistent discretisations can be derived e.g. from the ne grid discretisation Ln by taking the Galerkin approximation
Ln?e = Rn?e ;nLn Pn;n?e : If the three corrections cn?e were all restrictions of the (unknown) correction cn 2Vn indeed, then Rn?e+e ;n?e +1 cn?e +1 and Rn?e+e ;n?e ?1 cn?e ?1 would both have delivered the same result, viz., Rn?e+e ;n cn. This gives us the possibility to check how well such a function cn can be determined, by monitoring the quantities, j = 1; 2; 3, ? dn?e+e = 21 Rn?e+e ;n?e +1 cn?e +1 ? Rn?e+e ;n?e ?1 cn?e ?1 : (35) Summarising, our multigrid algorithm now reads: j
j
j
j
j
j
j
j
j
j
j
j
j
j
j
j
P
j
j
u(ni+1) = u(ni) + Pj=1;2;3 Pn;n?e cn?e (36) ? j=1;2;3 Pn;n?e+e cn?e+e + Pn;n?e cn;n?e ; where the corrections are given by (32), (33) and (34). This appears to be much similar to an multigrid algorithm for semi-coarsening, proposed by Mulder in [7]. The main di erence being that Mulder computes an approximation on the full grid Rn , whereas we compute the sparse grid approximation R^ m . The result of our algorithm is a solution on a sparse grid, i.e. a set of approximate solutions, viz. fun j jnj = mg, that are the solutions of the discrete equations Ln un = fn . All approximations un representing the same solution u of the continuous problem, we assume that they approximate the L2 ( )-projection of u in Vm = V(m;m;m) . To approximate this Rmu2Vm, we can construct a unique function um2Vm by means of the recursive interpolation formula that immediately follows from Theorem 2.1: j
j
j
uk =
X
j =1;2;3
uk?e ? uk?e+e + uk?e ; j
j
j
(37)
where 0 k m, jkj = m +1; :::; 3m; uk?e are the functions computed in the previous recursion cycle and uk?e+e and uk?e are approximations (possibly) derived as (33) and (34). In this way we nally obtain the unique representation j
j
um = v0 +
X
0jkjm
of
wk
um2Sm( ) V(m;m;m)( ) : This representation can be considered as the computed solution. 11
(38)
The same construction can be realised by a `decomposition and reconstruction' algorithm as used in wavelet theory [2]. Then the available approximate solutions fun g are decomposed into their components v0 and fwk g by
v0 =
P
jP nj=m R0 un jnj=m 1
and
P
nj=m Qk un wk = jP 1 jnj=m
and the reconstruction is performed by (38). This can conveniently be performed by a kind of a `pyramid algorithm'. This will be reported in a later paper. In practice, by the choice of the discrete operator our assumption that the L2 ( )projection of u was indeed consistently approximated on Vn , may not necessarily hold, and it can be checked by a monitor as (35). Now, e.g., the corresponding erroneous components might be removed from the sum (38). In the light of the treatment in Section 2 it is clear what restrictions and prolongations can be used between the di erent grids in the multigrid process. Because Vn L2 ( ), the obvious restriction Rn?e ;n is the L2 ( )-projection onto Vn?e : j
j
Rn?e ;n = Rn?e : j
j
This makes the diagram for the restrictions commutative: for any l m n we have Rl;mRm;n = Rl;n . An obvious prolongation can be the transposed restriction
Pn;n?e = RnT ?e ;n : j
j
However, this prolongation being of low order, it may be more appropriate to consider higher order prolongations. Such prolongations can always be represented by an additional operator Bn : Vn ! Vn so that we have
Pn;n?e = Bn RnT ?e : j
j
Here we will not elaborate on the di erent possibilities for Bn . The algorithm (36) shows that all relaxation processes for un on one and the same level m = jnj can be made in parallel. The cycling between the di erent (scalar) levels can be made as for the classical multigrid method: we can distinguish between V-, Wor F-cycles. However, in order to prove that the convergence of our multigrid-method is independent of the meshwidth, we now have to take into account that all aspect ratios will appear in the discretisations used.
4 Fourier convergence analysis In this section we rst summarise some results of Fourier analysis for more-dimensional discrete approximations and then we apply this to compute the convergence rate of our sparse-grid multiple-grid methods for the solution of the anisotropic Poisson equation. The approach is di erent from the usual treatment of Fourier analysis for multigrid for nite di erence methods for the following reasons. First, in view of the discretisation of conservation laws and divergence problems, we study nested box grids. This implies that mesh points in the coarse grids do not appear in the ne grids as well. The nesting of the (box) grids is di erent from the usual nesting of the ( nite di erence or nite element) grids. Second, we do not consider the usual sequences of ne and coarser meshes for multigrid methods, but we allow d-dimensional (d = 2; 3) sparse grids. Fourier analysis is one of the common tools to analyse linear constant coecient problems on regular grids, and it is particularly useful if the treatment of boundary conditions can be neglected. 12
In Section 4.1 we describe general tools that can be used for the Fourier analysis of functions de ned on more-dimensional box grids. The de nitions and theorems provide a useful machinery for the application of local mode analysis for the multigrid box-methods. For the technical proof and details related to this section we refer to [4]. In Section 4.2 we apply tools to analyse the multigrid algorithm introduced in section 3. The technical preparations in 4.1 allow us to be brief and clear in this treatment.
4.1 Fourier analysis for sparse box grids
For u2L2 (IR3 ) we introduce its Fourier Transform (FT) ub, scaled as Z e?ix! u(x) dx : ub(!) = (2)?3=2
(39)
Then we know that ub2L2 (IR3 ), and a back-transformation formula is available, Z e+ix! ub(!) d! ; u(x) = (2)?3=2
(40)
IR3
IR3
such that u(x) = u(x) almost everywhere on IR3 . Moreover, ub2L2 (IR3 ) and Parseval's equality holds: kukL2 (IR3) = k ubkL2(IR3) . We are interested in the Fourier transformation for an in nite set of equally spaced data. In this case the FT of such a \grid function" is a periodic function (a function de ned on a torus). Therefore we introduce a few de nitions. Let h2IR3 , h > 0, be given, then the h-periodisation of a function u : IR3 ! CI is de ned by X u(x ? kh); (41) u~(x) = 3 k2ZZ where kh = (k1h1 ; k2 h2 ; k3h3 ). We also introduce a notation for the three-dimensional torus Th = (?=h; =h] = (?=h1; =h1] ::: (?=h3 ; =h3]: (42) 3 Further we need the functions and Sinc , [1, pp.62,67] on IR , for jx j < 1=2 ; 1 i 3; i (x) = 10 otherwise (43) ; and iY =3 Sinc x = sinxxi : i=1
i
Using the relations mentioned in [1, p.98] we nd Y sinc ( 2!i ) = (2)?3=2 Sinc ( 2! ) : b (!) = (2)?3=2 i=1;2;3
(44)
For an h2IR3 , h > 0, we de ne the dilation operator Dh : L2 (IR3 ) ! L2 (IR3 ) by Dhf (x) = h?3=2 f (xh) ; (45) where h = (h1 h2 h3 )1=3 , and the convolution operator, ?, by (f ? g)(x) = (2)?3=2 We now know that Grid functions.
Z
IR3
f (y) g(x ? y) dy :
dh = D1=h fb and fd D ? g (!) = fb(!) bg(!) :
Here we introduce notations for the di erent types of grids and gridfunctions. 13
(46) (47)
De nition 4.1 For a xed mesh parameter h2IR3, h > 0, and for j 2ZZ 3 , we de ne an elementary cell by h;j = fx j jh < x < (j + e)hg, the volume of the cell is denoted by h3 = h1 h2 h3, and the box-grid is de ned by h = f h;j j j 2ZZ 3g : The regular in nite three-dimensional grid of vertices ZZh is de ned by ZZh = fjh j j 2ZZ 3g ; which should be well distinguished from the shifted grid which is de ned by ZZh? = f(j + e=2)h j j 2ZZ 3 g : Notice the relation with the grids as de ned in Section 2.2: n can be considered as a special case of h , and ?n as a special case of ZZh? . De nition 4.2 A complex or a real grid function uh is a mapping ZZh ! CI, or ZZh ! IR, and a shifted or box-grid function u?h is a mapping ZZh? ! CI or ZZh? ! IR. The vector space of such gridfunctions we denote by l(ZZh ) or l(ZZh? ), or brie y, by l. For any p 1 the space l(ZZh ) can be provided with a norm k kp
kuh kp = (h3
X
j 2ZZ 3
juh(jh)jp)1=p :
(48)
For a xed p, all grid functions with a nite norm k kp form a Banach space denoted by lp (ZZh ). For p = 2 we know that l2 (ZZh ) is a Hilbert space with the inner product,
< uh; vh >l2(ZZh) = h3
X
j 2ZZ 3
uh(jh)vh (jh) with uh; vh 2ZZh :
(49)
Similar de nitions are given for l(ZZh? ):
X ? ? < u?h ; vh? >l2(ZZ )= h3 uh (z)vh (z) with u?h ; vh? 2ZZh? : (50) h z2ZZh De nition 4.3 The direct restriction Rh : L2(IR3 ) ! l(ZZh ) is the operator that associates with a continuous7 u2L2 (IR3 ) the corresponding grid function on the grid ZZh , de ned by (Rh u)(jh) = u(jh) ; 8j2ZZ 3 ; (51) and the direct restriction Rh? : L2 (IR3 ) ! l(ZZh? ) on the shifted grid ZZh? is de ned by (Rh? u)((j + e=2)h) = u((j + e=2)h) ; 8j2ZZ 3 : (52) De nition 4.4 The box-restriction Rh : L2(IR3 ) ! L2(IR3 ) is the L2-projection on the piecewise constant functions on h , de ned by (cf. equation (6)) ?
?
(Rh u)(x) = h?3
Z
h j
u(z) dz 8x2 h;j :
(53)
;
The box-restriction RhB : L2 (IR3) ! l(ZZh? ) is de ned by RhB = Rh? Rh ; it associates the mean value of u on a cell h;j with the nodal value at the centre of h;j .
The box-restriction RhB u should be well distinguished from Rh? u. However, the L2 ( )projection Rh u in (53) and the restriction RhB u in (52) are conveniently related to each other by 2 3=2 Rh? ((Dh ) ? u) : (54) Rh u = h The Fourier transform of a grid function. Let u : ZZ ! CI be a grid function. We give the following
h
7
h
In case of a discontinuous function we can replace u by u as de ned in (40).
14
De nition 4.5 The Fourier transform uch 2L2(Th ) of a grid function uh 2l2 (ZZh) is a function Th ! CI, de ned by
h 3 X c e?ijh! uh (jh) : uh (!) = p 2 j 2ZZ 3
(55)
The inverse transformation is given by
u (jh) = h
1 3 Z p 2
!2Th
e+ijh! uch (!) d! :
(56)
Let u?h : ZZh? ! CI be a shifted grid function, then we have
De nition 4.6 The Fourier transform uc?h 2L2(Th ) of a shifted grid function u?h 2l2 (ZZh? ) is a function Th ! CI, de ned by
h 3 X ? c e?i(j +e=2)h! u?h((j + e=2)h) : uh (!) = p 2 j 2ZZ 3
(57)
Its inverse transformation is given by
u? ((j + e=2)h) = h
1 3 Z p 2
! 2 Th
e+i(j +e=2)h! uc?h(!) d! :
(58)
Remarks: We immediately see that uch (!) is [2=h]-periodic, whereas uc?h (!) is [2=h]antiperiodic, i.e. uc?h (! + 2=h) = (?)jej uc?h (!). We denote the Fourier transforms also by uch = F (uh) or uc?h = F (u?h); (59) i.e. we introduce the mapping F : l2 (ZZh ) ! L2 (Th) or F : l2 (ZZh? ) ! L2(Th ). At the end of this section we shall generalise this meaning of F . By the Parseval equality we have kuh k2 = k uch kL2(Th ) and ku?h k2 = k uc?h kL2(Th ) : (60)
Hence the Fourier transformation operators F : l2 (ZZh ) ! L2 (Th ) and F : l2 (ZZh? ) ! L2(Th ) are unitary operators. Because eh;! eh;!+2k=h or e?h;! (?)jekje?h;!+2k=h , for all k2ZZ 3, on a mesh of size h a frequency ! cannot be distinguished from a frequency ! + 2k=h. This phenomenon is called aliasing . The relation between FTs of a function restricted to different grids. In this section we rst present a few theorems associated with the di erent restrictions between two grids. We describe the relation between the FT of a continuous function de ned on IR3 and the FT of its restriction to the grid and then we show the relation between the FT of a ne grid function and the FT of its representation on a coarser grid. Next, we give the corresponding theorems for the prolongations. 15
Lemma 4.7 Let u2L2 (IR3) be a continuous function with FT ub. Its restriction uh to
the grid ZZh is de ned by (51). We have the following relation between ub and uch: X ub(! + 2k=h); (61) uch(!) = k2ZZ 3 i.e. uch is the [2=h]-periodisation of ub. Proof: For the proof we refer to [4]. 2 In the following lemmas q-restrictions are considered, with q2ZZ 3. Here q = (q1 ; q2 ; q3 ) is the coarsening factor, where usually 1 qj 2, j = 1; 2; 3. De nition 4.8 Let q2ZZ 3 with q > 0 and H = qh, then the canonical q-restriction Rq is the operator Rq : l(ZZh ) ! l(ZZqh ) = l(ZZH ), de ned by (Rq uh )(jH) = uH (jH) = uh (j qh) : (62) Theorem 4.9 We have the following relation between the FT of a grid function and that of its canonical q-restriction,
d X c Rq uh (!) = uh (! + 2p=h); 8!2TH ; H = qh: p2[0;q)
(63)
Proof: For the proof we refer to [4]. 2
Lemma 4.9 shows that, using the restriction Rq with q2ZZ 3 , q > 0, we get aliasing of q3 = q1:q2 :q3 frequencies onto one. Now we describe the relation between the Fourier transforms of a continuous function and its box restrictions. First we consider the direct restriction to the shifted grid, Rh? , and later the box-restriction, Rh and the q-restriction Rq? .
Lemma 4.10
? uc?h (!) = Rd h u(!) =
Proof: For the proof we refer to [4]. 2
X
k2ZZ 3
(?)jkj ub(! + 2k=h) :
(64)
For the Fourier transform of u2L2 (IR3 ) and Rh u2L2 (IR3 ) we have the following relation.
Theorem 4.11
uch (!) = Rd h u(!) =
X k2ZZ 3
(?)jkj Sinc ( h! 2 + k) ub(! + 2k=h) :
Proof: Using (54) we see uc?h (!) = F (Rh u)(!) ? = 2h 3=2 F (Rh? ((Dh ) ? u))(!) ? P = ? 2h 3=2 k (?)jkj F ((Dh ) ? u)(! + 2k=h) P = ? 2h 3=2 k (?)jkj (Dd h )(! + 2k=h) ub(! + 2k=h) 3=2 P 2 j k j = h k (?) D1=h b (! + 2k=h) ub(! + 2k=h) : P = h?3=2 k (?)jkj D1=h Sinc ( 2! + k=h) ub(! + 2k=h) : P = k (?)jkj Sinc ( h! 2 + k) ub(! + 2 k=h) : 2
16
(65)
(66)
De nition 4.12 Let q2ZZ 3 with q > 0 and H = qh, then, for s2[0; q), the s-frequency decomposition q-restriction is the operator Rqs : l(ZZh ) ! l(ZZqh ) = l(ZZH ) is de ned
by
(Rqs u?h)((j + e=2)H) = q ?3
X
(?)sk u?h ((qj + k + e=2)h) ;
k2[0;q)
(67)
where q 3 denotes q 3 = q1 q2 q3 .
Remarks: In the case s = 0 we call Rqs = Rq0 = Rq? simply the q-restriction. From the construction of the restriction operators RhB and Rq? it is clear that the following relation holds: B = R? RB : RH q h 3 , H = 2h, we have Theorem 4.13 Let q = 2e 2ZZ 3, then for all !2TH X m h! + (m + s) c? ( Rd (?) Cos( ) uh (! + m=h) (68) qs u?h )(!) = 2 q?3 is 2 m2[0;q) where
Cos(h!=2) =
Y
j =1;2;3
cos(hj !j =2) :
Proof: For the proof we refer to [4]. 2 De nition 4.14 The natural box-prolongation Ph? : l2 (ZZh? ) ! L2(IR3 ) is de ned by u(x) = (Ph? u?h )(x) = uh ((m + e=2)h); for all m2ZZ 3 and mh < x < (m + e)h. We also introduce a natural prolongation, Pq? , ? ) ! l2 (ZZ ? ) by P ? = R? P ? where from a coarse to a ne gridfunction Pq? : l2 (ZZH q h H h H = qh. The prolongation PhB : l(ZZh? ) ! L2 (IR3 ) is the operator dual to RhB in the sense that for all u?h 2l2 (ZZh? ) and v2L2 (IR3 ) we have
< PhB u?h ; v >L2(IR3) = < u?h; RhB v >l2 (ZZ
: (69) h) The following theorems show how we nd the FT of the prolongation if the FT of the source function is given.
Lemma 4.15
?
=2) c? ub(!) = F (Ph? u?h )(!) = Sin(h! h!=2 uh (!) :
(70)
Proof: For the proof we refer to [4]. 2 Theorem 4.16 With H = qh and q3 = q1 :q2 :q3, we have for the FT of the prolongation of a box gridfunction
=2) ud ? uc?h (!) = F (Pq? u?H )(!) = q?3 Sin(H! Sin(h!=2) H (!) :
Proof: For the proof we refer to [4]. 2 17
(71)
The Fourier transform of operators involving different grids. In (68) and (71) we see that, by the restriction and prolongation between functions on grids n and n+q , aliasing takes place and that q 3 frequencies on n+q correspond with a single frequency on n. This implies that, analysing a multigrid algorithm, we have to study the behaviour of the q 3 aliasing frequencies together. Collecting the q 3 corresponding amplitudes of the aliasing frequencies in a single q 3-vector, we extend the de nition (59) of F to the case q 3 > 1 and obtain F : l2 (ZZh ) ! [L2 (Tqh )]q3 or F : l2 (ZZh? ) ! [L2(Tqh)]q3 by ? F (uh)(!) = ubh (! + m=h) m2[0;q) ; !2Tqh : (72)
With these amplitude-vectors F (uh )(!) and F (u?h )(!) , we can introduce the linear operators F (Rq )(!) and F (Rq? )(!) by
and
F (Rq uh )(!) = F (Rq )(!) F (uh )(!) ; !2Tqh ;
(73)
F (Rq? u?h )(!) = F (Rq? )(!) F (u?h )(!) ; !2Tqh :
(74) We call the new operators, that depend on !, the Fourier transforms of the original operators. The new operators are ` q 3 ` matrices if ` aliasing frequencies are considered on the coarse grid. Similar to the restrictions, the prolongations can be associated with their Fourier transforms. (75) F (Pq uh )(!) = F (Pq )(!) F (uh)(!) ; !2Tqh ; and F (Pq? u?h )(!) = F (Pq? )(!) F (u?h)(!) ; !2Tqh : (76) These operators are q 3 ` ` matrices. For arbitrary linear constant coecient operators Ah : l2 (ZZh ) ! l2 (ZZh ), its Fourier transform F (Ah) : L2(Th ) ! L2 (Th ), can also be considered as a q 3 ` q 3 ` diagonal matrix F (Ahuh )(!) = F (Ah)(!) F (uh )(!) 8!2Tqh : Because of Parseval's equality we know that
( F (Ah) ) (!) k F (Ah)(!)k2 = !max kAhk2 = !max 2Tqh 2Tqh
(77)
with (A) the spectral norm (the largest singular value) of the matrix A, and
( F (Ah) ) (!) (78) (Ah ) = !max 2Tqh where denotes the spectral radius. This provides us with the means to study the norm and the spectral radius of the error-ampli cation operator of the multigrid iteration.
4.2 Fourier analysis convergence results
To get some insight in the behaviour of the more-dimensional multigrid algorithm introduced in Section 3, we use Fourier analysis to determine the convergence rate of the two-level algorithm for the two-dimensional anisotropic Poisson equation
uxx + "2 uyy = f ; discretised by the usual 5-point discretisation. 18
(79)
The error-ampli cation operator, Mn, of the two-level cycle (with pre- and postrelaxation steps) for the solution of (29) is described by
e(ni+1) = Mn e(ni) = Sn Cn Sn e(ni) ; where Sn denotes the smoothing, e.g. damped Jacobi iteration:
?
(80)
e(new) (81) n = Sn e(old) n = In ? Dn?1 Ln e(old) n ; with the damping parameter and Dn the main diagonal of the discrete operator Ln . The coarse grid correction is described by (30) and (31). For gridfunctions uh 2l2 (ZZh? ), i.e. neglecting boundary conditions, we nd, using (30) F (Mn )(!) = F (Sn ) (!) F (Cn )(!) F (Sn) (!) ; (82) with F (Sn ) = F (In) ? F (Dn)?1 F (Ln) ; (83) and F (Cn) = FX(In) (84) ? F (Pn;n?e ) F (Ln?e )?1 F (Rn?e ;n) F (Ln) +
j =1;2;3
X
j =1;2;3
j
j
j
F (Pn;n?e+e ) F (Ln?e+e )?1 F (Rn?e+e ;n) F (Ln) j
j
j
? F (Pn;n?e) F (Ln?e)?1 F (Rn?e;n) F (Ln) : To get an impression of the behaviour of the algorithm, keeping the explicit computation to reasonable size, we elaborate (82) for the equation (79), for the two-dimensional case, with q = (2; 2) and = = 1. Then F (Mn )(!) is a 4 4-matrix, which we derive from (82), (83) and (cf. (15)) F (Cn ) = F (In) ? F (Pn;n?e1 ) F (Ln?e1 )?1 F (Rn?e1;n ) F (Ln) ? F (Pn;n?e2 ) F (Ln?e2 )?1 F (Rn?e2;n ) F (Ln) + F (Pn;n?e) F (Ln?e )?1 F (Rn?e;n ) F (Ln) : From (68) and (71) we know F (Rn?e1;n ) = cos !0 2 h2 cos !0 2h2 sin !02 h2 sin !02h2 ;
cos ! h sin ! h 1 1 1 1 F (Rn?e2 ;n) =
and
0
0
0 0 cos !1 h1 sin !1 h1
0 cos(! h ) cos(! h ) 1T 1 1 2 2 B C : sin( ! h ) cos( ! 1 1 2 B F (Rn?e;n ) = @ cos(! h ) sin(! hh2 )) C A 1 1
2 2
sin(!1 h1 ) sin(!2 h2 ) So that with F (Pn;n?e1 ) = F (Rn?e1 ;n)T , F (Pn;n?e2 ) = F (Rn?e2 ;n)T and F (Pn;n?e) = F (Rn?e;n)T , the norm kMnk and the spectral radius (Mn ) can be computed by means of (77) and (78). To study the convergence behaviour of our algorithm, we consider the matrices (83), (84) and (82) as a function of ! 2Th = [?=h; =h]2, and for each ! we compute its 19
1
1
0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0
-0.2
-0.2
-0.4 0.8
0.8
0.6
0.6
0.4
0.4 0.2 0.2
0.6
0.4
0.8
0.2
(a) The eigenvalues of F (Sn)(!). 0 0
0.2
0.8
0.6
0.4
(b) The eigenvalues of F (Cn )(!). 0 0
0.1 0.08 0.3 0.06 0.25 0.04 0.2
0.02
0.15
0 -0.02
0.1
-0.04
0.05
0.8
0 0.8 0.6 0.6 0.4
0.8 0.6
0.4 0.2 0.2
0.4
0.6
(c) The eigenvalues of F (Mn)(!).
0 0
0.4
0.2
0.8
0.2
(d) The singular values of F (Mn )(!). 0
0
Figure 1: The frequency response of the operators Sn , Cn and Mn , for " = 1, q = (2; 2) and = 2=3. eigenvalues and singular values of these matrices. We show these values in Fig. 1 for the case = 2=3, " = 1. Without loss of generality we can take h = (1; 1), the parameter " taking care of the anisotropy. The damping parameter 2[0; 1] for the Jacobi relaxation, can be chosen freely. We select the value = 2=3 because it minimises
( F (Sn)(!)) : !=(0;=h)max ;( ;0);(=h;=h) This means that = 2=3 makes Sn a well balanced smoother for the di erent types of high frequencies (see Fig. 1.a). In Fig. 1.a we show the eigenvalues of the smoothing operator, and in Fig. 1.b of the coarse grid correction. In this gure we see that one eigenvalue of F (Cn ) is always equal to one. This eigenvalue corresponds with the highest frequencies, for which no correction can be obtained from any of the three coarser grids. The combined e ect of the smoother and the coarse grid correction is seen in Fig. 1.c, which shows that sup! (Mn (!)) 1=9, and also in Fig. 1.d, where we see sup! kMn (!)k 1=3. The rather low maximal values show that {at least for square ne-grid cells{ the multigrid algorithm has good convergence behaviour. Because it is important that the algorithm is e ective for an arbitrary cell aspect ratio, in the Figs 2 - 4 we show the singular values of F (Mn)(!) also for " = 1=8 and for the limit as " ! 0. Now it appears that for high aspect ratios the convergence behaviour deteriorates. We nd sup! lim"!0 ( F (Mn )(!)) 5. This implies that we cannot always guarantee error reduction for a single iteration sweep. Therefore we show in Fig. 4 also the behaviour of Mn2 (!). This shows that a couple of two consecutive 20
1.6 1.4 1.2
0.2
1
0.1
0.8
0
0.6
-0.1
0.4
0.8
-0.2 0.8
0.2 0.6 0 0.8
0.6 0.4
0.8 0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
(a) The eigenvalues of F (Mn)(!). 0
0.2
(b) The singular values of F (Mn)(!).
0
0
0
Figure 2: The frequency response of the operator Mn , for " = 1=8, q = (2; 2) and = 2=3. iteration steps does guarantee error reduction, and the convergence rate is signi cant: sup "lim ( F (Mn2 )(!)) 1=9 : ! !0 As a consequence we can expect that a W-type cycle of the multigrid algorithm will have good convergence properties.
From the computations of which the results are summarised in the Figs 2 - 4, we conclude that the eigenvalues of the iteration matrix are less that 1, uniformly in the parameter ". In fact, max" (Mn ) 0:33 and max" kMn k 5:0 and max" kMn2 k 0:11. The fact that kMn k > 1 and kMn2 k
1 Introduction In this paper we describe the approximation of a function on a nite-volume sparse grid, and a multigrid algorithm for the solution of partial di erential equations in three dimensions. The algorithm is intended for the solution of ow problems described by conservation laws, and therefore nite volumes are a natural choice for the discretisation. But to introduce the main principles, we will restrict the treatment here to second order elliptic equations, and in particular to the anisotropic Poisson equation. In contrast to the usual multigrid approach, we do not use a sequentially ordered set of discretisations on di erent meshes, but we use a partially ordered hierarchy of `semicoarsened' grids as proposed e.g. by Mulder [6, 7] and Naik-VanRosendale [8] or Zenger et al. [3, 9]. As indicated in [9], adaptive `sparse grid' discretisations can be constructed by taking a suitable subset of all possible discretisations in such a hierarchy. However, in contrast to the sparse grid approximation proposed in [3, 9], we base our approximation on nite volumes rather than on nite elements. The multigrid algorithm consists of damped Jacobi relaxation as a smoothing procedure and a coarse grid correction constructed by extrapolation from simultaneous corrections on several coarser grids levels. The algorithm is completely symmetric with respect to the three coordinate directions and it is suitable for combination with adaptive techniques. A description of a data structure to implement such adaptive three-dimensional algorithms is given elsewhere [5]. 1
2 Finite volume sparse grids In this section we introduce nite-volume sparse grids. We show the relation between the approximation by Haar wavelets (when this notion is extended to more dimensions) and the sparse-grid approximation. For the theory of wavelets, multiresolution analysis (MRA) etc. we refer to Daubechies [2].
2.1 The more-dimensional MRA A multidimensional multiresolution analysis of L2 ( ), = IR3 , is a partially ordered set of closed linear subspaces fVnjVn L2( )gn2ZZ3 with the properties:
T
S
(1) n Vn = f0g; n Vn dense L2 ( );3 m (2) f (x)2Vn , f (2 x)2Vn+m 8n2ZZ ; m2E; (3) f (x)2Vn , f (x ? 2?n k)2Vn 8k2ZZ 3 ; n2E; (4) 9 2V0 : f (x ? k)gk2ZZ 3 is a Riesz basis for V0 :
(1)
Here n = (n1 ; n2 ; n3)2ZZ 3 , and we denote jnj = n1 +n2 +n3 ; 2n = (2n1 ; 2n2 ; 2n3 ). We also use the notation 0 = (0; 0; 0)2IN 3 ; x = (x1 ; x2 ; x3 )2IR3 ; 2nx = (2n1 x1 ; 2n2 x2 ; 2n3 x3 ). Further we introduce in IN 3 the unit vectors ek , k = 1; 2; 3; as follows: e1 = (1; 0; 0); e2 = (0; 1; 0); e3 = (0; 0; 1), and we use e = (1; 1; 1). Finally we de ne E = fe1 ; e2; e3g. Although we are particularly interested in the three-dimensional case, generalisation to a di erent number of space dimensions is straightforward. The function (x) in (1.4) is called the the scaling function of the multiresolution analysis.
2.2 Piecewise constant function spaces Let either = IR3 be the three-dimensional Euclidean space, or let = (0; 1)3 IR3 be the open unit cube. For any n2ZZ 3 we introduce the function space Vn , i.e. the space of piecewise constant functions on a uniform grid with meshsize h = (2?n1 ; 2?n2 ; 2?n3 ). These grids are uniformly spaced in each of the three coordinate directions, but possibly with a di erent meshsize in the di erent directions. The volume of these cells is denoted by h3 = 2?jnj . The functions in Vn are all constant in each cell
n;k = [k12?n1 ; (k1 + 1)2?n1 ] [k22?n2 ; (k2 + 1)2?n2 ] [k3 2?n3 ; (k3 + 1)2?n3 ] ; and this family of cells forms the grid n = f n;k j n;k ; k2ZZ 3 g. The family of cell centres or cell nodes is denoted by n = fzn;k j zn;k = (k+ e=2)2?n ; k2ZZ 3 ; zn;k 2 g. The number of these nodes is equal to the dimension of Vn. Apparently, all grids are identi ed by a triple n; the number jnj is called the level of the grid n. Notice that {di erent from classical multigrid{ here and later in our multigrid algorithm, there is a clear distinction between the grid-identi cation index n and the level number jnj. Because for j = 1; 2; 3; (2) Vn Vn+e ; j
we see, by construction, that nesting relations exist between spaces Vn and that the 2
nesting provides a partial ordering1 : Vn Vm , n m: (3) Spaces Vn and Vm or grids n and m that satisfy this nesting relation n < m are called related. The construction of the spaces Vn shows that even a stronger relation holds than (2), namely Vn?e \ Vn?e = Vn?e ?e ; j; k = 1; 2; 3; j 6= k: (4) We see also that for = IR3 the spaces fVngn2ZZ 3 form a multiresolution analysis and that in this case the characteristic function on the unit cube, 2V0 , (x) = 01 ifif xx622
0;0 ;; (5) 0;0 serves as the scaling function. The set f nk j nk (x) = (2n x ? k); k2ZZ 3 g forms a basis in Vn , which corresponds with the usual Haar-basis for the one-dimensional case. In the case = (0; 1)3 we restrict ourselves to Vn with n1 ; n2 ; n3 0 and we see dim(Vn ) = 2jnj . Formally, for = (0; 1)3 and n1 , n2 or n3 negative we de ne Vn = Vn1;n2;n3 by Vn = f0g. For all spaces Vn we introduce the restriction operator Rn : L2 ( ) ! Vn, the L2 projection such that for u2L2 ( ) we have Rn u2Vn and j
j
k
k
(Rn u)(zn;k ) = 2jnj
Z
n k
u(z) dz :
(6)
;
2.3 More-dimensional wavelets
We introduce the wavelet space Wn Vn which consists of all functions in Vn that are not represented in any of the related function spaces on the next coarser level, i.e. they are in Vn but not in Span(Vn?e1 ; Vn?e2 ; Vn?e3 ), or (7) Vn = Wn Span(Vn?e1 ; Vn?e2 ; Vn?e3 ); (8) f0g = Wn \ Span(Vn?e1 ; Vn?e2 ; Vn?e3 ): This means that Wn contains the `di erence information' that is available in the ne grid Vn but not in the the span of the coarser grids Vn?e1 , Vn?e2 and Vn?e3 . In our case, where Vn contains the piecewise constant functions, it is simple to construct the spaces Wn such that (9) Wn ? Span(Vn?e1 ; Vn?e2 ; Vn?e3 ): This makes Wn the orthogonal complement of Span(Vn?(1;0;0) ; Vn?(0;1;0); Vn?(0;0;1)) in Vn , and (8) follows immediately from (9). For IR3 the relation (9) allows a straightforward decomposition of Vn. In the case = IR3 we have2 n1 n2 n3 M M M
Vn (IR3 ) =
j1 =?1 j2 =?1 j3 =?1 1 We de ne the inequalities between triples by
and
n<m,n nm,n
Wj ;
1
< m1 ; n2 < m2 ; n3 < m3 ;
1
m1 ; n2 m2 ; n3 m3 :
(10)
2 Notice that here, in the more-dimensional case, it is convenient to choose an indexing that is di erent from the usual indexing in the well-known one-dimensional case.
3
where all Wj are orthogonal to each other. To handle the bounded domain = (0; 1)3 we introduce the functions Vn0 ( ) Vn( ) which have a zero mean value on , i.e. Vn0 ( ) = fu2Vn( ) j R0 (u) = 0g, and we have3
Vn0 ( ) = and hence
n1 M n2 M n3 M
j1 =0 j2 =0 j3 =0
Vn ( ) = V0 Vn0 = V0
Wj ;
M
(11)
Wj :
(12)
0j n The `di erence information' between the approximations of a function f 2L2 ( ) on two successive levels, Rn f 2Vn on the one hand and Rn?e f 2Vn?e , j = 1; 2; 3; on the other hand, is given by the orthogonal projection Qnf of f onto the orthogonal complement Wn of Vn?e in Vn . This is described in the following theorem. Theorem 2.1 Let the operator Qn : L2 ( ) ! Wn be the orthogonal projection onto Wn , then it follows that Qnu = Rn u ? Rn?e1 u ? Rn?e2 u ? Rn?e3 u + Rn?e2 ?e3 u + Rn?e1 ?e3 u + Rn?e1 ?e2 u ? Rn?e u ; or, equivalently, X ? Qnu = Rnu ? Rn?eu + (13) Rn?e+e u ? Rn?e u : j
j
j
j
j
j =1;2;3
Proof: From (10) or (12) it follows that (possibly neglecting functions in V0 if = (0; 1)3 )
Vn =
so that
Rn u =
and
j n
X
j n
Wj ; Qj u ;
(14)
Rn ? Rn?e1 ? Rn?e2 ? Rn?e3 + Rn?e2 ?e3 + Rn?e1 ?e3 + Rn?e1 ?e2 ? Rn?e =
2
M
X
Qj ?
X
Qj ?
X
Qj ?
X
Qj
j n j n?e3 j n?e2 X j n?e1 X X X + Qj + Qj + Qj ? Qj j n?e2 ?e3 j n?e1 ?e3 j n?e1?e2 X X X Xj n?e = Qj + Qj + Qj + Qj n?e<j X n n?e1 <j n X n?e2<j n X n?e3<j n ? Qj ? Qj ? Qj n?e2?e3<j n n?e1?e3<j n n?e1?e2 <j n = Qn :
3 In the case of a bounded domain, ! = (0; 1)3 , W (0;0;0) ( ) is the zero-function and the functions W(j1 ;j2 ;j3 ) with some zero-index (e.g. W(j1 ;j2 ;0) or W(0;0;j3 ) ) have a shape di erent from those with j > 0.
4
Remarks: In the right-hand-side of equation (13) we recognise the information that can be represented on the levels jnj, jnj ? 3, jnj ? 2, jnj ? 1, respectively. In (13) the information on the level jnj ? 2 and jnj ? 3 can directly be derived from the information on level jnj ? 1, e.g. by Rn?e2 ?e3 u = Rn?e2 ?e3 (Rn?e3 u). Thus, expression (13) describes the di erence information between Rn u and its approximation on the related next coarser grids. Notice that in the two-dimensional case the relation (13) reduces to Qnu = Rnu ? Rn?e1 u ? Rn?e2 u + Rn?e u ; where e = (1; 1), and in the one-dimensional case we have
Qnu = Rn u ? Rn?e u :
(15) (16)
First, in the remainder of this section we restrict ourselves to the case of the unbounded domain = IR3 . The four relations (1.1) to (1.4) imply that also the spaces Wn are scaled versions of one space W0 , f (x)2Wn , f (2?n x)2W0 ; 8n2ZZ 3 ; (17) and, moreover, that they are translation invariant for the discrete translations 2?n ZZ 3 ,
f (x)2W0 , f (x ? k)2W0 ; 8n2ZZ 3 : (18) The relations (7) and (9) make that they are mutually orthogonal spaces, generating all functions of L2 (IR3 ), LWn ?3 WWm for nL6=2 (IRm3;): (19)
n2ZZ n
dense
Summarising, we obtain a nesting between the spaces fVn g and fWn g, 9 9 Vn?e1 > >> >
>
>> >> > = Vn > >> Vn?e3 > >> > >; = Wn Vn+e ; > > Vn+e ?e ?1 > > > > Vn+e ?e +1 > > > > Wn+e ; > Vn?e2 > =
j
j
j
j
j
j
that is essentially more complex than in the classical one-dimensional case, where there is a sequential ordering of the spaces fVn gand fWn g. As soon as we nd a function (x) with the property that (x ? k), k2ZZ 3 , is a basis of We , then by a simple rescaling, we see that (2n x ? k), yields a basis of Wn+e . Such a function is the more-dimensional generalisation Since L2 (IR3 ) is the directo n n+e of a wavelet. n + e sum of these Wn+e , the full collection k (x) j k (x) = (2nx ? k); n; k2ZZ 3 is a basis of L2 (IR3 ). 5
It is easy to check that the more-dimensional wavelet (x)2We , corresponding with the scaling function (x)2V0 , from the previous section, is the three-dimensional checkerboard basis function4 given by (5): 8 = 0 if x62 0;0 ; < (x) = : = +1 if x2 0;0 ; x2 e;k ; jkj even; = ?1 if x2 0;0 ; x2 e;k ; jkj odd: This function is the three-dimensional generalisation of the Haar-wavelet. In wavelet theory the spaces Wn are labelled channels, and the distinct channels are linearly independent. The rst decomposition of an arbitrary function from L2 ( ) P 3 with = IR consists in writing u(x) = n un (x), where each un belongs to the corresponding channel Wn with n2ZZ 3 . Similarly, we can write for functions de ned on = (0; 1)3 the relation (12) and make a similar decomposition in channels. Each subspace Wn+e , n 0, has its natural basis, the standard basis5 n n+e n+e(x) = (2nx ? k)o (x) j k k of functions with a minimal support. The basis function kn+e is a scaled, elementary checkerboard function, that may be characterised either by its support which is a single cell in n or by the centerpoint of this cell, zn;k . For = (0; 1)3 , the exceptions related with the boundary are found in the spaces Wn with a zero index (i.e. n1 n2 n3 = 0). These Wn have basis functions with di erent shapes. They are derived from the corresponding functions for the unbounded case, but T their support is restricted to n?e . Their corresponding nodal points zn?e;k are found on the boundary @ = n , n e, n 6= 0. Taking this into account, both for
= (0; 1)3 and for = IR we may write for each u2L2 ( ) a wavelet expansion X u(x) = an;k (2nx ? k): (20) n;k
2.4 Approximation results
The decompositions (10) or (12) clearly allow the approximation of a suciently smooth function in L2 ( ) by a series with elements in Wj . To obtain an impression of the quality of these expansions we derive some error estimates. As the case where boundaries are present is the more general one, we take = (0; 1)3 . To quantify the error of approximation on , we introduce for u2C 3 ( ) the seminorm6 @ 3 u(x) juj = sup @x @x @x + (21) x2 1 2 3 @ p @ q @ r max sup p;q;r=0;1 x2@ @x1 @x2 @x3 u(x) : Now we derive the following Theorem 2.2 If we consider an expansion of a C 3 ( )-function, u, in piecewise constant functions on the grid n, for an arbitrary n2ZZ 3 , n > 0, and if we write X Rn u = v0 + uj ; 0j n
Notice 2We Ve is a function piecewise constant on e . Notice that in more dimensions we use the indexing kn+e , whereas in the one-dimensional case one usually writes kn . 6 The necessity of the boundary terms in this seminorm is seen immediately if we want to approximate in L2 ( ) smooth functions u2C 3 ( ) that do not satisfy homogeneous boundary conditions. 4
5
6
with v0 2V0 and uj 2Wj , 0 j n, then
kuj kL2( ) 2?3jj j=2 juj ; and we get an estimate for the approximation error
ku ? Rn ukL2( ) 8 7?3=2 (2?3n1 + 2?3n2 + 2?3n3 )1=2 juj 8 7?3=2 (h31=2 + h32=2 + h33=2 ) juj :
(22)
Proof: We take the normalised f ~kj g = f2jj ?ej=2 kj g as a basis in Wj , 0 j n, j= 6 0. Clearly, all these functions are orthogonal to all functions v0 2V0 and mutually they form an orthonormal set (an orthonormal Haar basis) in Wj L2 ( ). We see this as follows
2We ; e kj 2We ; k 2 We ;
) = 0;0 ; e kj ) = 0;k ; k ) = j ?e;k ; or, in other words, kj 2Vj , but kj scales like a basis function in Vj ?e . Hence Z j d = 0 for k 6= m ; 2jj ?ej=2 kj 2jj ?ej=2 m and
Z
2jj ?ej=2 kj 2jj ?ej=2 kj d = 2jj ?ej
Thus, we nd with Now
support( support( support(
Z Z j j jj ?ej d
= 2 k k
Rn u = v0 + uj =
X
X
d = 1:
;
uj ;
0j n
X ajk ~kj = (u; ~kj ) ~kj : k k
ajk = (u; ~kj ) =
Z
u ~kj d =
By Taylor expansion around zj ?e;k , we have
Z
j ?e k
(23)
Z
j ?e k ;
u ~kj d :
u ~kj d 2?2jj jjuj :
j ?e k For j e the point zj ?e;k lies in the interior of and the estimate holds with ;
(24)
@ 3 u(x) juj = sup @x @x @x : x2 1 2 3
For j 6 e, i.e. for kj with a j-component equal to zero, the point zj ?e:k lies on the boundary and the function kj is constant in one direction over the whole domain , and 7
it is of Haar-wavelet type for the non-zero indices (or index). In this situation the same estimate (24) holds with, e.g. if j1 = 0, with
@ 2 u(x) juj = sup @x @x : 2 3
For j = 0 the relation (24) is trivially satis ed. Hence, the estimate (24) holds for j 0 if we use the seminorm (21), and we nd jaj ;k j 2?2jj j juj ; X X kuj k2 = jajk j2 2?4jj jjuj2 = 2?3jj jjuj2 ; k k so that kuj k 2?3jj j=2juj ; and
0 BB BB X X X X 2 2 2 + + ku ? Rn uk = kuj k juj B BB BB j1 > n1 j1 0 j 0 j1 > n1 @ j2 0 j21 > n2 j2 0 or j2 > n2 j3 0
or j3 > n3
And it follows that
j3 0
3 ? 78 2?3n1 + 2?3n2 + 2?3n3 juj2 :
j3 > n 3
1 CC CC CC 2?3jj j CC CA
3=2 2?3n1 =2 + 2?3n2 =2 + 2?3n3=2 juj ku ? Rn uk 78
2
3=2 h31=2 + h32=2 + h32=2 juj = 78
If we have no further a-priori knowledge about u, the most ecient approximation will be one with h1 = h2 = h3 because this equalises the three main terms in the error bound. We see that X Rn = Qj ; j n and the truncation error for u ? Rnu is neither particularly promising nor surprising: the major part of the error is produced by the largest meshwidth: (max(h1 ; h2 ; h3 ))3=2 , whereas the total number of degrees of freedom for an element in Vn is 2jnj . Following the idea of sparse grids, as introduced for nite elements in [3, 9], a better accuracy per degree of freedom is obtained for the approximation operator X R^m = Qj (25) jj jm with m2ZZ . 8
Theorem 2.3 For the approximation operator (25) we have the truncation error estimate ku ? R^muk < M (")2?(3?")m=2 juj = M (") (h1 h2 h3 )(3?")=2 juj : (26)
for some arbitrary small constant " and a constant M ("), depending on ". Proof: Following the same lines as in Theorem 2.2, and because kuj k2 2?3jj j juj, we get ku ? R^muk2 = Pjj j>m kuj k2 juj2 Pjj j>m 2?3jj j = juj2 Pl>m Pjj j=l 2?3jj j P (l+1)(l+2) 2?"l2?(3?")l (27) = juj2 l>m P 2 ?(3?")l 2 2 juj M (") l>m 2 = juj2 M (")2 2?m(3?") Hence ku ? R^m uk M (")2?m(3?")=2juj (28) ? m where 2 = h1 h2 h3 is the volume of the smallest cells in the sparse grid used for the approximation of u. 2 In (22) all hj need to be small and in (26) only their product. This means that for convergence in (22) all meshsizes should tend to zero, whereas in (26) only the area should vanish. Further, the estimate (26) is of the same order of accuracy as (22), except for a logarithmic small factor. However, the number of degrees of freedom for the approximation (26) is signi cantly less. Namely, in the unit cube, for Rn u the number of degrees of freedom is 2jnj , whereas for R^m u it is (m2 + m +2)2m ? 1. Because signi cantly less degrees of freedom are involved in the approximation R^m u than in the approximation of R(m;m;m) u, i.e. less coecients aj ;k and less gridpoints zj ;k , following [3], we call the approximation R^ m u the sparse grid approximation and
n
o
?m = zj ;k j zj ;k 2 ?j ; jjj m is the sparse grid or sparse box grid for this approximation on level m. In this paper we are interested in the approximate solution of PDEs discretised on a sparse grid, i.e. we are looking for an approximation of the solution of these equations in the space M M Sm (IR3) = Vn = Wn ; or, for = (0; 1)3 , in the space
Sm ( ) =
jnjm
M 0jnjm
jnjm
Vn = V0
M 0jnjm
Wn :
We call Sm ( ) the m-th level sparse-grid space.
3 The multigrid algorithm The basis principle of multigrid for the solution of the discrete equation Lh uh = fh is that the high frequencies in the error are reduced by relaxation on a ne grid, whereas the low frequencies are taken care of by coarse-grid discrete equations. The classical coarse grid correction (CGC) step is described by u(new) = u(old) + PhH L?H1 RHh (fh ? Lhu(old) ) ; 9
where LH is the coarse-grid discrete operator and PhH and RHh are the grid-transfer operators from the coarse-to-the- ne and ne-to-the-coarse grid respectively. Usually the coarse-grid mesh size is twice the mesh size on the next ner grid. The coarse grid problem is approximately solved by means of the recursive application of the same algorithm on the coarser level. In this classical procedure a linearly ordered sequence of ne and coarse discretisations is required. In the case of our sparse-grid nite-volume approximation, a discretisation should exist for all grids n , jnj m, ne and coarse. On each of these grids we can obtain approximations to un 2Vn, the solution of the discrete problem
Ln un = fn on n : (29) These discretisations, however, don't o er an ordered sequence. Nevertheless, the multidimensional wavelet decomposition of un2Vn, un = v0 +
X
wj ; with wj 2Wj ;
j n allows us to distinguish a high-frequency component, wn , that cannot be represented on coarser grids, and all other components, v0 and wj , j n, j 6= n, which can be present in coarser grid representations. Therefore we may consider the grid n to be solely responsible for the accurate (and ecient) representation of wn . This component is clearly a high-frequency function (in fact a checkerboard-type function), of which an error can be eciently reduced by a simple relaxation procedure as e.g. damped Jacobi. The decomposition (13) in Theorem 2.1 shows us how a CGC can be obtained from these coarser grids in n?e , j = 1; 2; 3, j
P
u(new) n + Pj=1;2;3 Pn;n?e L?n1??e1Rn?e ;n rn n = u(old) ? j=1;2;3 Pn;n?e+e Ln?e+e Rn?e+e ;n rn + Pn;n?e L?n1?e Rn?e;n rn ; j
j
j
j
with
j
j
(30)
(31) rn = fn ? Lnu(old) n : Here we denote by Rm;n : Vn ! Vm, m n, the restriction operator de ned by Rm;nun = Rmun for all un2Vn L2( ). The prolongation operator Pn;m : Vm ! Vn can be de ned e.g. as the adjoint of Rm;n. The third remark following Theorem 2.1 shows how the two- or one-dimensional case can be treated similarly and we see that {for the one-dimensional case{ our approach reduces to the classical scheme. The approach (30) would imply three coarser levels to be active for a CGC, and {as was shown in the remark after Theorem 2.1{ we can do with only one coarser level by deriving the information on the levels jnj ? 3 and jnj ? 2 from the information on level jnj ? 1. If we consider the corrections from level jnj ? 1, (32) cn?e = L?n1?e Rn?e ;n rn; j = 1; 2; 3; j
j
j
as approximating a single (but unknown) correction function cn 2Vn, the corrections from the levels jnj ? 2 and jnj ? 3 can be computed as the mean values ? cn?e+e = 21 Rn?e+e ;n?e +1 cn?e +1 + Rn?e+e ;n?e ?1 cn?e ?1 ; (33) j = 1; 2; 3, and X (34) Rn?e;n?e+e cn?e+e : cn?e = 31 j
j
j
j
j
j
j =1;2;3
10
j
j
j
This is justi ed by the fact that the restrictions are commutative, i.e.
m n l ) Rm;nRn;l = Rm;l and the following (trivial) lemma.
Lemma 3.1 If all discrete operators Ln are stable and relatively consistent, i.e. kRn;n+e Ln+e ? Ln Rn;n+e k O(2?jnjp) then
j
j
j
kcn?e ? Rn?e ;n cnk O(2?jnjp) : j
j
The consistent discretisations can be derived e.g. from the ne grid discretisation Ln by taking the Galerkin approximation
Ln?e = Rn?e ;nLn Pn;n?e : If the three corrections cn?e were all restrictions of the (unknown) correction cn 2Vn indeed, then Rn?e+e ;n?e +1 cn?e +1 and Rn?e+e ;n?e ?1 cn?e ?1 would both have delivered the same result, viz., Rn?e+e ;n cn. This gives us the possibility to check how well such a function cn can be determined, by monitoring the quantities, j = 1; 2; 3, ? dn?e+e = 21 Rn?e+e ;n?e +1 cn?e +1 ? Rn?e+e ;n?e ?1 cn?e ?1 : (35) Summarising, our multigrid algorithm now reads: j
j
j
j
j
j
j
j
j
j
j
j
j
j
j
j
P
j
j
u(ni+1) = u(ni) + Pj=1;2;3 Pn;n?e cn?e (36) ? j=1;2;3 Pn;n?e+e cn?e+e + Pn;n?e cn;n?e ; where the corrections are given by (32), (33) and (34). This appears to be much similar to an multigrid algorithm for semi-coarsening, proposed by Mulder in [7]. The main di erence being that Mulder computes an approximation on the full grid Rn , whereas we compute the sparse grid approximation R^ m . The result of our algorithm is a solution on a sparse grid, i.e. a set of approximate solutions, viz. fun j jnj = mg, that are the solutions of the discrete equations Ln un = fn . All approximations un representing the same solution u of the continuous problem, we assume that they approximate the L2 ( )-projection of u in Vm = V(m;m;m) . To approximate this Rmu2Vm, we can construct a unique function um2Vm by means of the recursive interpolation formula that immediately follows from Theorem 2.1: j
j
j
uk =
X
j =1;2;3
uk?e ? uk?e+e + uk?e ; j
j
j
(37)
where 0 k m, jkj = m +1; :::; 3m; uk?e are the functions computed in the previous recursion cycle and uk?e+e and uk?e are approximations (possibly) derived as (33) and (34). In this way we nally obtain the unique representation j
j
um = v0 +
X
0jkjm
of
wk
um2Sm( ) V(m;m;m)( ) : This representation can be considered as the computed solution. 11
(38)
The same construction can be realised by a `decomposition and reconstruction' algorithm as used in wavelet theory [2]. Then the available approximate solutions fun g are decomposed into their components v0 and fwk g by
v0 =
P
jP nj=m R0 un jnj=m 1
and
P
nj=m Qk un wk = jP 1 jnj=m
and the reconstruction is performed by (38). This can conveniently be performed by a kind of a `pyramid algorithm'. This will be reported in a later paper. In practice, by the choice of the discrete operator our assumption that the L2 ( )projection of u was indeed consistently approximated on Vn , may not necessarily hold, and it can be checked by a monitor as (35). Now, e.g., the corresponding erroneous components might be removed from the sum (38). In the light of the treatment in Section 2 it is clear what restrictions and prolongations can be used between the di erent grids in the multigrid process. Because Vn L2 ( ), the obvious restriction Rn?e ;n is the L2 ( )-projection onto Vn?e : j
j
Rn?e ;n = Rn?e : j
j
This makes the diagram for the restrictions commutative: for any l m n we have Rl;mRm;n = Rl;n . An obvious prolongation can be the transposed restriction
Pn;n?e = RnT ?e ;n : j
j
However, this prolongation being of low order, it may be more appropriate to consider higher order prolongations. Such prolongations can always be represented by an additional operator Bn : Vn ! Vn so that we have
Pn;n?e = Bn RnT ?e : j
j
Here we will not elaborate on the di erent possibilities for Bn . The algorithm (36) shows that all relaxation processes for un on one and the same level m = jnj can be made in parallel. The cycling between the di erent (scalar) levels can be made as for the classical multigrid method: we can distinguish between V-, Wor F-cycles. However, in order to prove that the convergence of our multigrid-method is independent of the meshwidth, we now have to take into account that all aspect ratios will appear in the discretisations used.
4 Fourier convergence analysis In this section we rst summarise some results of Fourier analysis for more-dimensional discrete approximations and then we apply this to compute the convergence rate of our sparse-grid multiple-grid methods for the solution of the anisotropic Poisson equation. The approach is di erent from the usual treatment of Fourier analysis for multigrid for nite di erence methods for the following reasons. First, in view of the discretisation of conservation laws and divergence problems, we study nested box grids. This implies that mesh points in the coarse grids do not appear in the ne grids as well. The nesting of the (box) grids is di erent from the usual nesting of the ( nite di erence or nite element) grids. Second, we do not consider the usual sequences of ne and coarser meshes for multigrid methods, but we allow d-dimensional (d = 2; 3) sparse grids. Fourier analysis is one of the common tools to analyse linear constant coecient problems on regular grids, and it is particularly useful if the treatment of boundary conditions can be neglected. 12
In Section 4.1 we describe general tools that can be used for the Fourier analysis of functions de ned on more-dimensional box grids. The de nitions and theorems provide a useful machinery for the application of local mode analysis for the multigrid box-methods. For the technical proof and details related to this section we refer to [4]. In Section 4.2 we apply tools to analyse the multigrid algorithm introduced in section 3. The technical preparations in 4.1 allow us to be brief and clear in this treatment.
4.1 Fourier analysis for sparse box grids
For u2L2 (IR3 ) we introduce its Fourier Transform (FT) ub, scaled as Z e?ix! u(x) dx : ub(!) = (2)?3=2
(39)
Then we know that ub2L2 (IR3 ), and a back-transformation formula is available, Z e+ix! ub(!) d! ; u(x) = (2)?3=2
(40)
IR3
IR3
such that u(x) = u(x) almost everywhere on IR3 . Moreover, ub2L2 (IR3 ) and Parseval's equality holds: kukL2 (IR3) = k ubkL2(IR3) . We are interested in the Fourier transformation for an in nite set of equally spaced data. In this case the FT of such a \grid function" is a periodic function (a function de ned on a torus). Therefore we introduce a few de nitions. Let h2IR3 , h > 0, be given, then the h-periodisation of a function u : IR3 ! CI is de ned by X u(x ? kh); (41) u~(x) = 3 k2ZZ where kh = (k1h1 ; k2 h2 ; k3h3 ). We also introduce a notation for the three-dimensional torus Th = (?=h; =h] = (?=h1; =h1] ::: (?=h3 ; =h3]: (42) 3 Further we need the functions and Sinc , [1, pp.62,67] on IR , for jx j < 1=2 ; 1 i 3; i (x) = 10 otherwise (43) ; and iY =3 Sinc x = sinxxi : i=1
i
Using the relations mentioned in [1, p.98] we nd Y sinc ( 2!i ) = (2)?3=2 Sinc ( 2! ) : b (!) = (2)?3=2 i=1;2;3
(44)
For an h2IR3 , h > 0, we de ne the dilation operator Dh : L2 (IR3 ) ! L2 (IR3 ) by Dhf (x) = h?3=2 f (xh) ; (45) where h = (h1 h2 h3 )1=3 , and the convolution operator, ?, by (f ? g)(x) = (2)?3=2 We now know that Grid functions.
Z
IR3
f (y) g(x ? y) dy :
dh = D1=h fb and fd D ? g (!) = fb(!) bg(!) :
Here we introduce notations for the di erent types of grids and gridfunctions. 13
(46) (47)
De nition 4.1 For a xed mesh parameter h2IR3, h > 0, and for j 2ZZ 3 , we de ne an elementary cell by h;j = fx j jh < x < (j + e)hg, the volume of the cell is denoted by h3 = h1 h2 h3, and the box-grid is de ned by h = f h;j j j 2ZZ 3g : The regular in nite three-dimensional grid of vertices ZZh is de ned by ZZh = fjh j j 2ZZ 3g ; which should be well distinguished from the shifted grid which is de ned by ZZh? = f(j + e=2)h j j 2ZZ 3 g : Notice the relation with the grids as de ned in Section 2.2: n can be considered as a special case of h , and ?n as a special case of ZZh? . De nition 4.2 A complex or a real grid function uh is a mapping ZZh ! CI, or ZZh ! IR, and a shifted or box-grid function u?h is a mapping ZZh? ! CI or ZZh? ! IR. The vector space of such gridfunctions we denote by l(ZZh ) or l(ZZh? ), or brie y, by l. For any p 1 the space l(ZZh ) can be provided with a norm k kp
kuh kp = (h3
X
j 2ZZ 3
juh(jh)jp)1=p :
(48)
For a xed p, all grid functions with a nite norm k kp form a Banach space denoted by lp (ZZh ). For p = 2 we know that l2 (ZZh ) is a Hilbert space with the inner product,
< uh; vh >l2(ZZh) = h3
X
j 2ZZ 3
uh(jh)vh (jh) with uh; vh 2ZZh :
(49)
Similar de nitions are given for l(ZZh? ):
X ? ? < u?h ; vh? >l2(ZZ )= h3 uh (z)vh (z) with u?h ; vh? 2ZZh? : (50) h z2ZZh De nition 4.3 The direct restriction Rh : L2(IR3 ) ! l(ZZh ) is the operator that associates with a continuous7 u2L2 (IR3 ) the corresponding grid function on the grid ZZh , de ned by (Rh u)(jh) = u(jh) ; 8j2ZZ 3 ; (51) and the direct restriction Rh? : L2 (IR3 ) ! l(ZZh? ) on the shifted grid ZZh? is de ned by (Rh? u)((j + e=2)h) = u((j + e=2)h) ; 8j2ZZ 3 : (52) De nition 4.4 The box-restriction Rh : L2(IR3 ) ! L2(IR3 ) is the L2-projection on the piecewise constant functions on h , de ned by (cf. equation (6)) ?
?
(Rh u)(x) = h?3
Z
h j
u(z) dz 8x2 h;j :
(53)
;
The box-restriction RhB : L2 (IR3) ! l(ZZh? ) is de ned by RhB = Rh? Rh ; it associates the mean value of u on a cell h;j with the nodal value at the centre of h;j .
The box-restriction RhB u should be well distinguished from Rh? u. However, the L2 ( )projection Rh u in (53) and the restriction RhB u in (52) are conveniently related to each other by 2 3=2 Rh? ((Dh ) ? u) : (54) Rh u = h The Fourier transform of a grid function. Let u : ZZ ! CI be a grid function. We give the following
h
7
h
In case of a discontinuous function we can replace u by u as de ned in (40).
14
De nition 4.5 The Fourier transform uch 2L2(Th ) of a grid function uh 2l2 (ZZh) is a function Th ! CI, de ned by
h 3 X c e?ijh! uh (jh) : uh (!) = p 2 j 2ZZ 3
(55)
The inverse transformation is given by
u (jh) = h
1 3 Z p 2
!2Th
e+ijh! uch (!) d! :
(56)
Let u?h : ZZh? ! CI be a shifted grid function, then we have
De nition 4.6 The Fourier transform uc?h 2L2(Th ) of a shifted grid function u?h 2l2 (ZZh? ) is a function Th ! CI, de ned by
h 3 X ? c e?i(j +e=2)h! u?h((j + e=2)h) : uh (!) = p 2 j 2ZZ 3
(57)
Its inverse transformation is given by
u? ((j + e=2)h) = h
1 3 Z p 2
! 2 Th
e+i(j +e=2)h! uc?h(!) d! :
(58)
Remarks: We immediately see that uch (!) is [2=h]-periodic, whereas uc?h (!) is [2=h]antiperiodic, i.e. uc?h (! + 2=h) = (?)jej uc?h (!). We denote the Fourier transforms also by uch = F (uh) or uc?h = F (u?h); (59) i.e. we introduce the mapping F : l2 (ZZh ) ! L2 (Th) or F : l2 (ZZh? ) ! L2(Th ). At the end of this section we shall generalise this meaning of F . By the Parseval equality we have kuh k2 = k uch kL2(Th ) and ku?h k2 = k uc?h kL2(Th ) : (60)
Hence the Fourier transformation operators F : l2 (ZZh ) ! L2 (Th ) and F : l2 (ZZh? ) ! L2(Th ) are unitary operators. Because eh;! eh;!+2k=h or e?h;! (?)jekje?h;!+2k=h , for all k2ZZ 3, on a mesh of size h a frequency ! cannot be distinguished from a frequency ! + 2k=h. This phenomenon is called aliasing . The relation between FTs of a function restricted to different grids. In this section we rst present a few theorems associated with the di erent restrictions between two grids. We describe the relation between the FT of a continuous function de ned on IR3 and the FT of its restriction to the grid and then we show the relation between the FT of a ne grid function and the FT of its representation on a coarser grid. Next, we give the corresponding theorems for the prolongations. 15
Lemma 4.7 Let u2L2 (IR3) be a continuous function with FT ub. Its restriction uh to
the grid ZZh is de ned by (51). We have the following relation between ub and uch: X ub(! + 2k=h); (61) uch(!) = k2ZZ 3 i.e. uch is the [2=h]-periodisation of ub. Proof: For the proof we refer to [4]. 2 In the following lemmas q-restrictions are considered, with q2ZZ 3. Here q = (q1 ; q2 ; q3 ) is the coarsening factor, where usually 1 qj 2, j = 1; 2; 3. De nition 4.8 Let q2ZZ 3 with q > 0 and H = qh, then the canonical q-restriction Rq is the operator Rq : l(ZZh ) ! l(ZZqh ) = l(ZZH ), de ned by (Rq uh )(jH) = uH (jH) = uh (j qh) : (62) Theorem 4.9 We have the following relation between the FT of a grid function and that of its canonical q-restriction,
d X c Rq uh (!) = uh (! + 2p=h); 8!2TH ; H = qh: p2[0;q)
(63)
Proof: For the proof we refer to [4]. 2
Lemma 4.9 shows that, using the restriction Rq with q2ZZ 3 , q > 0, we get aliasing of q3 = q1:q2 :q3 frequencies onto one. Now we describe the relation between the Fourier transforms of a continuous function and its box restrictions. First we consider the direct restriction to the shifted grid, Rh? , and later the box-restriction, Rh and the q-restriction Rq? .
Lemma 4.10
? uc?h (!) = Rd h u(!) =
Proof: For the proof we refer to [4]. 2
X
k2ZZ 3
(?)jkj ub(! + 2k=h) :
(64)
For the Fourier transform of u2L2 (IR3 ) and Rh u2L2 (IR3 ) we have the following relation.
Theorem 4.11
uch (!) = Rd h u(!) =
X k2ZZ 3
(?)jkj Sinc ( h! 2 + k) ub(! + 2k=h) :
Proof: Using (54) we see uc?h (!) = F (Rh u)(!) ? = 2h 3=2 F (Rh? ((Dh ) ? u))(!) ? P = ? 2h 3=2 k (?)jkj F ((Dh ) ? u)(! + 2k=h) P = ? 2h 3=2 k (?)jkj (Dd h )(! + 2k=h) ub(! + 2k=h) 3=2 P 2 j k j = h k (?) D1=h b (! + 2k=h) ub(! + 2k=h) : P = h?3=2 k (?)jkj D1=h Sinc ( 2! + k=h) ub(! + 2k=h) : P = k (?)jkj Sinc ( h! 2 + k) ub(! + 2 k=h) : 2
16
(65)
(66)
De nition 4.12 Let q2ZZ 3 with q > 0 and H = qh, then, for s2[0; q), the s-frequency decomposition q-restriction is the operator Rqs : l(ZZh ) ! l(ZZqh ) = l(ZZH ) is de ned
by
(Rqs u?h)((j + e=2)H) = q ?3
X
(?)sk u?h ((qj + k + e=2)h) ;
k2[0;q)
(67)
where q 3 denotes q 3 = q1 q2 q3 .
Remarks: In the case s = 0 we call Rqs = Rq0 = Rq? simply the q-restriction. From the construction of the restriction operators RhB and Rq? it is clear that the following relation holds: B = R? RB : RH q h 3 , H = 2h, we have Theorem 4.13 Let q = 2e 2ZZ 3, then for all !2TH X m h! + (m + s) c? ( Rd (?) Cos( ) uh (! + m=h) (68) qs u?h )(!) = 2 q?3 is 2 m2[0;q) where
Cos(h!=2) =
Y
j =1;2;3
cos(hj !j =2) :
Proof: For the proof we refer to [4]. 2 De nition 4.14 The natural box-prolongation Ph? : l2 (ZZh? ) ! L2(IR3 ) is de ned by u(x) = (Ph? u?h )(x) = uh ((m + e=2)h); for all m2ZZ 3 and mh < x < (m + e)h. We also introduce a natural prolongation, Pq? , ? ) ! l2 (ZZ ? ) by P ? = R? P ? where from a coarse to a ne gridfunction Pq? : l2 (ZZH q h H h H = qh. The prolongation PhB : l(ZZh? ) ! L2 (IR3 ) is the operator dual to RhB in the sense that for all u?h 2l2 (ZZh? ) and v2L2 (IR3 ) we have
< PhB u?h ; v >L2(IR3) = < u?h; RhB v >l2 (ZZ
: (69) h) The following theorems show how we nd the FT of the prolongation if the FT of the source function is given.
Lemma 4.15
?
=2) c? ub(!) = F (Ph? u?h )(!) = Sin(h! h!=2 uh (!) :
(70)
Proof: For the proof we refer to [4]. 2 Theorem 4.16 With H = qh and q3 = q1 :q2 :q3, we have for the FT of the prolongation of a box gridfunction
=2) ud ? uc?h (!) = F (Pq? u?H )(!) = q?3 Sin(H! Sin(h!=2) H (!) :
Proof: For the proof we refer to [4]. 2 17
(71)
The Fourier transform of operators involving different grids. In (68) and (71) we see that, by the restriction and prolongation between functions on grids n and n+q , aliasing takes place and that q 3 frequencies on n+q correspond with a single frequency on n. This implies that, analysing a multigrid algorithm, we have to study the behaviour of the q 3 aliasing frequencies together. Collecting the q 3 corresponding amplitudes of the aliasing frequencies in a single q 3-vector, we extend the de nition (59) of F to the case q 3 > 1 and obtain F : l2 (ZZh ) ! [L2 (Tqh )]q3 or F : l2 (ZZh? ) ! [L2(Tqh)]q3 by ? F (uh)(!) = ubh (! + m=h) m2[0;q) ; !2Tqh : (72)
With these amplitude-vectors F (uh )(!) and F (u?h )(!) , we can introduce the linear operators F (Rq )(!) and F (Rq? )(!) by
and
F (Rq uh )(!) = F (Rq )(!) F (uh )(!) ; !2Tqh ;
(73)
F (Rq? u?h )(!) = F (Rq? )(!) F (u?h )(!) ; !2Tqh :
(74) We call the new operators, that depend on !, the Fourier transforms of the original operators. The new operators are ` q 3 ` matrices if ` aliasing frequencies are considered on the coarse grid. Similar to the restrictions, the prolongations can be associated with their Fourier transforms. (75) F (Pq uh )(!) = F (Pq )(!) F (uh)(!) ; !2Tqh ; and F (Pq? u?h )(!) = F (Pq? )(!) F (u?h)(!) ; !2Tqh : (76) These operators are q 3 ` ` matrices. For arbitrary linear constant coecient operators Ah : l2 (ZZh ) ! l2 (ZZh ), its Fourier transform F (Ah) : L2(Th ) ! L2 (Th ), can also be considered as a q 3 ` q 3 ` diagonal matrix F (Ahuh )(!) = F (Ah)(!) F (uh )(!) 8!2Tqh : Because of Parseval's equality we know that
( F (Ah) ) (!) k F (Ah)(!)k2 = !max kAhk2 = !max 2Tqh 2Tqh
(77)
with (A) the spectral norm (the largest singular value) of the matrix A, and
( F (Ah) ) (!) (78) (Ah ) = !max 2Tqh where denotes the spectral radius. This provides us with the means to study the norm and the spectral radius of the error-ampli cation operator of the multigrid iteration.
4.2 Fourier analysis convergence results
To get some insight in the behaviour of the more-dimensional multigrid algorithm introduced in Section 3, we use Fourier analysis to determine the convergence rate of the two-level algorithm for the two-dimensional anisotropic Poisson equation
uxx + "2 uyy = f ; discretised by the usual 5-point discretisation. 18
(79)
The error-ampli cation operator, Mn, of the two-level cycle (with pre- and postrelaxation steps) for the solution of (29) is described by
e(ni+1) = Mn e(ni) = Sn Cn Sn e(ni) ; where Sn denotes the smoothing, e.g. damped Jacobi iteration:
?
(80)
e(new) (81) n = Sn e(old) n = In ? Dn?1 Ln e(old) n ; with the damping parameter and Dn the main diagonal of the discrete operator Ln . The coarse grid correction is described by (30) and (31). For gridfunctions uh 2l2 (ZZh? ), i.e. neglecting boundary conditions, we nd, using (30) F (Mn )(!) = F (Sn ) (!) F (Cn )(!) F (Sn) (!) ; (82) with F (Sn ) = F (In) ? F (Dn)?1 F (Ln) ; (83) and F (Cn) = FX(In) (84) ? F (Pn;n?e ) F (Ln?e )?1 F (Rn?e ;n) F (Ln) +
j =1;2;3
X
j =1;2;3
j
j
j
F (Pn;n?e+e ) F (Ln?e+e )?1 F (Rn?e+e ;n) F (Ln) j
j
j
? F (Pn;n?e) F (Ln?e)?1 F (Rn?e;n) F (Ln) : To get an impression of the behaviour of the algorithm, keeping the explicit computation to reasonable size, we elaborate (82) for the equation (79), for the two-dimensional case, with q = (2; 2) and = = 1. Then F (Mn )(!) is a 4 4-matrix, which we derive from (82), (83) and (cf. (15)) F (Cn ) = F (In) ? F (Pn;n?e1 ) F (Ln?e1 )?1 F (Rn?e1;n ) F (Ln) ? F (Pn;n?e2 ) F (Ln?e2 )?1 F (Rn?e2;n ) F (Ln) + F (Pn;n?e) F (Ln?e )?1 F (Rn?e;n ) F (Ln) : From (68) and (71) we know F (Rn?e1;n ) = cos !0 2 h2 cos !0 2h2 sin !02 h2 sin !02h2 ;
cos ! h sin ! h 1 1 1 1 F (Rn?e2 ;n) =
and
0
0
0 0 cos !1 h1 sin !1 h1
0 cos(! h ) cos(! h ) 1T 1 1 2 2 B C : sin( ! h ) cos( ! 1 1 2 B F (Rn?e;n ) = @ cos(! h ) sin(! hh2 )) C A 1 1
2 2
sin(!1 h1 ) sin(!2 h2 ) So that with F (Pn;n?e1 ) = F (Rn?e1 ;n)T , F (Pn;n?e2 ) = F (Rn?e2 ;n)T and F (Pn;n?e) = F (Rn?e;n)T , the norm kMnk and the spectral radius (Mn ) can be computed by means of (77) and (78). To study the convergence behaviour of our algorithm, we consider the matrices (83), (84) and (82) as a function of ! 2Th = [?=h; =h]2, and for each ! we compute its 19
1
1
0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0
-0.2
-0.2
-0.4 0.8
0.8
0.6
0.6
0.4
0.4 0.2 0.2
0.6
0.4
0.8
0.2
(a) The eigenvalues of F (Sn)(!). 0 0
0.2
0.8
0.6
0.4
(b) The eigenvalues of F (Cn )(!). 0 0
0.1 0.08 0.3 0.06 0.25 0.04 0.2
0.02
0.15
0 -0.02
0.1
-0.04
0.05
0.8
0 0.8 0.6 0.6 0.4
0.8 0.6
0.4 0.2 0.2
0.4
0.6
(c) The eigenvalues of F (Mn)(!).
0 0
0.4
0.2
0.8
0.2
(d) The singular values of F (Mn )(!). 0
0
Figure 1: The frequency response of the operators Sn , Cn and Mn , for " = 1, q = (2; 2) and = 2=3. eigenvalues and singular values of these matrices. We show these values in Fig. 1 for the case = 2=3, " = 1. Without loss of generality we can take h = (1; 1), the parameter " taking care of the anisotropy. The damping parameter 2[0; 1] for the Jacobi relaxation, can be chosen freely. We select the value = 2=3 because it minimises
( F (Sn)(!)) : !=(0;=h)max ;( ;0);(=h;=h) This means that = 2=3 makes Sn a well balanced smoother for the di erent types of high frequencies (see Fig. 1.a). In Fig. 1.a we show the eigenvalues of the smoothing operator, and in Fig. 1.b of the coarse grid correction. In this gure we see that one eigenvalue of F (Cn ) is always equal to one. This eigenvalue corresponds with the highest frequencies, for which no correction can be obtained from any of the three coarser grids. The combined e ect of the smoother and the coarse grid correction is seen in Fig. 1.c, which shows that sup! (Mn (!)) 1=9, and also in Fig. 1.d, where we see sup! kMn (!)k 1=3. The rather low maximal values show that {at least for square ne-grid cells{ the multigrid algorithm has good convergence behaviour. Because it is important that the algorithm is e ective for an arbitrary cell aspect ratio, in the Figs 2 - 4 we show the singular values of F (Mn)(!) also for " = 1=8 and for the limit as " ! 0. Now it appears that for high aspect ratios the convergence behaviour deteriorates. We nd sup! lim"!0 ( F (Mn )(!)) 5. This implies that we cannot always guarantee error reduction for a single iteration sweep. Therefore we show in Fig. 4 also the behaviour of Mn2 (!). This shows that a couple of two consecutive 20
1.6 1.4 1.2
0.2
1
0.1
0.8
0
0.6
-0.1
0.4
0.8
-0.2 0.8
0.2 0.6 0 0.8
0.6 0.4
0.8 0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
(a) The eigenvalues of F (Mn)(!). 0
0.2
(b) The singular values of F (Mn)(!).
0
0
0
Figure 2: The frequency response of the operator Mn , for " = 1=8, q = (2; 2) and = 2=3. iteration steps does guarantee error reduction, and the convergence rate is signi cant: sup "lim ( F (Mn2 )(!)) 1=9 : ! !0 As a consequence we can expect that a W-type cycle of the multigrid algorithm will have good convergence properties.
From the computations of which the results are summarised in the Figs 2 - 4, we conclude that the eigenvalues of the iteration matrix are less that 1, uniformly in the parameter ". In fact, max" (Mn ) 0:33 and max" kMn k 5:0 and max" kMn2 k 0:11. The fact that kMn k > 1 and kMn2 k
