Convergence of a Multigrid Method for Elliptic Equations with Highly
Convergence of a Multigrid Method for Elliptic Equations with Highly Oscillatory Coecients Bjorn Engquist
Erding Luo y
Abstract
Standard multigrid methods are not so e ective for equations with highly oscillatory coecients. New coarse grid operators based on homogenized operators are introduced to restore the fast convergence rate of multigrid methods. Finite di erence approximations are used for the discretization of the equations. The analysis of convergence is based on homogenization theory. Proofs are presented for a two level method with the homogenized coarse grid operator to solve two classes of elliptic equations with Dirichlet boundary conditions.
Key Words. elliptic equation, oscillation, nite di erence, multigrid method, homog-
enization theory, convergence
AMS subject classi cations. 65N06, 65N12, 65N55
Department of Mathematics, University of California at Los Angeles, LA, CA90024. Research supported by ARPA/ONR N00014-92-J-1890 and NSF DMS91-03104. y IMA, University of Minnesota, Mpls, MN 55455. Research supported by NSF through IMA.
1
1 Introduction Consider the multigrid method arising from the nite di erence approximations to elliptic equations with highly oscillatory coecients of the following type (1.1)
L u(x) =
X @ a (x) @ u (x) = f (x); @xi ij @xj i;j
where aij (x) = aij (x; x ); aij (x; y) = aji(x; y), strictly positive, continuous and 1-periodic in each component of y. Also, the operator L is uniformly elliptic. That is, there exist two positive constants q and Q independent of , such that (1.2)
0k02 ij M1 ij kh 1 1 1 1 2 2 2 kL;h2 (I L;h L;h1 L;h )khkL;h (1 L;h ) (i2+j2 >k02 ij ij )kh : 2 2 Since kL;h2 kh C and kI L;h L;h1 L;h kh C , 1
1
1
r
(2.22)
2 I2 C i2+j2 >k02 ij2 kL;h (1 L;h ) ij k2h 1 2 (1 C k2 max j ) j: 1
N
Assume 1 + 2 > k102 =
0
N
1N 1
k02
1N 1
( C hk022 ), we have
(2.23) I2 Ck0(1 k02) : For I2 12 , it is sucient to have 2 (2.24)
C h lnk ( 1 ): Combining (2.21) and (2.24), we have
0
k02
Ch
1
The proof is completed. Now, we are ready to show the main result. 11
k02
1 3
ln h:
2
Theorem 2.2 There exists constant C , the operator M de ned by (2.5) satis es (M ) 0 < 1; whenever h belongs to the set S (; h0) of Diophantine number, and
Ch
1 1=3 lnh:
Before we carry out the proof, we need the following lemma. Let
M2 = (L;h1 IHh LH1 IhH )L;hS : Lemma 2.5 For some constant C , kM2kh C : Proof. Since L;h and LH = L;H are the homogenized operators de ned respectively on
ne and coarser grid with constant coecients, they are well behaved [11, 12]. Furthermore, (2.25)
kL;h1 IHh LH1IhH kh Ch2:
Therefore,
(2.26)
kM2kh
kL;h1 IHh LH1IhH kh kL;hS kh Ch2h 2= C:
Proof of Theorem 2.2. Note that
2
M = M1 + M2: Therefore,
kM kh kM1kh + kM2kh:
Since (M ) kM kh , by Theorem 2.1 and Lemma 2.5, the rest of the proof can easily follow. 2 12
3 Oscillation Along the Diagonal Direction 3.1 Model Equation
Consider another special case of (1.1), a two-dimensional elliptic problem with coecients oscillatory diagonally,
@ a (x y) @ @ a (x y) @ = f (x; y); (x; y) 2 ; @x @x @y @y (x; y) = 0; (x; y) 2 @ ; where a(x) is a strictly positive continuous function, and a(x) = a(x=) = a(x= + 1). (3.1)
Also, the operator has the property of (1.2). From (1.3), the corresponding homogenized equation of (3.1) is: a + @ 2 + (a ) @ 2 a + @ 2 = f (x; y); (x; y) 2 ; (3.2) 2 @x2 @x@y 2 @y2 (x; y) = 0; (x; y) 2 @ ; R R where a = 01 a(x)dx, and = ( 01 1=a(x)dx) 1 . As goes to zero, we know that the solution of (3.1) converges to the solution of (3.2). Now, consider a corresponding discretized equation of (3.1), (3.3) D+i aij Di uhij D+j bij Dj uhij = fijh ; (i; j ) 2 h where aij = a( xi h= 2 yj ); bij = a( xi yj + h=2 ); (i; j ) 2 h : Here, we assume the discretized coecients have the following property, aj0 = b0N j+1; j = 0; ; N: Denote the discretization of the homogenized operator a + @ 2 + (a ) @ 2 a + @ 2 2 @x2 @x@y 2 @y2 in (3.2) by (3.4) L;h = h +2 ah (D+i Di + D+j Dj ) + h 2 ah (D+i D+j + Di Dj ); 13
where h = (h PNk=1 a10 ) 1, and ah = h PNk=1 ak0. The operator of the two level method can be expressed as k
M = (I IHh LH1 IhH L;h)S = (I IHh LH1IhH L;h )(I L;h ) ; where is the inverse of the largest eigenvalue of L;h, has order of h 2. And LH = L;H . (3.5)
3.2 Convergence Analysis
We still rst consider the simpli ed operator M1 de ned in (2.6). Theorem 3.1 If the ratio of h to is strictly irrational, h belongs to the set of Diophantine number S (; h0) de ned in (1.6), then there exist two constants C and 0 such that
kM1kh 0 < 1
(3.6) whenever Ch
1 2=3ln(h).
In order to prove Theorem 3.1, we introduce the following lemmate rst. Lemma 3.1 De ne two discrete functions on h , (3.7) (3.8)
Xi Xj Zij1 = 2h (i j h a1 + h b1 ); (i; j ) 2 h; k=1 kj i 1 h k=1 k0
X Zij2 = 2h (j i + a
k=1 0k j 1) h k=1 ik
X
b ; (i; j ) 2 h :
Then, Zij1 ; Zij2 are bounded, and satisfy
(3.9) for (i; j ) 2 h .
L;h Zij1 = 1 D+i aij ; L;h Zij2 = 1 D+j bij ;
Proof. Notice that by the assumption of the coecients, N h X
X h = 1; (i; j ) 2 : = h h k=1 akj k=1 bik N
14
Applying the operator L;h to Zij1 , (3.9) then follows. Now, we write Zij1 as follows, j 1 i 1 X X h h ) 2 (j h b ): ij = 2 (i h k=1 akj k=1 0k
Z1 By Lemma 1.1,
j 1 h (i Xi 1 ) = O(1); h (j X ) = O(1): h h 2 2 k=1 akj k=1 b0k Hence, Zij1 is bounded for (i; j ) 2 h , and satis es aij (1 Di Zij1 ) = aij +2 h ; bij Dj Zij1 = bij 2 h :
The result can be deduced similarly for Zij2 . We can also show that Zij2 = Zij1 , and satis es bij (1 Dj Zij2 ) = bij +2 h ; aij Di Zij2 = aij 2 h : This proves the lemma. 2 1 2 Remark. The explicit forms of Zij and Zij depend on a (x y) being a function of x y . For the general angular dependences a( x + y), these forms would not be possible.
Lemma 3.2 Assume discrete function ij de ned as
Xi Xj b ); (i; j ) 2 : ij = h ( akj iah) + h (jah 0k h k=1
k=1
Then, ij is bounded, and has the following properties
Di ij = aij ah; Dj ij = ah bij ; (i; j ) 2 h: Proof. Using the symmetric properties of the coecients, we can establish the results similarly as in Lemma 3.1. 2
Lemma 3.3 Assume Uij satis es (3.10)
L;h Uij = hij ; (i; j ) 2 h ; Uij = 0; (i; j ) 2 @ h; 15
where (; hij ) is a normalized eigenpair of L;h . Then, we have
X [(Di U h)2 + (Dj U h)2] = O( ); + ij + ij
N 1
kU k21 =
(3.11) and
i;j =0
X ((Di Di U
N 1
i;j =1
+
2 + (D i
ij )
Dj Uij )2 + (D+i D+j Uij )2 + (D+j Dj Uij )2)h2
(3.12) = O(2 ): Proof. First, we observe that
X [(Di h h)2 + (Dj h h)2] = O( ):
N 1
+ ij
+ ij
i;j =0
Applying Uij to the following equation and taking summation, we have
X L U U = NX1 L U : ;h ij ij ;h ij ij
N 1
i;j =1
i;j =1
(3.11) then follows. For (3.12), since Uij vanishes at the boundary, it satis es 1 NX1 [(Di Di U )2 + (Dj Dj U )2] NX1 (Di Di U )(Dj Dj U ) ij ij ij ij + + + 2 i;j=1 + i;j =1 (3.13)
=
N X (Di Dj U
i;j =1
ij
)2
=
X (Di Dj U )2: + + ij
N 1
i;j =0
Multiplying D+i Di Uij to (3.10), we get
(D+i Di Uij )2 + 2 aah + h (D+i D+j Uij )(D+i Di Uij ) + (D+i D+j Uij )2 + h
h
(D+i Di Uij )2 + 2 aah + h (Di Dj Uij )(D+i Di Uij ) + (Di Dj Uij )2 + h h j j i i ((D+ D Uij )(D+ D Uij ) (D+i D+j Uij )2) + ((D+j Dj Uij )(D+i Di Uij ) (Di Dj Uij )2) = a +2 hij D+i Di Uij ; (i; j ) 2 h : h
h
16
By the uniformly ellipticity property of the homogenized operator (3.4), we have (3.14)
X ((Di Di U
N 1
i;j =1
2 + (Di Dj U )2 + (Di + + ij
Dj Uij )2)h2 = O(2 ):
2 + (Di Dj U )2 + (Di + + ij
Dj Uij )2)h2 = O(2 ):
ij )
+
An similar argument gives (3.15)
X ((Dj Dj U
N 1
i;j =1
ij )
+
(3.13), (3.14) and (3.15) complete the proof of lemma. 2 Lemma 3.4 Assume that Uij de ned in Lemma 3.3 satis es the following boundary conditions, (3.16) D+i Uij + Di Uij = D+j Uij + Dj Uij = 0; (i; j ) 2 @ h: Then,
kU k22 = (3.17)
+
N X ((Di Di U
i;j =0 N 1
+
2 + (Dj Dj U )2 )h2 ij +
ij )
N X (Di Dj U h)2 + X (Di Dj Uij h)2 = O(2 =h): ij + +
i;j =0
i;j =1
Proof. Since
Di UNj = Di Uij + h we have Then,
(Di UNj )2 (Di Uij )2 + N X (Di U
j =1
2h
Nj )
k=i
X (Di Di U )2h: ij +
N 1 i=1
N N X X (Di Uij h)2 + (D+i Di Uij h)2 = O(2 ):
i;j =1
Combined (3.12) with the following relation, the rest of the proof follows.
X Di Di U ; kj +
N 1
i;j =1
D+i Di UNj = h2 D+i UN 1j ; 17
2
Theorem 3.2 Assume hij ; Uij are de ned in Lemma 3.3 and Lemma 3.4. Then, (3.18) k h U kh = O(p ): Proof. First, we introduce the following discrete functions
Gij = G1ij + G2ij ; G1ij = 21 Uij 2 (Zij1 D+i Uij + Zij2 D+j Uij ) 12 hij ; G2ij = 12 Uij 2 (Zij1 Di Uij + Zij2 Dj Uij ) 12 hij ;
for (i; j ) 2 h. By the assumption (3.16), Gij vanishes at boundary, i.e.,
Gij = 0; (i; j ) 2 @ h : For G1ij , we have 2Di G1ij = Di Uij (Di (Zij1 D+i Uij ) + Di (Zij2 D+j Uij )) Di hij = (Di Uij Di Zij1 D+i Ui 1j ) Zij1 Di D+i Uij
Di Zij2 D+j Uij Zi2 1j Di D+j Uij Di hij ;
2Dj G1ij = Dj Uij (Dj (Zij1 D+i Uij ) + Dj (Zij2 D+j Uij )) Dj hij = (Dj Uij Dj Zij2 D+j Uij 1 ) Zij2 Dj D+j Uij
Dj Zij1 D+i Uij Zij1 1Dj D+i Uij Dj hij :
Then,
(3.19)
2L;h G1ij = D+i aij +2 h Di Uij + D+i aij 2 h D+j Uij D+i (aij Zij1 D+i Di Uij + aij Zi2 1j Di D+j Uij ) +D+i bij 2 h D+j Uij + D+j bij +2 h Dj Uij D+j (bij Zij2 Dj D+j Uij + bij Zij1 1 Dj D+i Uij ) +L hij : 18
For G2ij , we establish similarly 2D+i G2ij = D+i Uij (D+i (Zij1 Di Uij ) + D+i (Zij2 Dj Uij )) D+i hij = (D+i Uij D+i Zij1 Di Ui+1j ) Zij1 Di D+i Uij
D+i Zij2 Dj Uij Zi2+1j D+i Dj Uij D+i hij ;
2D+j G2ij = D+j Uij (D+j (Zij1 Di Uij ) + D+j (Zij2 Dj Uij )) D+j hij = (D+j Uij D+j Zij2 Dj Uij+1 ) Zij2 D+j Dj Uij
D+j Zij1 Di Uij Zij1 +1D+j Di Uij D+j hij :
Then,
(3.20)
2L;h G2ij = Di ai+1j2+ h D+i Uij + Di ai+1j2 h D+j Uij +Dj bij+12 h Di Uij + Dj bij+12+ h D+j Uij Di (ai+1j Zij1 Di D+i Uij + ai+1j Zi2+1j D+i Dj Uij ) Dj (bij+1Zij2 Dj D+j Uij + bij+1Zij1 +1 D+j Di Uij ) +L hij :
For simplicity, we introduce another two operators, L1 and L by L1 = aij + bij2 2ah (D+i Di + D+i D+j + Di Dj + D+j Dj ); and L = 12 ( D+i aij +2 h Di + Di ai+1j2 h Dj + Dj bij+12 h Di + D+j bij +2 h Dj + Di ai+1j2+ h D+i + D+i aij 2 h D+j + D+j bij 2 h D+i + Dj bij+12+ h D+j ): Observe that (3.21)
L;h hij = L;h Uij = LUij L1Uij ; (i; j ) 2 h; 19
from (3.19), (3.20) and (3.21), we get 2L;h Gij = 2L;h G1ij + 2L;h G2ij = 2L1Uij D+i (aij Zij1 D+i Di Uij + aij Zi2 1j Di D+j Uij ) D+j (bij Zij2 Dj D+j Uij + bij Zij1 1 Dj D+i Uij ) Di (ai+1j Zij1 Di D+i Uij + ai+1j Zi2+1j D+i Dj Uij ) Dj (bij+1Zij2 Dj D+j Uij + bij+1Zij1 +1 D+j Di Uij ):
(3.22)
Meanwhile, by Lemma 3.2, taking summations by parts, we have
X (a
N 1
i;j =1
= =
N
i;j =1 N
i;j =0 N
=
X Di Dj Dj U G ij + ij ij
N 1
i;j =1
ij +
+ ij ij
i;j =1 N
i 1j +
ij
ij
X Dj ( G )Di Dj U X Dj Dj U Di G ij ij + + ij 1 i 1j + ij ij i;j =1 i;j =1 NX1 N X D+j (ij Gij )D+i D+j Uij i 1j D+j Dj Uij Di Gij i;j =0 i;j =1 NX1 NX1 Dj Di Dj U G + Dj G Di Dj U
=
=
ah)(D+j Dj Uij )Gij =
X Di (Dj Dj U G ) i 1j ij ij + i;j =1 N N X X j j i i 1j D D+ D Uij Gi 1j i 1j D+j Dj Uij Di Gij i;j =1 i;j =1 N N X X Di Dj Dj U G Dj Dj U Di G
=
=
ij
+ ij + + ij ij
X Dj Dj U Di G i 1j + ij ij
i;j =1 N 1 (h i;j =0 N
i;j =0
ij +1 + ij + + ij
X a b )Di Dj U G + NX1 Dj G Di Dj U ij +1 + + ij ij ij +1 + ij + + ij i;j =0 X i 1j D+j Dj Uij Di Gij : i;j =1
20
By the symmetric property of the coecients, aij = bij+1 , we get
X (a
N 1
(3.23)
ij
i;j =1 N 1
ah)(D+j Dj + D+i D+j )Uij Gij
N X Dj G Di Dj U X i 1j D+j Dj Uij Di Gij : ij +1 + ij + + ij
=
i;j =0
i;j =1
Proceeding in the same way as before, we obtain
X (b
N 1
i;j =1
(3.24)
=
ij
N X ij Di Gij D+j Dj Ui
i;j =1
X (a
N 1
(3.25)
ij
i;j =1 N 1
=
X Di Dj U Dj G : ij ij +1 + ij
N 1
i;j =0
ah)(D+i Di + Di Dj )Uij Gij
i;j =0
i;j =0
X (b
ij
i;j =1 N
=
1j +
X Dj G Di Di U NX1 Di Dj U Di G : ij +1 + ij + ij ij i+1j + ij
N 1
(3.26)
ah)(D+j Dj + Di Dj )Uij Gij
ah)(D+i Di + D+i D+j )Uij Gij
N X Dj G Di Di U X ij D+i D+j Ui 1j Di Gij : ij 1 ij + ij
i;j =1
i;j =1
Combining (3.23), (3.24), (3.25) and (3.26), we get N a +b X ij ij 2ah j j i Dj + Di Di + Di Dj )U G h2 ( D D + D + + + + ij ij 2 i;j =1
(3.27)
v u NX1 u u t [(D+i Gij h)2 + (D+j Gij h)2] CkU k2 i;j =0 v u u NX1 i p t [(D+Gij h)2 + (D+j Gij h)2]: C u i;j =0
21
Further,
X [Di (a Z 1 Di Di U + a Z 2 Di Dj U )]G h2 ij ij i 1j + ij ij + ij ij + i;j =0 NX1 = [ai+1j Zi1+1j D+i Di Ui+1j + ai+1j Zij2 Di D+j Ui+1j ]D+i Gij h2 i;j =0 v v v u u u NX1 NX1 u u NX1 i 2 2 u u u j i i i t (D+ Gij ) h C(t (D+ D Ui+1j )2h2 + t (D+ D+ Uij )2h2)u i;j =0 i;j =0 i;j =0 v u u NX1 i 2 2 p t (D+ Gij ) h : = C u
(3.28)
N 1
i;j =0
The exact same order for the last three terms in (3.22) can be established similarly. Consequently, from (3.27) and (3.28),
v u N u u t X (Di Gij h)2 + (Dj Gij h)2 = O(p ): i;j =1
By Poincare inequality, By Lemma 3.1 and Lemma 3.3,
kGkh = O(p ):
1 jkU k = O( ); kZ 1(D0i U D0j U )kh i;jmax j Z 1 ij 2
h
which implies, (3.29) k(I L;h1 L;h ) hkh C p : We hence complete the proof. 2 Remark. The result in (3.18) here is consistent with the result established in [8] for the continuous case. Proof of Theorem 3.1. The procedure is exactly the same as that in Theorem 2.1, except di erent inequalities estimated. By Theorem 3.2, we have (3.30)
p
I1 C hk03; 22
instead of (2.20). Therefore, in order to make I1 < 12 , we set (3.31) k0 Ch 1=6; instead of (2.21). For I2 12 , it is sucient to have (3.32)
C hk2 lnk0( k1 2 ): 2
0
0
Combining (3.31) and (3.32), we have
Ch
1
2 3
ln h:
2
As what we have done in previous section, we consequently establish the following main Theorem. Theorem 3.3 There exists constant C , the operator M de ned by (2.5) satis es (M ) 0 < 1; whenever the step size h belongs to the set S (; h0) of Diophantine numbers, and
Ch
1 2=3 lnh:
4 Conclusion
The analysis of the proof strongly indicates us the role of homogenization, which plays in the convergence process. If, for example, the coarse grid operator is replaced by its averaged operator in one dimensional problem [5], the direct estimate for multigrid convergence rate is not asymptotically better than just using the damped Jacobi smoothing operator. This follows from the e ect of the oscillations on the low eigenmodes. The homogenized coarse grid operator reduces the number of smoothing operation from O(h 2 ) to O(h 6=5 ln h), when the step size h belongs to the set S (; h0) of Diophantine numbers. In [9], it has also been shown that the number of smoothing iteration needed for the convergence of the multigird method with the average coarse grid operator guarantees the one with the homogenized coarse grid operator. The theoretical results established in this paper seem a little bit disappointing. From a number of numerical experiments [5, 6], we can get much faster convergence rate in practice 23
than that required in the theoretical results. However, numerical results do indicate that the convergent rate depends on the grid size h for these types of equations with oscillatory coecients [5, 6]. There are some inequalities in the implementation of the proof, which are potential to be improved so that a sharper convergent rate is possible. One of them is to enlarge the space of low eigenmodes, which can be approximated by the corresponding homogenized eigenmodes. Such as to improve (3.18) to (2.14), which we think is the sharpest inequality one can establish. We established the same inequality for the one dimensional case as (2.14) in the two dimensional case oscillatory along a coordinate direction [9]. However, the portion of the eigenmodes that can be approximated by the homogenized ones in later case is relatively much smaller than the previous one. That's why we obtain O(h 4=3 ln h) for the number of smoothing iterations instead of O(h 6=5 ln h) for the later case, although there have the same inequality (2.14) for the space of low eigenmodes. Nevertheless, from the analysis of homogenization, we understand that there always exists a boundary layer [2, 10], which makes it hard to get the rst lower order correction of the eigenfunctions. The case we discussed in chapter 2, which is equivalently to one dimensional problem, doesn't have such a boundary layer. We hence get an estimate as in (2.14). For the case in chapter 3, all we can establish is (3.18), which consists of the result established in [8] for the continuous case. And, it also de nes us a smaller low eigenspace. However, numerical examples tells us that there are also some di erence between these two cases. That a complete understanding of the rst lower order correction for the eigenfunctions is required to further improve the estimates.
References
[1] R.E. Alcou e, A. Brandt, J.E. Dendy, AND J.W. Painter, The Multi-Grid Method for the Di usion Equation with Strongly Discontinuous Coecients, SIAM J. Sci. Stat. Comput., Vol. 2, No. 4 (1981), pp. 430-454. [2] A. Bensoussan, J.L. Lions, AND G. Papanicolaou, Asymptotic Analysis for Periodic Structure, Studies in Mathematics and Its Applications, Vol. 5, North-Holland Publ., 1978. [3] A. Brandt, Multi-level Adaptive Solutions to Boundary-Value Problems, Mathematics of Computation, Vol. 31, No. 136 (1977), pp.333-390. 24
[4] B. Engquist, Computation of Oscillatory Solutions for Partial Di erential Equations, Lecture Notes in Mathematics 1270 (1989), pp. 10-22. [5] B. Engquist AND E. Luo: Multigrid Methods For Di erential Equations With Highly Oscillatory Coecients. Proceedings of the Sixth Copper Mountain Conference on Multigrid Methods, 1993, pp. 175-190. [6] B. Engquist AND E. Luo: New Coarse Grid Operators of Multigrid Methods for Highly Oscillatory Coecient Elliptic Problems. Preprint. [7] L. Giraud AND R.S. Tuminaro: Grid Transfer Operators for Highly Variable Coecient Problems. Preprint. [8] S. Kesavan, Homogenization of Elliptic Eigenvalue Problems: Part 1, Applied Mathematics and Optimization, Vol. 5, (1979), pp.153-167. [9] E. Luo: Multigrid Method for Elliptic Equation with Oscillatory Coecients. Ph.D. Thesis, UCLA, 1993. [10] F. Santosa AND M. Vogelius: First Order Corrections to the Homogenized Eigenvalues of a Periodic Composite Medium. SIAM J. Appl. Math., vol. 53, no. 6, 1993, pp. 16361668. [11] W. Hackbusch, Multi-Grid Methods and Applications, Berlin; New York: springerVerlag, 1985. [12] W. Hackbusch AND U. Trottenburg (eds), Multigrid Methods, Proceedings, KolnPolrz, Lecture Notes in Mathematics 960, Berlin; New York: springer-Verlag.
25
Erding Luo y
Abstract
Standard multigrid methods are not so e ective for equations with highly oscillatory coecients. New coarse grid operators based on homogenized operators are introduced to restore the fast convergence rate of multigrid methods. Finite di erence approximations are used for the discretization of the equations. The analysis of convergence is based on homogenization theory. Proofs are presented for a two level method with the homogenized coarse grid operator to solve two classes of elliptic equations with Dirichlet boundary conditions.
Key Words. elliptic equation, oscillation, nite di erence, multigrid method, homog-
enization theory, convergence
AMS subject classi cations. 65N06, 65N12, 65N55
Department of Mathematics, University of California at Los Angeles, LA, CA90024. Research supported by ARPA/ONR N00014-92-J-1890 and NSF DMS91-03104. y IMA, University of Minnesota, Mpls, MN 55455. Research supported by NSF through IMA.
1
1 Introduction Consider the multigrid method arising from the nite di erence approximations to elliptic equations with highly oscillatory coecients of the following type (1.1)
L u(x) =
X @ a (x) @ u (x) = f (x); @xi ij @xj i;j
where aij (x) = aij (x; x ); aij (x; y) = aji(x; y), strictly positive, continuous and 1-periodic in each component of y. Also, the operator L is uniformly elliptic. That is, there exist two positive constants q and Q independent of , such that (1.2)
0k02 ij M1 ij kh 1 1 1 1 2 2 2 kL;h2 (I L;h L;h1 L;h )khkL;h (1 L;h ) (i2+j2 >k02 ij ij )kh : 2 2 Since kL;h2 kh C and kI L;h L;h1 L;h kh C , 1
1
1
r
(2.22)
2 I2 C i2+j2 >k02 ij2 kL;h (1 L;h ) ij k2h 1 2 (1 C k2 max j ) j: 1
N
Assume 1 + 2 > k102 =
0
N
1N 1
k02
1N 1
( C hk022 ), we have
(2.23) I2 Ck0(1 k02) : For I2 12 , it is sucient to have 2 (2.24)
C h lnk ( 1 ): Combining (2.21) and (2.24), we have
0
k02
Ch
1
The proof is completed. Now, we are ready to show the main result. 11
k02
1 3
ln h:
2
Theorem 2.2 There exists constant C , the operator M de ned by (2.5) satis es (M ) 0 < 1; whenever h belongs to the set S (; h0) of Diophantine number, and
Ch
1 1=3 lnh:
Before we carry out the proof, we need the following lemma. Let
M2 = (L;h1 IHh LH1 IhH )L;hS : Lemma 2.5 For some constant C , kM2kh C : Proof. Since L;h and LH = L;H are the homogenized operators de ned respectively on
ne and coarser grid with constant coecients, they are well behaved [11, 12]. Furthermore, (2.25)
kL;h1 IHh LH1IhH kh Ch2:
Therefore,
(2.26)
kM2kh
kL;h1 IHh LH1IhH kh kL;hS kh Ch2h 2= C:
Proof of Theorem 2.2. Note that
2
M = M1 + M2: Therefore,
kM kh kM1kh + kM2kh:
Since (M ) kM kh , by Theorem 2.1 and Lemma 2.5, the rest of the proof can easily follow. 2 12
3 Oscillation Along the Diagonal Direction 3.1 Model Equation
Consider another special case of (1.1), a two-dimensional elliptic problem with coecients oscillatory diagonally,
@ a (x y) @ @ a (x y) @ = f (x; y); (x; y) 2 ; @x @x @y @y (x; y) = 0; (x; y) 2 @ ; where a(x) is a strictly positive continuous function, and a(x) = a(x=) = a(x= + 1). (3.1)
Also, the operator has the property of (1.2). From (1.3), the corresponding homogenized equation of (3.1) is: a + @ 2 + (a ) @ 2 a + @ 2 = f (x; y); (x; y) 2 ; (3.2) 2 @x2 @x@y 2 @y2 (x; y) = 0; (x; y) 2 @ ; R R where a = 01 a(x)dx, and = ( 01 1=a(x)dx) 1 . As goes to zero, we know that the solution of (3.1) converges to the solution of (3.2). Now, consider a corresponding discretized equation of (3.1), (3.3) D+i aij Di uhij D+j bij Dj uhij = fijh ; (i; j ) 2 h where aij = a( xi h= 2 yj ); bij = a( xi yj + h=2 ); (i; j ) 2 h : Here, we assume the discretized coecients have the following property, aj0 = b0N j+1; j = 0; ; N: Denote the discretization of the homogenized operator a + @ 2 + (a ) @ 2 a + @ 2 2 @x2 @x@y 2 @y2 in (3.2) by (3.4) L;h = h +2 ah (D+i Di + D+j Dj ) + h 2 ah (D+i D+j + Di Dj ); 13
where h = (h PNk=1 a10 ) 1, and ah = h PNk=1 ak0. The operator of the two level method can be expressed as k
M = (I IHh LH1 IhH L;h)S = (I IHh LH1IhH L;h )(I L;h ) ; where is the inverse of the largest eigenvalue of L;h, has order of h 2. And LH = L;H . (3.5)
3.2 Convergence Analysis
We still rst consider the simpli ed operator M1 de ned in (2.6). Theorem 3.1 If the ratio of h to is strictly irrational, h belongs to the set of Diophantine number S (; h0) de ned in (1.6), then there exist two constants C and 0 such that
kM1kh 0 < 1
(3.6) whenever Ch
1 2=3ln(h).
In order to prove Theorem 3.1, we introduce the following lemmate rst. Lemma 3.1 De ne two discrete functions on h , (3.7) (3.8)
Xi Xj Zij1 = 2h (i j h a1 + h b1 ); (i; j ) 2 h; k=1 kj i 1 h k=1 k0
X Zij2 = 2h (j i + a
k=1 0k j 1) h k=1 ik
X
b ; (i; j ) 2 h :
Then, Zij1 ; Zij2 are bounded, and satisfy
(3.9) for (i; j ) 2 h .
L;h Zij1 = 1 D+i aij ; L;h Zij2 = 1 D+j bij ;
Proof. Notice that by the assumption of the coecients, N h X
X h = 1; (i; j ) 2 : = h h k=1 akj k=1 bik N
14
Applying the operator L;h to Zij1 , (3.9) then follows. Now, we write Zij1 as follows, j 1 i 1 X X h h ) 2 (j h b ): ij = 2 (i h k=1 akj k=1 0k
Z1 By Lemma 1.1,
j 1 h (i Xi 1 ) = O(1); h (j X ) = O(1): h h 2 2 k=1 akj k=1 b0k Hence, Zij1 is bounded for (i; j ) 2 h , and satis es aij (1 Di Zij1 ) = aij +2 h ; bij Dj Zij1 = bij 2 h :
The result can be deduced similarly for Zij2 . We can also show that Zij2 = Zij1 , and satis es bij (1 Dj Zij2 ) = bij +2 h ; aij Di Zij2 = aij 2 h : This proves the lemma. 2 1 2 Remark. The explicit forms of Zij and Zij depend on a (x y) being a function of x y . For the general angular dependences a( x + y), these forms would not be possible.
Lemma 3.2 Assume discrete function ij de ned as
Xi Xj b ); (i; j ) 2 : ij = h ( akj iah) + h (jah 0k h k=1
k=1
Then, ij is bounded, and has the following properties
Di ij = aij ah; Dj ij = ah bij ; (i; j ) 2 h: Proof. Using the symmetric properties of the coecients, we can establish the results similarly as in Lemma 3.1. 2
Lemma 3.3 Assume Uij satis es (3.10)
L;h Uij = hij ; (i; j ) 2 h ; Uij = 0; (i; j ) 2 @ h; 15
where (; hij ) is a normalized eigenpair of L;h . Then, we have
X [(Di U h)2 + (Dj U h)2] = O( ); + ij + ij
N 1
kU k21 =
(3.11) and
i;j =0
X ((Di Di U
N 1
i;j =1
+
2 + (D i
ij )
Dj Uij )2 + (D+i D+j Uij )2 + (D+j Dj Uij )2)h2
(3.12) = O(2 ): Proof. First, we observe that
X [(Di h h)2 + (Dj h h)2] = O( ):
N 1
+ ij
+ ij
i;j =0
Applying Uij to the following equation and taking summation, we have
X L U U = NX1 L U : ;h ij ij ;h ij ij
N 1
i;j =1
i;j =1
(3.11) then follows. For (3.12), since Uij vanishes at the boundary, it satis es 1 NX1 [(Di Di U )2 + (Dj Dj U )2] NX1 (Di Di U )(Dj Dj U ) ij ij ij ij + + + 2 i;j=1 + i;j =1 (3.13)
=
N X (Di Dj U
i;j =1
ij
)2
=
X (Di Dj U )2: + + ij
N 1
i;j =0
Multiplying D+i Di Uij to (3.10), we get
(D+i Di Uij )2 + 2 aah + h (D+i D+j Uij )(D+i Di Uij ) + (D+i D+j Uij )2 + h
h
(D+i Di Uij )2 + 2 aah + h (Di Dj Uij )(D+i Di Uij ) + (Di Dj Uij )2 + h h j j i i ((D+ D Uij )(D+ D Uij ) (D+i D+j Uij )2) + ((D+j Dj Uij )(D+i Di Uij ) (Di Dj Uij )2) = a +2 hij D+i Di Uij ; (i; j ) 2 h : h
h
16
By the uniformly ellipticity property of the homogenized operator (3.4), we have (3.14)
X ((Di Di U
N 1
i;j =1
2 + (Di Dj U )2 + (Di + + ij
Dj Uij )2)h2 = O(2 ):
2 + (Di Dj U )2 + (Di + + ij
Dj Uij )2)h2 = O(2 ):
ij )
+
An similar argument gives (3.15)
X ((Dj Dj U
N 1
i;j =1
ij )
+
(3.13), (3.14) and (3.15) complete the proof of lemma. 2 Lemma 3.4 Assume that Uij de ned in Lemma 3.3 satis es the following boundary conditions, (3.16) D+i Uij + Di Uij = D+j Uij + Dj Uij = 0; (i; j ) 2 @ h: Then,
kU k22 = (3.17)
+
N X ((Di Di U
i;j =0 N 1
+
2 + (Dj Dj U )2 )h2 ij +
ij )
N X (Di Dj U h)2 + X (Di Dj Uij h)2 = O(2 =h): ij + +
i;j =0
i;j =1
Proof. Since
Di UNj = Di Uij + h we have Then,
(Di UNj )2 (Di Uij )2 + N X (Di U
j =1
2h
Nj )
k=i
X (Di Di U )2h: ij +
N 1 i=1
N N X X (Di Uij h)2 + (D+i Di Uij h)2 = O(2 ):
i;j =1
Combined (3.12) with the following relation, the rest of the proof follows.
X Di Di U ; kj +
N 1
i;j =1
D+i Di UNj = h2 D+i UN 1j ; 17
2
Theorem 3.2 Assume hij ; Uij are de ned in Lemma 3.3 and Lemma 3.4. Then, (3.18) k h U kh = O(p ): Proof. First, we introduce the following discrete functions
Gij = G1ij + G2ij ; G1ij = 21 Uij 2 (Zij1 D+i Uij + Zij2 D+j Uij ) 12 hij ; G2ij = 12 Uij 2 (Zij1 Di Uij + Zij2 Dj Uij ) 12 hij ;
for (i; j ) 2 h. By the assumption (3.16), Gij vanishes at boundary, i.e.,
Gij = 0; (i; j ) 2 @ h : For G1ij , we have 2Di G1ij = Di Uij (Di (Zij1 D+i Uij ) + Di (Zij2 D+j Uij )) Di hij = (Di Uij Di Zij1 D+i Ui 1j ) Zij1 Di D+i Uij
Di Zij2 D+j Uij Zi2 1j Di D+j Uij Di hij ;
2Dj G1ij = Dj Uij (Dj (Zij1 D+i Uij ) + Dj (Zij2 D+j Uij )) Dj hij = (Dj Uij Dj Zij2 D+j Uij 1 ) Zij2 Dj D+j Uij
Dj Zij1 D+i Uij Zij1 1Dj D+i Uij Dj hij :
Then,
(3.19)
2L;h G1ij = D+i aij +2 h Di Uij + D+i aij 2 h D+j Uij D+i (aij Zij1 D+i Di Uij + aij Zi2 1j Di D+j Uij ) +D+i bij 2 h D+j Uij + D+j bij +2 h Dj Uij D+j (bij Zij2 Dj D+j Uij + bij Zij1 1 Dj D+i Uij ) +L hij : 18
For G2ij , we establish similarly 2D+i G2ij = D+i Uij (D+i (Zij1 Di Uij ) + D+i (Zij2 Dj Uij )) D+i hij = (D+i Uij D+i Zij1 Di Ui+1j ) Zij1 Di D+i Uij
D+i Zij2 Dj Uij Zi2+1j D+i Dj Uij D+i hij ;
2D+j G2ij = D+j Uij (D+j (Zij1 Di Uij ) + D+j (Zij2 Dj Uij )) D+j hij = (D+j Uij D+j Zij2 Dj Uij+1 ) Zij2 D+j Dj Uij
D+j Zij1 Di Uij Zij1 +1D+j Di Uij D+j hij :
Then,
(3.20)
2L;h G2ij = Di ai+1j2+ h D+i Uij + Di ai+1j2 h D+j Uij +Dj bij+12 h Di Uij + Dj bij+12+ h D+j Uij Di (ai+1j Zij1 Di D+i Uij + ai+1j Zi2+1j D+i Dj Uij ) Dj (bij+1Zij2 Dj D+j Uij + bij+1Zij1 +1 D+j Di Uij ) +L hij :
For simplicity, we introduce another two operators, L1 and L by L1 = aij + bij2 2ah (D+i Di + D+i D+j + Di Dj + D+j Dj ); and L = 12 ( D+i aij +2 h Di + Di ai+1j2 h Dj + Dj bij+12 h Di + D+j bij +2 h Dj + Di ai+1j2+ h D+i + D+i aij 2 h D+j + D+j bij 2 h D+i + Dj bij+12+ h D+j ): Observe that (3.21)
L;h hij = L;h Uij = LUij L1Uij ; (i; j ) 2 h; 19
from (3.19), (3.20) and (3.21), we get 2L;h Gij = 2L;h G1ij + 2L;h G2ij = 2L1Uij D+i (aij Zij1 D+i Di Uij + aij Zi2 1j Di D+j Uij ) D+j (bij Zij2 Dj D+j Uij + bij Zij1 1 Dj D+i Uij ) Di (ai+1j Zij1 Di D+i Uij + ai+1j Zi2+1j D+i Dj Uij ) Dj (bij+1Zij2 Dj D+j Uij + bij+1Zij1 +1 D+j Di Uij ):
(3.22)
Meanwhile, by Lemma 3.2, taking summations by parts, we have
X (a
N 1
i;j =1
= =
N
i;j =1 N
i;j =0 N
=
X Di Dj Dj U G ij + ij ij
N 1
i;j =1
ij +
+ ij ij
i;j =1 N
i 1j +
ij
ij
X Dj ( G )Di Dj U X Dj Dj U Di G ij ij + + ij 1 i 1j + ij ij i;j =1 i;j =1 NX1 N X D+j (ij Gij )D+i D+j Uij i 1j D+j Dj Uij Di Gij i;j =0 i;j =1 NX1 NX1 Dj Di Dj U G + Dj G Di Dj U
=
=
ah)(D+j Dj Uij )Gij =
X Di (Dj Dj U G ) i 1j ij ij + i;j =1 N N X X j j i i 1j D D+ D Uij Gi 1j i 1j D+j Dj Uij Di Gij i;j =1 i;j =1 N N X X Di Dj Dj U G Dj Dj U Di G
=
=
ij
+ ij + + ij ij
X Dj Dj U Di G i 1j + ij ij
i;j =1 N 1 (h i;j =0 N
i;j =0
ij +1 + ij + + ij
X a b )Di Dj U G + NX1 Dj G Di Dj U ij +1 + + ij ij ij +1 + ij + + ij i;j =0 X i 1j D+j Dj Uij Di Gij : i;j =1
20
By the symmetric property of the coecients, aij = bij+1 , we get
X (a
N 1
(3.23)
ij
i;j =1 N 1
ah)(D+j Dj + D+i D+j )Uij Gij
N X Dj G Di Dj U X i 1j D+j Dj Uij Di Gij : ij +1 + ij + + ij
=
i;j =0
i;j =1
Proceeding in the same way as before, we obtain
X (b
N 1
i;j =1
(3.24)
=
ij
N X ij Di Gij D+j Dj Ui
i;j =1
X (a
N 1
(3.25)
ij
i;j =1 N 1
=
X Di Dj U Dj G : ij ij +1 + ij
N 1
i;j =0
ah)(D+i Di + Di Dj )Uij Gij
i;j =0
i;j =0
X (b
ij
i;j =1 N
=
1j +
X Dj G Di Di U NX1 Di Dj U Di G : ij +1 + ij + ij ij i+1j + ij
N 1
(3.26)
ah)(D+j Dj + Di Dj )Uij Gij
ah)(D+i Di + D+i D+j )Uij Gij
N X Dj G Di Di U X ij D+i D+j Ui 1j Di Gij : ij 1 ij + ij
i;j =1
i;j =1
Combining (3.23), (3.24), (3.25) and (3.26), we get N a +b X ij ij 2ah j j i Dj + Di Di + Di Dj )U G h2 ( D D + D + + + + ij ij 2 i;j =1
(3.27)
v u NX1 u u t [(D+i Gij h)2 + (D+j Gij h)2] CkU k2 i;j =0 v u u NX1 i p t [(D+Gij h)2 + (D+j Gij h)2]: C u i;j =0
21
Further,
X [Di (a Z 1 Di Di U + a Z 2 Di Dj U )]G h2 ij ij i 1j + ij ij + ij ij + i;j =0 NX1 = [ai+1j Zi1+1j D+i Di Ui+1j + ai+1j Zij2 Di D+j Ui+1j ]D+i Gij h2 i;j =0 v v v u u u NX1 NX1 u u NX1 i 2 2 u u u j i i i t (D+ Gij ) h C(t (D+ D Ui+1j )2h2 + t (D+ D+ Uij )2h2)u i;j =0 i;j =0 i;j =0 v u u NX1 i 2 2 p t (D+ Gij ) h : = C u
(3.28)
N 1
i;j =0
The exact same order for the last three terms in (3.22) can be established similarly. Consequently, from (3.27) and (3.28),
v u N u u t X (Di Gij h)2 + (Dj Gij h)2 = O(p ): i;j =1
By Poincare inequality, By Lemma 3.1 and Lemma 3.3,
kGkh = O(p ):
1 jkU k = O( ); kZ 1(D0i U D0j U )kh i;jmax j Z 1 ij 2
h
which implies, (3.29) k(I L;h1 L;h ) hkh C p : We hence complete the proof. 2 Remark. The result in (3.18) here is consistent with the result established in [8] for the continuous case. Proof of Theorem 3.1. The procedure is exactly the same as that in Theorem 2.1, except di erent inequalities estimated. By Theorem 3.2, we have (3.30)
p
I1 C hk03; 22
instead of (2.20). Therefore, in order to make I1 < 12 , we set (3.31) k0 Ch 1=6; instead of (2.21). For I2 12 , it is sucient to have (3.32)
C hk2 lnk0( k1 2 ): 2
0
0
Combining (3.31) and (3.32), we have
Ch
1
2 3
ln h:
2
As what we have done in previous section, we consequently establish the following main Theorem. Theorem 3.3 There exists constant C , the operator M de ned by (2.5) satis es (M ) 0 < 1; whenever the step size h belongs to the set S (; h0) of Diophantine numbers, and
Ch
1 2=3 lnh:
4 Conclusion
The analysis of the proof strongly indicates us the role of homogenization, which plays in the convergence process. If, for example, the coarse grid operator is replaced by its averaged operator in one dimensional problem [5], the direct estimate for multigrid convergence rate is not asymptotically better than just using the damped Jacobi smoothing operator. This follows from the e ect of the oscillations on the low eigenmodes. The homogenized coarse grid operator reduces the number of smoothing operation from O(h 2 ) to O(h 6=5 ln h), when the step size h belongs to the set S (; h0) of Diophantine numbers. In [9], it has also been shown that the number of smoothing iteration needed for the convergence of the multigird method with the average coarse grid operator guarantees the one with the homogenized coarse grid operator. The theoretical results established in this paper seem a little bit disappointing. From a number of numerical experiments [5, 6], we can get much faster convergence rate in practice 23
than that required in the theoretical results. However, numerical results do indicate that the convergent rate depends on the grid size h for these types of equations with oscillatory coecients [5, 6]. There are some inequalities in the implementation of the proof, which are potential to be improved so that a sharper convergent rate is possible. One of them is to enlarge the space of low eigenmodes, which can be approximated by the corresponding homogenized eigenmodes. Such as to improve (3.18) to (2.14), which we think is the sharpest inequality one can establish. We established the same inequality for the one dimensional case as (2.14) in the two dimensional case oscillatory along a coordinate direction [9]. However, the portion of the eigenmodes that can be approximated by the homogenized ones in later case is relatively much smaller than the previous one. That's why we obtain O(h 4=3 ln h) for the number of smoothing iterations instead of O(h 6=5 ln h) for the later case, although there have the same inequality (2.14) for the space of low eigenmodes. Nevertheless, from the analysis of homogenization, we understand that there always exists a boundary layer [2, 10], which makes it hard to get the rst lower order correction of the eigenfunctions. The case we discussed in chapter 2, which is equivalently to one dimensional problem, doesn't have such a boundary layer. We hence get an estimate as in (2.14). For the case in chapter 3, all we can establish is (3.18), which consists of the result established in [8] for the continuous case. And, it also de nes us a smaller low eigenspace. However, numerical examples tells us that there are also some di erence between these two cases. That a complete understanding of the rst lower order correction for the eigenfunctions is required to further improve the estimates.
References
[1] R.E. Alcou e, A. Brandt, J.E. Dendy, AND J.W. Painter, The Multi-Grid Method for the Di usion Equation with Strongly Discontinuous Coecients, SIAM J. Sci. Stat. Comput., Vol. 2, No. 4 (1981), pp. 430-454. [2] A. Bensoussan, J.L. Lions, AND G. Papanicolaou, Asymptotic Analysis for Periodic Structure, Studies in Mathematics and Its Applications, Vol. 5, North-Holland Publ., 1978. [3] A. Brandt, Multi-level Adaptive Solutions to Boundary-Value Problems, Mathematics of Computation, Vol. 31, No. 136 (1977), pp.333-390. 24
[4] B. Engquist, Computation of Oscillatory Solutions for Partial Di erential Equations, Lecture Notes in Mathematics 1270 (1989), pp. 10-22. [5] B. Engquist AND E. Luo: Multigrid Methods For Di erential Equations With Highly Oscillatory Coecients. Proceedings of the Sixth Copper Mountain Conference on Multigrid Methods, 1993, pp. 175-190. [6] B. Engquist AND E. Luo: New Coarse Grid Operators of Multigrid Methods for Highly Oscillatory Coecient Elliptic Problems. Preprint. [7] L. Giraud AND R.S. Tuminaro: Grid Transfer Operators for Highly Variable Coecient Problems. Preprint. [8] S. Kesavan, Homogenization of Elliptic Eigenvalue Problems: Part 1, Applied Mathematics and Optimization, Vol. 5, (1979), pp.153-167. [9] E. Luo: Multigrid Method for Elliptic Equation with Oscillatory Coecients. Ph.D. Thesis, UCLA, 1993. [10] F. Santosa AND M. Vogelius: First Order Corrections to the Homogenized Eigenvalues of a Periodic Composite Medium. SIAM J. Appl. Math., vol. 53, no. 6, 1993, pp. 16361668. [11] W. Hackbusch, Multi-Grid Methods and Applications, Berlin; New York: springerVerlag, 1985. [12] W. Hackbusch AND U. Trottenburg (eds), Multigrid Methods, Proceedings, KolnPolrz, Lecture Notes in Mathematics 960, Berlin; New York: springer-Verlag.
25
