- Email: [email protected]
A study of some multigrid ideas
A Study of Some Multigrid Ideas* David Kamowitz and Seymour V. Parter Computer Sciences Department University of Wisconsin - Mad&m Madison, Wisconsin 53706
Transmitted by S. McCormick
ABSTRACT In an effort to understand certain ideas and concepts associated with multigrid iterations, we give an indepth study of a particular simple problem. We consider a standard finite-difference system associated with the twqoint boundary-value problem
-(pd)‘+bd+qu=o,
u(0) = u(1) = 0
The operators Zth, Zt,, are “operator” based interpolation and projection operators, while the smoothers are the damped Jacobi iterations with parameter a > 0. We determine the exact rates of convergence for the “twc@d” scheme and upper bounds ( < 1) for the multigrid schemes. Experimental results are discussed.
1.
INTRODUCTION
The multigrid approach for the numerical solution of boundary-value problems for elliptic partial differential equations is proving itself as one of the fastest and most efficient methods-see [l], [3], [4], [5], [6], [17]. Moreover, there are a large number of theoretical papers on this subject-see [2], [3], [6], [7], [8], [9], [lo], [ll], [13], [16]. Nevertheless, it seems (at least to these authors) that we are just beginning to understand this powerful idea. In particular, there are the questions: how do we choose the interpolation and *Supportedby the U.S. Air Force Office of Scientific Research under Contract No. AFOSRSZ 0275. APPLIED
MATHEMATICS
AND COMPUTATION
0 Elsevier Science Publishing Co., Inc., 1985 52 Vanderbilt Ave., New York, NY 10017
17:153-184 (1985)
153 -/85/$03.30
DAVID KAMOWI’IZ
154
AND SEYMOUR V. PARTER
projection operators? How do we choose the smoothing operators? What do we mean when we say smoothing? etc. This report is a reflection of our effects to understand and appreciate the theoretical insights of Frederickson [7], McCormick, and Ruge [13], McCormick [14, 151, and Greenbaum [8] and apply those ideas to extend the explicit convergence rates given by Hackbusch [12, (2.21); 111 for the very simplest problem uf’=
u(0) = u( 1) = 0.
f,
Specifically, we consider the twopoint boundary-value problem
Lu := -
(put)‘+ b(x)u’+ qu = f, u(0) = u(l) = 0
O
(1.1)
(1.2)
where p(x), b(x), Q(X) are smooth functions and
P(X)2
PO ’
0,
q(x) > 0.
(1.3)
In Section 2 we describe a basic approach to multigrid which is based on the ideas of Frederickson [7], McCormick and Ruge [ 131, and Greenbaum [8]. In Section 3 we describe a discretization (finite-difference) of the problem (l.l), (1.2) and a specific twogrid iterative procedure for its solution. In Section 4 we describe the class of damped Jacobi “smoothers” and use our knowledge of these schemes and their eigenvectors to describe the basic spaces: Range: Z& = 9 and Nullspace I,2hLh = q. In Section 5 we obtain estimates for the norm decay of a single step in a twogrid scheme for two different norms. In addition we obtain a better estimate for the norm decay for all iterative steps beyond the first in both these norms. Interestingly enough, this estimate is the same in both norms. This latter result is an improvement over the estimates of Hackbusch [ll, 121. It is well known that the problem (l.l), (1.2) is equivalent to a self-adjoint problem. Moreover, the discretization (3.4) is also equivalent to a symmetric problem. In fact, our multigrid treatment of this problem is equivalent to the “same” multigrid treatment of this symmetric problem. This equivalence is not needed for the discussion in Sections l-5. However, as we turn to the extension of a twogrid scheme to a true multigrid scheme, we require this information. In Section 6 we demonstrate this equivalence. In Section 7 we describe the n-grid “sawtooth” multigrid schemes and give a theory (closely
155
A Study of Some Multigrid Ideas
related to a theorem of McCormick [14, 151) which describes the rates of convergence of this scheme. In addition Section 7 contains some experimental results.
2.
A BASIC THEORY
The theory presented in this section is based on the work of Frederickson [7], McCormick and Ruge 1131, and Greenbaum [8]. We consider a finite dimensional linear space S, and a problem
L/p= f,
U,f
ESh,
(24
where Lh:S,-+Sh
is a linear, nonsingular operator. Multigrid is an iterative method for the solution of this problem. The basic idea is to utilize another finite-dimensional space S,, with dim Ss,, < dim S,.
(2.2)
Hence we require operators I,2h, Z& which enable us to effect communication between these spaces. In particular, we have
z,zh:s, + s,,
(projection),
(2.3a)
I,:, : s,, + s,
(interpolation),
(2.3b)
where Zih, Z&, are linear operators. We also require a “smoothing” operator y,, and a “COaW2 grid” Operator &,’ The SImOthing OperatOr y,, iS an affine operator which has U, the unique solution of (2.1), as its only fixed point. That is, (2.4a) where G, : Sh + S, is a linear operator, and if u is the solution of (2.1), then YhU =
u.
(2.4b)
DAVID KAMOWITZ
156
AND SEYMOUR V. PARTER
Finally, the “coarse grid” operator
(2.5) is a linear, nonsingular operator taking Ssh onto itself. Let U” E S, be a guess for the solution U of (2.1). Set (2.6a) (2.6b) (26c)
r=f-L,cl=LL,(U-O)=L,E, R=12hrh
(26d) (26e)
*
Solve
L,,P=
R,
i.e.,
2= L,;;‘R.
(2.7)
Set
(2.8) uo := u’, and return to (2.6a).
REMARK. While one might say that we have merely described a twogrid scheme, the iterative scheme described above does, in fact, describe “multigrid’ schemes. The point is that the “coarse-grid correction” operator L,, may be a complicated procedure involving more grids. The work of Frederickson [7], McCormick and Ruge [13], and Greenbaum [8] suggests we study two fundamental subspaces. These are (2.9a) n = n&pace
of IihLh.
(2.9b)
157
A Study of Some Multigrid Ideas In addition we consider a special two-grid “coarse grid” operator e
2h
= 12hL h Ih2h*
(2.10)
h
This particular operator is the “Gale&n choice” and is optimal in a certain sense. This fact is emphasized by the following
LEMMA 2.1. with
Consti
one iteration
as &scribed
L,, = i,, suppose
EE 9.
above
by (2.6a)-(2.8)
.
(2.11)
That is, suppose there is a w(2h) E s2,, so that E= a(h) = z;hW(2h).
(2.12)
C= w(2h)
(2.13a)
U’=U.
(2.13b)
Then
and
Hence, the problem is solved.
PROOF. From (2.6d) and (2.12) we have r = L,I = L,,Zthw(2h). Thus R=
z,fhr= ( z,2hL,z;h)w(2h).
That is, R = i2hW(2h).
(2.14)
Hence, from (2.7) we have (2.13a). Finally, (2.13b) follows from (2.8) and (2.12). n
158
DAVID KAMOWITZ
AND SEYMOUR V. PARTER
Now suppose we can write S, as the direct sum (not necessarily orthogonal) of Range Z& and Nullspace ZlhL,. That is, every grid function w(h) E S, can be uniquely written as w(h) = Z,h,w’(2h) + w”(h),
(2.15a)
where z,29+?(
h) = 0.
(2.15b)
Let us apply this decomposition to the intermediate error I = q(h). Then E(h) = Z&w’(2h) + w”(h).
(2.16)
We will now show that &1= w2(h). We have r = L,z;~w1(2h)
+ L,w2( h)
and R = Z;% = [ &$!d,z,h,] w’(2h). Using Lemma 2.1 we see that
u’ = 0 + z;,$(2h) and hence
That is, &I= W”(h). Thus, the convergence of the process can be tested by
where
(2.17)
159
A Study of Some Multigrid Ideas
and the norm involved is any norm. The authors mentioned above all dealt with self-adjoint L, and assumed that
zlh = c(Zih)T In this case it is an easy matter to see that Range Zih and NuBspace ZihLh are Lb-orthogonal complements of each other in S,. In the general case we must assume the decomposition (2.15a). However, in Section 4 we will establish the validity of (2.15a) from the following simple criteria: (2.15a) holds if L,, is nonsingular and (i) rank Z,:, = dim Sa,,, (ii) dim Nukpace ZihLh Z dim Sh - dim S,,. This follows from a simple counting argument and the fact that Lemma 2.1 implies that the zero vector is the only vector common to both subspaces. 3.
THE DISCRETE
PROBLEM
Let an integer N > 1 be chosen, set 1
’
h=
2(&7+1)
(3.1)
= (2N+1)+1’
and let the fine “grid points” be given by
xj( h)
j=O,1,2
= jh,
,..., 2N+2.
(3.2)
Define the difference operator Lh by
[Lhu]k=
(3.3a)
-akU,-,+PkUk--YkUk+l~
where
1
b,
pk-l/2
ak=
-
h2
I
pk+l/2
Bk =
Pk+l/2 Yk =
7-s
+2h
1’
+ pk-l/2
+ok
h2 bk I
*
1 9
(3.3b)
Then L, is a consistent approximation to the operator L described in (1.1).
160
DAVID KAMOWITZ
AND SEYMOUR V. PARTER
We assume h is so small that lyk> 0, yk > 0 for all k. We are concerned with the system of linear equations
P%~lk=f,~
k=1,2
,..., 2N+l,
u, = u,,,, = 0.
(3.4)
We shall solve this system with a particular multigrid iterative scheme. Consider a coarse grid on which we have a mesh spacing of 2h and the coarse grid points are given by X,(2h)
= 2hk,
k=0,1,2
,..., N+l.
Let S,bethespaceofgridfunctions{U,(h); k=O,1,...,2N+2}definedon the fine grid, and S,, the space of grid functions { U,(2h); k = 0, 1,. . . , N + l} defined on the coarse grid. Our first step is to construct the interpolation operator I& which maps Ss,, into S,, i.e.,
Following the experience of Dendy [4, 51 we define I!& according to [ Z$,U(2h)],,
= U,(2h)
(common points)
(3.5a)
and (new points).
This choice of “operator” interpolation may be described in the following way: If the physical point x j( h) of the fine grid is also a physical point of the coarse grid, i.e., if j is even, say j = 2k, we set Uj(h)= U,(2h) equal to the same value as the coarse grid function assumed at that point; ‘hat is, (3.5a) holds. If the physical point xl(h) of the fine grid is not a point of the coarse grid, i.e., j is odd, say j = 2k - 1, we require that
~~&i’S@h)]},,_, =o.
(3.5c)
161
A Study of Some Multigrid Ideas
We formalize this remark as
LEMMA 3.1. Let the interpolation operator l& be defined by (3.5a), (3Sb). Then a fin&on U(h) E S, is in the range of I.$ if and only if k=1,2
[L,U(h)],,_,=O, The projection
,...,
w
N+l.
operator Iih: S, + S,, is defined by
REMARK. If b(x) = 0 and the operator L is self-adjoint, then ffk=Yk-l
and we see that z;“=C[z;h]T.
(3.7)
As we have said in Section 2, the relationship (3.7) between Zih and Zlh is the “ variational choice” and is frequently recommended for selfadjoint problems-see [13], [7], [4]. For the purposes of this exposition we take the optimal choice of “coarsegrid correction,” i.e., L 2h :=e
2h
=12hL h
h
Ih
2h’
REMARK. A direct computation shows that
[L,,U@h)],
= - (Y~2)uk-1+fi~2’uk-
yi2)u,+,,
(3.8a)
where a2ka2k- 1 OLk (2, = 12 [ 82k-1
(3.8b) 1 Y2ka2k+
1 (2) _
fi2k-,
Y2k-t2k+l
yk
-
ii
[
p2k+l
1
b2k+l
I
1
(3.&)
(3.8d)
DAVID KAMOWITZ AND SEYMOUR V. PARTER
162
Hence, (3.3a).
L,, is again a diagonally dominant three-term operator of the form
From the theory developed in Section 2 we see that, if we were to choose red-black Gauss-Seidel (indeed, black Gauss-Seidel only) as our smoother, the twogrid scheme would converge in one iteration. Thus we would have an explicit direct solver. However, since we are interested in multigrid, we will choose other smoothing operators. 4.
JACOBI
ITERATIVE
SCHEMES
A stationary iterative method for the solution of (3.4) is described by an operator splitting of the form
L,=M-N. Then, given a first guess Ua, we define successive iterates by the formula MUi+‘=
NUj+
f.
(4.1)
The convergence of this scheme is determined by the eigenvalue problem
XMU=NU.
(4.2)
As is well known, the scheme (4.1) is convergent iff maxIhI< 1. It is easy to verify that if (h, +) is an eigenpair of (4.2), then
L,&= (l-
A)M+
(4.3)
In this section we are concerned with a particularly simple class of such iterative methods, namely the Jacobi methods. We set M=(l+a)B,
a 2 0,
(4.4a)
where B = diag(p,).
(4.4b)
A Study of Some Multigrid Ideas
163
In this case we may rewrite (4.1) as
&B-Q-
ui+l=ui+
L,uj).
When a = 0 we call this scheme the Jacobi method. When a > 0 we call this scheme a damped Jacobi method. In these cases we are able to give a relatively complete discussionof the eigenvalue problem (4.2). We have the following facts: (i) The method is convergent for a > 0. (ii) Let Q = 0. Let (p, Cp)be an eigenpair, i.e.,
Let 6 be defined by (4.6) Then ( - /.L, 4) is also an eigenpair. The eigenvalues p are real and distinct. Furthermore, as h + 0 the { p } fill out the interval [ - 1,11.For completeness we repeat the basic relationshipbetween + and 6. Namely, if
k isodd,
+k = ‘?‘k;
(4.7a)
if
kiseven,
+k =
(433)
-
6,.
Since dim S, = 2N + 1, there is a single eigenvector 4 associated with the eigenvalue p = 0. This eigenvector satisfies i,,
(4.8)
= 0.
For u > 0 there are corresponding eigenpairs (A, +>. (x,6>, are the eigenfunctions described in (4.7) and
h_Cl+a
where 9 and 6
X,a-cl
1+u'
1+u'
The eigenpair (0,Cp) corresponds to an eigenpair (A, (p) with
j+a
lfu’
(4.9)
DAVID KAMOWITZ
164
AND SEYMOUR V. PARTER
We now turn to the determination of the three important subspaces
RangeZth, Nullspace If",
Nullspace ZLhL,,
where the operators Z&, Zih are given by (3.5a), (3.5b), and (3.6).
LEMMA4.1. Let (p, $J), ( - p, 4) be the two eigenpairs described above with p # 0. Let @=
N1+d+(1-d61.
in (ii)
(4.10)
Then Q E Range Z&. Further, since the vectors Cp corresponding to different are linearly independent, this construction provides N elements of Range Zl,,.
pairs (P,+),(--II~~) linearly independent
PROOF. Using (4.3) we have
Lh@=
LO+do - CL)W- Cl- cL)(1+PM4
=(1+P)(w4m41. Thus, the lemma follows from Lemma 3.1, which characterizes Range Zih, n and from (4.7a), (4.7b).
LEMMA4.2. (i) For all choices of a > 0, the vector B? (associated with p = 0) is an element of NullspaceZth and hence (4.3) implies that 6 is an element of Nullspace ZfhL,. (ii) Let (CL,$I), ( - p, 6) with p # 0 be the two eigenpairs described above. Let
‘~=B[(~-cL)~+(~+P)+].
(4.11a)
Then \kE NullspaceZih and (4.3) implies that C#I + 4 E Nullspace ZihL,.
(4.11b)
165
A Study of Some MultigridZdeas
Further, since the vectors \k associated with different pairs (CL,+), ( - p, 6) are linearly independent, we have N linearly independent elements of this nullspace and N linearly independent vectors of Nullspace IthLh. The vectors Bc&,6 provide one more independent vector of each of these subspaces.
PROOF. The result follows from a direct computation using (3.6) and the n defining eigenvalue problem.
COROLLARY.
N&pace
There is a unique decomposition
of sh into RangeI.& and
zfhL,.
PROOF. Since L, is nonsing.ular, we have shown that dim Range Zih > N, dim Nullspace zE”Lh = dim N&pace
Zih > N + 1,
dimSh=2N+1. Thus the corollary follows from the observation that Lemma 3.1 implies that these two subspaces have only the zero vector in common. n
5.
SOME ESTIMATES
Let a > 0. We take as our smoother m applications of the corresponding Jacobi iteration. That is, given U ’ = U’(h), we obtain 6 [as in (2.6b)] from the formula Mui+‘=
NUj+
f,
j=O,l,...,
m-l,
G=Urn*
(5.lb)
First let us consider the special vector 6 with its associated eigenvalue h=L Suppose
(5.la)
l+a’
DAVID KAMOWITZ
166
AND SEYMOUR V. PARTER
Then
L,E = CA”‘B+, and, using Lemma 4.2, we see that ZihL hE -= 0 . Hence, in this case, for any norm,
(5.2)
We now consider two special norms defined on S,. Let w, v E S, be of the form
w =
5 A,c#B~+&+ 5 dj$j,
j=l
v =
f
(5.3a)
j=l
cjgBj +Gj
j=l
+ f ej&j.
(5.3b)
j=l
Define
(w,v)~=
f
A,C,+dc+
j=l
(w,u)~=
z
5 djcj, j=l
AjC,(l-~j)+AC+,~l~jd,(l+pi),
j=l
Ilwlli = (WI w>o, llwllf=(w,~)~=
; j=l
A;(l-pj)+ii2+
5 d2j(l+pj). j=l
A Study of Some Multigrid Z&as
167
REMARK. In the simplest case where b = 0 and a and c are constant, the vectors $j are orthogonal and Ilwll; a
lbllr2,~
Ibaa (LhW7w>= lb&f,. LEMMA5.1.
Suppose e’=tJ-U’=ccp+d;,
(5.4)
and let .9 be defined by the two-grid iteration scheme. Then
(5.5a)
IIM max-=L[(l+(l)2m(l-p)+(~)2m(l+p)]. (5Sb) 2 l+a c,d
ll~“ll~
PROOF. From (5.4) we see that I = CA’?++ dk”& Following the theory of Section 2, we write (see 2.16) (5.6)
where Z;hL,e’(
h) = 0.
We claim that Z&w’(2h)
= “”
; d^xm[(l+
81 = A”( 1 - cl) + dk( 2
/.L)$- (1 - /A)+] 1 + ,u)
[~+~I.
(5.7a)
(5.7b)
DAVID KAMOWITZ
168
AND SEYMOUR V. PARTER
To verify this we need merely verify that the sum of the right-hand sides of (5.7a) and (5.7b) is E, and use (4.10) and (4.11b). Having verified (5.6), (5.7a), (5.7b) we proceed as follows: 11&111~=~[cXm(l-c,)+dlXm(l+CL)]2,
(5.8a)
lleO1l;= c2 + d2,
(5.8b)
ll~‘ll;=&A~(l-~)+d~~(l+~)]~,
(5.&)
ll~‘ll: = c2(1 - Jo)+ d2(1+
(5.8d)
p).
Thus, a simple argument shows that sllp~=~[A~“(l-p)2+P~(l+p)2],
(5.10a)
0
sup~=I[hZ”(1-p)+~2*(I+lr)]. (5.1Ob) 1
Using the basic formulae (4.9) we obtain
Il4G sup ,l$l,;
-=q(q2m(l-p)2+(~)2m(l+p)2], 2 1+a
Il~‘ll: sup 11~011:
-=+[(~)2m(l-~)+(~)2m(l+~)].
Thus, the lemma is proven.
THEOREM
5.1.
In the general case we have
(5.11a)
(5.11b)
169
A Study of Some Multigrid Ideas
Moreover, there is a choice of e” which makes these estimates sharp (with an error of order h).
PROOF. The theorem follows immediately from the previous lemma.
n
We observe that (5.11a) is precisely the formula obtained by Hackbusch [12, (2.21)] in the special case p(x) = 1, b(x) = q(x) = 0, and a = 1. To make a complete identification we merely set
1-p
u=2’
1_,=+.
(5.12)
However, while (S.lla), (5,llb) describe the worst-case decay in one multigrid iteration in a twegrid scheme, it does not give the estimate of real interest. From the discussion in the proof of Lemma 5.1-and (5.7) in particular-we realize that, even though the constants c and d of (5.4) may be arbitrary for E’, that is not true for Ed, k > 1. We have [from (5.7b)]
El=.[++~]
(5.13)
with (J=
~h~(l--~)+d~~(l+~) 2
Therefore, following the argument of Lemma 5.1,
Hence [l&q; = Il&(‘)ll~ = 2a2 , ll&‘2’~~~=ll&‘2’11~=~[~~(1-p)+X”(l+~)]2, and for j = 0,l we have
~=i[h-(l-p)ix_(l+p)l”. I
(5.14)
DAVID KAMOWITZ
170
AND SEYMOUR V. PARTER
Thus we have proven
THEOREM 5.2.
forj = 0,l and all k > 1 we have
In the general case,
II&(k+q;
=
lle’k’ll;
sup
(5.15)
{+[Xm(l-/&)+1~(l+r_l)]“}.
-l
REMARK. The distinction between Hackbusch’s result (5.11a) and (5.15) is nontrivial for large m. We have, as m + co,
ll4lO -_--
l+a
1
11~0110fi
(5.16a)
fT7n’
while
II&(k+‘)llj _-- It-a llEkllj
2
1
(5.16b)
m *
Thus for k >, 1 we have
(5.17)
6.
SYMMETRIZATION Consider the difference equation (3.4) described by (3.3a), (3.3b). Let (6.1)
& = dkvk, where the coefficients
d, are computed recursively by (6.2a)
do = 1, ^.
d k+l- -d
Uk+l
k
k=O,l,...
d
-fk
’
.
(6.2b)
171
A Study of SOme Multigrid ~deos
Then (3.4) becomes
or
(6.3)
that is, (6.3a) where (6.4a)
(6.4b)
Hence iik = ?k- 1,
(6.5)
so that E, is symmetric. We now turn our attention to the Jacobi iterative schemes of Section 4. We have ui+l=ui+
&rl(
f
- L,uj).
(6.6)
The change of variables (6.1) is conveniently described by U=DV
(6.7a)
where D=diag(d,,...,d,,+,).
(6-W
DAVID KAMOWITZ
172
AND SEYMOUR V. PARTER
With this change of variables the iteration (6.6) becomes Vi+r=Vi+
(f-
L,DVj).
But, since D and B are both diagonal matrices, we have vi+l=vi+
&3-‘(D-‘f-
E,vj),
(6.8a)
where
D-‘f=(~,~,...,fZiN+l)‘.
(6.8b)
Since B is also the diagonal of i,, (6.8) is precisely the same Jacobi iterative scheme for modified symmetric equations. Now let us study the effect of this transformation on the twogrid iterative scheme. We compute
Imagine U(2h), = d,,V(2h),
is given in S,,. Then, from (35a),
Gf [w@w]2k = Vcwk~ d,:-z[ZzhhU(2h)],,-z=V(2h)k-l.
(6.9a) (6.9b)
Thus (3.5b) yields
d,;-,[ZkWh)],,_, = ~[a2,_,d,:_,d,,_,V(2~)k-,+y,t-,d~~-,d,,V(2~)k] 2k 1 =
jj+[h2k-lv(2h)k-l+ 2k 1
?2k-1vt2h)k]*
Thus, with this change of variables the mapping Z$, of our original unsymmetric problem becomes i&, the appropriate mapping associated with the new symmetric problem. Finally, let us consider
A Study of Some Multigrid Ideas
173
A straightforward calculation verifies that
TV,,-,+V,,I,,+~V(h),,,, Zk+l
. I
Thus, following the remarks of Section 3 [see (3.7)], we see that
the appropriate projection operator. For our purposes, the major significance of these calculations is that the “l-norm” introduced in Section 5 is the “operator norm” for the symmetric problem. Hence, we have a norm which is well defined on alI spaces Shi. 7.
MULTIGRID
AND EXPERIMENTAL
RESULTS
The results of the previous sections, and Theorem 5.1 in particular, provide exact estimates of the decay of the error (in two norms) in one iteration of a twogrid scheme-in the worst case. Since L,, is again a three-term (diagonally dominant) operator of the form (3.3a)-and given specificaIly by (3.8)-we may apply our multigrid approach inductively as follows: Assume that the n-grid multigrid scheme based on “smoothing” with m applications of the damped Jacobi iteration with parameter a is defined. Suppose h=$
n > 2,
(7.la)
where H is of the form Hz-
1 p+l’
pal.
(7.lb)
We wish to solve (3.4) on the h-grid. The iterative scheme is given by the following inductive description.
(1) OfI the
h-grid (h = 2-“H):
(4 Let U” be chosen. 03)Apply the damped (cl
obtain 0. Form r(h)=f-
Jacobi (with parameter a) iteration m times to
L,,fi.
DAVID KAMOWITZ
174
AND SEYMOUR V. PARTER
(2) Transfer information:
(a) Set r(2h) = Z:%(h) (3) On the 2h grid: (a) Consider the problem &#(2h)
= r(2h).
(b) If 2h = H, solve exactly. (c) If 2h < H, set 6’(2h) = 0 and apply the n-grid iterative scheme (based on smoothing with m applications of the damped Jacobi iteration with parameter a). Let U1(2h) be the result of this step. (4) Transfer infmtion: (a) U1 = 6+ (b) u’+UO.
Z,h,U1(2h).
Return to l(b). In the multigrid jargon this is the so-called slash or sawtooth cycle, which we indicate schematically as follows: h o~nmc&ng.5~~ \
/Oh 0 2h
m “smiwthings” 4hO ‘0
/ o/04h
‘. 2&o
0 H H
/’
/’
Solve
NOTE. There are no smoothing steps during the transfers from coarse to fine grids. McCormick [ 14, 1.51 calls such a multigrid cycle an M \ ,, cycle. When the smoothing occurs only on the way “up” the cycle and the errors are merely restricted on the fin&o-coarse transfers, he calls this an M/,, cycle. For the symmetric case, using Richardson’s iteration (see [14]), he shows that IIM/hlli = IIM\hlll*
(7.2)
In discussing the symmetric M/,, cycle he obtains the following estimate. Let &O= ?I0+ I$#,&
(7.3a)
A Study of Some Multigrid Ideas
175
Suppose LY,0 < (r < 1, satisfies IPoll;
d 41~“ll;
+ ll~2hh~“II:.
v&O-
(7.3b)
Then IW/hlll
d a1’2>
(7.4a)
that is,
llE1lll= JIU- Ullll Q a”211e0111.
(7.4b)
Since our Jacobi iteration is not all that different than Richardson’s iteration, it is not surprising that a similar result holds in our case. Indeed, if one applies his argument to our multigrid cycle, i.e. the M \h cycle, one gets the following result.
LEMMA7.1.
Consider the symmetric case, and suppose II%
Q 1.
(7.5)
Let Ge”=E=q+Z$p.
(7.6a)
Suppose 2, 0 < & < 1, satisfies 11~11;+ 4lGPll~
Q WOIl~*
(7.6b)
z%f?n lIM\,Ji Q 2’2.
(7.7)
PROOF. The proof proceeds by induction. Since we use an exact solver on the coarsest grid,
(IM\,(J,
= 0 < ~3”~.
(7.8a)
Q h”2,
(7.8b)
Assume that llM\2,Jl1
176
DAVID KAMOWITZ
AND SEYMOUR V. PARTER
that is,
IP4\,,42~)- 4w IL =s~“21142~) IL.
(7.8~)
Then
Il~‘llf= 11~11~ + (/Go4\2h@
-
0)
11:.
(7.9)
Because this is a symmetric problem we know that (7.1Oa) and that
IlZ&ullf= (Z&u,L,Z,h,u)= (wz;hLhz,“,u> =
Y(U9
E,,u)
=
Yll~ll~.
(7.1Ob)
Therefore, Il&‘llf = llsll~ + YllM\z+J
- 4lE
By the inductive hypotheses (7.8b) we have ll~lll;
G 111711: + Yw4I;
= llsll;
+ ~llGPllf*
(7.11)
NOTE. In (7.1Ob) and in this calculation the symbols ]]w]]r and ](l&+~]]r refer to the designated norms on the spaces S,,, and S,. By the basic inequality (7.6b) we have ll~‘ll? d Ilsll~ + 4GPll~
which proves the lemma.
G 41~“ll:~ n
Since the proof of this lemma is immediate once one understands the proof of McCormick’s Lemma 2.2 of [15], one would expect that a=
&.
(7.12)
177
A Study of SomeMulti&d Ideas
Indeed, this is the case. Direct but messy calculations based on the results of Section 5 yield a=&
~f(~)zm(l+p)+~(~)zm(l-P,-(~)zm(~)zm\
=
w,n
1-1 (-)P-a 2
1+a
I
2m (l-P)-;(~)2m(I+P)
(7.13) Moreover, for all choices of a and m, the supremum is attained at p=l.The corresponding values of CY”~are displayed in the following table: BOUNDSONTHECONVERGENCEBATE
In view of the results of Section 6 which demonstrate the complete equivalence of our problem to a related symmetric problem, these upper bounds apply in our case. However, the estimate a’/2 = &li2 is only an upper bound for the rate of convergence of the multigrid iterative scheme. In order to complete our investigation we have undertaken an experimental project. A computer program was written with the following capabilities: The user supplies
P(X).b(x), Q(X),f(x), m, (1, n, and M = i - 1 = number of points on the finest grid, where p(x), b(x), g(x), f(x) and
are the coefficients of the problems (l.l),
m = number of applications of the damped Jacobi iteration, a = parameter of the damped Jacobi iteration, n = number of grid levels.
(1.2)
178
DAVID KAMOWITZ
AND SEYMOUR V. PARTER
The user also supplies an initial guess U” and a tolerance E. The program then executes multigrid iterations until the I, norm [see (7.14a)l of the residual is below the indicated tolerance E. The program is run in an interactive fashion which allows the user to change the parameters M, m, a, and n. The experiments reported here were run on the VAX 786 in both singleand double-precision arithmetic (approximately eight and sixteen decimal digits of accuracy respectively). The single-precision results were qualitatively similar to the double-precision results; for increased accuracy the double-precision results are reported here. For our present purposes the basic program was modified to enable us to estimate the “rates of convergence” of the multigrid iteration. For each test problem we used a known solution u(x) of the boundary-value problem (l.l), (1.2). Then we computed the exact solution U(X, h) of the algebraic system (2.4). Then using the two norms (7.14a)
ll4ll= (w uw2,
(7.14b)
we computed the norms of the error u - ui at each iteration. The rate of convergence was measured by computing
-=Ikill IIEi-llj
Pi
(7.15)
at each iteration i = 1,2,3,. . . . To check that the program was working correctly, a number of measures were taken. The most simplistic was to carry out some of the iterations by hand and to compare the hand computed calculations with the iterates generated by the machine. In addition, since the discretization error is 0( h2), it is not unreasonable to expect that halving the step size should reduce the final error in u by a factor of four. This property was checked and found to be true. One of the requirements for Z22:h+lh is that (7.16) (by Lemma 2.1). After each coarsetofine grid transfer, (7.16) was computed and checked. Finally, from (5.7b) one sees that the error, .sZk(h), on the even
179
A Study of Some Multi&d Ideas
points of the coarsest grid should be zero. This requirement was also verified after each iteration. The test problems are best described by giving the choices of p(r), b(x), 9(x), and u(x), the true solution of the differential equation (l.l), (1.2) [which determines f( r )] . As a basic case we took P(X) = I,
9(x)=b(x)=O
and
u(x)=O.
(7.17a)
This test was merely to be sure the program worked on this simple case. Other tests involved six different problems based on two additional sets of coefficients p(x), b(x), 9(x) and three “solutions” u(x). The coefficients used were p(x)=l++sin47rr,
b(z)=l+x,
q(r) = (sin57rx)2, (7.17b)
P(X) = ex,
9(x) = (1-
+)=1+x2,
x)ex/2, (7.17c)
and the “solutions” were
ul(x) = x(e
- e’),
u2(x) = P(l-
x),
u,(r)=sin147rx.
(7.18a) (7.18b) (7.3
In principle, unless one is terribly unlucky in the choice of initial guess, almost any initial guess will lead to the same (and the correct) asymptotic error decay. Indeed, the conventional wisdom asserts that even if the initial guess lacks a component of the worst error vector, roundoff error will introduce such a component. However, in practice the multigrid convergence is so fast that it is essential that the initial error contain a reasonably large component of the worst error vector. Therefore, for each problem test runs were made with a variety of initial guesses. Each initial guess consisted of a smooth component u;=20sinM+l
ka
(where M is the number of points on the finest grid) and a rough component.
180
DAVID KAMOWITZ AND SEYMOUR V. PARTER
The rough component was chosen in various ways in order to have different compositions on the coarser grids. The rough components of the initial guesses are best described schematically by setting Uk = us, +406,,
where IS,1 = 1, and the sign of Sk follows the patterns Initial guess
Pattern for 6,
A
+-+-+-+-+-+-+-+
B C D E
++--++--++--+++++---+++---+++ ++++----++++---I---+++----+++++
(7.19) Runs were made with a, the damped Jacobi parameter, equal to 0.333,0.5,0.667,0.75,1.0,1.333, while m, the number of smoothing iterations, ran from 1 to 4 and the number of grid layers varied from 2 to 5. For each test problem, the program stopped when the discrete 1, norm of the residual vector was less than 0.00005. The most recently computed rate of convergence
Il%nallll, Il%&lll,, was saved and is recorded in Tables 5-8. The theoretical rate for a twogrid iteration scheme was computed from Theorem 5.1 and Theorem 5.2 by solving for the maximum of
and
using Newton’s method. Table 1 exhibits [max F(p)] ‘I2 (a predicted rate of
A Study of Some Multi&
181
Ideas
TABLE 1 PBEDICI’ED
BATE
BASED
F(p)
ON
TABLE 2 DAMPED
JACOBI
PREDICTED
BATE
PABAMETEB
/l
TABLE 3 m
a=.333
1
500
2
.250
3 4
.125 a68
BASED
ON
FI( fi)
.500
.667
.750
1.000
1.333
333 .lll .078 .062
A00 ml .088 .066
A29
.500
.572
.184 ,093 .072
.250 .125 .083
.326 .187 ,109
convergence) as a function of m and u. The value of p at which the maximum of F(p) occurred can be found in Table 2. Table 3 exhibits [max F1(p)] ‘I2 (a better rate of convergence) as a function of m and a. The value p at which the maximum of FI( p) occurred can be found in Table 4. Tables 5 through 8 contain the worst rate of convergence found experimentally as measured in the 2, and 11)I1 norms. The letters in the tables TABLE 4 DAMPED m
1 2 3 4
a=.333
1 1 1 .883
JACOBI
PARhMETER
F
.500
.667
.750
1.000
1.333
0 0 .612 .707
0 0 .577 .667
0 0 xi0 .650
0 0 0 .577
0 0 0 .370
DAVID KAMOWITZ AND SEYMOUR V. PARTER
182
WORST m
1 2 3 4
TABLE 5 CASE, 2 GRIDS
.5
667 .400(b)
.aQ(b)
.5w3)
.if@(b)
JWJ)
.124(c)
.333(b) .111(b) .075(b)
.063(a)
.062(r)
.068(g)
.093(e) .071(h)
.250(b) .125(b)
(1=X33
.499(a) .250(c)
.087(d)
.75
WORST
TABLE 6 CASE, 3 GRIDS
WORST
TABLE 7 CASE, 4 GRIDS
WORST
TABLE 8 CASE, 5 GRIDS
1.0
.082(e)
1.333 .571(b) .327(b)
.187(b) .107(i)
A Study of Some Multigrid Ideas
Problem ; : e f hg i : 1 m n 0
P 9
183
Coefficients
“Solution”
Pattern for (11~in initial guess
(7.17b) (7.17c) (7.17t.I) (7.17c) (7.17b) (7.17c) (7.17a) (7.17c) (7.17c) (7.17a) (7.17b) (7.17c) (7.17a) (7.17c) (7.17b) (7.17c) (7.17b)
r(e - 8) sin 14 rrr r5’2(1 - x) sin 14 rrx r( e - ex) r( e - e*) 0 r( e - e*) sin 14 rrr 0 x( e - e’) r5’2(1 - r) 0 r5’2(1 - x) x(e - ex) r( e - e’) sin 14 rrx
B B C C C D D D E E E A E C A A A
correspond to the choices of coefficients, “solutions,” and patterns for rough components in the initial guess [see (7.17), (7.18), (7.19)] displayed in
Table 9. CONCLUDING
REMARKS
As can be seen from the computational results, no particular choice of problem or initial guess always resulted in the worst case. Moreover, it appears that alI2 = $I2 is an upper bound on the rate of convergence of the multigrid scheme but does not yield the exact rate of convergence. Notice that there seems to be no degradation at all for m = 1. However, as m increases we find some degradation in the rate of convergence, though it appears to be quite less than orI2 = &‘12. REFERENCES 1
Ft. E. Alcouffe, A. Brandt, J. E. Dandy Jr. and J. W. Painter, The multigrid methods for the diffusion equation with strongly discontinuous coefficients, SIAM I. Sci. Statist. Comput. 2439-454 (1981).
184 2
6
7 8 9 10 11 12
13 14 15 16 17
DAVID KAMOWITZ AND SEYMOUR V. PARTER N. S. Bahvalov, On the convergence of a relaxation method with natural constraints on the elliptic operator. Zh. Vychisl. Mat. i Mat. Fiz. 6:861-883 (1966); U.S.S.R. Comput. Math. and Math. Phys. 6:101-135. A. Brandt, Multilevel adaptive solutions to boundary-value problems, Math. camp. 31333-390 (1977). J. E. Dendy, Jr., Black box multigrid, J. Cumput. Phys. 48:366-386 (1982). J. E. Dendy, Jr., Black box multigrid for nonsymmetric problems, Los Alamos Report LA-UR-83-625, presented at International Multigrid Conference, Copper Mountain, Colo., 6-8 Apr. 1983. R. P. Fedorenko, A relaxation method for solving elliptic difference equations, Zh. Vychisl. Mat. i Mat. Fiz. 1:922-927 (1961); U.S.S.R. Comput. Math. and Math. Phys. 1:1092-1096 (1962). P. 0. Frederickson, Fast approzimate inversion of large sparse linear systems, Math. Report 7-75, Lakehead Univ., Ontario, Canada, 1975. A. Greenbaum, Analysis of a multigrid method as an iterative technique for solving linear systems, SIAM 1. Nuw. Anal., to appear. W. Hackbusch, On the convergence of multi-grid iterations, Be&i&e Numer. Math. 9:213-239 (1981). W. Hackbusch, Convergence of multi-grid iterations applied to difference equations, Math. Cump. 34:425-440 (1986). W. Hackbusch, On the multi-grid method applied to difference equations. Ccnnputing 20:291-366 (1978). W. Hackbusch, Introduction to multi-grid methods for the solution of boundary value problems, Seminaire d’Analyse Numerique No. 358, Inst. National Polytechnique de Grenoble, France, 1981. S. McCormick and J. Ruge, Multigrid methods for variational problems, SIAM J. Nuw. Anal. 19:924929 (1982). S. McCormick, Multigrid methods for variational problems: Further results, SIAM J. Numer. Anal. 212552263 (1984). S. McCormick, Multigrid methods for variational problems: General theory, to appear. R. A. Nicolaides, On the Z2 convergence of an algorithm for solving finite-element equation, Math. Camp. 31:892-966 (1977). J. C. South (Jr) and A. Brandt, Application of a multi-level grid method to transonic flow calculations, in Trunsonic Flow ProHems in Turbumachinery (T. C. Adamson and M. F. Platzer, Eds.), Hemisphere, Washington.
Transmitted by S. McCormick
ABSTRACT In an effort to understand certain ideas and concepts associated with multigrid iterations, we give an indepth study of a particular simple problem. We consider a standard finite-difference system associated with the twqoint boundary-value problem
-(pd)‘+bd+qu=o,
u(0) = u(1) = 0
The operators Zth, Zt,, are “operator” based interpolation and projection operators, while the smoothers are the damped Jacobi iterations with parameter a > 0. We determine the exact rates of convergence for the “twc@d” scheme and upper bounds ( < 1) for the multigrid schemes. Experimental results are discussed.
1.
INTRODUCTION
The multigrid approach for the numerical solution of boundary-value problems for elliptic partial differential equations is proving itself as one of the fastest and most efficient methods-see [l], [3], [4], [5], [6], [17]. Moreover, there are a large number of theoretical papers on this subject-see [2], [3], [6], [7], [8], [9], [lo], [ll], [13], [16]. Nevertheless, it seems (at least to these authors) that we are just beginning to understand this powerful idea. In particular, there are the questions: how do we choose the interpolation and *Supportedby the U.S. Air Force Office of Scientific Research under Contract No. AFOSRSZ 0275. APPLIED
MATHEMATICS
AND COMPUTATION
0 Elsevier Science Publishing Co., Inc., 1985 52 Vanderbilt Ave., New York, NY 10017
17:153-184 (1985)
153 -/85/$03.30
DAVID KAMOWI’IZ
154
AND SEYMOUR V. PARTER
projection operators? How do we choose the smoothing operators? What do we mean when we say smoothing? etc. This report is a reflection of our effects to understand and appreciate the theoretical insights of Frederickson [7], McCormick, and Ruge [13], McCormick [14, 151, and Greenbaum [8] and apply those ideas to extend the explicit convergence rates given by Hackbusch [12, (2.21); 111 for the very simplest problem uf’=
u(0) = u( 1) = 0.
f,
Specifically, we consider the twopoint boundary-value problem
Lu := -
(put)‘+ b(x)u’+ qu = f, u(0) = u(l) = 0
O
(1.1)
(1.2)
where p(x), b(x), Q(X) are smooth functions and
P(X)2
PO ’
0,
q(x) > 0.
(1.3)
In Section 2 we describe a basic approach to multigrid which is based on the ideas of Frederickson [7], McCormick and Ruge [ 131, and Greenbaum [8]. In Section 3 we describe a discretization (finite-difference) of the problem (l.l), (1.2) and a specific twogrid iterative procedure for its solution. In Section 4 we describe the class of damped Jacobi “smoothers” and use our knowledge of these schemes and their eigenvectors to describe the basic spaces: Range: Z& = 9 and Nullspace I,2hLh = q. In Section 5 we obtain estimates for the norm decay of a single step in a twogrid scheme for two different norms. In addition we obtain a better estimate for the norm decay for all iterative steps beyond the first in both these norms. Interestingly enough, this estimate is the same in both norms. This latter result is an improvement over the estimates of Hackbusch [ll, 121. It is well known that the problem (l.l), (1.2) is equivalent to a self-adjoint problem. Moreover, the discretization (3.4) is also equivalent to a symmetric problem. In fact, our multigrid treatment of this problem is equivalent to the “same” multigrid treatment of this symmetric problem. This equivalence is not needed for the discussion in Sections l-5. However, as we turn to the extension of a twogrid scheme to a true multigrid scheme, we require this information. In Section 6 we demonstrate this equivalence. In Section 7 we describe the n-grid “sawtooth” multigrid schemes and give a theory (closely
155
A Study of Some Multigrid Ideas
related to a theorem of McCormick [14, 151) which describes the rates of convergence of this scheme. In addition Section 7 contains some experimental results.
2.
A BASIC THEORY
The theory presented in this section is based on the work of Frederickson [7], McCormick and Ruge 1131, and Greenbaum [8]. We consider a finite dimensional linear space S, and a problem
L/p= f,
U,f
ESh,
(24
where Lh:S,-+Sh
is a linear, nonsingular operator. Multigrid is an iterative method for the solution of this problem. The basic idea is to utilize another finite-dimensional space S,, with dim Ss,, < dim S,.
(2.2)
Hence we require operators I,2h, Z& which enable us to effect communication between these spaces. In particular, we have
z,zh:s, + s,,
(projection),
(2.3a)
I,:, : s,, + s,
(interpolation),
(2.3b)
where Zih, Z&, are linear operators. We also require a “smoothing” operator y,, and a “COaW2 grid” Operator &,’ The SImOthing OperatOr y,, iS an affine operator which has U, the unique solution of (2.1), as its only fixed point. That is, (2.4a) where G, : Sh + S, is a linear operator, and if u is the solution of (2.1), then YhU =
u.
(2.4b)
DAVID KAMOWITZ
156
AND SEYMOUR V. PARTER
Finally, the “coarse grid” operator
(2.5) is a linear, nonsingular operator taking Ssh onto itself. Let U” E S, be a guess for the solution U of (2.1). Set (2.6a) (2.6b) (26c)
r=f-L,cl=LL,(U-O)=L,E, R=12hrh
(26d) (26e)
*
Solve
L,,P=
R,
i.e.,
2= L,;;‘R.
(2.7)
Set
(2.8) uo := u’, and return to (2.6a).
REMARK. While one might say that we have merely described a twogrid scheme, the iterative scheme described above does, in fact, describe “multigrid’ schemes. The point is that the “coarse-grid correction” operator L,, may be a complicated procedure involving more grids. The work of Frederickson [7], McCormick and Ruge [13], and Greenbaum [8] suggests we study two fundamental subspaces. These are (2.9a) n = n&pace
of IihLh.
(2.9b)
157
A Study of Some Multigrid Ideas In addition we consider a special two-grid “coarse grid” operator e
2h
= 12hL h Ih2h*
(2.10)
h
This particular operator is the “Gale&n choice” and is optimal in a certain sense. This fact is emphasized by the following
LEMMA 2.1. with
Consti
one iteration
as &scribed
L,, = i,, suppose
EE 9.
above
by (2.6a)-(2.8)
.
(2.11)
That is, suppose there is a w(2h) E s2,, so that E= a(h) = z;hW(2h).
(2.12)
C= w(2h)
(2.13a)
U’=U.
(2.13b)
Then
and
Hence, the problem is solved.
PROOF. From (2.6d) and (2.12) we have r = L,I = L,,Zthw(2h). Thus R=
z,fhr= ( z,2hL,z;h)w(2h).
That is, R = i2hW(2h).
(2.14)
Hence, from (2.7) we have (2.13a). Finally, (2.13b) follows from (2.8) and (2.12). n
158
DAVID KAMOWITZ
AND SEYMOUR V. PARTER
Now suppose we can write S, as the direct sum (not necessarily orthogonal) of Range Z& and Nullspace ZlhL,. That is, every grid function w(h) E S, can be uniquely written as w(h) = Z,h,w’(2h) + w”(h),
(2.15a)
where z,29+?(
h) = 0.
(2.15b)
Let us apply this decomposition to the intermediate error I = q(h). Then E(h) = Z&w’(2h) + w”(h).
(2.16)
We will now show that &1= w2(h). We have r = L,z;~w1(2h)
+ L,w2( h)
and R = Z;% = [ &$!d,z,h,] w’(2h). Using Lemma 2.1 we see that
u’ = 0 + z;,$(2h) and hence
That is, &I= W”(h). Thus, the convergence of the process can be tested by
where
(2.17)
159
A Study of Some Multigrid Ideas
and the norm involved is any norm. The authors mentioned above all dealt with self-adjoint L, and assumed that
zlh = c(Zih)T In this case it is an easy matter to see that Range Zih and NuBspace ZihLh are Lb-orthogonal complements of each other in S,. In the general case we must assume the decomposition (2.15a). However, in Section 4 we will establish the validity of (2.15a) from the following simple criteria: (2.15a) holds if L,, is nonsingular and (i) rank Z,:, = dim Sa,,, (ii) dim Nukpace ZihLh Z dim Sh - dim S,,. This follows from a simple counting argument and the fact that Lemma 2.1 implies that the zero vector is the only vector common to both subspaces. 3.
THE DISCRETE
PROBLEM
Let an integer N > 1 be chosen, set 1
’
h=
2(&7+1)
(3.1)
= (2N+1)+1’
and let the fine “grid points” be given by
xj( h)
j=O,1,2
= jh,
,..., 2N+2.
(3.2)
Define the difference operator Lh by
[Lhu]k=
(3.3a)
-akU,-,+PkUk--YkUk+l~
where
1
b,
pk-l/2
ak=
-
h2
I
pk+l/2
Bk =
Pk+l/2 Yk =
7-s
+2h
1’
+ pk-l/2
+ok
h2 bk I
*
1 9
(3.3b)
Then L, is a consistent approximation to the operator L described in (1.1).
160
DAVID KAMOWITZ
AND SEYMOUR V. PARTER
We assume h is so small that lyk> 0, yk > 0 for all k. We are concerned with the system of linear equations
P%~lk=f,~
k=1,2
,..., 2N+l,
u, = u,,,, = 0.
(3.4)
We shall solve this system with a particular multigrid iterative scheme. Consider a coarse grid on which we have a mesh spacing of 2h and the coarse grid points are given by X,(2h)
= 2hk,
k=0,1,2
,..., N+l.
Let S,bethespaceofgridfunctions{U,(h); k=O,1,...,2N+2}definedon the fine grid, and S,, the space of grid functions { U,(2h); k = 0, 1,. . . , N + l} defined on the coarse grid. Our first step is to construct the interpolation operator I& which maps Ss,, into S,, i.e.,
Following the experience of Dendy [4, 51 we define I!& according to [ Z$,U(2h)],,
= U,(2h)
(common points)
(3.5a)
and (new points).
This choice of “operator” interpolation may be described in the following way: If the physical point x j( h) of the fine grid is also a physical point of the coarse grid, i.e., if j is even, say j = 2k, we set Uj(h)= U,(2h) equal to the same value as the coarse grid function assumed at that point; ‘hat is, (3.5a) holds. If the physical point xl(h) of the fine grid is not a point of the coarse grid, i.e., j is odd, say j = 2k - 1, we require that
~~&i’S@h)]},,_, =o.
(3.5c)
161
A Study of Some Multigrid Ideas
We formalize this remark as
LEMMA 3.1. Let the interpolation operator l& be defined by (3.5a), (3Sb). Then a fin&on U(h) E S, is in the range of I.$ if and only if k=1,2
[L,U(h)],,_,=O, The projection
,...,
w
N+l.
operator Iih: S, + S,, is defined by
REMARK. If b(x) = 0 and the operator L is self-adjoint, then ffk=Yk-l
and we see that z;“=C[z;h]T.
(3.7)
As we have said in Section 2, the relationship (3.7) between Zih and Zlh is the “ variational choice” and is frequently recommended for selfadjoint problems-see [13], [7], [4]. For the purposes of this exposition we take the optimal choice of “coarsegrid correction,” i.e., L 2h :=e
2h
=12hL h
h
Ih
2h’
REMARK. A direct computation shows that
[L,,U@h)],
= - (Y~2)uk-1+fi~2’uk-
yi2)u,+,,
(3.8a)
where a2ka2k- 1 OLk (2, = 12 [ 82k-1
(3.8b) 1 Y2ka2k+
1 (2) _
fi2k-,
Y2k-t2k+l
yk
-
ii
[
p2k+l
1
b2k+l
I
1
(3.&)
(3.8d)
DAVID KAMOWITZ AND SEYMOUR V. PARTER
162
Hence, (3.3a).
L,, is again a diagonally dominant three-term operator of the form
From the theory developed in Section 2 we see that, if we were to choose red-black Gauss-Seidel (indeed, black Gauss-Seidel only) as our smoother, the twogrid scheme would converge in one iteration. Thus we would have an explicit direct solver. However, since we are interested in multigrid, we will choose other smoothing operators. 4.
JACOBI
ITERATIVE
SCHEMES
A stationary iterative method for the solution of (3.4) is described by an operator splitting of the form
L,=M-N. Then, given a first guess Ua, we define successive iterates by the formula MUi+‘=
NUj+
f.
(4.1)
The convergence of this scheme is determined by the eigenvalue problem
XMU=NU.
(4.2)
As is well known, the scheme (4.1) is convergent iff maxIhI< 1. It is easy to verify that if (h, +) is an eigenpair of (4.2), then
L,&= (l-
A)M+
(4.3)
In this section we are concerned with a particularly simple class of such iterative methods, namely the Jacobi methods. We set M=(l+a)B,
a 2 0,
(4.4a)
where B = diag(p,).
(4.4b)
A Study of Some Multigrid Ideas
163
In this case we may rewrite (4.1) as
&B-Q-
ui+l=ui+
L,uj).
When a = 0 we call this scheme the Jacobi method. When a > 0 we call this scheme a damped Jacobi method. In these cases we are able to give a relatively complete discussionof the eigenvalue problem (4.2). We have the following facts: (i) The method is convergent for a > 0. (ii) Let Q = 0. Let (p, Cp)be an eigenpair, i.e.,
Let 6 be defined by (4.6) Then ( - /.L, 4) is also an eigenpair. The eigenvalues p are real and distinct. Furthermore, as h + 0 the { p } fill out the interval [ - 1,11.For completeness we repeat the basic relationshipbetween + and 6. Namely, if
k isodd,
+k = ‘?‘k;
(4.7a)
if
kiseven,
+k =
(433)
-
6,.
Since dim S, = 2N + 1, there is a single eigenvector 4 associated with the eigenvalue p = 0. This eigenvector satisfies i,,
(4.8)
= 0.
For u > 0 there are corresponding eigenpairs (A, +>. (x,6>, are the eigenfunctions described in (4.7) and
h_Cl+a
where 9 and 6
X,a-cl
1+u'
1+u'
The eigenpair (0,Cp) corresponds to an eigenpair (A, (p) with
j+a
lfu’
(4.9)
DAVID KAMOWITZ
164
AND SEYMOUR V. PARTER
We now turn to the determination of the three important subspaces
RangeZth, Nullspace If",
Nullspace ZLhL,,
where the operators Z&, Zih are given by (3.5a), (3.5b), and (3.6).
LEMMA4.1. Let (p, $J), ( - p, 4) be the two eigenpairs described above with p # 0. Let @=
N1+d+(1-d61.
in (ii)
(4.10)
Then Q E Range Z&. Further, since the vectors Cp corresponding to different are linearly independent, this construction provides N elements of Range Zl,,.
pairs (P,+),(--II~~) linearly independent
PROOF. Using (4.3) we have
Lh@=
LO+do - CL)W- Cl- cL)(1+PM4
=(1+P)(w4m41. Thus, the lemma follows from Lemma 3.1, which characterizes Range Zih, n and from (4.7a), (4.7b).
LEMMA4.2. (i) For all choices of a > 0, the vector B? (associated with p = 0) is an element of NullspaceZth and hence (4.3) implies that 6 is an element of Nullspace ZfhL,. (ii) Let (CL,$I), ( - p, 6) with p # 0 be the two eigenpairs described above. Let
‘~=B[(~-cL)~+(~+P)+].
(4.11a)
Then \kE NullspaceZih and (4.3) implies that C#I + 4 E Nullspace ZihL,.
(4.11b)
165
A Study of Some MultigridZdeas
Further, since the vectors \k associated with different pairs (CL,+), ( - p, 6) are linearly independent, we have N linearly independent elements of this nullspace and N linearly independent vectors of Nullspace IthLh. The vectors Bc&,6 provide one more independent vector of each of these subspaces.
PROOF. The result follows from a direct computation using (3.6) and the n defining eigenvalue problem.
COROLLARY.
N&pace
There is a unique decomposition
of sh into RangeI.& and
zfhL,.
PROOF. Since L, is nonsing.ular, we have shown that dim Range Zih > N, dim Nullspace zE”Lh = dim N&pace
Zih > N + 1,
dimSh=2N+1. Thus the corollary follows from the observation that Lemma 3.1 implies that these two subspaces have only the zero vector in common. n
5.
SOME ESTIMATES
Let a > 0. We take as our smoother m applications of the corresponding Jacobi iteration. That is, given U ’ = U’(h), we obtain 6 [as in (2.6b)] from the formula Mui+‘=
NUj+
f,
j=O,l,...,
m-l,
G=Urn*
(5.lb)
First let us consider the special vector 6 with its associated eigenvalue h=L Suppose
(5.la)
l+a’
DAVID KAMOWITZ
166
AND SEYMOUR V. PARTER
Then
L,E = CA”‘B+, and, using Lemma 4.2, we see that ZihL hE -= 0 . Hence, in this case, for any norm,
(5.2)
We now consider two special norms defined on S,. Let w, v E S, be of the form
w =
5 A,c#B~+&+ 5 dj$j,
j=l
v =
f
(5.3a)
j=l
cjgBj +Gj
j=l
+ f ej&j.
(5.3b)
j=l
Define
(w,v)~=
f
A,C,+dc+
j=l
(w,u)~=
z
5 djcj, j=l
AjC,(l-~j)+AC+,~l~jd,(l+pi),
j=l
Ilwlli = (WI w>o, llwllf=(w,~)~=
; j=l
A;(l-pj)+ii2+
5 d2j(l+pj). j=l
A Study of Some Multigrid Z&as
167
REMARK. In the simplest case where b = 0 and a and c are constant, the vectors $j are orthogonal and Ilwll; a
lbllr2,~
Ibaa (LhW7w>= lb&f,. LEMMA5.1.
Suppose e’=tJ-U’=ccp+d;,
(5.4)
and let .9 be defined by the two-grid iteration scheme. Then
(5.5a)
IIM max-=L[(l+(l)2m(l-p)+(~)2m(l+p)]. (5Sb) 2 l+a c,d
ll~“ll~
PROOF. From (5.4) we see that I = CA’?++ dk”& Following the theory of Section 2, we write (see 2.16) (5.6)
where Z;hL,e’(
h) = 0.
We claim that Z&w’(2h)
= “”
; d^xm[(l+
81 = A”( 1 - cl) + dk( 2
/.L)$- (1 - /A)+] 1 + ,u)
[~+~I.
(5.7a)
(5.7b)
DAVID KAMOWITZ
168
AND SEYMOUR V. PARTER
To verify this we need merely verify that the sum of the right-hand sides of (5.7a) and (5.7b) is E, and use (4.10) and (4.11b). Having verified (5.6), (5.7a), (5.7b) we proceed as follows: 11&111~=~[cXm(l-c,)+dlXm(l+CL)]2,
(5.8a)
lleO1l;= c2 + d2,
(5.8b)
ll~‘ll;=&A~(l-~)+d~~(l+~)]~,
(5.&)
ll~‘ll: = c2(1 - Jo)+ d2(1+
(5.8d)
p).
Thus, a simple argument shows that sllp~=~[A~“(l-p)2+P~(l+p)2],
(5.10a)
0
sup~=I[hZ”(1-p)+~2*(I+lr)]. (5.1Ob) 1
Using the basic formulae (4.9) we obtain
Il4G sup ,l$l,;
-=q(q2m(l-p)2+(~)2m(l+p)2], 2 1+a
Il~‘ll: sup 11~011:
-=+[(~)2m(l-~)+(~)2m(l+~)].
Thus, the lemma is proven.
THEOREM
5.1.
In the general case we have
(5.11a)
(5.11b)
169
A Study of Some Multigrid Ideas
Moreover, there is a choice of e” which makes these estimates sharp (with an error of order h).
PROOF. The theorem follows immediately from the previous lemma.
n
We observe that (5.11a) is precisely the formula obtained by Hackbusch [12, (2.21)] in the special case p(x) = 1, b(x) = q(x) = 0, and a = 1. To make a complete identification we merely set
1-p
u=2’
1_,=+.
(5.12)
However, while (S.lla), (5,llb) describe the worst-case decay in one multigrid iteration in a twegrid scheme, it does not give the estimate of real interest. From the discussion in the proof of Lemma 5.1-and (5.7) in particular-we realize that, even though the constants c and d of (5.4) may be arbitrary for E’, that is not true for Ed, k > 1. We have [from (5.7b)]
El=.[++~]
(5.13)
with (J=
~h~(l--~)+d~~(l+~) 2
Therefore, following the argument of Lemma 5.1,
Hence [l&q; = Il&(‘)ll~ = 2a2 , ll&‘2’~~~=ll&‘2’11~=~[~~(1-p)+X”(l+~)]2, and for j = 0,l we have
~=i[h-(l-p)ix_(l+p)l”. I
(5.14)
DAVID KAMOWITZ
170
AND SEYMOUR V. PARTER
Thus we have proven
THEOREM 5.2.
forj = 0,l and all k > 1 we have
In the general case,
II&(k+q;
=
lle’k’ll;
sup
(5.15)
{+[Xm(l-/&)+1~(l+r_l)]“}.
-l
REMARK. The distinction between Hackbusch’s result (5.11a) and (5.15) is nontrivial for large m. We have, as m + co,
ll4lO -_--
l+a
1
11~0110fi
(5.16a)
fT7n’
while
II&(k+‘)llj _-- It-a llEkllj
2
1
(5.16b)
m *
Thus for k >, 1 we have
(5.17)
6.
SYMMETRIZATION Consider the difference equation (3.4) described by (3.3a), (3.3b). Let (6.1)
& = dkvk, where the coefficients
d, are computed recursively by (6.2a)
do = 1, ^.
d k+l- -d
Uk+l
k
k=O,l,...
d
-fk
’
.
(6.2b)
171
A Study of SOme Multigrid ~deos
Then (3.4) becomes
or
(6.3)
that is, (6.3a) where (6.4a)
(6.4b)
Hence iik = ?k- 1,
(6.5)
so that E, is symmetric. We now turn our attention to the Jacobi iterative schemes of Section 4. We have ui+l=ui+
&rl(
f
- L,uj).
(6.6)
The change of variables (6.1) is conveniently described by U=DV
(6.7a)
where D=diag(d,,...,d,,+,).
(6-W
DAVID KAMOWITZ
172
AND SEYMOUR V. PARTER
With this change of variables the iteration (6.6) becomes Vi+r=Vi+
(f-
L,DVj).
But, since D and B are both diagonal matrices, we have vi+l=vi+
&3-‘(D-‘f-
E,vj),
(6.8a)
where
D-‘f=(~,~,...,fZiN+l)‘.
(6.8b)
Since B is also the diagonal of i,, (6.8) is precisely the same Jacobi iterative scheme for modified symmetric equations. Now let us study the effect of this transformation on the twogrid iterative scheme. We compute
Imagine U(2h), = d,,V(2h),
is given in S,,. Then, from (35a),
Gf [w@w]2k = Vcwk~ d,:-z[ZzhhU(2h)],,-z=V(2h)k-l.
(6.9a) (6.9b)
Thus (3.5b) yields
d,;-,[ZkWh)],,_, = ~[a2,_,d,:_,d,,_,V(2~)k-,+y,t-,d~~-,d,,V(2~)k] 2k 1 =
jj+[h2k-lv(2h)k-l+ 2k 1
?2k-1vt2h)k]*
Thus, with this change of variables the mapping Z$, of our original unsymmetric problem becomes i&, the appropriate mapping associated with the new symmetric problem. Finally, let us consider
A Study of Some Multigrid Ideas
173
A straightforward calculation verifies that
TV,,-,+V,,I,,+~V(h),,,, Zk+l
. I
Thus, following the remarks of Section 3 [see (3.7)], we see that
the appropriate projection operator. For our purposes, the major significance of these calculations is that the “l-norm” introduced in Section 5 is the “operator norm” for the symmetric problem. Hence, we have a norm which is well defined on alI spaces Shi. 7.
MULTIGRID
AND EXPERIMENTAL
RESULTS
The results of the previous sections, and Theorem 5.1 in particular, provide exact estimates of the decay of the error (in two norms) in one iteration of a twogrid scheme-in the worst case. Since L,, is again a three-term (diagonally dominant) operator of the form (3.3a)-and given specificaIly by (3.8)-we may apply our multigrid approach inductively as follows: Assume that the n-grid multigrid scheme based on “smoothing” with m applications of the damped Jacobi iteration with parameter a is defined. Suppose h=$
n > 2,
(7.la)
where H is of the form Hz-
1 p+l’
pal.
(7.lb)
We wish to solve (3.4) on the h-grid. The iterative scheme is given by the following inductive description.
(1) OfI the
h-grid (h = 2-“H):
(4 Let U” be chosen. 03)Apply the damped (cl
obtain 0. Form r(h)=f-
Jacobi (with parameter a) iteration m times to
L,,fi.
DAVID KAMOWITZ
174
AND SEYMOUR V. PARTER
(2) Transfer information:
(a) Set r(2h) = Z:%(h) (3) On the 2h grid: (a) Consider the problem &#(2h)
= r(2h).
(b) If 2h = H, solve exactly. (c) If 2h < H, set 6’(2h) = 0 and apply the n-grid iterative scheme (based on smoothing with m applications of the damped Jacobi iteration with parameter a). Let U1(2h) be the result of this step. (4) Transfer infmtion: (a) U1 = 6+ (b) u’+UO.
Z,h,U1(2h).
Return to l(b). In the multigrid jargon this is the so-called slash or sawtooth cycle, which we indicate schematically as follows: h o~nmc&ng.5~~ \
/Oh 0 2h
m “smiwthings” 4hO ‘0
/ o/04h
‘. 2&o
0 H H
/’
/’
Solve
NOTE. There are no smoothing steps during the transfers from coarse to fine grids. McCormick [ 14, 1.51 calls such a multigrid cycle an M \ ,, cycle. When the smoothing occurs only on the way “up” the cycle and the errors are merely restricted on the fin&o-coarse transfers, he calls this an M/,, cycle. For the symmetric case, using Richardson’s iteration (see [14]), he shows that IIM/hlli = IIM\hlll*
(7.2)
In discussing the symmetric M/,, cycle he obtains the following estimate. Let &O= ?I0+ I$#,&
(7.3a)
A Study of Some Multigrid Ideas
175
Suppose LY,0 < (r < 1, satisfies IPoll;
d 41~“ll;
+ ll~2hh~“II:.
v&O-
(7.3b)
Then IW/hlll
d a1’2>
(7.4a)
that is,
llE1lll= JIU- Ullll Q a”211e0111.
(7.4b)
Since our Jacobi iteration is not all that different than Richardson’s iteration, it is not surprising that a similar result holds in our case. Indeed, if one applies his argument to our multigrid cycle, i.e. the M \h cycle, one gets the following result.
LEMMA7.1.
Consider the symmetric case, and suppose II%
Q 1.
(7.5)
Let Ge”=E=q+Z$p.
(7.6a)
Suppose 2, 0 < & < 1, satisfies 11~11;+ 4lGPll~
Q WOIl~*
(7.6b)
z%f?n lIM\,Ji Q 2’2.
(7.7)
PROOF. The proof proceeds by induction. Since we use an exact solver on the coarsest grid,
(IM\,(J,
= 0 < ~3”~.
(7.8a)
Q h”2,
(7.8b)
Assume that llM\2,Jl1
176
DAVID KAMOWITZ
AND SEYMOUR V. PARTER
that is,
IP4\,,42~)- 4w IL =s~“21142~) IL.
(7.8~)
Then
Il~‘llf= 11~11~ + (/Go4\2h@
-
0)
11:.
(7.9)
Because this is a symmetric problem we know that (7.1Oa) and that
IlZ&ullf= (Z&u,L,Z,h,u)= (wz;hLhz,“,u> =
Y(U9
E,,u)
=
Yll~ll~.
(7.1Ob)
Therefore, Il&‘llf = llsll~ + YllM\z+J
- 4lE
By the inductive hypotheses (7.8b) we have ll~lll;
G 111711: + Yw4I;
= llsll;
+ ~llGPllf*
(7.11)
NOTE. In (7.1Ob) and in this calculation the symbols ]]w]]r and ](l&+~]]r refer to the designated norms on the spaces S,,, and S,. By the basic inequality (7.6b) we have ll~‘ll? d Ilsll~ + 4GPll~
which proves the lemma.
G 41~“ll:~ n
Since the proof of this lemma is immediate once one understands the proof of McCormick’s Lemma 2.2 of [15], one would expect that a=
&.
(7.12)
177
A Study of SomeMulti&d Ideas
Indeed, this is the case. Direct but messy calculations based on the results of Section 5 yield a=&
~f(~)zm(l+p)+~(~)zm(l-P,-(~)zm(~)zm\
=
w,n
1-1 (-)P-a 2
1+a
I
2m (l-P)-;(~)2m(I+P)
(7.13) Moreover, for all choices of a and m, the supremum is attained at p=l.The corresponding values of CY”~are displayed in the following table: BOUNDSONTHECONVERGENCEBATE
In view of the results of Section 6 which demonstrate the complete equivalence of our problem to a related symmetric problem, these upper bounds apply in our case. However, the estimate a’/2 = &li2 is only an upper bound for the rate of convergence of the multigrid iterative scheme. In order to complete our investigation we have undertaken an experimental project. A computer program was written with the following capabilities: The user supplies
P(X).b(x), Q(X),f(x), m, (1, n, and M = i - 1 = number of points on the finest grid, where p(x), b(x), g(x), f(x) and
are the coefficients of the problems (l.l),
m = number of applications of the damped Jacobi iteration, a = parameter of the damped Jacobi iteration, n = number of grid levels.
(1.2)
178
DAVID KAMOWITZ
AND SEYMOUR V. PARTER
The user also supplies an initial guess U” and a tolerance E. The program then executes multigrid iterations until the I, norm [see (7.14a)l of the residual is below the indicated tolerance E. The program is run in an interactive fashion which allows the user to change the parameters M, m, a, and n. The experiments reported here were run on the VAX 786 in both singleand double-precision arithmetic (approximately eight and sixteen decimal digits of accuracy respectively). The single-precision results were qualitatively similar to the double-precision results; for increased accuracy the double-precision results are reported here. For our present purposes the basic program was modified to enable us to estimate the “rates of convergence” of the multigrid iteration. For each test problem we used a known solution u(x) of the boundary-value problem (l.l), (1.2). Then we computed the exact solution U(X, h) of the algebraic system (2.4). Then using the two norms (7.14a)
ll4ll= (w uw2,
(7.14b)
we computed the norms of the error u - ui at each iteration. The rate of convergence was measured by computing
-=Ikill IIEi-llj
Pi
(7.15)
at each iteration i = 1,2,3,. . . . To check that the program was working correctly, a number of measures were taken. The most simplistic was to carry out some of the iterations by hand and to compare the hand computed calculations with the iterates generated by the machine. In addition, since the discretization error is 0( h2), it is not unreasonable to expect that halving the step size should reduce the final error in u by a factor of four. This property was checked and found to be true. One of the requirements for Z22:h+lh is that (7.16) (by Lemma 2.1). After each coarsetofine grid transfer, (7.16) was computed and checked. Finally, from (5.7b) one sees that the error, .sZk(h), on the even
179
A Study of Some Multi&d Ideas
points of the coarsest grid should be zero. This requirement was also verified after each iteration. The test problems are best described by giving the choices of p(r), b(x), 9(x), and u(x), the true solution of the differential equation (l.l), (1.2) [which determines f( r )] . As a basic case we took P(X) = I,
9(x)=b(x)=O
and
u(x)=O.
(7.17a)
This test was merely to be sure the program worked on this simple case. Other tests involved six different problems based on two additional sets of coefficients p(x), b(x), 9(x) and three “solutions” u(x). The coefficients used were p(x)=l++sin47rr,
b(z)=l+x,
q(r) = (sin57rx)2, (7.17b)
P(X) = ex,
9(x) = (1-
+)=1+x2,
x)ex/2, (7.17c)
and the “solutions” were
ul(x) = x(e
- e’),
u2(x) = P(l-
x),
u,(r)=sin147rx.
(7.18a) (7.18b) (7.3
In principle, unless one is terribly unlucky in the choice of initial guess, almost any initial guess will lead to the same (and the correct) asymptotic error decay. Indeed, the conventional wisdom asserts that even if the initial guess lacks a component of the worst error vector, roundoff error will introduce such a component. However, in practice the multigrid convergence is so fast that it is essential that the initial error contain a reasonably large component of the worst error vector. Therefore, for each problem test runs were made with a variety of initial guesses. Each initial guess consisted of a smooth component u;=20sinM+l
ka
(where M is the number of points on the finest grid) and a rough component.
180
DAVID KAMOWITZ AND SEYMOUR V. PARTER
The rough component was chosen in various ways in order to have different compositions on the coarser grids. The rough components of the initial guesses are best described schematically by setting Uk = us, +406,,
where IS,1 = 1, and the sign of Sk follows the patterns Initial guess
Pattern for 6,
A
+-+-+-+-+-+-+-+
B C D E
++--++--++--+++++---+++---+++ ++++----++++---I---+++----+++++
(7.19) Runs were made with a, the damped Jacobi parameter, equal to 0.333,0.5,0.667,0.75,1.0,1.333, while m, the number of smoothing iterations, ran from 1 to 4 and the number of grid layers varied from 2 to 5. For each test problem, the program stopped when the discrete 1, norm of the residual vector was less than 0.00005. The most recently computed rate of convergence
Il%nallll, Il%&lll,, was saved and is recorded in Tables 5-8. The theoretical rate for a twogrid iteration scheme was computed from Theorem 5.1 and Theorem 5.2 by solving for the maximum of
and
using Newton’s method. Table 1 exhibits [max F(p)] ‘I2 (a predicted rate of
A Study of Some Multi&
181
Ideas
TABLE 1 PBEDICI’ED
BATE
BASED
F(p)
ON
TABLE 2 DAMPED
JACOBI
PREDICTED
BATE
PABAMETEB
/l
TABLE 3 m
a=.333
1
500
2
.250
3 4
.125 a68
BASED
ON
FI( fi)
.500
.667
.750
1.000
1.333
333 .lll .078 .062
A00 ml .088 .066
A29
.500
.572
.184 ,093 .072
.250 .125 .083
.326 .187 ,109
convergence) as a function of m and u. The value of p at which the maximum of F(p) occurred can be found in Table 2. Table 3 exhibits [max F1(p)] ‘I2 (a better rate of convergence) as a function of m and a. The value p at which the maximum of FI( p) occurred can be found in Table 4. Tables 5 through 8 contain the worst rate of convergence found experimentally as measured in the 2, and 11)I1 norms. The letters in the tables TABLE 4 DAMPED m
1 2 3 4
a=.333
1 1 1 .883
JACOBI
PARhMETER
F
.500
.667
.750
1.000
1.333
0 0 .612 .707
0 0 .577 .667
0 0 xi0 .650
0 0 0 .577
0 0 0 .370
DAVID KAMOWITZ AND SEYMOUR V. PARTER
182
WORST m
1 2 3 4
TABLE 5 CASE, 2 GRIDS
.5
667 .400(b)
.aQ(b)
.5w3)
.if@(b)
JWJ)
.124(c)
.333(b) .111(b) .075(b)
.063(a)
.062(r)
.068(g)
.093(e) .071(h)
.250(b) .125(b)
(1=X33
.499(a) .250(c)
.087(d)
.75
WORST
TABLE 6 CASE, 3 GRIDS
WORST
TABLE 7 CASE, 4 GRIDS
WORST
TABLE 8 CASE, 5 GRIDS
1.0
.082(e)
1.333 .571(b) .327(b)
.187(b) .107(i)
A Study of Some Multigrid Ideas
Problem ; : e f hg i : 1 m n 0
P 9
183
Coefficients
“Solution”
Pattern for (11~in initial guess
(7.17b) (7.17c) (7.17t.I) (7.17c) (7.17b) (7.17c) (7.17a) (7.17c) (7.17c) (7.17a) (7.17b) (7.17c) (7.17a) (7.17c) (7.17b) (7.17c) (7.17b)
r(e - 8) sin 14 rrr r5’2(1 - x) sin 14 rrx r( e - ex) r( e - e*) 0 r( e - e*) sin 14 rrr 0 x( e - e’) r5’2(1 - r) 0 r5’2(1 - x) x(e - ex) r( e - e’) sin 14 rrx
B B C C C D D D E E E A E C A A A
correspond to the choices of coefficients, “solutions,” and patterns for rough components in the initial guess [see (7.17), (7.18), (7.19)] displayed in
Table 9. CONCLUDING
REMARKS
As can be seen from the computational results, no particular choice of problem or initial guess always resulted in the worst case. Moreover, it appears that alI2 = $I2 is an upper bound on the rate of convergence of the multigrid scheme but does not yield the exact rate of convergence. Notice that there seems to be no degradation at all for m = 1. However, as m increases we find some degradation in the rate of convergence, though it appears to be quite less than orI2 = &‘12. REFERENCES 1
Ft. E. Alcouffe, A. Brandt, J. E. Dandy Jr. and J. W. Painter, The multigrid methods for the diffusion equation with strongly discontinuous coefficients, SIAM I. Sci. Statist. Comput. 2439-454 (1981).
184 2
6
7 8 9 10 11 12
13 14 15 16 17
DAVID KAMOWITZ AND SEYMOUR V. PARTER N. S. Bahvalov, On the convergence of a relaxation method with natural constraints on the elliptic operator. Zh. Vychisl. Mat. i Mat. Fiz. 6:861-883 (1966); U.S.S.R. Comput. Math. and Math. Phys. 6:101-135. A. Brandt, Multilevel adaptive solutions to boundary-value problems, Math. camp. 31333-390 (1977). J. E. Dendy, Jr., Black box multigrid, J. Cumput. Phys. 48:366-386 (1982). J. E. Dendy, Jr., Black box multigrid for nonsymmetric problems, Los Alamos Report LA-UR-83-625, presented at International Multigrid Conference, Copper Mountain, Colo., 6-8 Apr. 1983. R. P. Fedorenko, A relaxation method for solving elliptic difference equations, Zh. Vychisl. Mat. i Mat. Fiz. 1:922-927 (1961); U.S.S.R. Comput. Math. and Math. Phys. 1:1092-1096 (1962). P. 0. Frederickson, Fast approzimate inversion of large sparse linear systems, Math. Report 7-75, Lakehead Univ., Ontario, Canada, 1975. A. Greenbaum, Analysis of a multigrid method as an iterative technique for solving linear systems, SIAM 1. Nuw. Anal., to appear. W. Hackbusch, On the convergence of multi-grid iterations, Be&i&e Numer. Math. 9:213-239 (1981). W. Hackbusch, Convergence of multi-grid iterations applied to difference equations, Math. Cump. 34:425-440 (1986). W. Hackbusch, On the multi-grid method applied to difference equations. Ccnnputing 20:291-366 (1978). W. Hackbusch, Introduction to multi-grid methods for the solution of boundary value problems, Seminaire d’Analyse Numerique No. 358, Inst. National Polytechnique de Grenoble, France, 1981. S. McCormick and J. Ruge, Multigrid methods for variational problems, SIAM J. Nuw. Anal. 19:924929 (1982). S. McCormick, Multigrid methods for variational problems: Further results, SIAM J. Numer. Anal. 212552263 (1984). S. McCormick, Multigrid methods for variational problems: General theory, to appear. R. A. Nicolaides, On the Z2 convergence of an algorithm for solving finite-element equation, Math. Camp. 31:892-966 (1977). J. C. South (Jr) and A. Brandt, Application of a multi-level grid method to transonic flow calculations, in Trunsonic Flow ProHems in Turbumachinery (T. C. Adamson and M. F. Platzer, Eds.), Hemisphere, Washington.