LIBRARY OF THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN Gop- 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 LSEP21 19 DEC 9 -7 Reel L161 — O-1096 6 io. &<+ %JH>^ uiucdcs-r-7 1-U92 #T ' AN ECONOMICAL ALGORITHM FOR THE SOLUTION OF FINITE DIFFERENCE EQUATIONS Martin A. Diamond December 1971 DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN URBANA, ILLINOIS ^TOEUSRARIOg UNIVERSITY C ATV AN ECONOMICAL ALGORITHM FOR THE SOLUTION OF ELLIPTIC DIFFERENCE EQUATIONS INDEPENDENT OF USER- SUPPLIED PARAMETERS BY MARTIN ALLEN DIAMOND B.A. ; Arizona State University, 19&7 M.A. , Arizona State University, 1968 THESIS Submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Computer Science in the Graduate College of the University of Illinois at Urbana-Champaign, 1972 Urbana, Illinois Digitized by the Internet Archive in 2013 http://archive.org/details/economicalalgori492diam AN ECONOMICAL ALGORITHM FOR THE SOLUTION OF ELLIPTIC DIFFERENCE EQUATIONS INDEPENDENT OF USER-SUPPLIED PARAMETERS hy Martin Diamond December 5 , 1971 Department of Computer Science University of Illinois at Urbana-Champaign Urbana, Illinois 61801 This work was supported in part by the National Science Foundation under Grant No. GJ812 and was submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Computer Science, February, 1972. Ill ACKNOWLEDGEMENT The author wishes to express his appreciation and gratitude to Professor P. E. Saylor for his help in editting which has greatly improved the readability of this thesis. The author would also like to thank Professor D. B. Gillies for his financial support. This work was supported in part by the National Science Foundation under Grant No. GJ812. Finally, the author wishes to thank Mrs. Connie Slovak for her aid in preparing the manuscript. IV TABLE OF CONTENTS Page 1. INTRODUCTION 1 2. PRELIMINARIES 6 2.1 Origin Of The Problem 6 2.2 Known Results 10 3. CONVERGENCE CRITERIA 13 3.1 Convergence Criterion When A And A+B Are Arbitrary 13 3.2 The Positive Definite Case 15 3.3 The Use Of A Sequence of t's When A+B Is Positive Definite . . 17 3. If Obtaining A Sequence of t's When A+B Is Not Positive Definite. 20 K. USING A SEQUENCE OF t's WHEN LITTLE IS KNOWN ABOUT THE EIGENVALUES OF (A+B)-^ 31 4.1 Introduction 31 k-2. Calculating An Eigenvector and Eigenvalue of (A+B) _1 A 3^- 4.3 Application Of The Method . . . . 1+1 k.k Using The Recursive Property Of The Tchebychev Polynomials . . 43 4.5 The Algorithm And Its Convergence 53 4.6 Experimental Results 6l 5. AN ESTIMATE OF X 78 min 5.1 The Limiting Values Of ft 78 5. 2 An Estimate Of K 91 min LIST OF REFERENCES 96 VITA 98 1. INTRODUCTION The principal result of this thesis is the development of an algorithm to solve the large sets of linear equations that arise in the solution of elliptic partial differential equations by implicit numerical techniques. The algorithm is an important development for the following reasons: 1. The algorithm is always convergent. 2. The algorithm requires no set of parameters. 3. The algorithm is not dependent upon special properties of the partial differential equation or the matrix of coefficients of the resulting linear equations. !+. The algorithm has proved to be rapidly convergent on a variety of test problems.* 5. The use of this algorithm in conjunction with a factorization due to Y. J. Kim [11] promises to be the first effective iterative method to solve the systems of equations generated in the application of finite element techniques. The algorithm developed in this thesis applied in conjunction with a factorization is the only existing technique with all of these properties. ^Although a priori estimates of the rate of convergence can be established, they depend upon the extreme eigenvalues of the iteration matrix. If accurate approximations of these eigenvalues are available and are used in the algorithm, then the algorithm reduces to the method of Dupont, Kendall and Rachford [5] with optimal parameters. Thus in this case the rate of convergence can be computed as described in [5]. The application of the algorithm, however, assumes no knowledge of the eigenvalues of the iteration matrix. Therefore no a priori bound for the rate of convergence is computable. 2 Since the application of direct methods to the large systems we are considering is either impossible or impractical, iterative methods are usually used. Unfortunately the standard iterative techniques are also deficient. The Jacobi method is quite slow even on some simple problems. The Gauss-Seidel [19] method can be shown to be simply a factor of two better than the Jacobi iteration regardless of the number of equations. Successive-Over-Relaxation [19 ] is dependent on the problem and may become quite slow. The Alternating Direction Implicit iteration [20] is more efficient than either Gauss-Seidel or Successive-Over-Relaxation. However, it requires the selection of a set of iteration parameters to be applied cyclically during the iteration. In general, there is no theoretical basis for selecting a set that will give rapid convergence. Moreover, the method loses its fast convergence when the region of definition of the boundary value problem is irregular. The algorithm developed in this thesis is based on factorization techniques. To elucidate the difference between this algorithm and other algorithms based on factorizations, we shall give a brief explanation of factorization methods. The solution of a system of equations Ax = q by factorization techniques can be separated into two segments. The first step is the definition of an auxiliary matrix, B, which is chosen so that A+B = LU where L and U are sparse lower and upper triangular matrices respectively. In practical computer programs the elements of B are not computed; instead the elements of L and U are computed. Then by using the factors L and U, systems of equations whose matrix of coefficients is A+B can be easily solved. 3 The second part of the method consists of an algorithm which uses an iteration of the form (A+B)x n+1 = (A+B)x n - (Ax n -q), n = 0, 1, 2, ... to generate a sequence of vectors {x } which converges to the solution of Ax = q. The convergence of the iteration is dependent upon the definition of the auxiliary matrix (through the definition of L and U) and a sequence of parameters. One of two types of parameters is employed in the iteration. The first type may be used to alter the auxiliary matrix for each n. There has been no analysis of the resulting iteration, (A+B )x , = (A+B )x - (Ax -q). n n+1 n n n Therefore there is no theoretical basis for the selection of such a set of parameters. The second set of parameters, denoted by x , multiplies the term (Ax -q) in the iteration as follows n (A+B)x _ = (A+B)x -x (Ax -q). ' n+1 n n n ^ In previous algorithms the definition of the parameters x is based upon a %Ay x^ *\Ax x ^ lower bound for min yjrrhn ^ and an upper bound for Max prrfrrl <> I x ^ <{A+B)x,x> ** x £ <{A+B)x^x> ' equivalently they are based on lower and upper bounds for the extreme eigenvalues of (A+B) A. These bounds must be calculated every time the matrix of coefficients A, is altered. Furthermore, the rates of convergence of the previous algorithms are dependent upon the accuracy of the bounds. The new algorithm requires no knowledge of the eigenvalues of (A+B) A. Moreover the rate of convergence is improved by quantities calculated during the execution of the algorithm in such a way that the rate of convergence approaches the optimal rate. k The new algorithm is also superior to algorithms which vary the auxiliary matrix for each iteration. Its primary advantage is that there is a mathematical basis which guarantees convergence. A second advantage is that only one auxiliary matrix is used. Thus only one pair of factors L and U is defined and consequently less work is required. Although no knowledge of the eigenvalues of (A+B)~ A is necessary, information about the eigenvalues is useful. Furthermore, the rate of convergence is governed by the extreme eigenvalues. Therefore the selection of the auxiliary matrix is still important. For matrices which arise in the solution of elliptic partial differential equations by finite difference methods, the symmetric auxiliary matrix defined by Stone [15] (see section 2.2) seems to be the best choice of the known factorizations. We shall now summarize the results of this thesis. In Chapter 3 it is shown that for arbitrary A and A+B the iteration (A+B)x , . = (A+B)x -t (Ax -q) n+1 n n n ^ is convergent with t a constant t iff (if and only if) the real part of each of -1 t the eigenvalues of A B is greater than p - 1. By modifying this result for positive definite A and Hermitian B, it is shown that the iteration is convergent with t =1 iff A+2B is positive definite. This strengthens the long established result that when A is positive definite and B is Hermitian, the iteration is convergent with % = 1 if A+2B is positive definite, see Wachspress [20]. Next it is shown that if A and A+B are positive definite then the iteration is convergent with a constant value t iff 2 < T < Max X. ((A+B) _1 A) where \ ((A+B)" A) is an eigenvalue of (A+B) -1 A. 5 In section 3.3 the analysis by Dupont, Kendall and Rachford [5] concerning the selection of a best sequence of x ' s to use when the eigenvalues of (A+B) A are real is presented. The results are extended in section 3-h establishing a best sequence of x 's to use when the n eigenvalues of (A+B) A are complex. In Chapter k an algorithm based on the iteration (A+B)x . * (A+B)x -x (Ax -q) n+1 n n n *' is developed. In Chapter 5 an approximation of the smallest eigenvalue of (A+B Q )~ A, where B is the symmetric auxiliary matrix defined by Stone , is b b derived for some special cases. This approximation can be used in starting the algorithm of Chapter >+ when B Q is the auxiliary matrix. b 2. PEELIMINARIES. Notation: Throughout this thesis, X. (M) will denote an eigenvalue of the matrix M, X . (M) will denote the eigenvalue of minimum algebraic value, \^ (M) will denote the eigenvalue of largest algebraic value, and p(M) will be used to denote the spectral radius of M, i.e., p(M) = Max \. (M) . i - 1 I All matrices used in this thesis, even those specified as general, are real matrices. 2.1 Origin Of The Problem . Consider the Dirchlet problem j - ^ (a i (x) 1^ - ox; (a 2 (x) si: ) = q(x) > x eD ' u(x) = {x) > for ( 2.f) D where x = (x ,x p ),D is the interior of a compact region with boundary c^D and where the differential operator in (2.1) is elliptic, that is a, (x) and a p (x) are strictly positive on D. A common method of solving (2.1) is to approximate the differential operator with a difference operator. A system of linear equations results and the problem is then to find a solution to the system of equations. The system is derived as follows. Assume D lies in the first quadrant; cover D with two families of parallel lines, x 1 = jh j = 1, . . . , n ± and x 2 = kh k = 1, . . . , n p . Denote the points of intersection (jh, kh) by (j,k). The points (j,k) are called the mesh points or gridpoints . Let bj), be the polygon formed by 7 joining those mesh points (j,k) for which one of (j+l,k) or (j,k±l) is not in D. Similarly let D be the open and connected set hounded by dD . At each gridpoint of D the differential operator is approximated by a difference operator and the equations are in one-to-one correspondence with the gridpoints. For simplicity of notation we shall consider the region D to be the unit square bounded by the x-. and x p axes and the lines x =1 and x p =l. Let x = (x 1}1 , r^g, ..., x^, x 2 A , ..., x n ^) T where n ± = n? _ = n. We obtain the matrix equation Ax = q (2.2) defined by (Ax) . * q. . , where the subscript corresponds to gridpoint (j^k) and q is derived from q(jh,kh) and the boundary conditions. The difference operator can be chosen so that matrix A is diagonally dominant and positive definite with a... > and a. , . < for i / j. One method for deriving such an A is to approximate 3 — (a. (x) 3 — u(x)) by i i h" {a i (x+ g-e ) [u(x)-u(x+he i ) ] + a ± (x - ^e ± ) [u(x)-u(x-he i ) ] } where e and e are unit vectors along the x n and x p axes respectively. This defines A by (Ax)., n B. , X. , , + D. . x. , , + E . . + D. ,, . x.,_ , + B. , .x. ,,. , y j,k j,k j,k-l j,k j-l,k j,k o+l, k j+l,k j,k+l j,k+l j^k ss 1, , . . t n, where B • k = -a 2 (jh,(k- g-)h), Dj k = -a x ((j- §)h,kh), (2.5) and E. = a 2 (jh,(k- |)h) + a 2 (jh,(k + |)h + a^j- §)h,kh) + a x {U + |)h,kh). When 1 = 1 or k=l we define D. . = or B. = respectively. The following properties of A and the quantities defined in (2.3) will be used below. ^,k - E d,K*a,k + B j,k x j,*.i + V x d-i,* + ViaV,* + B j,k + i x j, k+ i , and E. > -2(B. . + D ) with equality holding if B. , ^ and D. . / J^k ~~ J^k j,k J^-K J^-K for j=l, 2, . . . , n and k=l, 2, . . . , n. The subscript («j,k) refers to the gridpoint (jh,kh). In this thesis we shall study the solution of equations Ax = q by factorization techniques. As sescribed in the introduction, the idea of a factorization method is to add an auxiliary matrix, B, to A so that A+B factors into known sparse matrices L and U. An algorithm based on the iteration (A+B)x . = (A+B)x -t (Ax -q) (2.5) n+1 ' 11 n n is then used to compute a sequence of vectors {x ) which converges to the solution x. The parameters 1 are defined by the algorithm. 9 A factorization is defined by the auxiliary matrix or the factors L and U. We shall give particular attention to the symmetric factorization due to Stone [15 ]. This factorization which can be applied to any matrix A with the properties listed in {2.1+) , is defined as follows. where j,k ' j,k X j,k-l c j,k x j-l,k H j,k x j,k ' (Ux) . , a x . , + e. ,x. r , , + t. ,X. , -, > y j,k j,k j,k j+l,k j,k j,k+l b . , = B . , -ac . . _ f . . , , , j,k j,k j,k-l j-l,k-l (2.6) c . 1 = D . , -ab e J,k j,k j-l,k j-l,k-l j,k j,k d,k-l j,k j-l,k ' j,k j,k-l J-l,k-l +ab e j-l,k j-l,k-l , d. e . s D., n . -ab . . e . . . , j,k j,k j+l,k j,k j,k-l (2.7a) and j,k j,k "" j,k+l" °j,k j-l,k The quantities B. . , D. , and E. . are the elements of A defined in (2.4), j,k j,k j,k y ' and the parameter a is used to vary the auxiliary matrix if Stone's algorithm [16] based on the iteration (2.5) is used. We are interested in algorithms that don 't vary the auxiliary matrix and will assume a = 1. In addition to the quantities defined in (2.7a) we define C i k " b i k e i k 1 and G i k = C i k f i 1 k ' (2 ' Tb) Since there should be no confusion as to whether the auxiliary matrix is defined by (2.6) and (2.7) or is arbitrary,, B will be used to denote either. 10 With L and U of the form (2.6), A+B is given by ((a+bW^ = ((MWj,* - V^A-i + t J , k e j,k.i x j + i,k-i + c j,k x j-l,k + j,k H 3,k j,k-l "* c j,k e j-l,k' ,x j,k + d j,k e j,k x j+l,k + c j,k B j-l,k x j-l,k+l + d j,k f j,k x j,k+l * (2.8) 2.2 Known Results The following definitions are standard. T \T Definition 1.1 : If x=(x 1 , . . . , X ) and y~(y,, ..., y ) are n-component n * * vectors then ^c,y^ = T\ x.y. where y. denotes the complex conjugate of L=l X X X y. . -NX,y> is called the innerproduct of x and y. jL, 2 Definition 1.2 : ||x||, the L 2 norm of x, is defined by ||x|| = ||x || _ = . 2 If A is positive definite then the A- norm of x, ||x||., is defined as 2 . Dupont, Kendall and Rachford [Jj,] have proved the following lemma. Lemma 2.1 : Let A and A+B be positive definite. Suppose there exist positive numbers e and e p such that for all nonzero x we have <(A + B)x,x> e [e l' e 2 ] (2 ' 9) Then for k <0, M d .i,* >0 > (a) -1< e < 0, (e) -1< f < 0, J >& and (f) 1+e. ,+f.. . > 0. Bracha has also shown in [1] that for the factorization defined in (2.6) and (2.7) A+B is positive definite when A satisfies (2.4). A more concise proof is included here. Lemma 2.2 : If A satisfies (2.k) and B is defined by (2.6) and (2.7) then A+B is positive definite. Proof: From (2.7) it is seen that b . , =d.. , f . , n and c. n = d. - . e. - . . v " j,k j,k-l j,k-l j,k j-l,k j-l,k Thus (2.6) implies (L) = d • •, • ( u ) where row r of A+B corresponds to r,s j ,k. s,r 12 T the gridpoint (j,k) and s =1, . . . , r. We have L = U diag(d ), where j )£■ T diagfd. . ) is the diagonal matrix with d. , along the main diagonal and U J,k ot* is the transpose of U. Using property (c) from above, we have A+B = LU = U T diag(d 1 / 2 ) diag(d 1 / 2 )U which is clearly positive definite. In addition Bracha proved the following theorem [2]. Theorem 2.1 ; Let A and A+B he the matrices defined in (1.6) and (l. 10) respectively. Let e. and f be defined hy (1.10). Let p ,j,k = -e. ; -f. ; J** - x > z > - •• > n > J,k j,k and p = min p . Then ^ 1 Max ttttsR — r < r* x? fo <(A+B)x,a> - 1 _ J^ ««j • "^Ax.x)" . 1 and min yl . *S g > x/0 <(A+B)x,x> " l + ki where k is a positive constant. 3- CONVERGENCE CRITERIA In this chapter necessary and sufficient conditions are derived for the convergence of (A+B)x ,,= (A+B)x -t (Ax -q). The conditions we n+1 n n n obtain do not seem to he know in the literature. The conditions are first developed for arbitrary A and B. For this case it is shown that the iteration is convergent with a constant x iff i is selected from a certain region in the complex plane. This result is modified for the case when A is positive definite and further modifed for the case when both A and A+B are positive definite. When A and A+B are positive definite, a best constant value of r is established. Section ^.k establishes a best sequence of t 's to use when the eigenvalues of (A+B)" A are complex. The derivation follows that of the known results presented in Section 3-3;> which establish a best sequence of t 's when the eigenvalues of (A+B) A are real. n 3.1 Convergence Criterion When A and A+B Are Arbitrary In the following two lemmas convergence criteria are derived in terms of t for (A+B)x =(A+B)x -t(Ax -q). (3-1) n+1 n n Lemma 3.1: Iteration (3.1) is convergent iff Re(X. (A B)) > _ - 1. th Proof: Let x be the exact solution of Ax=q. Let the error in the n iterate be given by e = x-x . Subtracting (3-l) from n n (A+B)x= (A+B)x-T(Ax-q) yields (A+B) (x-X n+1 ) = (A+B) (x-X Q )-T(A(x-x n )) or e = (A+B)" (A+B-rA)e . n+1 n 13 Ik The matrix (A+B)" (A+B-tA) is the iteration matrix of the iteration (3.1)- It is well known that the iteration is convergent iff p( (A+B)" 1 (A+B-tA))<1. (3.2) The iteration matrix can he rewritten as (A+B)" 1 (A+B- xA) = (I-t(A+B)" 1 A) = (I-t(A- 1 (A+B))- 1 ) = (I-Td+A^B)- 1 ), where I is the identity matrix. Thus (3.2) can he rewritten as ^.(I-Td+A^B)" 1 ) eS (0,0; 1) where 5(u>v;R) = {(u,v):(u-x) + (v-y) (5-3) ^((I+A^B)" 1 ) e 5 (-i,0; ±) , or \. (l+A~ B) has real part greater than x/2. It now follows immediately that (3«l) is convergent iff Re(\ (A _1 B))> 1- -1. (3. k) Corollary 3. 1.1: If A is positive definite and B is Hermitian, then (3-1) is convergent iff A+B- ^-A is positive definite. Proof: If A and B are Hermitian, each \.(A~ B) is real and (3.4) becomes ?,.(A _1 B)> -L-l . This is equivalent to K. (l(l- i)+A B )X> ■ But 1(1- 1. )+A _1 B = A" 1 (A(1- 1 ) +B) = A _1 (A+B- -£-A). 15 Thus (3. if) is equivalent to X. (A~ (A+B- -*-A)) > 0. (3.5) Since A -1 (A+B- ^-A) is similar to A" 1 / 2 (A+B ^-A)A _1 / 2 , (3.5) may "be restated as K A ' (A+B- -*-A)A~ x,x) > >0, for all nonzero vectors x. Therefore < (A+B J-A)A" ' x,A~ ' x > > 0, for all nonzero x. Letting y=A~ x we obtain K (A+B- -—-A)y,y ^ > 0, for all nonzero y. completes the proof. Hence (3. l) is convergent iff A+B- -^A is positive definite, and this For x = 1 the corollary implies (3«l) converges iff A+2B is positive definite. This strenghtens the previously known fact that a sufficient condition for the convergence of (3.1) is that A+2B he positive definite. 3.2 The Positive Definite Case We now modify the above results for an A+B which is positive definite, as it is in the factorization due to Dupont, Kendall and Rachford [ 5] and the symmetric factorization due to Stone [ 15] . Lemma 3.2: If A and A+B are positive definite then (3-1) is convergent 2 iff < x < . Proof: (3«l) is convergent iff I ?,.(I-t(A+B) _1 A) I < 1. Since \.((A+B)~ A) > , this condition is equivalent to 1 " t Niax ((a+b) " 1a)> - 1 2 \ , and this completes the proof. 3--J-A 1 W (A+B A > 16 We can make a further refinement by deriving a "best value for x as seen in the next lemma. L emma 3.3: If A and A+B are positive definite then 2 X, = * " A m ((A + B)- 1 A) + X mln ((A + B)- 1 A) is the best choice for t in (3.1). -1, Proof: The iteration matrix of (3.1) is T = (I-t(A+B)~ A). Thus, since all the eigenvalues of (A+B)~ A are positive, the eigenvalues of T are decreasing as a function of x. Clearly the spectral radius of T is minimize' with respect to 1, when | ^-ji-C^) I = I ^mAX^t^ I ' ± ' e ' > I l-T\. /rA ..((A+B)' 1 A) I = 1-t\ . ((A+B) _1 A) , 1 MAX. ' ' mm ' or 2 = t[ \ MAX ((A+B)- 1 A)+\ min ((A+B)" 1 A)]. Therefore t, = 2/[\,, / .„((A+B)" A) ] is the best choice of x as claimed. b m 2X . ((A+B)"^) Moreover for x = x ,p(T ) = 1- min 1 , '" , '" r " T^.MA+Byh^X . ((A+B)"^) MAX v ' ' mm There is a strong relation between Lemmas 3. 2 and 3-3 an( i Lemma 2.1. Lemma 3.2 states a necessary and sufficient condition while Lemma 2.1 states only a sufficient condition. The hypothesis of Lemma 2.1, however, includes a supposition that positive numbers e and e exist such that for all nonzero x < Ax,x > e [e, ,e ] . < (A+B)x,x > X d This is certanily valid when both A and A+B are positive definite. In fact e, = K»- '(A) mm 1 W (A+B) 17 and X MAX (A ) e 2 ~ [a+bJ iran are sufficient. The next lemma shows precisely the connection of Lemmas 2.1 and J. 2. Lemma 3. k: Let A and A+B be positive definite. Then < AX.X > r n < (A+B)x,x > € Le l ;6 2 J for all nonzero x iff \.((A+B _1 A)) e [e ± ,e 2 ] . Proof: The conclusion may be restated as 1 < ay,y > X . ((A+B) A) = Min j ,. . <— . c . min vv , < (A+B)y,y > (3.6) and y v 2_ < ay,y > \ ((A+B)' A) - Max . ™ y^o < (A+B)y,y > Now since (A+B)~ A is similar to A •'< (A+B) A '^ which is Hermitian (3.7) -1a ^ _, ^V2 r A o.^-l 4 V2 > _ m,„ \_. m ( (A+B) A) - \_. „ (A ' (A+B) A ' ) = Min x/o mm man x . 1/2 , ,_i 1/2 - -1/2 1/2 Letting x = LA (A+B) A J Ay and us i n g "the selfadjointness of [A 1 / 2 (A+B)"^ 2 j" 1 / 2 and A 1 / 2 (3-7) becomes \ . ((A+B)"^) = min Mnn < Ay,y > $0 < (A+B)y,y > ' Similarly ^((A+B)'^) = Max < (A+B^y > This completes the proof. 3.3 The Use Of A Sequence Of t's When A+B Is Positive Definite There are no new results presented in this section. However, Section J.k and Chapter k are developed from the results described here. 18 In this section a "best" sequence of t 's to be used in the iteration (A+B)x _ = (A+B)x -T n (Ax n -q) , is derived for the case when both A and A+B are positive definite. A "best" sequence is one that minimizes the spectral radius of the iteration matrix. As in Lemma 3.1 the error vector e is multiplied at each step by (I-t (A+B)" A). Hence N-l , 6 N = n io ^- T n (A+B) ' A) e ' (5 ' 8) To make the iteration converge as rapidly as possible, the t 's should be chosen to minimize the spectral radius of N-l 1 "n ■ nSo [I -Vi ( < A+B »" A ' ! and since \.((A+B)~ A) is real, • \.(^L) is real. The problem of minimizing the \. (HL ) thus reduces to selecting a sequence of t 's that minimizes the maximum absolute value of N-l P N (X) =n& (l " T n X) (5 ' 9) on the interval X . ((A+B)" A) ^ X ^ \,,._-((A+B)" A). No normalization mm ' J MAX is necessary since it is clear that all P„ of the form (3.9) satisfy P (o)=l. By classical Tchebychev analysis it is seen that the desired polynomial is given by p (y) _ V (? W +: SlAX - 2X) / (X MAX- X min )) , (3.IO) N v mm MAX ' TdAX mm' where T is the Tchebychev polynomial of degree N and \ . and K^jr ^ en0 ^ e \ da ((A«)" 1 A>ana ^((A+l of P H ( X ) on ft^,^] is ^• min ((A4B)" A) and ?\^^((A+B)~ A) respectively. The maximum absolute value 19 |~ / Nnin +X MAX \ ~~| L \ MAX mm/ J . (3-H) Note that if \ . ^ < \. flY then the maximum of P must he greater than or equal to 1, since P (0)=1. However a and A+B positive definite implies that X^ and \^^ are positive and therefore assures the reduction of the error vector when {t } is defined by (3.9) and (3.10) There are two ways of implementing a sequence of x ' s so that (3.9) is the polynomial given in (3. 10). In the first method an N is selected and the sequence x n=0,l. . . ,N-1 is computed as the reciprocal of the roots of T / min MAX \ . The sequence is then used repeatedly ^MAX'Nnin until the norm of the residual vector, r = Ax -q, is sufficiently small. ' n n J A sequence of x's derived as ahove is called a Tchebychev sequence of length N. The other method of implementing the sequence, which was first suggested by Stiefel[l|], uses the recursive property of the Tchebychev polynomials. If this technique is used then every x is the result of a Tchebychev sequence of length n. if the first method is used then x is the result of a Tchebychev sequence of length n only when n=N. The Stiefel sequence is defined as follows. Suppose x is the result of using a Tchebychev sequence of length k. Define x as where ( \lAX " WW*' \^> " ' with 2 (A+B)" 1 (q-AxJ 5x =- ~ V ; V4 y Nnin + ""MAX 20 for Y _ HlAX +X min *"MAX ~^min then x. is the result of using a Tchebychev sequence of length k+1. K+l In general \ . and \.,.„ are not known. Tchebychev minimization ° mm MAX is typically implemented by finding an interval (e ,e g ) containing (>, min ,> mx ) and using e ± and e 2 instead of X^ and k^ respectively. In place of (3»9) we have, e.+e^-2X N-l p„(x) = n _ (i-T x) = And the Stiefel sequence is defined by H e 2" e l J t r 6 l +e 2 1 JM e„-e„ Vi = \ + 5 \ > W k (y) T k 1 (y) 5 *k - (e 2 - ei )T k+1 (y) (A+B) ( ^V + T^y) Vl , S x^ = (A+B) _1 (q-Ax n ) e ! +e 2 °* where y = (e +e )/(e 2 -e ) . 3.1+ Obtaining A Sequence Of t's When A+B Is Not Positive Definite In this section we shall derive a best sequence of t's to be used in (A+B)x _ = (A+B)x -t (Ax -q) n+1 ' n n n ^ 21 -1, when (A+B) A has complex eigenvalues. Before proceeding, it is necessary to define what "best" means. The iteration (A+B)x , , - (A+B)x -t (Ax -q) n+1 n n n N-l n= is most rapidly convergent when p( r| II n (l-T (A+B)~ A)) is minimized. This is accomplished when N-l -1, Max n (1-t X. (A+B) A) ' n i i n= (3-13) is minimized. Since not all of the X. are known, the t could he selected 1 n to minimize Max zeR V z > N-l where P, T (z) = II (l-x z) N a».0 b-ik) and R is some region containing all of the eigenvalues of (A+B This concept is formalized with the following definition. -1 A. Definition 3- !• A sequence of parameters 1 n = 0,1,«««,N-1 used in iteration (2.5) is called a best sequence of length N in a region R containing { \. ( (A+B) _1 A)} if N-l Max zeR n (1-t z) n= = min lT n J n=0 Max zeR N-l n (i-x z) n n n= More conveniently, we will say the best sequence in R. In the remainder of this section a best sequence of parameters will be derived for a circular region and then for an elliptic region. 22 To begin the derivation the following definitions are required. Definition 3.2: Let g (z) be any complex polynomial of degree n, then M \z (r )1 denotes the maximum modulus function of g (z) for |z|g . r , i.e., M [g (r)]= Max n izU r g n (z) Definition 3. 3: S is the set of all polynomials g of degree n such that g (1) = 1. to n Definition 3. k: M [g (r,d)] = Max d- z I g r g n (z) Definition 3. 3: S is the set of all polynomials g of degree n such that g (0) = 1. n * The polynomials which are required in (3-1^) belong to S . However n we will need the following theorem, due to E.H. Zarantonello (see Varga [19]), which is stated for polynomials in S . Therefore, after stating the Theorem it will be modified in Theorem 3.2 to apply to polynomials in S . n Theorem 3.1: For all r such that g r ^ 1 , min {M[g (r)]] g eS n to n n = r for all positive integers n. 23 Theorem 3.2: For all real r and complex d such that ^ r < |d|, mi g eS & n n » (M[g n (r,d)], « -£, for all positive integers n. Proof: Let r and d be such that ^r < |d| and suppose there exists some g eS such that M fg (r,d)l< n n L n J Let h (z)=g (d-dz). Then n ' n r T n > Max d- z I ^ r g n (z) Max d- z I ^ I r M) z g Max |h (z) I = M[h ( -£■ )]. 1 n n d r Since h (l) = g (o) = 1, the above inequality contradicts Theorem 3.1, and M [g (r,d)]sL«j- for all g is S . On the other hand, if n n S n (z) = ("^r) n thenM^[g n (r,d)] = d-z \ n theorem follows and the conclusion of the n Corollary 3-2.1: The polynomial P (z) = ( — ■— — j is the polynomial in "X- -& -#- S satisfying M [P (r,d)] = min„ M [g (r,d)] for all integers n> and n n _"*" n g eS n n all r and d such that ^ r ^ I d I . Proof : This follows immediately from the proof of the theorem. Corollary 3. 2.1 can be used to determine a best sequence of t 's to use on a disc D = {z:|z-d| < r, where < r < |d|] . That the 2k \.((A+B)~ A lie in such a disc is seen as follows 1 The idea of choosing a sequence of t ' s is to speed convergence of iteration (2.5). It is reasonable to assume there is a constant value t for which the iteration is convergent. From this and (3-3) it is seen that for some t > 0, X ± ((A+B)"^) e B(-~ ,05-i-)- {z: ~~ -1 * -i-} - (3-15) Consequently Corollary 3.2.1 can he used to determine a best sequence of t ' s in a disc as follows, n Theorem 3.3: If the eigenvalues of (A+B)" A satisfy (3*15) then the best sequence of t 's to use in (2.5) is t = x a constant. Moreover t = — — where d is the center of the smallest disc containing all X. ((A+B)"^) . Proof: Since no confusion will arise, X. will he used to denote 1 X.((A+B)~ A) in this proof. Let s(d,0;r) he the smallest disc containing all the X. . That the smallest disc is of this form is seen as follows. 1 Since the complex eigenvalues occur in complex conjugate pairs, the smallest disc containing them has a real center, (d,0). In addition (3.15) implies g r < d. Now from Theorem 3.2 and Corollary 3. 2.1 it is seen that of all polynomials in S , the minimax polynomial on S(d,0;rj 13 .H v*> = -¥ and M*[P N (r,d)] = N r d It is apparent that either at least k of the X. lie on the 25 circumference as in Figure 3- la, at least 2 as in Figure 3.1b or 3 with one on the real axis and two as in Figure 2.1a. Figure 3.1 Let s(&+€, ,6pjr') "be any other disc containing all the X.. The polynomial g^S which is the minimax polynomial on this disc satisfies M*[q N (r',(d+e i; e 2 ))] [(a^r+Cegn "2-j N (3-17) The Theorem will be established if it is shown that the equantity in (3.I7) is greater than the quantity in (3.16). To establish this we shall assume, without loss of generality, that e, and e p are non-negative. If this were not the case, the following argument would be changed only in the selection of the X. of Figure 3.1. Since 5(d+e ,e p ;r ' ) contains all the X. it in particular contains one which is as X of Figure 3.1a or Figure 3. lb, i.e.;, X = (u,-v) where u > and v ^ 0. Therefore, 26 (r , ) 2 ^(d.+e 1 -u) 2 +(e 2 +v) = (d+€ 1 ) 2 +u 2 +e 2 2 +v 2 +2e 2 v-2u(d+€ 1 ) (d+€ 1 ) 2 +u 2 +€ 2 2 +v 2 -2u[(d+e 1 ) 2 +e 2 2 ]~2" = [((d +ei ) 2 +e 2 2 )"2".u] 2 +v 2 . Hence the quantity in (3. IT) is greater than {[(d+ £l ) 2 +e 2 2 ] 2 -u} 2 +v 2 2" (**!>' + e. N 2 Thus M*[q N (r',(d+e 1 +G 2 ))] > 2 2 where e_=[(d+e_ ) +e ] 2 -d > 0. 3 J- <= /* \2 2 (d+e -u) +v 3 (a+e 3 )< And M*[q N (rS(d+e 1 ,e 2 ))] > (d-u) 2 +v 2 N = M [P N (d,r)]. The last inequality is derived from d > u and the Law of Sines. We have proved that the polynomial P M (z) of (3.16) has the minimum maximum on the smallest disc containing each "K- > where the minimum is taken over all polynomials in S ; and this maximum is less than the minimum maximum of all polynomials in S on any other disc containing the X.. The t should, therefore be chosen so that n N-l . .N N i.e. ; t = — — . This completes the proof. 27 If more is known about the X. ((A+B)" A) then some more useful results concerning a best sequence of parameters can be obtained. In particular, if the \. ((A+B)" A) are contained in a ellipse then the best sequence is a Tchebychev sequence as in the real case. Before proving this we need to recall some terminology 2 2 Definition 3.6: In the ellipse ^ X "^ + — «Z|J — = 1 with a 2 s b 2 , a fe 2 JL_ a is called the semi -major axis , e = (l- — - ) 2 is the eccentricity b and (p,q) is the center . A modification of the following result due to Clayton, see Wrigley [23], will be used in determining a best sequence of x 's when the X.((A+B) A) are contained in an ellipse. Theorem 3-k: The polynomial in S„ which, over the ellipse of semi -major axis a < 1, eccentricity e and center at the origin, has the smallest maximum modulous value is given by P N (z) = T N (z/ae)/T N (l/ae) . (3-18) If the center of the ellipse, still assumed to lie within the unit circle, is translated to (d,0) then the minimax polynomial is P N (z) = T N ((z-d)/ae)/T N ((l-d)/ae). (3-19) Remark : As e -> the ellipse becomes a circle and, as expected, the 28 polynomial P^ of (3. 18) approaches z N . This is easily seen, for when e is small , T ( z / ae ) = -|- [exp(Ncosh" 1 (z/ae))+exp(-Ncosh' (z/ae))] = -i- exp(Ncosh" (z/ae)) 1 -I N = 4" [IT +(Wae) 2 -D-] . 2 N - 1 (z/ae) N . Hence lim P(z) - ^"fo"^ ■ z N . In the first equality the fact e-*Q iN 2 (l/ae) that T (z) - cosh(Ncosh' 1 z) is valid throughout the complex plane was used (see for example Wrigley [23] ). * Since we are interested in polynomials in S^ rather than S^ , Theorem J>,k is rephrased as follows. The orem 3.5: Let E be the ellipse with semi -major axis a, eccentricity N e and center (d,0) where a < d. The polynomial in S which has the smallest maximum modulous value on E is given by m I d-Z \ T N (-^~ ) P N (z) = T (d/ae) ' (3- 20 ) N Proof: First note that P (z) of (3-20) is in S* . Now suppose there is some q N (z) in S with smaller maximum modulous than P^ on E. Then Max zeE' q N (d-dz) Max < zeE' 29 P N (d-dz) (3.21) where E* is the ellipse with semi -major axis a/d , eccentricity e and center (0,0). But p H (d.dz) . Jalzkil ■ e d Thus (3.21) contradicts Theorem 3-h and the conclusion of the theorem follows . The following corollary shows that polynomial P N (z) in (3-20) can be put in the form of the polynomial (3. 12). Thus the best sequence of t 's in an ellipse can be implemented in the same way that the best sequence on an interval (e ,e ) is implemented. Corollary 3. 5. 1: Let E be as in Theorem 3» 5 and suppose each \j((A+B) A) is in E. The best sequence of parameters to use in (2.5) on E is the sequence {x } for which N-l T N ((e 1 +e 2 -2z)/(e 2 -e 1 )) n i d-V } =P N (2 > = T N ((e 1+ e 2 )/(e 2 - ei )) (3-22; where e = d-ae and e = d+ae. Proof: Let E be as in Theorem 3-5- As noted when considering a disc, the existence of such an E which contains each X. is guaranteed by (3.I5). Let e = d-ae and e = d+ae. The best sequence of parameters on E yeilds a polynomial 30 V>-JSo (1 -v) = rrx- * (3.23) IT ae ; But e +e = 2d and e -e = 2ae so that (3.23) is g^V^f ) /e 2+ e r 2z N ' / V e l \ . / e 2 +e l e 2" e i J N V 6 2 " e i as claimed. Remark : As the ellipse degnerates to a line, e -* 1 so that e -* d-a and e -* d+a which is exactly the result of the previous section. k. USING A SEQUENCE OF t's WHEN LITTLE IS KNOWN ABOUT THE EIGENVALUES OF (A+B) _1 A lf.1 Introduction In this chapter we will develop an always convergent algorithm to solve Ax=q. The algorithm is similar to the method presented in section 3- 3 hut is superior since it does not require any knowledge of the eigenvalues of (A+B) _1 A. The theoretical basis is derived in sections k.2, k.J> and k.k. The algorithm is stated in 4.5 and some experimental results are given in 4.6. The sequence of t 's in the iteration (A+B)x n+1 = (A + B)x n - T n (Ax n -q) {kml) is defined during the execution of the algorithm. As in the methods presented in Chapter 3, the sequence is dependent upon an interval. Instead of using a constant interval as the previous algorithms do; the algorithm of this chapter generates a sequence of intervals (a., b.) which converge to the optimal interval. The following is a brief description of the algorithm. It is shown that 5 n = V X n-l' which is automatically calculated in the practical implementation of {k.l) , can he made to approach an eigenvector of (A+B) _1 A corresponding to an extreme eigenvalue. This approximate eigenvector is used to compute an approximation to the corresponding eigenvalue. The interval (a^ Td ± ) is then improved by replacing one endpoint with the approximate eigenvalue. 31 32 The algorithm to be developed can he applied to any iteration of the form (l+.l) in which A and A+B are positive definite. If A is derived from the discretization of an elliptic partial differential equation as described in Chapter 2 and A+B is defined by one of the factorizations due to Kim [8], Stone p_5 ], or Dupont, Kendall and Rachford [5] then A and A+B will be positive definite. Thus the algorithm can be applied to these factorizations. It is assumed in the sequel that A and A+B are positive definite. Additional Notation ; Throughout this chapter the subscripts and super- scripts will denote the following: i corresponds to the interval (a., b.). n and N correspond to the number of iterations. k and K correspond to the dimension of the matrix A. Thus the vector v calculated on the n iteration using the interval (a., b.) will be denoted by v . This vector will be denoted by v when the interval J n n used is not being stressed. The algorithm will use vectors which are calculated during execution of computer programs written to solve the set of equations Ax=q by using a factorization and iteration (4.1). Therefore, we shall begin with a description of the practical utilization of the methods presented in Chapter 3. Two methods are presented in Section 3.3. In the first a vector Xq, an interval (a,b) and an integer N are selected. The iterates x . n-1, . . . , N, are then defined by (A+B)x n+1 = (A+B)x n - -r n (Ax n -q) (4.1) 33 where the T are selected so that p w (x) = (i-T x) = Nl b -» ; . N nio n T (£*) In the second method the iterates are defined by x , , = x +5x . n+1 n n 5x „ )Tn(y) , ( A+ B)- 1 (q-^ ) + ^44s^ „ (4.2) 85 o=Jf( a+b )"VaS ), where y = 4; r and x_ is the initial vector. (b-aj Throughout the remainder of this chapter x will denote a vector defined by (4.1) and x will denote a vector defined by (k.2). The vectors in the following definitions are calculated in the practical use of (4.1). Definition 4.1: The n residual is the vector r = q-Ax . n ^ n Definition 4.2 : The n incremental vector is 8 = x -x , . Using these definitions the iteration (A+B)x ,_ = (A+B)x - t (Ax -q) x J n+1 v ' n n n can be rewritten as (A+B)6 =t ,r .. ' n n-1 n-1 This equation is solved by replacing A+B by LU and calculating 5 from Lz = t ,r . n-1 n-1 and U5 = z. n The next iterate is then x = x -, + P> n n-1 u n The vectors x are similarly calculated. Instead of 6 and r n n n however,, we define A „ A Definition k.J> : ? n = q - ^ n - Iteration (U-. 2) is performed hy computing -1a z = L r , n ' *n + l - U " lz > and , A ^Tn(y) a Vl^) a ^.5) 6x n = (b-a)T n+1 (y) Vl + T^yJ 5 Vl n > 1. The vector Sx» is computed from t -1a z =L r Q and 5x n = — -=- U z . a+b For every n > x , is calculated from A A A x , , = x + 8x n+1 n n Note that 6 defined in (h-.~*>) satisfies n ,-1a The importance of this vector will he shown in Section 4.4. 4.2 Calculating An Eigenvector And Eigenvalue Of (A+B)~ A In this section it is proved that if y is any vector and if % =P H ((A+B)- 1 A)y Wh6re T I a + b-2X v« - ,*:;■' N Ib-a 35 and a> X . or b < \ ..._, then y. T approaches an eigenvector of (A+B)~ A. mm MAX N It will also be proved that W N <(A+B)y N ,y > approaches the corresponding eigenvalue. (4.6) -1, In the next section 6„ , is shown to satisfy 5 . = P ((A+B)~ A)S,. From this it follows that 8 . is approaching an eigenvector of (A+B)~ as N increases. We shall begin with some motivation. Suppose that y n is any initial vector and that N-l 'N (l-x*(A+B)" 1 A)y = P*((A+B) _1 A)y , where P*(x) = II (1-t*X). (k.l) n=0 N-l n n=0 The object is to select the x such that y approaches one of the extreme eigenvalues of (A+B)~ A. If K y = Y a v w v * where the w , k=l, ..., K, form a complete set k=l of eigenvectors of (A+B)" A, then y is given by K 'N L °k ] k=l N ^((A+B)"^) w. = P. N ^((A+B)"^) I a ^[^((A^B)-^ )] k=l To make y most rapidly approach w , {t )in (k-7) should be selected so that Max P;[x k ((A + B)- 1 AJ P;[X 1 ((A + B)" 1 A)] 36 is minimized. Since \ , . .., \ K are unknown, we minimize Max \ K < x<*. 2 *J« F* [V^Br'A)] (^.8) where \, denotes \ k ((A+B) _1 A). -1, Note that y converges to an eigenvector of (A+B)" A provided n lim Max — : — £-p r — r N^oo k=2,...,K P N I ^((A+B)" A)J = , a.9) or lim N-*oo Max \ < x< \ 2 4 q__(y) *B' m IV 'T *W for all y /[-1,1]. 37 Corollary U.l-1 : The polynomial P*(X) = I a+b-2X ] Nl b+a / N I b-a I satisfies Max a< X X > . . . > X v i > K = x • > °- MAX 1 2 — — K-l K min K Let y n = V a,w be any vector such that c* a /: and define the sequence k=l of t . n=0, .... N-l. such that a+b-2X N-l b-a p (x) = n d-x n x) = n=0 t(— ] Nlb-a/ (U-.10) where a> X . or b < X...„ and min MAX then either N v mm 7 ^ VHlAX^ If % = V (A+B) " lA)y o lim 1 N v MAX' y N = Vl or (4.11) S! r ;' . y N = VK 38 The convergence of y N to w 1 or w K is most rapid when a = \^ n and b = L or when a = X^ and b = X^, respectively. And y^ converges to w only if b< K m ^ and converges to w K only if a > X^^. Proof: Suppose a> X . or b < \ ,t ** mm MAX' Suppose vw < "VSlAr This implies b< \,,.^ and MAX 'H (X MAX } > ? N ( V k = 1, . .. , K. From the definition of y Q and y^ we have i * w — t % = L. " ^T^MAX' i" k VW " k a+b K -»* ■E a, "Nl b-a j k=l T. N a+b -2 :> *MAX = a i W l + . J b-a a. a+b -2^ •N b-a w. v ra,»i*TJ a + b - a W l " I b-a a+b-2\. 1 k ) T Nl b-a j St? .-, k m /a+b-2\ >^[a,b] ^ T N SlAx] w. b-a (4.12) Applying Corollary if. 1.1 to the terms of the first sum we see that as N increases these terms are decreasing at an optimal rate. For the terms of the second sum we have from (if. 12) that 1 < a+b-2\. k |b-a| < a+b-2\. MAX b-al 59 Thus cosh -1 a+b-2\, I b-a - cosh -1 a+b-2\. b-a MAX < 0. Using the inequalities and the relation T (X) = [exp(N cosh _1 X) + exp(-N cosh _1 X) ]/2 W for X > 1 yields lim N a+b-2\ k b-a = Tj a+b -^MAx' b-a Thus lim Tr-77 r y* T = a - 1 w n which establishes the first half of (k.ll). The second half of (k.ll) is similarly derived assuming P, T (\ . ) > P, T (X„„„) N min J N MAX That the convergence is most rapid when a = \ . and b = \_ or min 2 a = \ , and b = ^ MAY follows directly from Corollary 4.1.1. Finally we note that y, , y p , ... converges to w only if V^mat > V\iin } to w_. only if a > X . K J mm , i.e., only if b < X This completes the proof. MAX . Also y 1 , j 2 , ... converges Any approximation to an eigenvector of (A+B)" A yields a corresponding eigenvalue approximation. This is derived from a theorem in [20], presented here for convenience. Theorem L. 5 [20]: If A and B are positive definite then the eigenvectors, x, , and the corresponding eigenvalues, \, f of the generalized eigenvalue problem, satisfy the following: 1, \ > 0; 2, x k V. =0 j /k. A^ = X kBXk , ko Proof: Since B is positive definite there exists an orthogonal R such that R BR = diag(& ) = D , where P k is an eigenvalue of B. (A-XB) = RD(D" 1 R T ARD" 1 -XI)DR T . Hence det(A-XB) = (det R) 2 (det D) det(P-Xl), where P = D -1 R ARD - . Since (det R) =1 and (det D) = H(p, ) > 0, the zeros of det(A-XB) are those of det(P-Xl). These are the eigenvalues of a symmetric matrix and are therefore real. The matrix P has a complete set of eigenvectors z, which may be taken to be orthogonal. Hence D _1 R T ARD _1 z =\ z ; and A(RD _1 z k )=?^RDz k = \ k RD(DR T RD" 1 )z k = X, B(BD" z. ). Thus x. = RD~ z, is an eigenvector of the generalized eigenvalue problem, corresponding to X,. Since {z } is orthogonal and RD" is nonsingular, {x, } forms a complete set of real eigenvectors. Furthermore m m m rn mm m = z k z. = (DEx,) DBx. =x k RBDRx. = x^ Bx . , so that {x } is orthogonal with respect to B. Finally, X, > since < Ax, , x, > = <\, Bx, , x> This completes the proof. Observe that this theorem applied to A and A+B implies that the eigenvectors , w , of (A+B)~ A are orthogonal with respect to (A+B). The next lemma defines an approximate eigenvalue. K Lemma l+.l : Suppose A+B and A are positive definite. Let y = *£, e \r\ k=d where w k is an eigenvector of (A+B)~ A corresponding to X , . Let U = NYjAy > . If y is approximately equal to w with errors € in the J K kl direction w , , and e » e ^ 2 k j= J _, then \i = X with error of order J Proof : Assume J = 1. The remaining cases are similar. Assume {w } is orthonormal with respect to (A+B). Let M ~ U.13) K Substituting ]T e w for y and (A+B)(A+B)" A for A in (4.13) yields k=l "k k W = K K 2_ € v w v T (A+B)(A+B)' 1 A ; , e,w. k=i "k k I J-l J J K K k=d j=l € .W . 1 T Now using the facts that (A+B)" Aw. = \.w. and w (A+B)w- = 5 , where J J J .k J j ^k 5. . is Kronecker delta, we obtain Al k k 1 k^2 k k X £ 2 " £ 2/ 2 Z \ 1+ L v € i k=l k k=2 k X This completes the proof. In the case A+B=I the above lemma reduces to a result of Hestenes [21]. k*3 Application Of The Method To compute an eigenvalue using the above method it is necessary to determine a sequence of vectors y , ..., y such that y = P((A+B)~ A)y where V X) = T H a+b-2X b-a "N a+b b-a k2 The next theorem shows that 8 , n. = 1> . .., N + l is the desired sequence of vectors. Since s is calculated when iteration (4.1) is used to solve n Ax = q, no additional work is required to obtain the approximate eigenvector. Theorem 4.4 ; Let x b e an y vector and define x , n = 1, 2, ... , N+l by- iteration (4.1). If the sequence of fr )used in the iteration is such that NJ . T A^m) P ^ (X) = nQo (1 " ^ = N I b-a i where a > \ . or b < \„ AV and mm MAX P N (X min ) / P N^MAX' ' then Vi = Vr x N approaches an eigenvector of (A+B)~ A as N increases. Proof: Since 8 n = x n - x^ = (i - ^(M^A)*^ + Vl^'W > ,-1 -1, 5 -t ,(A+B)"q=-T _(A+B) Ax . n n-l v ' H n-l v ' n-1 (4.U) and 5 ^.-t (A+B)" 1 q = -t (A+B) _1 Ax . n+l n v / m, n ^ / n (4.15) Subtracting (t /t ) times (4.15) from (4.14) yields n-1' n 5-(x ,/t )8 ^ n =t . (A+B)' 1 A(x -x n ) =t n (A+B) _1 AS , n n-1' n' n+l n-l v ' v n n-1' n-l v ' n ' k -l. which implies (I - x^A+B)" A)8 n = ( T n _l/ T n^ 5 n+1 * Hence (T n / T n-l )(l "Vl (A+B) " lA)5 n = 5 n+l' And N-1 V T IL (I " T n (A+B )" A > 5 1 = Vl ' Since the sequence of t , n =0, . . . , N-1, is used cyclically, T = t +3 and thus N-l 5 N+1 = (I-t (A+B)- A) 5l . (4.l6) n=0 The conclusion follows from Theorem 4-2. An approximation to the corresponding eigenvalue can be calculated, according to Lemma 4.1, by Mm " ' ( ' 4.4 Using The Recursive Property Of The Tchebychev Polynomials From the theory developed in Sections 4-2 and 4.3 it is possible to use iteration (4.1) to solve Ax=q and simultaneously compute an approximate eigenvector of (A+B)~ A. The process is started by selecting N and using (4.1) to compute the sequence of vectors x , n=l, . . . , N+l. The resulting approximate eigenvector is 5^ _ = x , -x_, . In the next theorem it is shown that if x is computed from (4.2) then % n+1 = (A + B)- 1 (q-A^ n ) (4.18) is an approximate eigenvector for every n =1, 2, ... . Recall that § ,, /x , , -x . Note that § n is computed when (4.2) is used to solve n+l n+l n n+l Ax=q and thus the approximate eigenvector is again obtained with no additional work. Theorem 4. 5 : Let x be any vector and define the sequence of vectors x , n =1, 2, ..., by iteration (4.2). Assume the interval (a,b) used in the kk iteration satisfies (a/b) i> (X^, >- mx )- ^t r n = q-Ax n - It follows that 5 = (A+B)" r approaches an eigenvector of (A+B)~ A as n increases. Furthermore 5 approaches the eigenvector corresponding to X^^ only if b < X and approaches the eigenvector corresponding to X . only if JXlA-X. iiij.ii a > X min Proof: Define t n = 0, ... N - 1, such that a + b-2x iio n t. r a+b [ a+b I » L b-a J and let t = t / , __v for n > N. Let x_ = x n and define n n(mod NJ — W N W N N / \-l/ \ x = x n + 5 where 8 = r _ (A+B) (q-Ax ,), n n-1 n n n-l v ' VM - n-1' The vectors £ defined be (4.2) have the property that n a n x = x for every n. n n ° Therefore ^ = (A+B) _1 (q-Ax n ) = (A+B) _1 (q-AxJJ) = -^(A+B) _1 (q-Ax^). Ttf-1 But x n = t^ and T n (A+B) _1 (q-Ax n ) = 5 n , . Hence n n v ' VH zr n+1 'n From (ij..l6) we obtain n J n+1 n u n+l ' : o n-1 l.\JB J n+1 n" IL (l - T d (A+B) A K (4.19) 3 = Since n w l 8?= -4 T* (A + B)" 1 (q-Ax^) ,-1 = (A + B)- ± (q-Ax ) = 5\ , *5 (4.19) can be rewritten as % n+1 = H (I-t^A+B)- 1 ^ . (1..20) j = J The conclusion follows from Theorem 4.2. The approximate eigenvalue is now given by We now show that no matrix- vector multiplications are necessary to compute the quantities (A+B)s _ and Ag -, which are needed in the calculation of ju . The vector (A+B)^ = ^ n and therefore is available. That A5 - can be calculated without a matrix-vector multiplication is proved in the next theorem. More precisely, it is shown that AS , may be calculated from known quantities with only two vector additions, at the expense of additional storage to hold (£ , - £ ). Theorem 4.6: If £ is the n iterate defined by (4.2) and $ = q-Ax , then - — n n ^ n' T n+2 (y)(b-a) A K+1 = 4T n+1 (y) L^n+l"^n+2^ q n+l^ r n" r n+l^J ' where y = (a+b)/(b-a) as in (4.2) and q^ = (T n (y))/(T n+2 (y)). 4T (y) Proof: From (4.2) we have £ ,_ = x , n + % ,. / w , — r + q ^-, 5x , where v ' n+2 n+1 n+1 T pKYiVo-a.) n n+l n' y = (a+b)/(b-a), q . =(T (y))/(T. p (y)) and & is given in (4.2). Therefore ■n+2 + T_ n (y) ^n+2 = » -^n + 2 = Vl " A Vl V "(y)(b-a) " A vA 14-6 Replacing q A§ n by (y)(b-a) r A Vl= T^73 [ \ . and approaches an eigenvector corresponding to X„„^ only 1 mm jr -o MAX if b. < \,, AV . Thus ju. can approximate \ . only if a. > X . and jti. can 1 MAX 1 r * mm J 1 man 1 approximate T^y only if b. < Tw™-. Furthermore we are developing a method which will generate a sequence of intervals (a. ,b. ) that converge to (X . A„ AV .). 1' 1 mm' MAX Therefore we must be able to improve the left endpoint of some intervals and the right endpoint of others. This implies that the condition a. > X . and 1 mm hi b. < \j, AV must be satisfied for each i. If this condition holds it can be concluded that ju. approximates \ . if ju. < a. and approximates X if ju. > b.. l—i Note that the condition a. > X . and b. < \.... r can be replaced l mm i MAX * by a. > X . and b. < \.. AV , i.e., (a.,b.)c: (\ . \„.„). With this relaxed J i — mm i — MAX i' i mm' MAX condition, the endpoints of each (a. , b. ) which are not already optimal can be improved. The next lemma shows that if (a. ,b. ) cz (\ X.„„) and i' l mm' MAX' (a. , ,b. , ) is defined by replacing either a. or b. with a. , then l+l ' l+l 11 i' (a. ,_ ,b. , _ ) c (\ . , \., AV ). Therefore if (a-, ,b n ) c (\ . , \,„ AV ) we may v i+I ' i+I mirr MAX v 1' l y v mm' MAX ^ assume (a. ,b. ) c (\ . , \., A „) for each i. How to select (a n ,b, ) c (\ . .\ rliV ) v i' l mm' MAX 1' 1' v min' MAX is discussed in section k. 5* Lemma 1+.2 ; If ju = ^vTlpS 5. where A and A+B are positive definite, then HW (a+b >~ 1a >>W (a+b >" 1a >] • Proof : This is a restatement of (3.6) Throughout the remainder of this chapter P (X) will denote the polynomial ., —„. ,, T ( a :T b - 2X )/T (^L.), a 1
  • k8 3. If n ± < a ± define a ±+1 = H ± and b ±+1 = b ± It u. > b i define a i+1 = a ± and b ±+1 = fi ± In either case let x^ 1 = x^J and return to step 1. We now consider the case a. < ju. < b.. If the approximate eigenvalue is repeatedly in the interval (a.,b.) then one of the following conditions holds: Either 1. (a.,b.) ^ (^ Xmx ) or 2. P_(\ . ) = IT mm' P ]T X MAX' Remark: We can choose (a, ,b, ) c (\ . _, \ MAV ) and thus assure (a. ,b. ) c: (\ . , \,„__). However we will present a heuristic method of v x f i nmr MAX selecting the initial interval that does not guarantee (a., ,b ) c (\ . ,\ ) (see Section 4-5). For this reason it is necessary to consider case 1 above. In either case the vectors x will converge to the solution of Ax=q. In the first case the new method reduces to the method of Chapter 3« The initial interval can be chosen so that if this case occurs, convergence is rapid (see section k-5)- In the second case convergence holds since < P (X . ) < 1 for all n (see Fig. 4.1). Nevertheless, there is a n man v ° ' ' complication. We shall consider two subcases. 2a. a. = X . and b. = X,,,,. . l min l MAX 2b. a . £ X . and b. ^ \., A __ l ' mm i ' MAX Consider 2a. If a. = X . and b. = \„.„ then K T (\ . ) = l mm l MAX N v nmr P N^ X MAX^ Thus if the interval is improved so that (a.,b.) = (X . , \„ AV ), then v i' i' v min' MAX" P N (X min ) V^MAX^ In this case the interval is nearly optimal and should be used in the subsequent iterations to compute AjL-., ^+2* * * * k-9 In case 2b the situation may be as pictured below. p„(«) Convergence of fx_ T l will be slow since the interval (a. ,b. ) is a crude approximation to (X . , X ). Furthermore the vectors 5 do not approach an eigenvector of (A+B)" A. Instead they approach a linear combination of Note that merely increasing the degree of P, i.e., performing more iterations, will not guarantee the eigenvectors corresponding to X . and X r ,.„ to *• ~o mm MAX convergence of { § ) to an eigenvector of (A+B)~ A. It follows that the interval cannot be improved using the method described above. We shall now describe a technique which will force convergence of Performing one iteration of (4.1) with [&-) when Pjx . ) = 1 n J I\P mm' P N^MAX' t = l/b. will diminish the components in the error vector and the approximate eigenvector corresponding to the larger eigenvalues. Therefore it will increase the relative size of the components corresponding to the smallest eigenvalues. This idea can be stated more rigorously as follows. If N+l iterations have been done using the interval (a. ,b. ) and & =P ((A+B)~ A)s where P„(X . ) N min P N^MAX^ , then one iteration with x = l/b. yields Va = W ((A+B ^ A ^ * ( \ 50 where P^X) = (l -£) PjCX). Furthermore P^O^) # We could next compute r^ +1 = q-Ax^.^ V2 ■ (A+B )"\i • i < V2, A V2> and improve the interval if u^ / (a i ,b ). However, if u p e (a. ,T>. ) we assume that case 2a has occurred and thus want to continue iterating using (a. ,b.). There are two methods of using the interval for the next W steps. The first is to define xj +1 =^ +1 and compute x^ +1 n=l, ... N using (a i+1 , b ±+1 ) = (a^K). The second is to continue using (a. ,b.) and thus define £LpJ ^i\r-V """ ^PN+l We take the second alternative. Furthermore the calculation of jli__, and the corresponding attempt to improve the interval is omitted. Instead ^L., is defined using T=l/b. and ^.p; • • • .> ^pivr+p a 2,6 immediately- computed. Then i Nb 2N+l ,Ab 2N+1 is computed and used to improve the interval if ^otvu-i e ( a -.»b.) . This procedure is chosen for the following reasons. An approximate eigenvalue u is computed only once every N iterations instead of twice. Next, we know that case 2a will eventually occur in every application of this method and this case is not benefitted by the extra calculation of uL . Finally case 2h happens only infrequently. (This can be established by fixing X . and considering the range of values of \,,._ r that cause case 2b to occur. ) MAX 51 A technical problem now remains. The vectors x , ... x^ cannot be computed by simply restarting the iteration using (4.2) since SxV. is not defined. The following will elucidate the problem. Computing x p requires £ , and 6X . First,, 5x is computed and then x 1 „ = & , + 5x , , . If iteration (4.2) is used to calculate x^", ... x^l n+2 n+1 n+1 1' TT then the last vector 5x to be calculated is 8$ , . If £ is computed by using t = l/b. in iteration (4.1) then it is necessary to also define 5x7, before resuming iteration (4.2). The vector S^L. must be defined so that the iterates & ? , ^vr,,^ •••> converge to x and so that oVr +2 ^ %4-v '•'> converge to an eigenvector of (A+B)~ A. Convergence of [zl } is established by proving, in the following theorem, that the error vectors e = x - £ n > N+3 satisfy e =P ((A+B)~ A)e_, ' n n — ^ J n n 0' where again P (X) = (l -y—)? -, (X). Then using the fact that the same i polynomial PUA+B)-^) = n (I-t (A+B)"^) 2=0 J appears in and e = P ((A+B) 1 A)e. n n U=V< A+B ) A )^i (compare (3.8) and (4.20)) we conclude U = p > +b >~ 1a >\ > and thus 8 approaches an eigenvector of (A+B) A. Theorem 4.8 : Let x n be any initial vector and let x n=l, . . . , N and 5x n n=0, ..., N-l be defined by (4.2). Let x^ +1 = ^ N + £^+1 > 52 where and let Vi = ^b)"Va\) e£ N = (i-tCa+b)- 1 ^^^ . If x , n > N+2 is generated by a a , A x . = x +5x , n+1 n n ^n - (b-a)T n ( y ) (A.B)"^^) + ^- ^ , ft.*) where y = (a+b)/(b-a) . Then • a+ - b _ 2 x 'n 1 "b-a C n for n > N+2. e n = (l-T(A+B)" 1 A)P n ((A+B)- 1 A)e , where P q (x) = n l g ^ Proof : The proof will he by induction one , j=l, 2, ... . The error vectors e = x-x n=l, .... N satisfy e = P ((A+B)" A)e^ where n n ' ' ^ n n T p (x) =. -S I a+b-2X b-a T n I b-a Since T +1 (X) = 2XT (X) -T (x) for all n, and y = (a+b)/(b-a) it follows that T n (y)U T (y) K 4-1 ( X ) = - rp / \/v \ p (x) + P (x) + m , . (P (X)-P . (X)). n+1' T (y)(b-a) n v n v ' T (y) v n v y n-l v " ni n+1 (If. 23) Consequently the polynomial P (X) = (1-tX)P 1 (X) satisfies * . . T n (y)^ * * T (y) ^ 1 (y)(b-a) W X > + W X > + T^ £* W = - T° (y)(b-a) 4i»Kkl« + T^UJ^tl^K^' ^ 53 -1, Let Vl = \ + ^N+l ' Where Vl = T ( A+B )~ (ci-A^f) and let -1. gftjj - (I-tCCa+B)-^)^^. Obviously e N+1 = P^ +1 ((A+B) _1 A)e . (4.25) Now suppose x n = N+2, . . . , is defined by (if.. 22). We have A _ A A e N+2 ~ X " X N+2 ~ X " ^+1 " 6X N+1 h-T (y) T (y) But (q-Ax N+1 ) =A(x-x N+1 ) = A(l - T(A+B) _± A)e N and = (l-r(A + B)- 1 A)(e N _ 1 -e N ) Therefore (4.26) can be rewritten as e N+2 = (I- T (A+B) _1 A) L W N (y) x T N1 (y) e N" (b-ajT^y) (A+B) Ae N + T^T?) <•*-*-!> With (4.23) and e = P ((A+B) X A)e this becomes e N+2 = (l- T (A+B)" 1 A)P N+1 ((A+B)" 1 A)e () = P* +2 ( (A+B) _1 A)e . (l|.27) From (I4.25) and (4.27) e = P* ((A+B) _1 A)e for j=l and j =2 . The general inductive step follows trivially from (4.24) and 5x = e . -e 4.5 The Algorithm And Its Convergence In this section a method of selecting the initial interval (a. ,b. ) will be developed. A precise definition of the algorithm to use (4.2) and simultaneously improve the endpoints a. and t>. is given. Convergence is also shown. In the next section some empirical results -will be given. Effective application of the methods of the previous section requires the selection of an initial interval (a^t^) c (^j^Wy). Lemma 4-3 (below) implies that the only interval guaranteed to satisfy this condition for every problem is the point interval a, =b. =1. Although the methods of the previous section work when this point interval is used as the initial interval, the resulting iteration may start inefficiently. Indeed, the error vectors e may even grow until the interval is improved. Therefore the iterations performed after the interval is improved must annihilate an error vector that is larger than the initial error vector. We shall now consider an alternate method of selection to overcome this inefficiency. Our goal is to develop a method that requires no estimate of \ . or \„., r . Since \ . and \...„ may be arbitrarily close to unity (see mm MAX mm MAX * Lemma k-3) it is impossible to develop a general a priori method of selecting an initial interval (a n ,b, ) c (\ . ,\,,.„) which avoids the inefficiency of 1' 1 mm' MAX a point interval. We can, however, present a heuristic method of selecting a nondegenerate initial interval that yields a "good" rate of convergence. The interval may fail to satisfy (a-, ,b, ) C (\ . ,X,,.__). However, it will in * l 1 oan MAX ' this case yield a rate of convergence defined by the user to be "good". Our method of selection is based on the following facts. First, Nnin - 1 — \iAX* An<3 " secon< ^-y> the L -norm of the error vector is reduced by a factor o < — ; — -—- , after N iterations when (a,b) -> (\ . . \ t .J) is ■" T ( a+b ] v ' ' x mm' MAX^ Nl b-a ) used as a constant interval in iteration (J*. 2). We will show that if 55 (l) a, and b are chosen from small intervals to the left and right of unity respectively. and ri- a J where the user decided that a reduction of the error by a factor less than or equal to € is "good", then the rate of convergence will be either the "good" rate selected by the user or a rate which is optimal or nearly optimal. There are three cases: (i) (a-, ,b, ) cz (\ . X.....) , \ / ^1*1' v min' MAX ' (ii) (a, ,b n ) 3 (\ . *.,,.„) and v ' x 1' 1 / v min' MAX (iii) eu, < \ . and b, < \,. AV or a, > \ . and b n > \„ n „ . v ' 1 mm 1 MAX 1 mm 1 MAX If (a, ,b n ) cz (x . , \ rir .„) then the intervals are improved by the methods of v 1' l y v min' MAX the previous section, and approach the optimal interval. Thus the rate of convergence approaches its optimal value. It should be noted that since (a-. ,b.. ) c (x . X...-U-) we will have V T.' l y v min' MAX' < Tj a l +b ll WNnin^MAX I b n -a, / I X^.-^-K . I \ 1 1 1 \ MAX min / Thus the optimal rate of convergence for this problem is slower than the "good" rate selected by the programmer. Also note that since X . Z 1 < X.,,.,,, J * & min - - MAX' selecting a, and b from small intervals to the left and right of unity is a good heuristic for satisfying (a, ,b_, ) a (\ X r ,. v ). Comments on exact J to v 1' 1 / min' MAX numbers to use as a, and b are given below. 56 If fa i"b, ) =3 (\ . , X„ A v) then the method never changes the v 1' 1 min MAX interval and the rate of convergence is the "good" rate selected by the programmer. Finally, if a^X^ and b 1 >\ m or a^X^ and \ X . and b > \* A y The left endpoints, a., will be decreased to approximately X . . The right endpoint s, however , will never be improved. The rate of convergence will approach the rate obtained using the constant interval (X . , b_ ). But b, was chosen in a small interval to the right of mm 11 unity, say b 1 = 1 + 5. In addition X^ >1 so that | ^ -^ = h ± -X^^, <6. Thus, intuitively, although b, >X^ b, is not a gross overestimate of X and convergence is not significantly slower than the optimal rate. Remarks : The values of a, and b, are selected to yield a "good" rate of convergence as a precaution for cases in which (a, ,b. ) d (X . , X r ,. v ) . The v 1' 1' T v min-' MAX' only purpose in selecting this "good" rate is to prevent the selection of too large an interval with an unacceptable rate. It is impossible to create a rate of convergence that is better than the optimal rate of the problem. That is, selecting (a, ,b_ ) such that a N (a 1 ,b 1 ) = * V a J is very small cannot improve the rate of convergence. In fact, as o' N (a 1 ,b 1 ) goes to zero the interval becomes a point, which is the exact condition we are attempting to avoid. 57 Our experience is that a., e (l/k, 8/l0) and b e (1.25, 3) yield "reasonable" initial intervals. The interval (1/3, 5/2) has proved to be effective for problems in which nothing is known about X . or X mm MAX In Chapter 5 a value for a, in some special cases is derived. In these cases the interval (a,, 5/2) has been effective. The following lemma establishes the basis for the heuristic method of selecting the initial interval. Lemma 1+ . 3 : Let A be as in (2.6) and let B be given by (2.7). Then 1 is an eigenvalue of (A+B)~ A. Furthermore for every e > there is some A such that 1-X . ((A+B) - ^) < e and \ MAV ( (A+Br^) -1 < e . min MAX Proof : If B is given by (2.7) then the first row and column are null. Hence B has a zero eigenvalue. Let w be the corresponding eigenvector. Then (A+B)w = Aw, i.e., (A+B)~ Aw ~ w. To see that X . and X.... r can be arbitrarily close to unity min MAX J J consider the following. As the elements B. in matrix A (2. If) tend to J >*■ zero, elements b. and f in L and U (2.6) also tend to zero. In the limiting case B. , = 0, b =0 and f . , = 0. Moreover LU = A, that is, B is the null matrix. In this case (A+B)~ A = A~ A =1. Hence every eigenvalue of (A+B) A is unity. Since the eigenvalues vary continuously with the elements of the matrix, for a non-null B with arbitrarily small elements X . and \. r . v may be arbitrarily close to unity, min MAX Remark : This lemma is valid for the factorizations due to Kim [11] and Bupont, Kendall and Rachford [5]. We shall now give a statement of the algorithm. 58 The Algorithm . 1. Select an initial interval (sUjlO according to one of the following rules. (a) Define (a.. ,b ) according to the method presented in the first part of this section. For example a, € (.25, .75) . and b. e (1.5; 3.0)- In practice, a 1 = l/j5, b_ = 5/2 have been effective. (t>) If possible use the approximation of a, developed in Chapter 5, Theorem 5-3* Select b.. as in (a). (c) Let a, = b = 1. Comment : Although (c) yields a convergent and practical- it eration, it is recommended that (a) or (b) be used to define the initial interval. 2. Select an initial vector x and some e>0 which is a bound on the maximum acceptable error, measured by the norm of the residual. 3. Select an integer N which is the number of iterations to perform before attempting to correct the current interval. Comment: As N is increased the convergence of a. and b. to \ . 11 min and ^- MA y is improved but the convergence of x to x is retarded if the current values of a. and b. are crude approximations to \ . 11 ^ min and a.^.. In experiments using the algorithm, values of N = 5, 6, 7, 8, 9, and 10 were equally effective. k. Execute N steps of iteration (k-2) generating Ai /\i , Ai X N ' 8 N-1 and ^-1 ' where i corresponds to (a.,b.). 59 5. Iflig All q-Ax II < e then stop, otherwise compute 11, in. 1 <^*K.i$-i> u. = 6. i. If ii. < a. , define a . _ = min (min a. , /i. ) and b = b. i+1 i+1 j h. , do one iteration of (k.l) with x = l/ju. . Define a. , = b. and "b. , = ju. . Add one to i and go to k. l+l i l+l i ° Comment: The definition of the next interval when u.. > b. i i is different from the definition when ju. < a. for the l i following reason. Components of the error vector corresponding to eigenvalues X, > b. will have been amplified if 3 X> > 1, where the i indicates correspondence to the interval (a.,b.). When \i. > b. , the next interval is defined to emphasize the annihilation of the dominant components of the error. The components which became dominant when ju. < a. , on the other hand, can ii become dominant only if they are not diminished as rapidly as other components. They are never amplified since P H (x) < 1 for all x e (0,a.). Therefore when p.. < a. ' l ii it is unnecessary to concentrate on the reduction of the dominant components. If u. e [a. .b. ] , do one iteration with t = l/b. . Add one to i i' l ' ' l i and do N iterations according to ('4.22), then go to 5- Comment : This requires essentially no new code since (>4.22) is identical to (1+.2) except that the subscripts on the 6o right side of the equations are decreased. Furthermore the code required to form (A+B) _1 (q-A^) and (i, t(A+B)~ A)a£J before proceeding with (4.22) already exists in the code to perform (k.2). The purpose of this step is to correct the condition as described preceeding Theorem if. 8. vw P N^MAX^ Convergence : The convergence of the algorithm is established by- considering three types of initial intervals (a,,b-). 1. If (a., ,b, ) ^ Ck . , \,,., r ) then the algorithm reduces to the method v 1' 1' mm mAX of Dupont, Kendall and Rachford E?] and is therefore convergent. 2. If (a,,bj i> (X . , \., AV ) and P.A . ) = P_ T (\. AV ) then Theorem If. 5 guarantees the convergence of the intervals (a. ,b. ) to an interval 1' 1 (a.b) such that a =X . or b = X.,. v . Furthermore if (a, ,K')c(\ . A,,.,,) v ' ' mm MAX 1 1 mm' MAX' then the (a. ,b. ) approach (\ . , X,,.,,) and the rate of convergence ^ V \ * mm' MAX of x to x approaches the optimal rate. 3. The case P,_(\ . ) = x N mm' < P (\ . ) < 1. N v min 7 P N^MAX^ is trivially disposed of since Work : The following considerations show that if N" is the number defined in step 3 of the algorithm the amount of additional computation per iteration required by the new method as compared to the method of Dupont, Kendall and Rachford, (if. 2), is approximately y^ times the number of computations in one iteration of their method. There are two differences in computation of the methods. The first is that our algorithm computes an approximate eigenvalue on every th N iteration. This requires 2 vector additions and the accumulation of 2 6l inner products. Therefore the number of extra calculations required per iteration to compute the approximate eigenvalue is (3k additions + 2K multiplications )/N, where K is the dimension of A. The second difference is that after every N iteration the new algorithm employs a new interval. This creates a difference in the work performed because the first iteration with each interval is a special case. In general 4 T (y) , T ,(y) where § = (A+B)~ £ , and y = (a.+b. )/(b. -a. ). When n=0, however, n n-1 1 1 ' 1 1 ' ' 1 1 is computed. In every case x" , is defined as x" + S& . Thus the first * n+1 n n iteration with each new interval requires K fewer additions and K fewer multiplications. The additional calculations due to these differences is (2K additions + K multiplications )/N The number of computations in one iteration of (J+.l) or (if. 2) without the calculation of the approximate eigenvalue is (I3K additions + UK multiplications). Considering only the multiplications we see that the additional calculations per iteration of the new algorithm 1 UN calculations per iteration of 0.1) or (if. 2) if. 6 Experimental Results Four distinct problems will be studied in this section. Stone [l6] and Taranto [17 ] have also made studies using the type of problems 62 considered here. They compared the amount of work required to solve the system of equations Ax=q using first Stone's factorizations and then using AD I. Taranto also compared Stone's symmetric factorization to Stone's nonsymmetric factorization and studied the effect of changing the parameters a and x (see Chapter 2 for definitions of a and t). The following tests make two comparisons. First the number of iterations and the amount of work required to solve systems of equations using Stone's symmetric factorization and the algorithm developed in this chapter is compared with the number of iterations and the work required when the problems are solved using Stone's factorization with a single best parameter. Next the comparisons are made for problems solved by Stone's symmetric factorization with the new algorithm and both of Stone's factorizations used with Stone's algorithm to vary the auxiliary matrix, [15, 16]. Problems a and c described below are identical to problems studied by Stone and Taranto. Problem b is similar to a problem they studied but differs in the ratio of the functions a, (x) and a~(x). (see Chapter 2, where a, and a. are defined. ) This change was made because the problem studied by Stone and Taranto converges so rapidly using the methods under consideration that no significant difference in the relative rates of convergence results. The last problem studied involves random numbers in the matrix A. Since the random numbers generated in the problem below are not the same as those generated in the previous studies, no comparison of results is attempted. In fact the form of the problem is changed slightly since the comparison is impossible. In the problem studied below the functions a-, and a 2 of (2.1) are random throughout the region. In the problem studied 63 by Stone and Taranto the functions were constant on subregions but were otherwise random. Each problem is solved on the unit square, {(x, ,x p ) : 0< x. < 1 i=l,2), with zero boundary condition, i.e., Y(x) of (2.1) is identically zero. The solution x is chosen to be x. = cos(jHh).cos(kHh) j,k=l, ..., 30, where h =1/(30+1). Similar tests were run with cosine replaced by sine. The results were almost identical and therefore are not presented. The vector q of (2.2) is computed from The initial vector, x n , is chosen to be the zero vector. Problem a . This problem is the model problem, defined by letting a, (x) and a p (x) be identically -1. The by the analysis of Chapter 5 is \. and a (x) be identically -1. The initial estimate of X . , suggested 2 rain Problem b . For this problem, called the generalized model problem, a, (x) = -I/9 and a~(x) = -1. The initial estimate of X . as suggested 1 ' 2 mm in Chapter 5 is 5/8 • Problem c . This problem is defined by a heterogeneous region containing homogeneous subregions see Figure 1+.2. In region A, the ratio of a (x) to a p (x) is unity with a, (x) = -1. In region B, this ratio is 10 with a, (x) = -1; in C it is 0.1 with a, (x) = -l/lO and in D it is 1 with a..(x) = -l/lO. The theory of Chapter 5 does not yield an estimate of X . in this case. 6k Problem d. In this problem random values of a. (x) are chosen from the set {z : 0m. 3. Stone's symmetric factorization implemented with a sequence of a' s as in 2. The results for methods 2 and 3 are taken from the study done by Taranto [17]. The graphs in Figures 4. 7 and 4.8 have the. reduction of the error is plotted against the number of iterations. The reduction of the error is calculated as l- « / HeJI where e~ is the error in the initial vector and e is the error in the n n iterate. The graphs in Figures 4-9 and 4.10 have the reduction of the error plotted against the work performed. One unit of work is taken to he the amount required to do one iteration of (4.2). It is assumed that cc all the auxiliary matrices, B _, are calculated only once and that the amount of work needed to calculate the elements of the factor L and Tj is equivalent to the amount required to do one iteration. 73 _ -1 -2 -3 -4 -5 t- -6 -7 -8 -9 -10 -11 -12 1- New Algorithm with Initial Interval (0.5,2.5) 2- Stones Nonsymmetric Factorization with a sequence of 3- Stones Symmetric Factorization with a sequence of — I 1 14 21 -+- 28 35 NUMBER OF ITERATIONS Figure 4.7. Model Problem Ik O IxJ O O I— o 10 10 10 10 10 10 10 10 10 _ -1 -2 -3 -4 -5 -6 -7 E- -8 — — -in » 10 10 10 10 -10 -11 -12 14 21 28 NUMBER OF ITERATIONS Figure 4.8. Heterogeneous with Homogeneous Subreglons 75 10 10 10 10 10 10 _ -1 -2 -3 -4 -5 S10" 6 on oz -7 uj-JO ' -8 o10 iio o 10 £l0" 11 10 -12 -9 tr- -10 14 21 28 35 WORK REQUIRED Figure 4.9. Model Problem 76 10 10 10 10 10 _ -1 10 -2 -3 -4 -5 at 10 o Q& 7 -6 — LU 10 10 10 -8 -9 F- 5 io- 10 a 10 10 V -11 - ^ ^ - — i 1 h 7 14 21 28 35 WORK REQUIRED Figure 4.10. Heterogeneous with Homogeneous Subregions 77 Conclusions : The tests comparing the algorithm developed in this chapter to the use of a single best f clearly indicate that fewer iterations are required when the new algorithm is used. Since the tests were performed with N of step 3 in the algorithm equal to 6, the work estimate of the last section implies that the additional work required by the new algorithm is approximately l/66 of the work performed by the iteration which doesn't change the interval. Thus the same graphs may be used to conclude that the new algorithm requires less work than the iteration using a single best parameter. The graphs of Figures k>l and 1+.8 indicate that in some cases the required number of iterations of the algorithm developed is greater than the number required by Stone's algorithm. However, Figures k>9 and J+.IO indicate less work is required by the new algorithm. 5. AN ESTIMATE OF \ . mm Throughout this chapter the following notation will be used. B will denote the auxiliary matrix defined using (2.6) and (2.7) with a . 1. The quantities ^ = ^((A+B^A) , b^ fe , c.^ , d j%k , e. , , f , , C. . and G. . will be as in (2.7). B. . , D. . and j,k ' j,k ' j,k j,k j,k j,k E will be as in (2. If). When B. . , D and E. , have constant j,k J^-K j, K J,.K values, except B. = and D. , = , they will be denoted by B, D and E respectively. Whether B is the auxiliary matrix of (2.6) or the scalar constant B. of (2. If) will be obvious from the use. Finally, J >k A the elements of the matrix B will be referred to as b. .to avoid any confusion with b. , of (2.6). A 5.1 The Limiting Values Of b. . iii In this section it will be shown that if B. , = B and D. , = D, J,k j,k A except for D, , = and B = , then the elements b. . , of B, approach J-^k J? 1 ijj a limit as both subscripts increase. In particular G . , -* G =\IED as j and k J,k increase. This will be used in the next section to derive an estimate of X . . mm The following lemma will establish some useful properties of A+B. Although we are only concerned with the case a = 1, the lemma is valid for all ae[0,l]. Lemma 5.1; The quantities defined in (2.6), (2.7) exist and satisfy the following relations: 78 79 (a) V SB j,k*° (b) c o,k*V so (0) d J,k <0 (d) -1 s e.^ k t= (e) -1 * f Jjfc £ (h > c j,k = d j-l,k e j-l,k j,k j,k j,k " j,k e j,k-l ~ c j,k j-l,k £ Proof; Parts (a) through (h) were proved by Bracha [1]. To prove (i) we have, j,k j,k j,k " j,k e j,k-l " c j,k J-l,k j,k j,k J " j,k e j,k-l " C j,k j-l,k j,k " j,k + D - c. - b. . f . . n - c. . e J,k j,k j,k j,k-l j,k j-l,k - -V {1 + e a,k.i + f j; k-i) - "j.kfr + e j.i,k + 'j.i, k ) • Using (a), (b) and (f) yields d. , + B. , + D. , - b. . e. . , - c. . f . n . § j,k j,k j,k o,k j,k-l j,k j-l,k which completes the proof. Before demonstrating that G . , = c . . f . , . tends to sT BD as j,k j,k j-l,k j and k increase, it is shown that for fixed k the quantities of (2.7) approach limits as j increases. 80 Lemma 5-2: Let B. . = B and D = D, except B = and D = . r . J,.K J,.K. J, J- -L,K For fixed k the elements b.^ , c.^ , d.^ , e.^ , f .^ , C.^ = b.^ e.^ and G = c f . , , tend to limits as j increases. These limits will j,k j,k 3-1, k be called b fc , c fc , d^. , e k , f fc , C fc and G k respectively. Proof: The proof will be done by an induction on k. For k = 1 we have the following. b . _ = B . . =0 and c . _ = D . . = D . 3A jji oA jA d n = _ 2(B+D)- De. _ _ and e. _ = D/d. . hence d. A= -2(B + D)- D2/d._ 1A (5 - l} )^I which tends to a limit d = -(B+D)+Wb +2BD Consequently e . _ -* e n = D/d. . 3,1 1 ' 1 Finally f = (B-Df. -,)/d. . For large j this becomes 3,1 J ""1,1 3,1 B-Df. . . f _ J- 1 ^ 3,1 \ which approaches B/(d_+D) as j increases. Therefore G. -» G = c f Now assume the truth of the lemma for k < k. For k = K we have the following. b . „ -> B-G^ . = b^ . 3,K K-l K c . Jr - D-G . = c . 3,K K-l K d. jr -» -2(B+D)+2G t ^ . - b f . + c 2/d . 3,K K-l K K-l K ' K This is in the same form as (5.1) with only the constants changed. Thus d, „ tends to 3,£ 81 a K" -( M H-i- T \ B > D. Lemma 5-k then shows that the b, are increasing in magnitude for any B, D < 0. This result is used in Lemma 5-5 to prove that as j and k increase the quantities t> . , and d. are approaching limits h and d respectively. Lemma 5.3: Suppose > B > D . Let b and cl be as in Lemma 5-2. Then the b, are decreasing and the d are increasing as k increases. Proof; The proof will again be by an induction. b. = 0. c = D. d = -(B+D)+\yB^+2BD . f = B/(d + D) = (B +VB 2 +2BD)/2D . and G, = c f = (B +\/ B 2 +2BD)/2. To show that b < b and d > d we have the following relations. 82 b 2 = B - G„ C 2 = (B^\jB 2 +2BD)/2 < = 1d 1 ? o D - G = (2D-B'^\/B 2 +2bD)/2 , and d 2 = -(B+D)+G 1 - -|> b 2 f x + ~ {[-2(B+D)+2G 1 -b 2 f 1 f -l+c^} 2 Substituting the above expressions for f , G , b and C yields d 2 =—--]> -|^\/b 2 +2BD + -|- (- -2- B 2 + 6BD-3B~\/B 2 +2BD)"2" . Since B > D, -3ByB 2 +2BD > -3b\/B 2 +2BD + ^p- + B 2 + 2BD-6BD. (5.2) Collecting terms taking the square root and multiplying by l/2 we obtain -i(- -|- B 2 +6bD-3BVB 2 +2BD) 1 / 2 > - -|- B+|"WB 2 +2BD . Using this inequality and the last expression for d yields d 2 > -B-I>Vb 2 +2BD = d . The following relations, which come from substituting the limiting values of b and d. into (2.7) will be required in the general step of 3 )*■ J )K- the induction. b. , -B+D -1 k-1 \ =B(1 + \~T ] * =2 >3>'" (5-3) and p b (b -B+D)^ Note that b < B < implies k b -B+D and thus - k 4^ — < 0. (5.5) d k. 1 83 Before starting the inductive step, it is necessary to show that b, < b . This inequality holds if 3 2 or equivalently if b -B-D < "4 11 - • (5.6) d 2 d l Replacing b^ and d by (B^\/b +2BD)/2 and -(B+D) +WB +2BD (5.6) becomes -3BD-2D 2 +3D B 2 +2BD 2d^ < -B + D . Substituting for d and collecting terms we get (B+2D) (- -|- ^yB 2 +2BD)- -i BD < (D-B) (- -2- B 2 +6BJ)-jS\J B 2 +2BI)) 1 ' 2 . (5-7) This will be true if B 2 (B 2 +l+B]>lfD 2 ) ( -2- -ByB 2 +2BD + B 2 +2BD) - 7BD(B+2D)*Vb 2 +2BD is greater than (D 2 -2BB+B 2 ) (- -£■ B 2 +6BD-3ByB 2 +2BD) . 4- (That D-B is negative was used in obtaining this condition. ) Multiplying out the factors of the above inequality leads to 2B^-l/2 B 5 D +25-2- B 2 D 2 +2BD 5 -ByB 2 +2bd (-2B 2 +15D 2 +17BD) > which is obviously true. Thus (5-7) is valid and therefore (5. 6) is valid and b^ < b„ . 3 ^ We now start the inductive step. Assume b < b and 8k \-l > \-2 ■ ^ \ = Vl " \ > ° ' aM let e k = ^ "^-1 • From (5.U) b 2 (b. -B+D) 2 k v k Replacing b fc with b fe _ 1 - ^ and d^ with € k + d^ introduces - 2 (l> f "b k _ 1 )' Then using (5.0 and the definition of e fc yields Since d. 2 < d, , by assumption, ^_;L +2 Wl'^ * B 2 +I) g -2Bb k _ 1 +2I>b k _ 1 -2BJ (yyM) 8 e * >28 * \7i \.x ^ 2& e (b -B+D) 2 ■^ K-A-*^' + a^ • (5 - 9) Now from part (i) of Lemma 5*1 we have d, ,+B+D-b, , e p -c, ,f . 1 > 0. But -Vi e k _i = " c k- 2 f k-2 = \-r B > and -c. _f. _ = b -B k-1 k-1 k Hence -1 ■"""" "k-l~~'"k d i: and t^VVi >0 • (5 - 10) Now (5.9) and (5. 10) together with 5, > and d. , > yield e (b -B+D) 2 k > K V 4 1 k 85 Furthermore (5-10) implies | (b k -B+D)/d k _ 1 | <1 and (5.5) implies I Ob. -BKD)/d. I < 1. Hence e > and therefore d> d It remains necessary to show that b, _ < b . We will show that e < e , where e is defined as in Lemma 5.2 (This in turn requires f k-l d. . and b < b . If in addition f < f k-1 k-2 (5-12) 86 Then ^f w < -^-\_ ± \_ 2 and _ -2B-b, ,f, p p (1+e/ < k - X X - 2 - (1+e^) 2 . (5.13) But l+e n and l+e n , are positive so that (5.13) implies k k-1 1+e. < 1+e _ k k-1 or e, < e. n . (5.H) k k-1 To see that f. . < f , n we write k-1 k-2 k k-1 k-1 k-1 k-2 k-2 Since c, -e, , = b, -b, _ , the last equation implies k k-1 k k-1 •WWA-i - c k-i + (° k .i +8 k.i )f k- e • "w'Vl^ " 5 k +8 k-l f k-2 • Vi (f k-r f k-2» +( - M)(f k.r f w» - VV-Ajs • f k-2<-\- 2 +B - D ' > f k .l ( -Vl +M) • But c. =b n -B+D so k k -c , f > f £Si_ > f k-2 k-1 -c, _ k-1 k-2 Hence (5-12) holds and implies (5-llj.) which implies b < b . This completes the proof. If we remove the condition B > D we can still prove b < b 87 for all k although it is no longer true that d. , > d, . The restriction B > D was used only in (5-2) to show d > d and in (5.7) to show b < b . If D ^ B then (5*7) is trivially true and thus b, < b_ . Furthermore if for some k, d, ^ cL, and b < b , then from (5.3) b -B+D \ _1 / b, -B+D \ _1 That is, b < b k whether d^^d k orcL. Thus the previous lemma handled the more difficult case of the following lemma. Lemma 5>k: The b as defined in Lemma 5»2 are decreasing. Lemma 5.U can now be used to prove all the elements of the matrix B approach limits. Lemma 5 . g : Let B = B and D. = D. The quantities b , c. , 3 >& 3 >*■ 3 >&■ 3 >x- d. . , e. . and f . . as defined in (2.1) approach limiting values as J,k j,k j,k ** *° j and k increase. We will call these values b,c,d,e, and f respectively. Proof: Since all elements can be defined in terms of b . , and d. . , J,k j,k we will prove only that b . , -» b and d . , -*• d. J ,& 3 >*■ By Lemma 5.1+, the b are increasing in magnitude as k increases. It is sufficient to show that the b are bounded below for then they must approach a limit, b. Moreover the d^ must approach a limit since 88 b. .-B+D , bk = B(1+ ^_)- Let d . = rain {&. } and let b, = , b = b and mm k j 1 d d b = b ( Ssia ) k d . +b\ ,-B+D mm k-1 Now suppose b ^ b then d . k+1 d . +b ? -B+D mm k d . mm d . +b, -B+D mm k since the denominator is positive and increased. The fraction is greater than one so that adding the same positive number to the numerator and denominator decreases it. Hence \ n = B \ = b, n . k+1 vv» k+1 Thus b\ ^ b n for all k. Furthermore the b n decrease to k k k -d . + B-D+V (d . -B+D) 2 +l^Bd . limh = 212 V min n 21HL . k-* oo The b , therefore, are bounded below and the proof is complete. The next theorem establishes values for the limits of b, and Theorem 5.1: If b and d are as in Lemma 5.2 then \ 89 lim b = B- VbD and lim cL = -B-D+2 V"bD . k -* oo k «* oo Proof: From Lemma 5-5 ^ have that b, and cL approach limits b and d respectively. In addition these limits satisfy (5.3) and (5.k) since the b and cL satisfy them. Hence Bd J^ b-B+ D \ - 1 . , b = d+b-B+D =B 1+ — (5.15) and -1 o 2 d = -2(D+b) - i (b + (b-B+D) ) . Therefore d 2 +2d(D+b)+b 2 + (b-B+D) 2 = . (5.16) From (5.15) b 2 m B 2 (l+ 2(b-B>D) + (b-B+D) 2 f l _ d 2 2 Using (5. 16) to replace (b-B+D) /d yields b 2 = B 2 (i + 2(b-B+D) _ x _ 2(D+b) _ _bf_ ) -1 d d ,2 d 2 2 B d -2Bd-b 2 Consequently k P ? P -b -b (2Bd)-B d = 0. and d 2 2 2 2 ,2 2Bd - V kB d -UB d b = — ■ -Bd . Since b < , the last equality implies b = -V^id . (5- IT) 90 Substituting (5-17) into (5.16) yields d 2 +2d(D-V^Bd)-2Bd -2"V-Bd (-B+D)+(-B+D) 2 = 0. Collecting terms we get (d-B+D) 2 -2V^Bd~(a-B+D) = 0. Hence -2b = 2"V -Bd = d-B+D . (5-18) Squaring both sides of the second equality and again collecting terms = d 2 +B 2 +D +2Bd+2dD-2BD. Therefore d +\ 2 2 2 2 -2B-2D - V 1|3 +4D +8BD+8BD-^B -W> d = 2 = _B-D -2 VBD Note that if d = -B-D-2 VED then (5. 18) implies b = B+V BD . But b ^ b_ < B. k Therefore we must have d = -B-D+2 VBD and b = B-VBD . This completes the proof. Corollary 5. 1.1; If B. . = B and D. , = D then G. , of (2.7) tends to 0,k J,k j,k =V^ G = yBD as j and k increase. Proof: This is obvious from the theorem. 91 5.2 An Estimate Of X. rain In this section an estimate of X . to be used, in the algorithm mm to of Chapter k is developed. We "begin "by examining the case B. , m I and D. . a 1 except B. ■ D, . =0. J,-K o,k o,l l,K: Then some more general cases are considered. -1, It will he shown that for B = D =1, X . ((A+B)~ A) approaches mm l/2 from ahove as the meshsize, h, of the grid used in deriving (2.2) from (2.1), tends to zero. We begin with the following lemma which will be used to show X . ( (A+B) A)s l/2 for all h > 0. mm ' Lemma ^.6: Let A be as in (2.k) with B. , = D =1 except B. = D = 3 >& J^-K J,l 1,£ j,k = \,2,'" , n . Let B = ($. .) be the auxiliary matrix derived from ^- > 3 (2.7) using a = 1. Then for all x £ 0, > . Proof: Bracha has shown in [1] that if A = (a. .) is an N by N matrix in the form of (2.4) and B = (b\ .) is given by (2.7) then 1M-J. = J -a_ i+1 |2 ¥ I x.,n-x. + y -a. . x. -x. 1+1 1 .*• 1,1+n i+n 1 1 1=1 ' N > (a. .+a. . ,-,+a. ., +a. . .+a. , £*. 1,1 1,1+1 i,i + n i-l,i i-n,l" 1=1 ' x. and V- A . ,2 V A = - Z. b. . ,-. x. , -x. + ) b. . , f*n 1,1+1 j i+n i| f*p 1,1+n x. 1 N-l N-l V b. . ,_ x. , ,-x. + Y b. A* , 1,1+1 1+1 i A» ,, 1,1+1 i=n+l ' ' i=n+l ' x. 1 (5.19) 92 N-n-1 k *"> i+n x. , , -x. , -l+l i+n 2 ^Ji- 1 + £ b i+l,i+n ( l X i+lf i=l + x. i+n ■), where N = n The second and fourth sums may be rewritten as N-n-1 I e, i=l i+l,i+n+l X i+1 (5-20) and N-^" 1 a 2 > h. . , x. *-• i+n,i+n+l i=l i+n (5-21) Now from the definition of B i+l,i+n+l A l+l, i+n and A A A "b. . . = - ah. . . = - ah. _ . i+n,i+n+l i+n,i+l i+l,i+n Hence for a = 1, (5.20) and (5.21) become N-n-1 i=l i+1, i+n x i+1 2 N "^- 1 a and - ) D.,.. . A- l+l, i+n i=l x. i+n respectively. Therefore, N-l <(A-B)x,x>= Y -a. . _ X. , - -X. 1+1 1 2 r 1 ft i=n+l ' x _ x i+1 i N-n i=l ^ 1+n x. , -x. i+n i 2 N ~ n / + ft L ' 1+n x. -x. i+n i (5.23! 93 N + 5^ (a. ,+a. ,.,+a. ,, +a. . +a )!x N-n-1 i=l i+l,i+n l+l i+n N-l f^ i,i+l i,i+l' N-n A + T. (-a. . +b. . ) £■*, v 1,1+ n 1,1+rr N-n t ltd N I i=l N-n-1 x -x. l+l l X. , -X. i+n i + > (a. .+a. , ,,+a. ., +a. -. . + a. .) «-* 1,1 i,i+l i,i+n i-l,i i-n,i x. A T b - -i • i+1 X i+l" X i+l (5.2 5 ) From Lemma 5-5 and Theorem 5-1 it is seen that b 1;J < 1 and from Lemma 5-1 it is seen that b. ., . > 0. Hence all the coefficients H-l, i+n in (5-23) are positive since a. . = -1 for i ^ j. Therefore <(A-B)x,x> > proving the lemma. Theorem $.2; If h, A and B are as in Lemma 5-6, X . ((A+B)" A) > l/2 mm and X . ((A+B) _1 A) -* l/2 as h -» 0. mm" ' Proof: From Lemma 5-6 we have <(A+B)x,x> <2 . ^„, . ThUS <(A+B)x,x> > 2~4x^ = 1/2 . Using (3.6) yields A . ((A+B) -1 A) > 1/2 mm ' X 9h To establish the remainder of the theorem consider the vector i \T , ,, , _/l if j+k=n j > J and k > K / E OJ % = (x •••x ) such that x. , =(~ ,, . (.5. 2k. ; k 1,1' n,n y j,k 10 otherwise where J and K are chosen so large that the elements of B have nearly constant value. That is, if row i of B corresponds to the point (j,k) with j > J and k > K then b ± ± _^ = -1 ; \> ± i _ n+1 s 1 I 1 ° ± ^j- - 1 5 b = 2 ; b. . .= -1 ; b. . .= 1 and b. . = -1. ii ' 1,1+1 ' 1,1+ n-1 1,1+ n = ( n-K-1 )><. and = (n-K-3)^+6. _ ( n-K-1 K 1 ft, pe -\ <(A+B)x,x> " 2(n-K-l)4-2 ~ I 2 ' ^'^ 2 " (n-K-l) as h goes to 0, (i.e., as n goes to infinity) the quantity in (5.25) tends to l/2 from the above. Using (3-6) again completes the proof. The argument used to prove the second half of the above theorem can be used to prove a similar result for a more general matrix A. Suppose A is derived from a differential operator in (2.1) with constant a, (x) and a (x). Then A will have constant diagonals a. . = B. a. . _= D 2 D 1,1-n i,i-l and a. . = -2(B+D) except for some zero elements on rows corresponding to boundary points. Furthermore the elements of the auxiliary matrix B tend to the following values: b. .+ -> -V BD : b. .+ , „ N -A/BD • 1,1-n ' i,i-(n-l) b + n -» -V BD and b. . -> 2VBD . Thus if X is defined as in (5-21+) then 1 ^i-J- 1 ,1 (n-K-1) (-2) (B+D) <(A+B)x,x> (n-K-1) (~2) (B+D)+ (n-K-3)My^ )^^~ -2 (n-K-1) (B+D) 2 (n-K-1 J (-B-D+yBD) -2yBD" 95 -(B+P) 1 -(B+D) +2VbF p + 2VBT (3.6) and the above yield the following theorem. Theorem 5.3: If A is in the form of (2.6) with B. . = B and D. , = D except B . . = and D = and if B is the auxiliary matrix defined in (2.10) using a = 1, then lim^aA+B)- 1 ^ t - \ m , (5.26) where h is the meshsize of the grid used in the discretization of the differential problem. There is no lemma coresponding to Lemma 5*6 which implies the inequality in (5*26) can be reversed. That is, the quantity in (5.26) is not lim X . ((A+B)~ A). However this value has proved to be effective h-> mn as the initial left endpoint, a, , in the algorithm of Chapter k- This value is suggested. The result of the previous theorem can be generalized to a martix A which is block tridiagonal with constant diagonals within each block. The result is of little value however since again there is no lemma which like Lemma 5-6 establises a lower bound for the eigenvalues. Furthermore, the value, which arises from this modification of the argument leading to Theorem 5«3> is so crude an estimate of X . that it is not worth calculating ^' mm for use in the algorithm. As mentioned in Chapter k the initial left endpoint, a, , to use in the algorithm should be calculated according to Theorem 5-3 when possible. If the Theorem is not applicable then a e[l/l,l/2] is sufficient. LIST OF REFERENCES [I] Bracha, A. , "An Implicit Method for the Solution of Elliptic Partial Differential Equations", Department of Computer Science, University of Illinois, Urbana, Illinois, 61801, Report No. 392, April, 1970. [2] Bracha, A. "A Symmetric Procedure for the Solution of Ellipic Boundary Value Problems", Department of Computer Science, University of Illinois, Urbana, Illinois, Report No. kk-O, April, 1971. [3] Carre, B. A. , "The Determination of the Optimum Accelerating Factor for Successive Over-Relavation", The Computer Journal , Vol. k, April, 1961. [4] Dupont, T. , "A Factorization Procedure for the Solution of Elliptic Difference Equations", SIAM Journal On Numerical Analysis , December, 1968. [5] Dupont, T. , R. P. Kendall, and H. H. Rachford, Jr., "An Approximate Factorization Procedure for Solving Self -Adjoint Elliptic Difference Equations", SIAM Journal on Numerical Analysis , September, 1968. [6] D'Yakonov, E. G. , "On the Solution of Some Elliptic Difference Equations", Moscow State University, U.S.S.R. , 1969. [7] Engeli, M. , "Over-Relaxation and Related Methods "(Chapter IV of Refined Iterative Methods for Computation of the Solution and the Eigenvalues of Self -Adjoint Boundary Value Problems, by M. Engeli, Mitteilungen aus dem Ins ti tut fur angewandte Mathematik. an der Eidgenossischen Technischen Hochechule Zurich, Nr. 8). Birkhause Verlag, Basel/Stuttgart, 1959. [8] Frank, ¥. L. , "Solution of Lemar Systems by Richardson' s Method", J. Assoc . Compu . March , Vol. 7, July, i960. [9] Frankel, S. P., "Convergence Rates of Iterative Treatments of Partial Differential Equations", MTAC , Vol. k, 1950. [10] Hageman, L. A. and R. B. Kellogg, "Estimating Optimum Overrdavation Parameters", Math . Comp . , Vol. 22, 1968. [II] Kim, Y. J. , Private Communication, 1971. [12] Rigler, A. K. , "Estimation of the Successive Over-Relaxation Factor", Math . Comp . , Vol. 19, 1965. [13] Saulyev, V. K. , "integration of Equations of Parabolic Type by the Method of Nets', Pergamon Press, I96I4.. 96 97 [Ik] Sheldon, J. W. , "On the Numerical Solutions of Elliptic Difference Equations", MTAC , Vol. 9, 1955- [15] Stone, H. L. , Private Communication, April, 1969. [16] Stone, H. L. , "iterative Solution of Implicit Approximations of Multidimensional Partial Differential Equations", SIAM Journal On Numerical Analysis , September, 1968. [17] Taranto, R. , "Numerical Studies of Stones' Symmetric Factorization and the Iteration Parameters a and t", Department of Computer Science, University of Illinois, Urbana, Illinois, 61801, Report No. U23, January, 1971- [18] Varga, R. S. , "A Comparison of the Successive Overrelaxation Method and Semi-Iterative Methods Using Thebychev Polynomials", T. Soc . Indust . Appl . Math . , Vol. 5, June, 1957. [19] Varga, R. S. , Matrix Iterative Analysis , Prentice -Hall, Englewood Cliffs, New Jersey, 1962. [20] Wackspress, E. L. , Iterative Solution of Elliptic Systems , Prentice- Hall, Englewood Cliffs, New Jersey, I9E0". [21] Wilkenson, J. H. , The Algebraic Eigenvalue Problem , Clarendon Press, Oxford, 1965. [22] Wilkanson, J. H. , "The Calculation of the Latent Roots and Vectors of Matrix on the Pilot Model of the ACE", Proc. Camb . Phil . Soc . 50, 195!+. [23] Wrigley, H. E. "Accelecating the Jacobi Method for Solving Simultanous Equations by Thebyshev Extrapolation when then Eigenvalues of the Iteration Matrix are Complex", The Computer Journal, Vol, 6, July, 1963. 98 VITA Martin Allen Diamond was born in Los Angeles, California in September 19^-5 • After receiving a B.A. in Mathematics from Arizona State University in 1967, he was granted an NDEA fellowship and continued studying at that university until August 1968 when he received an M.A. in Mathematics. Since September 1968 he has been a teaching assistant and research assistant in the Department of Computer Science of the University of Illinois at Urbana-Champaign. ^ \ST*