LIBRARY OF THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAICN cop.'Z. The person charging this material is re- sponsible for its return to the library from which it was withdrawn on or before the Latest Date stamped below. Theft, mutilation, and underlining of books are reasons for disciplinary action and may result in dismissal from the University. UNIVERSITY OF ILLINOIS LIBRARY AT URBANA-CHAMPAIGN OCT TWTI SEP! 5^ c * L161 — O-1096 Digitized by the Internet Archive in 2013 http://archive.org/details/implicitmethodfo392brac ' / / l4-s r *-'V—' REPORT NO. 392 y 7 :L - 7 AN IMPLICIT METHOD FOR THE SOLUTION OF ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS April 15, 1970 by Aran on Brae ha NOV 91972 UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN URBANA, ILLINOIS Report No. 392 AN IMPLICIT METHOD FOR THE SOLUTION ^ OF ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS by Amnon Bracha April 15, 1970 Department of Computer Science University of Illinois Urbana, Illinois 6l801 This work was supported in part by the National Science Foundation under Grant No. US NSF-GJ-328 and was submitted in partial fulfillment for the degree of Master of Science in Computer Science, June, 1970. Ill ACKNOWLEDGMENT I would like to express my appreciation to Professor P. E. Saylor, who suggested the thesis topic and supervised its development. This thesis is dedicated to the memory of my mother. IV TABLE OF CONTENTS Page 1. INTRODUCTION 1 2. PRELIMINARIES 3 3. THE FACTORIZATION 7 h. ANALYSIS AND CONVERGENCE CONDITIONS Ik 5. COMPUTATIONAL NOTES 26 LIST OF REFERENCES 28 1 . INTRODUCTION Let Ax = b be a linear system of algebraic equation resulting from the discretization of an elliptic boundary value problem. A direct method for computing the solution typically is equivalent to factoring A as the product A = LU of a lower triangular matrix L by an upper triangular matrix U. Such methods are computationally inefficient because the factors L and U are not sparse matrices. For this reason, an iterative procedure is generally used rather than a direct method. Iterative procedures have recently been proposed by several people, notably H. L. Stone [5], and Dupont, Kendall and Rachford [2], based on the idea of finding a matrix A + B that can be factored as the product of sparse matrices L and U, A + B = LU. The efficient solution of (A + B)u=b is then used in the iterative scheme defined by (A + B)u n+1 = (A + B)u n - x(Au n - b), where t is a parameter. Stone refers to such procedures as strongly implicit (Si) procedures. The subject of this thesis is an SI procedure due to Stone [k] . The factorization in the procedure depends on a parameter a, and we indicate this dependence below with a subscript. Stone designed the factorization to make L U = A + B symmetric. The main result of this thesis is that A + B is defined and positive definite for a satisfying < a < 1, when A is derived from a discretization of the Dirichlet problem. We also show that, for A and A + B positive definite, the procedure (A + Bju . = (A + B )u - t(Au - b) v a n+1 an n 2 converges for values of i satisfying < t < A II ' 112 2 . PRELIMINARIES Consider the Dirichlet problem - £- ^ ^ > - 4 (a 2 (x) h ] = Q(x) > xeD (2.1) u(x) = 0, xedD where x = (x n ,x ), D is the interior of a compact region, and where the differential operator in (2.1) is elliptic, that is, a i (x)S l + a 2 (x)£ 2 > ° for (\,Z 2 ) + (0,0) and x € D. Cover the region D with a rectangular grid system, replace the differential operator in (2.1) with a difference operator, and replace u with a vector whose components correspond to the grid points of the system. The result is a set of simultaneous linear equations with a one-to-one correspondence between the gridpoints and the equations. We shall study a method for the iterative solution of this set of equations. In order to make the problem more specific, let D be the unit square, D = {(x.,x 2 ) : < X;l ,x 2 < 1} . Let h = — - , n a positive integer. Let D be the grid system over D, with uniform spacing h , defined by D h = {(jh,kh) : 1 < j,k < n}. Denote the grid point (jh,kh) by (j,k) The boundary, ^D , of D is defined by £D h = {(jh,0)} U {(l,kh)} U {(jh,l)3 U {(0,kh)} for < j,k < n+1. Order the points of D in a left -to-right, down-to-up fashion. h That is, if (j ,k ) and (j ,k ) are two gridpoints, then (j ,k ) precedes («5 2 >k ) if k x < k 2 or if ^ = k 2 and ^ < o 2 - 2 Let u be an n -vector whose components are ordered according to the order of D . Denote the components of u by u. , (j,k) a gridpoint. h j,k Let the matrix A defined by (Au) . , = B. . u. . _ + D. . u. . . + E. . u. . j,k j,k j,k-l j,k j-l,k j,k j,k (2.2) + D. u. + B. . , u. 0+1, k o+l, k o,k+l o,k+l : version of the differential operator defined in (2.1), and the equat; (Au) 4,k = V he boundary value problem in (2.1). The important properties required of A are that (1) V' V* ' (ii) D = for J = 1, k = 1, 2, ..., n; (iii) B. . = for k = 1, J = 1, 2, ..., n; (2.3) (IT) E. )k =-(B. )k + D . ;k + D. +1;k + B. >M )>0 ) if D . ;k ^O and B.^0; (v) V > " ( V + V + W + V^' if D. . = or B. . = 0. These properties result from any of the common discretizations of the quantity r — (a.(x) r — u(x)). dx. l dx. l l For example, §xT (a i U) §x7 U(X)) " h " 2[a i (x + I *JM*) - u(x + he.) + a t (x - - e i )[u(x) - u(x - he^]}, vhere e , e are the unit vectors along the x and x axes This particular discretization yields, (2.10 E j,k = a 2 ( J^ (k " | )h) + a 2 ^ h ^ k + l^ + a x ((j - -)h,kh) + a 1 ((j + -)h,kh). When j = 1 or k = 1, we adopt the convention that B = or D =0 3 >& 3 >k respectively. Relations (2.3) 'follow. 3. THE FACTORIZATION The difficulty in using Gaussian elimination to solve Ax = b is that forward reduction transforms A with only five non-zero diagonals to an upper triangular matrix with n non-zero diagonals. Since Gaussian elimination is equivalent to factoring A as A = L»U where L is lower triangular and U is upper triangular with its main diagonal consisting of l's, and since the factorization is unique, it follows that A cannot be factored as the product of sparse lower and upper triangular matrices. The idea of a strongly implicit procedure is to replace the matrix A with a matrix of the form A + B, where A + B can be written as the product of sparse lower and upper triangular matrices A + B = L-U. The SI procedure of Stone constructs lower and upper triangular matrices, each with three non-zero diagonals in the same positions as the non-zero diagonals of A. That is, (3.1) (Lu) . , = b. , u. , , + c. . u. n , + d. . u. v y j,k j,k j,k-l j,k j-l,k j,k 3, (Uu) . . = vl. . + e. . u. . . + f . , u. v , 1 v y j,k j,k j,k j+l,k j,k j,k+i k As a preliminary to the formal statement of the algorithm for computing L and U, observe that the product L*U is given by [(A + B)u] j>k = KL-nh] = b jjk » j,k j,k-l j+l,k-l j,k j-l,k (3 - 2) + (a j,k + b j,k ' f j,k-i + c j,k •j-i.k) V + d. . * e . _ u. + c. • f. u. j,k j,k j+l,k j,k j-l,k o-l, k+1 + d. • f . u. . , . 0,k j,k j,k+l Let < CC < 1. Stone suggests in [k] the following algorithm for computing L and U: j,k " j,k j,k-l j-l,k-l . c- v = D. - cft>. . . • e. 0,k j^k 0-1, k j-l,k-l d. + b. • f + c. . * e. = E. . j,k o,k j,k-l j,k j-l,k j,k + olc . . . • f . _ . . + ab . . . • e . . . . J,k-1 j-l,k-l j-l,k j-l,k-l d. • e. = D. _ . - eft). , • e. . . J>k j,k o+l, k j,k j,k-l d 0, • =B 0,K + l- aC 0,k' f o-l,k J^- 1 ' 2 ' — n -he lefl Jide of (3.3) forms the elements of the (j,k) row of 9 Stone proved in [h] that L-U is symmetric and his proof is reproduced here. First, observe that in the actual computational algorithm it is not necessary to compute b. and c. since it can be seen by adjust- ing subscripts in (3.3) that and (3.5) c. = d • e. . j,k j-l,k o-l,k Changing the subscripts in (3.5) gives (3.6) C n i T= d -T -, * e - ■, -i • 0+1, k-1 j,k-l j,k-l Multiplying the right hand side of (3.6) by the left side of (3.*0 gives b. . (d. . _ • e. ) = c. _ . n (d. • f. .) 0,k v j,k-l a,k-l' j+l,k-l x j,k-l j,k-l or (3 7) b. , • e. , _ = c. .. , , * f . . , Ki ' u o,k j,k-l j+l,k-l j,k-l Equation (3.7), along with (3.U) and (3.5), establishes the symmetry of A + B as defined in (3.2). We now prove a lemma which establishes a useful set of relations between the elements of A and the elements of A + B. Lemma 1 : The quantities defined in equation (3-3) exist and satisfy the following relations: 10 (3.8) Proof: < a > \ k < B j^ (b) c J)k ^ d j,^° (c) V >0 (a) -l0 ' D 2 1 -l 1 E l.l B ] - 1 ^ f l,l = ^-< ^ 1,1 .1 B l D 2 .1 + B l .2 + E l.l . ti i,r i -E 1 1 -1,1 -1,1 11 Now assume that the relations (a) - (f) hold for all points preceding (j,k). To prove (3.8) (a) and (b), we have c . . . < o, f ; o, b. - . < 0, e . , , _ < and < a < 1. j-l,k-l - - - It follows that j,k " j,k j,k-l j-l,k-l - j,k - and c. ,=D. -ab._. • e . _ . . < D . . < 0. J,k j,k j-l,k j-l,k-l - j,k - To prove (3.8) (c) we have, by assumption, 1 + f _ > and 1 + e. -, , > 0. Therefore, j-l,k - > d. . = E. . - b. . ■ f. . _ - c. . • e. J,k j,k j,k j,k-l j,k j-l,k + ccc . . _ . f . _ . _ + ab . _ . • e . _ j,k-l j-l,k-l j-l,k j-l,k-l = E. . + B. . +D. n -b. . (l + f . . _) 0,k j,k j,k j,k v j,k-l Next we prove (3.8) (d) and (e). 12 Certainly, D. . . - Ob. e. . .. e = , l+ 1 » k ,],k j,k-l < ^ k d j,k and f = ■i» k + 1 ^ k J" 1 ^ < o j' k d j,k To prove that -1 is a lower bound for these quantities, observe that e, k • d + d Since (3.8) (d) and (e) are valid for (j,k-l), it follows from (3.3) that 0,k j,k j,k " j+l,k j,k j,k-l j,k j,k j+l,k j,k j,k-l j,k o,k-l c . . • e . . . + ac . . . • f . _ . _ + ab . _ . • e . j,k j-l,k j,k-l j-l,k-l o-l, k o-l,k-l = E . . + D . . . + B . . + D . . - b . . ( 1 + ae , . _ + f . . , ) 0,k j+l,k j,k j,k j,k' 0,k-l o,k-l' . . (1 + e 4 . . ) > Therefore, -1 < e. . . , k • _1 < f • v < °- ,k To prove the remaining assertion, we have from (3.3) D, . , - Ob ' e. . . e f + 1 m J±Lk j^k jJszl a. 0",k B, . , - ac, . • f. . . + J» k+1 .). k .1-1^ +1 13 d+D+B - Ob • e . - ac , . • f . _ . ,i.k .1+1. k .i.k+1 ,i.k .i.k-1 ,i.k .i-l. k d. . . .i,k + ,1+1, k + ,1,k+l " ,1,k e ,j,k-l c ,1,k ,1-l,k~,1,k* j,k-l -c. • e. _ + ac. • f. + ab. • e. . . , >^ k .Tl. k ,i,k-l .1-1, k-1 ,t-1, k ,i-l, k-1 d - V = ,i,k + ,j,k + ,j,k + ,j,k+l + ,1+1, k _ ,1 ? k f 1+0Le + f ) c . - -i*£ (1 + af . _ . + e. ) d j-l,k j-l,k' By the induction hypothesis, and from the facts that 1 + Qte . . . + f . . _ > 1 + e . . n + f . . n > j, k-1 j, k-1 - o,k-l J, k-1 - and 1 + af . . . + e. . . > 1 + e. + f > 0, j-l,k j-l,k - J-l,k j-l,k - it follows that e. . + f . . + 1 > 0. Ik h. ANALYSIS AND CONVERGENCE CONDITIONS The iteration ve consider is defined as follows: (h.l) (A + B)u - = (A + B)u - t(Au - q). n+1 n n Our objectives in this section are to prove that A + B is positive definite and that (4.1) converges for % e (0,2/||A|| ). First, we establish the notation and conventions that we will use. Let n denote the number of grid 2 points on a horizontal or vertical line. Let N = n . In the remainder of this thesis, it is convenient to rely on standard matrix notation. Thus, if C = (c . . ) is a matrix, 1 < i,j < n, then c. . is the element at the i-th row i.j '—'"—■' i j and the j-th column of the matrix, and the notation c. . no longer refers to an element of the row of the matrix corresponding to the gridpoint defined by the intersection of the i-th horizontal line with the j-th vertical line. The next lemma provides a convenient representation of the quadratic forms and . Lemma 2 : For the matrix A, whose properties are listed in (2.3), and for B defined by B = LU - A, we have N-l 2 N-n 2 = Z-a. . n |u. ., -11.'"+ E-a. . I u. - u. ' . n l.i+l' l+l i' . . 1,1+n 1 l+n i' i=l ' i=l ' N Z( + L (a. . + a. . . + a. . + a. n . + a. .)|u. . , 1,1 i,i+l i,i+n i-l,i i-n,i ' l 15 = a Z b. . lu. - u. I* - a E b. . lu I 2 1,1+n 1 l+n 1 1 1,1+n 1 i 1 i=2 ' i=2 ' N-n N-n N-l N-l x E b. . _ lu. _ - u. I - a Z , l.i+i 1 l+i i 1 i=n+l ' i=n+l + a Z b. . _ lu. - - u. I" - a Z b. . Ju.r 1,1+1' i+I i' . _ , i,i+l i' N-n-1 Z b. , . lu. - - u. . n i+I, i+n' i+I i+n' i=l ' N-n-1 2 2 + Z b. . ( lu. _ I" + lu. ) . n i+I, i+n ' i+I' ' i+n' i=l ' Proof: We only prove the second part of the assertion. The first part is a familiar fact. We have N = - OL Z b. . • u. * u. ' . _ 1,1 -n l-n l i=n+2 ' N-l + Z b. . _ • u. , ' u. n i,i+l-n i+1-n l i=n+l ' - a Z b. . , • u. _ • u. ~ i,i-l i-l i i=n+2 ' N + a Z b. . lu. i=n+2 ' N-l - a Z b. . n • u . . • u. n i,i+l i+l i i=n+l ' 16 N-n + Z b. . . • u. , • u. . ^ i,i+n-l i+n-1 1 i=2 ' N-n - a Z b. . * u. • u . , . i,i+n l+n l i=2 Since B is symmetric, and since b. . = b. . + b. . ,, 1,1 i,i-n i,i-l' N u. = - a Z b. . * u. „ i-n,i l-n i=n+2 ' N-l + Z b. , . • 11. . • u. , i+l-n,i i+1-n i i=n+l ' N - a Z b. . . • u. _ • u. i=n+2 1 " 1 ' 1 1_1 X N p + a Z (b, i=n+2 >. . + b. _ . ) In. i-n,i i-l,i ' i N-l - a Z b. . _ • u. _ • u. , i,i+l i+I i i=n+l ' N-n + Z b. . . • u. -, * u. , p i,i+n-l i+n-1 i N-n - a I b . . • u . • u . P i,i+n i+n i IT In the expression on the right, change the variable of summation in the first sum and combine with the last sum. Also, change the variable of summation in the third sum and combine with the fifth sum. Finally, change the variable of summation in the second and the sixth sums. We have, N-n = - a Z b. . (u. • u. + u. * u. ) . _ 1,1+n 1 l+n l+n 1 1=2 ' N-l a Z b. . n (u. • u. . + u. , • u.) _ 1,1+1 1 1+1 1+1 1 i=n+l ' N-n-1 + Z b. _ . (u. ' u. + u. • u , ) . n i+l,i+n v l+l l+n l+n l+l i=l ' and N z i=n+2 2 + a Z (b. . + b. . ) lu. . From the identities -u. u. - u. u. = lu. - XL. I " " - (lu |'" + |u | ), i+n i i l+n l+n l ' ' i+n ' i -u. _ u. - u. u. n = hi. - - u. - u J + u ) i+1 i l l+l ' l+l i 1 ' i+l 1 i u. n u. + u. u. ., = - u. n - u. t \ + U. - " + u. , i+l i+n i+n i+l ' i+l i+n 1 ' i+l 1 ' i+n 1 ' we have 18 N-n = a Z b. . lu. - u. ._ p 1,1+n 1 l+n 1 N-n 2 2 - a Z b. . ( lu. I " + lu. I ) . „ i,i+n v ' i+n' ' i' i=2 ' N-l 2 + a Z b. . _ lu. , - u. | _ l.i+l 1 i+I i' i=n+l ' N-l 2 g - a Z b. . _ ( lu. _ I + lu. I ) n i,i+l V| i+l ' ' i 1 i=n+l ' N-n-1 2 Z b. . . lu. _ - u. . . i+l, i+n' i+l i+n' i=l ' N-n-1 + Z b. . . ( u. n + u. ) . _ i+l, i+n ' i+l ' ' i+n' i=l ' N-n 2 + a Z b. . lu. . . i,i+n' i+n' i=2 ' N-l 2 + a Z b. lu. ,, I . , 1,1+1' 1+1 ' i=n+l ' It follows that N-n = a Z b. . |u. - u. | . _P 1,1+n 1 i+n i' p i=2 I, i+n 1 i 1 19 N-l + a E b. . |u. . - u. i=n + l 1 ' 1+1 ' 1+1 x N-l - a E b. . Ju. n l.i+l 1 i i=n+l ' N-n-1 Z b. _ . lu. _ - u. . , l+l, i+n' i+I i+n N-n-1 I 2 i i^ + u. ) , ._, i+l,i+n V| i+l' ' i+n 1 " which was to be proved. The following lemma of Dupont, Kendall and Rachford [1,2] will be used in the proof of theorem 1. Lemma 3 n Let c. be nonnegative for i = 1, . ,.,n. Let E c . be positive. i=l Let e and a., i = 1, 2, . . ., n be complex. Then n E c. i=l " E c.c la.-a I < E c.la.-eT Theorem 1: For < a < 1 the matrix A+B, defined by (3.2), is positive definite, Proof: Let A = (a. ). 20 We have from lemma 2, N-l 2 N-n <(A+B)u,u> = Z - a. . . lu. n -u. I + Z - a. . lu. -u. I ' . , l.i+l 1 l+l i' . -. i,i+n' l+n i' i=l ' i=l ' N + Z (a. .+a. . .,+a. . +a. n .+a. . ) I u.. 1 ._, 1,1 i,i+l i,i+n i-l,i i-n,i'' i 1 N-l N-l 2 + a Z b. . _|u. _-u.| " - a Z b. . -lu.l' _ i,i+i' i+i i 1 . .. 1,1+1' i' i=n+l ' i=n+l ' N-n 2 N-n 2 + a Z b. . lu. -u. I - a Z b. . lu. I . _ 1,1+n 1 l+n i' . _ 1,1+n 1 i' 1=2 ' i=2 ' N-n-1 + Z b. - . (|u. J 2 + lu. | 2 ) . n i+i, l+n ' i+l' ' i+n' i=l ' N-n-1 2 Z b . I u . -u . I . ._, i+l, i+n 1 i+l i+n 1 Since b. . _ = for i = 1, 2,..., n, and b, n = 0, i,i+l ' ' ' ' 1,1+n ' N-l <(A+B)u,u> = Z (a b. -a. . )|u, . n -u, ' 1=1 i,i+l i,i+l ' i+l i N-n + Z (a b. . -a. . )|u. -u. , - i,i+n 1,1+n' 1 i+n i N Z (a. ,+a. . _+a J . +a . _ .+a .)|u. i,i i,i+l i,i+n i-l,i i-n,i ' i 21 N-n-1 E b. . . lu. -u. . . 1+1,1+n 1 l+l i+n 1 i=l " N-n + E (b. . .-Of b. . ) . _ i.i+n-1 i,i+n i=2 " " u. l N-l Z (b. . _ -oc b. . Jlu. i=n+l 1 ^ 1+1 ~ n i,i+l" i By lemma 3, applied to the first two sums, and since b = 0, N-n <(A+B)u,u> > Z " i=l (a b. . -a. . , )(a b. . -a. . ) i,i+l i,i+l i,i+n i,i+n b. OL b. . -a. . _ + a b. . -a. . i+l,i+n i,i+l i,i+l i,i+n i,i+n ' N-l 2 > + E (ab. . _ - a. . -,)|u. . - u. I . „ } i,i+l i,i+l ' i+I i 1 i=N-n+l ' ' i+1 i+n 1 (h.2) N E (a. ,)k + Li (a. ,+a. . n +a. . +a. _ . +a. . ; u. ._ n 1,1 i,i+l i,i+n i-l,i i-n,i ' i N-n + E (b. . .-a b. . )|u. . _ i,i+n-l i,i+n ' i i=2 ' N-l 2 + Z (b. . _ -a b. . . )|u. I . n v i,i+l-n i,i+l ' i 1 i=n+l To prove that <(A+B)u,u> > 0, consider each sum in the last expression separately. The coefficient of term i in the first sum is (cCb. . -,-a. . ^)(cfo. . -a. . % . v i,i+l i,i+l M i,i+n i,i+n) _ fe i ab. . n -a. . , + CCb. . -a. . i+1, i+n* i,i+l i,i+l i,i+n i,i+n * 22 It is necessary to replace the quantities in r. by the expressions which define them. To do so, let matrix row i correspond to point (j,k) of the grid. Then it follows that r _ ,j,k e ,j,k-l ,j+l,k A C ,j,k ,j-l,k ,j,k+l j,k j,k-l j+l,k j,k j-l,k j,k+l J jk j ,k j ,k By (3.3), (_d i k * e i k )(_d i k " f i k } 1 -d . e - d. • f . j,k j,k j,k J,k j,k j,k j,k d ' °' d ' d, e . f . . ( 1+e . . +f . . ) J,k j,k Therefore, by lemma 1, I\ > 0. In order to see that the arithmetic quantities and operations in ve reduction leading to (^.3) are defined, we have for denominator of a b. . • e. - D. . . + a c. . • f . . , j,k j,k-l j+l,k j,k j-l,k hi - j+l,k j,k+l Lies that the sum on the right is positive. the ser- 0. j,k j,k-l v Therefore the fourth and fifth sums are nonnegative. It follows that <(A+B)u,u> > 0. The inequality is strict, since A+B is nonsingular. For, det(A+B) = det(L) = H d. f 0. This completes the proof of Theorem 1. We can now show there exist values of the parameter i which assure convergence of the iteration scheme,. (A+B)u n+1 = (A+B)u n - x(Au n -q). Let u satisfy (A+B)u = (A+B)u - r(Au-q). 2U We have (A+B)e . = (A+B-xA)e v n+1 n where e = u-u . n n Let w , w , ..., w be a complete set of eigenvectors of (A+B) (A+B-tA), orthonormal with respect to the inner product ( v -i> v p)(A t>\ defined by >)( A+B ) = (( A + B ) v 1? v 2 )- Let ^ ^-2'""' X n be the corres P° ndin g eigenvalues . Let N ( s e = E d< n) w. n . _ 1 1 i=l and e , = Zj d; w. . n+1 . _ i i i=l bince (A+B-tA)v. =X.(A+B)w. , . lows that X. = X.((A+B)w.,w. ) = ((A+B-tA)w.,w.) = 1-t(Aw.,w.). i i vv V l vv ' V 1 J v i> i' M and m be the largest and smallest eigenvalues of A respectively, ' hat 1 - tM < X. = 1 - t(Aw. ,./. ) < 1 - -nn. — i I' i — I - |1" • )| < 1 25 if 1 - xM > -1 and 1 - xm < 1 i.e., If 0. < t < n. For this range of x , li U "Vll | A + B