nil 
 
L I B RAHY 
 
 OF THE 
 
 U N IVERSITY 
 
 OF ILLI NOI5 
 
 S\0-8A- 
 
 r\o. 2A3-2A9 
 
 cot 
 
 >. fc 
 
The person charging this material is re- 
 sponsible for its return 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 
 
 ►W 13 1372 
 
 g£B 1 6 Kecd 
 
Digitized by the Internet Archive 
 in 2013 
 
 http://archive.org/details/variablesuccessi244mcdo 
 

 •7 
 
 - \ t u 
 
 
 ' ' E ER-RELAXA ! t 
 
 b.\ 
 
 md K . Mc 
 
 THE UBRARY OF THE 
 
 u£C * 4 1967 
 UNIVERSITY OF ILLINOIS 
 
 i 
 
Report No<, 24 4 
 
 •AETA&LE SUCCESSIVE OvER-FELAXAT 1 ' DN 
 
 by 
 
 Leland K„ McDowell 
 
 September 18, 196c 
 
 Department of Computer Science 
 iversity of .Illinois 
 Jrbana, Cllinois 618OI 
 
 j in partial fulfillment of the Doctor of 
 1 Mathematics and was supported 1 n part by 
 
 al Science Foundation under Grant No. NSF-GP-H6360 
 
Ill 
 ACKNOWLEDGMENT 
 
 The author acknowledges with pleasure the assistance and 
 encouragement of Professor Donald B, Gillies, in whose researches the 
 problem treated here originated and whose many helpful comments and 
 suggestions contributed greatly to the present work. 
 
 The author would also like to express his gratitude to the 
 Department of Computer Science of the University of Illinois for 
 providing him with financial support during the preparation of this 
 work and for making available to him the facilities of the ILLLAC II 
 computer system, This work was supported in part by the National 
 Science Foundation under Grant No, jP-U636„ 
 
 Many thanks are also due to Miss Bonnie J, Malcor for patient 
 and careful typing of the manuscript. 
 
IV 
 
 TABLE OF CONTENTS 
 
 Page 
 
 1. INTRODUCTION . . . , . . . . . . „ . . . . . „ . Q . . . . . . . 1 
 
 1.1 Relaxation Methods, . . . . . „ „ , „ . . . . „ „ . . . 1 
 
 1.2 Variable Extrapolation Factors. ..,...,.„.„„,. 3 
 
 2. VARIABLE SUCCESSIVE OVER -RELAXATION FOR LINEAR SYSTEMS WHOSE 
 COEFFICIENT MATRICES ARE CONSISTENTLY ORDERED S-MATRICES . . . . 6 
 
 2.1 Basic Concepts. ........ . „ ... ... . . . . . „ 6 
 
 2.2 Successive Over-Relaxation (SOR). ............. 10 
 
 2.3 VSOR With Two-Extrapolation Factors ............ 12 
 
 2.U General VSOR; Extrapolation Matrices ........... 27 
 
 2,5 Extrapolation Matrices When the Block Matrices Have a 
 
 "jnon Basis of Eigenvectors, ..,..„„.„..„... 31 
 
 3. SEQUENTIAL EXTRAPOLATED IMPLICIT METHODS ..«'".. . . 36 
 
 3.1 A Model Problem ....................... 36 
 
 3.2 The Sequential Extrapolated Implicit Method (SEl) ..... 37 
 
 3.3 Cyclic Chebyshev SEl. ................... k2 
 
 3.U Numerical Results . . . . .......... ...... . ^5 
 
 3.5 The Variable Sequential Extrapolated Implicit Method 
 
 (VSEl) 49 
 
 3.6 VSEI with o^ =* p fi =» ................. . 
 
 DN ........................ o . . 55 
 
 1 Summary .......................... 55 
 
 2 Extensions. ........................ 56 
 
 ST OF REFERENCES ....................... 57 
 
 "^A s8 
 
 7 •*■ J-* 1 000090000009000000000000000000 V*- / 
 
VARIABLE SUCCESSIVE OVER -RELAXATION 
 
 by 
 
 Leland Kitchin McDowell 
 
 Department of Mathematics 
 
 University of Illinois, 1961 
 
 This thesis investigates the solution of linear systems by 
 extrapolated Gauss-Seidel iteration using a multiplicity of extrapolation 
 parameters, Phis is a generalization of the method of successive over- 
 relaxation., which uses a single scalar extrapolation parameter. The linear 
 systems considered are those which arise in the numerical solution of 
 boundary value problems for self-adjoint,, elliptic partial differential 
 equations , 
 
 In Chapter 2 it is shown that the use of two appropriately 
 chosen scalar extrapolation factors yields an iteration having a higher 
 rate of convergence than SOR, and formulas are derived for choosing the 
 factors optimally,, Also, it is shown that by the use of extrapolation 
 matrices, an iteration can be constructed whose matrix is nilpotent, i„e„ 5 
 whose rate of convergence is infinite. 
 
 Chapter 3 considers a more limited class of linear systems for 
 which a certain sort of spectral decomposition is possible. For the 
 solution of such systems the SE" and VSEJ- methods are introduced and shown 
 to be equivalent to several simultaneous extrapolated lauss-Seidel iterations 
 on certain sabspaces^ each with a different extrapolation factor or set of 
 factors o Theoretical and experimental results are presented which show that 
 SEIj which requires less work per iteration than SOR ? has the same 
 asymptotic rate of convergence and. for most starting vectors, an improved 
 actual rate of convergence. Experimental evidence is presented showing 
 that certain versions of VSE7 have a higher asymptotic rate of convergence 
 for the problems considered than 30R„ 
 
1. INTRODUCTION 
 
 1 „ 1 Relaxat i on Metho ds 
 
 The use of relaxation for the solution of linear systems is 
 documented as early as 1823;, when Gauss [3] wrote favorably of his 
 experience with such a procedure,, which he called indirec t elimination 
 The l a cob i Iteration for the solution of a linear system 
 
 AZ -- K (l.l.l) 
 
 where A is symmetric was discussed in 1846 by Jacobi [7], who used plane 
 rotations to increase the diagonal dominance of A and hence to improve 
 convergence, 
 
 In 1862, Seidel [10] considered the method now known as 
 Gauss-Seidel iteration for the case in which A is m x n with m > n„ 
 Seidel proved that his iteration applied to 
 
 A T A Z = A T K L.1.2) 
 
 converges to a point which best satisfies (1,1,1) in the sense of least 
 squares and hence to the unique solution if A is square and non-singular, 
 
 Seidel noted that if there exists a subsystem of k unknowns 
 and k equations in which the unknowns within the subsystem are coupled to 
 only a few outside unknowns, then it is profitable to relax repeatedly 
 the residuals of the subsystem until they become small while treating the 
 outside unknowns as constants, Thus he anticipated modern block relaxation 
 methods^ in which a number of residuals are relaxed simultaneously, 
 
 * Numbers in square brackets refer to items in the List of Referent 
 
2 
 
 Both the Jacobi and the Gauss-Seidel iterations are known to 
 be convergent if, for example,, the coefficient matrix is either strictly 
 or jrreducibly di agonally dominant [ 5 ] .- 
 
 Experience with hand computation showed that it is often 
 profitable to over-relax, and intuition indicates that the greatest 
 advantage lies in relaxing the largest residual at each stage. If this 
 is done, then the judgment of the relaxer intervenes to alter the 
 algorithm from one iteration to the next, making the iteration non-cyclic . 
 
 The use of electronic digital computers permits the application 
 of relaxation methods to linear systems of very large order, but such 
 machines are best suited to cyclic iterative methods , in which the 
 algorithm once defined is repeated without alteration at each iteration,, 
 For linear systems arising in the discretization of boundary value 
 problems for self-adjoint, elliptic partial differential equations, 
 Frankel [4] and Young [12] in 1950 described such a cyclic procedure 
 for e xtra p olated Gauss-Seidel iteration , often called successive over- 
 relaxation (SOR) o Young showed that for certain orderings of the 
 unknowns, which he termed consistent orderings, a functional relationship 
 exists between the extrapolation factor w and the rate of convergence of 
 SOR, Using this relationship. Young derived formulas for the optimum 
 value of w and for the rate of convergence of the resulting iteration, 
 which he showed to be substantially faster than Gauss-Seidel iteration, 
 especially when the latter is slow 
 
 Ybung ; s results, which applied to point iterations, were 
 extended to block iterations in 195& by Arms, Gates, and Zondek [l]. 
 
3 
 
 Few results are known concerning the use of extrapolation in 
 case no consistent ordering of the unknowns exists, but Ostrowski [9] 
 proved in 195^ that if A is Hermitian„ then SOR converges if and only if 
 A is positive definite and < w < 2„ 
 
 In 1959 Golub [6] introduced the Chebyshev semi-iterative 
 method and the modified (or cyclic ) Chebyshev semi-iterative method . 
 The latter is a variation of SOR asing a particular consistent ordering 
 of the unknowns and a sequence of extrapolation factors,, a different 
 pair of factors being used for each iteration. The cyclic Chebyshev 
 semi-iterative method provides an improved average rate of convergence 
 over SOR, although the asymptotic rates of convergence of the two 
 iterations are the same, 
 
 1,2 Variable Extrapolation Factors 
 
 Whereas SOR uses the same scalar extrapolation factor for each 
 unknown of the linear system being solved, this thesis investigates 
 extrapolated Gauss -Seidel iteration using extrapolation factors which 
 vary from one unknown to another or from one block of unknowns to 
 another. The use of matrices as extrapolation parameters is also 
 investigated. We call such iterations variable successive over -relaxation 
 (VSOR) 
 
 The linear systems considered are those for which the use of a 
 single extrapolation factor is known to be profitable ; e g„_, linear 
 systems arising from the discretization of boundary value problems for 
 self-adjoint elliptic partial differential equations 
 
k 
 
 Chapter 2 begins by introducing 2-factor VSOR, an extrapolated 
 Gauss-Seidel iteration using two scalar extrapolation factors. Formulas 
 for the optimum pair of factors are obtained, and it is shown that the 
 rate of convergence of 2-factor VSOR is greater than that of SOR, 
 
 In Section 2 of Chapter 2 the concept of extrapolation matrices 
 is introduced, and it is shown that by their use an extrapolated 
 Gauss-Seidel iteration can be constructed whose iteration matrix is 
 nilpotent, Hence, in the absence of rounding errors, the exact solution 
 is obtained after a finite number of iterations. Moreover, it is shown 
 tnat if the non-null blocks of the associated block Jacobi matrix are all 
 square of order r and have a common basis of eigenvectors, then VSOR with 
 extrapolation matrices is equivalent to r simultaneous extrapolated Gauss- 
 Seidel iterations each using a sequence of scalar extrapolation factors , 
 Thus, a decomposition theorem for VSOR is obtained. 
 
 In Chapter 3 we consider the linear system resulting from the 
 discretization of the Dirichlet problem for Poisson's equation on a 
 rectangle, The resulting rq x rq linear system is the model problem for 
 whose solution the sequential extrapolated implicit method (SET. ) is 
 introduced. This method accelerates Gauss-Seidel iteration by means of a 
 shift of origin, which requires less work per iteration than extrapolation 
 by a scalar, It is shown, however, that SEI is actually equivalent to the 
 use of certain extrapolation matrices, and the decomposition theorem of 
 Chapter 2 is then applied to show that SEI is equivalent to performing 
 SOR simultaneously in each of r subspaces, a different scalar extrapolation 
 factor being used in each subspace. It is then shown that the asymptotic 
 rates of convergence of SEI and SOR are identical when the optimum 
 acceleration parameter is used for each, but that the actual rate of 
 
5 
 convergence of SEI is always at least as high as that of SOR and for 
 certain starting vectors is substantially higher, 
 
 In Section 3»3 cyclic Chebyshev SEI is introduced. The 
 relationship of this iteration to SEI is analogous to the relationship of the 
 cyclic Chebyshev semi-iterative method to SOP. : 
 
 Section 3.^ presents numerical results comparing the actual 
 rates of convergence using various starting vectors of SOR,, SEI , the 
 Chebyshev semi-iterative method., and cyclic Chebyshev SEI , 
 
 Finally, the variable sequential extrapolated implicit method 
 ( VSEI ) for the solution of the model problem is introduced. Several 
 versions of this iteration are investigated, each of which uses the criteria 
 of Chapter 2 for the construction of nilpotent iteration matrices 
 together with the decomposition theorem to make the VSEI iteration matrix 
 nilpotent on certain q-dimensional subspaces of rq-space Experimental 
 evidence is presented which shows that VSEI is asymptotically faster than 
 SOR for the model problem- 
 
VARIABLE SUCCESSIVE OVER -RELAXATION FOR LINEAR SYSTEMS 
 WHOSE COEFFICIENT MATRICES ARE CONSISTENTLY ORDERED S-MATRICES 
 
 2,1 Basic Concepts 
 
 Throughout this thesis we will be concerned with the iterative 
 solution of the linear system 
 
 AZ ^ K 
 
 (2.1.1) 
 
 where A is a matrix usually of large order and sparse; i.e., having only 
 few non-zero entries. A detailed development of the ideas outlined in 
 this section may be found in [11, Chap. 1,3] or [13, Chap. 1] . 
 
 Our results will apply in particular in case A is an S-matrix; 
 i.e., A satisfies 
 
 (i) A is real and symmetric 
 'ill a . . < for i ^ j 
 
 id - 
 
 (iii) A is irreducible; i.e., there exists no permutation 
 matrix P such that 
 
 T 
 PAP = 
 
 A 
 
 11 
 
 
 
 A 21 A 22 
 
 where the diagonal submatrices are square 
 (iv) A is diagonally dominant; i.e., 
 
 li 
 
 > £ 
 d# 
 
 ij 
 
 (v) The inequality in (iv) is strict for at least one i 
 S-matrices arise, for example, in the discretization of self-adjoint, 
 elliptic partial differential equations. 
 
7 
 
 A stationary, linear iteration for the solution of (2„l.l) Ls 
 an aff ine transformation 
 
 ? m4l .HZ m +F (2.1.2) 
 
 with the property that 
 
 Z - HZ + F ['2.1.3) 
 
 The error in the m-th iterate is E " -- Z - Z , and subtraction of (2.1,3) 
 from (2.1.2) shows that 
 
 m+1 TT Z*' m rT nn-JJf o 
 
 E = HE =: H 
 
 1 fT° 1.10 
 
 Hence ; the convergence properties of the iteration (2,1.2) depend 
 
 only upon the matrix H„ In particular, lim E = for arbitrary E 
 
 ' m-»<» 
 
 m 
 if and only if lim H is the null matrix. This is the case if and only 
 m->°° 
 
 if all eigenvalues of H have modulus less than one. The spe ctral radius 
 
 of H, denoted by p(H), is the maximum of the moduli of the eigenvalues of 
 
 H. and the ( asympto tic) rate of convergence of the iteration (2,1,2) is 
 
 R = -In p. 
 
 A partitioning of A in which the i-th diagonal submatrix A 
 
 is r, X r. induces a partitioning Z ' = [Z . \ Z _ J — |Z of the unknown 
 11 ^ to 1 1 2 ' ' q' 
 
 vector Z such that the block Z. contains r, of the unknowns, Tne constant 
 
 1 1 
 
 vector K is similarly partitioned. 
 
 If A is an S-matnx, then A is positive definite, Eact: A, . is 
 
 * 11 
 
 therefore also positive definite and hence non-singular. The block 
 Jacobi iteration corresponding to the particular partitioning of A is 
 
 Z nifl = - Z A" 1 A . Z m + A" 1 K. ; 1< i < q, (2.1.U) 
 
 J+1 ii lj j il i " " 
 
where the partitioning is assumed to be such that A. . is r. x r.. 
 
 ij 1 J 
 
 Throughout this thesis M will denote the matrix for the iteration (2*1.^); 
 
 i e„ 3 M denotes the block Jacobi matrix derived from A, 
 
 The block Gauss - Seidel iteration differs from the block Jacobi 
 
 in that as soon as a block Z . of the next vector iterate is calculated, 
 
 1 ' 
 
 it and not Z . is used in the computation of later blocks Z . _ , 
 l l+l 
 
 Z . ' , etc. Thus, 
 
 14-2 
 
 Z m+1 = - Z. AT?" A.. Z m+1 - .X. A" 1 A . Z m 
 
 l j<i li ij j j >i li ij j 
 
 (2.1.5) 
 
 + A _1 K. . 1 < i < q. 
 
 If P is a permutation matrix, then the coordinate 
 transformation 
 
 P A P T Z = P K (2.1.6) 
 
 
 
 corresponds to a reordering of the unknowns. The ordering o is 
 
 T 
 consistent if A = P A P is block tri-diagonal (i.e., A. . is null if 
 ooo ' ij 
 
 j|>l) with square diagonal blocks, and if so, then A is said to be 
 consistently ordered. The block Jacobi matrix M derived from A is said 
 to be consistently ordered if A is. Throughout this thesis A and M 
 will be assumed consistently ordered already and the subscript o will be 
 omitted, 
 
 If M is consistently ordered, it is of the form 
 
q-lq-1 q-lq 
 
 a" 1 A _ 
 qq q^q-i 
 
 o B 
 
 q-lq 
 
 so B. . = for ji-j| =f= 1° Kence M = L + U ; where L is strictly lower 
 triangular and U is strictly upper triangular „ Moveover,, if M is 
 consistently ordered and |i is an eigenvalue, then -u is an eigenvalue 
 of the same multiplicity. 
 
10 
 
 If 
 
 D = 
 
 -1 -1/2 -1/2 
 
 then M = I - D A. Hence M is similar to I - D ' AD ' , and it 
 
 follows that this is symmetric because the square root of a positive 
 
 definite symmetric matrix is symmetric . Hence all of the eigenvalues 
 
 of M are real, 
 
 2 , 2 Su ccessive Over-Relaxation (SOR) 
 
 Define the intermediate value 
 
 _ m+- _ . 
 
 Z . 2 = B. . _ Z m+ ] + B. . 
 
 1 i,i-l l-l 1,1+1 l+l 11 1 
 
 z m . + a: 1 k 
 
 (2.2.1) 
 
 and extrapolate by a fixed scalar W; 
 
 1 
 
 Z m+1 = "[Z ^ - Z m ] + Z ?, 1 < i < q 
 
 i j_ i i 
 
 (2.2.2) 
 
 Then 
 
 Z m+1 - u>[B. - , Z ^ + B. . . Z . 
 1 i,i-l l-l i,i+l l+l 
 
 H WAT. K., 1 < i < q< 
 11 1 — - 
 
 m ] + (i-«) z m 
 
 (2.2.3) 
 
 The iteration (2.2.3) is the successive over-relaxation method with 
 extrapolation factor <£. The matrix for the iteration (2.2.3) is 
 
11 
 
 £ =» [I - coL]" 1 [uU + (l - co)l] 
 
 CO 
 
 Let < Li < i-t < < Li = (J. be the positive eigenvalues of 
 
 M, and let p be the multiplicity of zero as an eigenvalue. For 1 < i < l y 
 let 
 
 CP.(\) = \ - [2(1-W) + CO Li ] X + (l-w) . 
 Then tne characteristic polynomial of £ is 
 
 CO 
 
 $;\) ^ [(i - co) - \]* n cp.(\) . 
 
 isl X 
 
 The maximum of the moduli of the roots of cp (A.) is minimized as 
 
 JO 
 
 a function of lo when the roots are real and equal; that is,, when the 
 discriminant of 9/A) is zero„ This is the case when 
 
 co = co s 
 b 
 
 + \|l - |I 
 
 For co - co the modulus of the roots is 
 b 7 
 
 p SOR 
 
 1 - 1-L - Li 
 
 i + \fT7^ 
 
 If the discriminant of cp (\) is negative or zero for a 
 particular value of co,, then the discriminant of cp.(A,)., 1 < i < l is 
 also negative or zero. It follows that for co ^ co the roots of <E>(A.) a ll 
 lie on a circle of radius p ariD p and hence 
 
 min pUJ ^ p(j^ ) =p S0R . 
 co b 
 
 This completes the summary of known results. 
 
12 
 2.3 VSOR With Two Extrapolation Factors 
 
 In this section we consider a linear, stationary iteration using 
 two distinct extrapolation factors . It will be shown that for certain 
 consistently ordered S-matrices an optimally chosen pair of factors used 
 in alternation on consecutive blocks yield an iteration having a higher 
 rate of convergence than SOR. 
 
 For any two non-zero scalars w and w consider the iteration 
 
 T 1 " 1 = u> [B ,Z^ + B » m 
 1 l 1,1-1 i-l l 
 
 (2.3.1) 
 
 Z. 1+X = u.[B. . , Z n + B. . n Z m I + (1 - U> ) Z m 
 
 -.,1+1 r+l J s i' l 
 
 + U). A? 1 K., 1 < i < q, 
 l li i' - - ^ 
 
 where 
 
 w is odd 
 
 CO. 
 
 1 w if i is even 
 
 Let Ju denote the iteration matrix for (2.3. l), and suppose that X 
 
 1' 2 _ 
 partitioned like Z is an eigenvector of £ w belonging to the eigenvalue 
 
 1' 2 
 It follows from (2.3. l) that 
 
 AX. = w.[\B. . n X. n + B. . _ X. J + (1 - w. ) X. , 1 < i < q, 
 l l i,i-l i-l i,i+l l+l l l' - — ' 
 
 or, after dividing by w. , 
 
 (l - w . ) - A 
 
 - [X.B. . X + B. X. 1 + [ —1 ] X. . 
 
 1,1-1 i-l 1,1+1 l+l w. i 
 
 ' ' l 
 
 — » \^ — " — ■> 
 
 This linear transformation on X we write as G X = 0. If 
 
 7, = 
 
 (1 - u.) - \ 
 
 i w. 
 
 l 
 
 then 
 
13 
 
 G 
 
 q-i ; q 
 
 \B . y I 
 
 q ; q-l 7 q q 
 
 Here I is the identity matrix of the same order as A..„ 
 1 11 
 
 \s — * — * 
 
 The equation G X = has a non-trivial solution if and only if 
 
 det (Gj = 0, and since det (G.) is a polynomial of degree n in A. ; its 
 
 roots are the eigenvalues of £ „ Assuming A. f (an assumption which 
 
 V 2 
 will later be justified;, we premultipy G by 
 
 V 
 
 1 
 . "2 
 
 and postmultiply by 
 
 A -- 
 r 
 
 \ * 
 
 X 
 
lit 
 
 This transformation multiplies each subdiagonal block by A. , each 
 
 1 
 
 "2 
 diagonal block by \ } and each superdiagonal block by 1. Hence, 
 
 1 
 
 ~2 n/ 
 letting 7. = A. 7 . , we have 
 
 G = A. G A 
 
 q-i,q 
 
 B 7 I 
 q,q-l q q 
 
 But det (A.) det (A ) = X. , and so the following lemma has been proved, 
 
 Lemma 2,3^1 : X. ./• is an eigenvalue of Jo if and only if it is a 
 
 1- 2 
 zero of det (G). 
 
 Since A. . is of order r. the vector Z. contains r 
 11 111 
 
 components, and it follows that w is used with r = Z r. components 
 
 i odd 
 
 and w is used with r = Z r. components. 
 2 e . 1 
 
 1 even 
 
 Lemma 2.3 >2 ; Let the matrix G' be obtained by replacing each diagonal 
 
 \ 1/2 
 element 7. of G by (7 7 J . Then 
 
 r - r 
 o e 
 
 r - r 
 e o 
 
 det(G) = 7 
 
 1 
 
 det(G' ) 
 
15 
 
 Proof 
 
 Let 
 
 1 
 
 y l 1 
 
 r = 
 
 r 
 
 1 
 r 2 2 
 
 7 2 I 
 
 and 
 
 r 
 
 i 
 
 r" 2 I 
 7 2 1 
 
 1 
 
 7 _2 I 
 r l 2 
 
 1 
 7 _2 I 
 
 \ 
 
 Then G = T G" r - and so det (g) = det (r.) det (G' ) det (r ). But 
 
 det(r^) det(r r ) - 7 1 
 
 r - r r - r 
 
 o e e o 
 
 2 2 
 
 7 , -which proves the lemma 
 
16 
 Lemma 2,3^3 ; Let U , u. , ..,, u be the positive eigenvalues of M, and 
 let p denote the multiplicity of zero as an eigenvalue. For 1 < i < £., 
 define 
 
 <P (X) = X - [(1 - w i ) + (1 - o> 2 ) + ^uy^] \ + (i _ Wi )(i _ w 2 ) 
 
 Then the zeroes of det (G* ) are precisely the zeroes of 
 
 £ 
 
 E E 
 
 [(1 - ^) - ^] 2 [(1 - " 2 ) - ^f 
 
 i=l i 
 
 Proof : 
 
 G' = 
 
 l/2 ]/2 
 is, except for the diagonal blocks y ' y ' l. } identical to M. Hence 
 
 1/2 1/2 
 
 det(G ! ) = if and only if 7 7 = \i., where \x. is an eigenvalue of M, 
 
 l/? l/? 
 Moreover, the multiplicity of (7 7 - u.) as a factor of det (G' ) is 
 
 equal to the multiplicity of u. as an eigenvalue of M, If |J.. / 0, then 
 
 / 1/2 1/2 N , N 
 
 (7, 7„ + M.J is also a factor of det (G' ) and therefore so is 
 
17 
 
 2 2 2 
 
 17! y k But ^ 7 2 - --— 
 
 from wh •■■ , s that corresponding to the 21 non- zero eigenvali 
 
 t 
 of M are 2# zeroes of di I ,^mei> the 21 roots of II 
 
 The remaining zeroes of det are the ze: 
 
 1 1 f, 
 
 genvalues 
 
 • 2 2 p 
 r, / - 1 
 
 l-w ] I - m[ l-a> 2 ) - M 
 
 uj a) A. 
 
 from whicn the conclusion of the lemma follows „ 
 
 p + (r - r p + (r - r ) 
 
 o e e o 
 
 Lemma 2 - 3 and p a — are both non- 
 negative integers. 
 
 Proof' r - r = n = p + 2£„ Hence, r + r ) has the same parity as 
 o e ' o e ' 
 
 p, and so, therefore, does r - r Then p + (r - r and 
 
 o e e 
 
 - r i are both even, and so p^ and p^ are integers, 
 ■ o e' 12 
 
 By lemmas 2 3*2 and 2 „ 3 - 3 - the zeroes of det 5) are precisely 
 
 P x P 2 I 
 
 the zeroes of [{l - w ) - Vj [(i - w K] n cp ' \ • , rh< are 
 
 exactly n such zeroes because each one is an eigenvalue of £ by 
 
 2 
 lemma 2 Hence p < p and p < p„ But p + p p. so both p and p 
 
 are non-negaiiv- 
 
 orera_2 The cnaracteristic polynomial of £, is 
 
 2 
 
 P] ^ 
 
 $..■ ~ [(l - II cp 
 
 i 
 
 f: B;y lemmas 2 2 3 2- and 2 of $>(> ) are 
 
 :isely the eigenvalues of £. , But by lemma 2-3-' 4 <£(. V 1 is a 
 
 2 
 
 nomialj which proves the theorem. 
 
18 
 In practical problems the eigenvalues of M are usually not 
 Known, but an interval can be determined in which the eigenvalues are 
 known to lie. We therefore consider choosing to and to to minimize the 
 maximum of the moduli of the roots of 
 
 cp (\) = \ 2 - [(1 - u^) + (l - u> 2 ) + ^ 2 [i 2 ]\ + (1 - w i )(l - w ) 
 
 where it is known only that 0<u_<u-<u<l 8 
 
 The discriminant of cp (X) is 
 
 u 
 
 D(n) - w 1 w 2 ki 2 [w i w 2 H 2 + 2(1 - w ) + 2(1 
 
 to ) 1 + (to 
 2 ;j i- v 2 
 
 V 
 
 Lemma 
 
 2,3.5: If ^ > l, to > i, D(u_) < 0, and D(u) < 0, then D(u) < 
 
 for \i < |i < \i, 
 
 P roof : This is evident from the graph of D(u). 
 
 D(n) 
 
 J*» H 
 
 Figure 1 
 

 note the maximum of the moduli of the roots of 
 
 Then max ITS minimal if and only if to and to are such 
 H J 1 2 
 
 that D(u/ = D(ff : 
 
 Proof The equation cp (\) = can be written 
 
 - U - (1 - u> )][\ - (1 - u) )] = u 2 \ 
 12 X 
 
 Completing the square oh the left hand side yields 
 
 (i - k) - (i - <*> tp 
 
 CO to 
 
 1?L 
 
 \ 
 
 (to - to V 
 2 J 2 l' 
 a U A. + 
 
 i|CO U) 
 
 12 
 
 Some algebraic manipulation converts this to 
 
 \pl 
 
 i_ (X . 1) + i 
 
 u 
 
 \ co \ CO 
 \ 2 N .1 
 
 LO 
 
 2 1 
 = u \. + £ 
 
 
 =1 2 
 
 CO \ CO 
 
 o \ i 
 
 It follows that the roots of 9 'M = are the abscissas of the points of 
 intersection in the (a. a) plane of the parabola 
 
 2 2, 2.1 
 
 
 to 
 
 , 
 
 CO 
 
 
 \ Cu cu 
 
 \| 1 2 
 
 
 to 
 
 
 ■ 
 
 
 
 
 
 
 
 
 
20 
 
 Figure 2 
 
 (jj 
 
 2 2 
 
 77- = I, the vertices of the parabolas a , (p.) are at the origin. For 
 
 Ui 
 
 1 
 
 U) 
 
 0) 
 
 7^ 1, the parabolas intersect on the a-axis, the vertices lie on the 
 
 1 
 
 negative \-axis, and the abscissas of the vertices increase toward zero 
 
 2 
 as u increases, For fixed — > 1 and ^-.^p sufficiently small, the line 
 
 2 ~~ 
 a intersects the parabola a (u) in two distinct points. As w -, w p 
 
 increases, the slope of the line decreases, and so p— decreases until the 
 
 line becomes tangent to the parabola. At this point the roots of cp— (\) 
 
 are real and equal to \| ( 1 - w )(l - w ). If w > 1 and w > 1 (as will 
 be shown to be the case of interest), then further increase of w u 
 
21 
 
 causes the roots of cpH\) to become complex and to increase in modulus c 
 
 u 
 
 2 
 
 It follows then., tnat for — fixed, the value of ^-.^ which minimizes p— 
 
 w, ' 12 u. 
 
 o .— 
 
 is the value such that the line g is tangent to the parabola a, (n)<> 
 
 co 
 2 
 For — > 1 but sufficiently small, the line a„ does not intersect 
 
 the parabola a (u_) when- it is tangent to the parabola a,(|^)j and so the 
 
 CO 
 
 roots of <?..(*■ ) are complex of modulus Wl - w ,)(l - w )o As — — increases 
 
 2, 
 
 the parabolas a-, [n) and a, (u) move to the left and the distance between 
 
 their vertices increases, 
 
 CO CO CO 
 
 2 2 2\ 
 
 For -— sufficiently large, say — = (— v, , there exists a value 
 
 CO •'W CO h 
 
 1 11 
 
 of co w -which makes the line n simultaneously tangent to the parabolas 
 
 2 - 2 / ■— 
 
 7 and a, >n). Hence, D(n) = D(^) --■ 0, and the roots of <p—(?0 a re both 
 I *"~ 1 
 
 equal to U,l - w -.)(l - w „) while the roots of cp \; are botn equal to 
 
 - \| 1 - w )(l - w ) „ Moreover., it follows from lemma 2,3-5 that the roots 
 
 of cp (X) for \x < [i < |i are complex of modulus \|'l - w -.)(l - ^ p 
 
 cu 
 2 
 
 ..creases further and w ,^ is adjusted so that the line 
 "J 12 
 
 p , 
 remains tangent to the parabola a-,'-.'; the point of tangency moves to the 
 
 left because the parabola moves to the left Hence min p.. regarded as a 
 
 function of U) -, w increases, and the line ceases to intersect the parabola 
 
 2- - 
 a, ,- ■) , It follows that max p |p is minimal when to and w are such 
 
 b. _ ~ : 
 
 that the line is simultaneously tangent to both parabolas , which is the 
 
 case if and only if D 0, 
 
22 
 
 o 
 This argument shows that if — > 1, then the pair (oj w ) is 
 
 1 
 
 optimal if and only if D(u) = D(u_) = 0. This latter condition also 
 
 00 
 
 2 
 characterizes an optimal pair such that — < 1, as can be shown by 
 
 OJ 
 
 1 
 
 letting — decrease from 1 in the above argument. Indeed, since cp (X.) is 
 
 tc 
 
 a symmetric function of w and co it follows that (oo oo ) is an optimal 
 
 pa 
 
 ir if and only if (w oo ) is alsc 
 
 2 
 Theorem 2 C 3.3 : Let " - , 
 
 b 1+ u_u + \J(l-u_ 2 )(l-JI 2 ) 
 
 w 
 
 1 - u u + \|(l-u. 2 )(l-Ti 2 ) 
 
 , and p 
 
 W - \ 
 
 1-J 2 
 
 2S0R 
 
 \Jl^? + \|l - ? 
 
 Then 
 
 p u = P 2S0R ^ 
 
 max _ p 
 
 u_<u<^ w ]_> u 2 | li < ^ < I 1 
 
 M- 
 
 00 = OJ j 00 = 00 
 
 1 and only ' either < - > or < , , , ) holds. 
 
 » OJ = 00 / \ 0J = 00 
 
 2 b I 2 -b 
 
 S0R 
 
 , ) If w denotes the optimum extrapolation factor for SOR and 
 p(£ ), thenl < to < w < ST and p os „ < p 
 
 U> 
 
 b - b - b 
 
 "2S0R - M S0R' 
 
 (iii) uo b = oo b = oo b and p^ = p g0R if and only if u_ = 0, 
 
 Proof: 
 
 (i) By Theorem 2,3.2, max _ p is minimal if and only if 
 
 u_ < u < u 
 
 Dl ,, ) - D(u) = 0. Since 
 
 D(|i) - CdjUyi [oo^u 2 + 2(1-^) + 2(l-w 2 )] + (oj^oj^ 
 
 2 2-2 
 
 is quadratic in u- with roots u_ and |i } we have 
 

 (cu I U) - ) 
 
 1 
 
 CO to 
 1 2 
 
 
 
 (w - to V 
 
 £ u = 2 p 
 
 CO CO 
 
 1 2 
 
 (2.3.3) 
 
 CO - CO 
 
 Adding plus and minus the positive square root Li u_ -• — — — — of (2,3 = 3) 
 
 to [2 3-2' yields 
 
 CO CO 
 
 1 2 
 
 (u + u_) co^ - Uco 2 + h 
 
 and 
 
 - u_) co^ 2 - to^ + h = 0, 
 
 whose solutions are 
 
 and 
 
 CO ^ CO 
 
 1 -b 
 
 CO = CO 
 
 2 b 
 
 1 
 
 "\ 
 
 1 t J J - \|(l - U )(l - ji ) 
 
 U) 
 
 fc» f- 
 
 \J(1 - ^)(1 - M 2 )j 
 
 (l-<*> )(l-u> )| is smaller for the solution 
 
2h 
 
 CO = CO 
 
 1 -b 
 
 CO = CO 
 
 2 b 
 
 An analogous development using the negative square root (i (i = 
 
 of (2.3.3) yields the soluti 
 
 on 
 
 co = CO 
 
 1 b 
 
 CO = CO 
 
 2 ~b 
 
 Hence, 
 
 min / max _ p > = ^(l - w ) (l - cq ) } 
 
 CO - CO 
 
 _1 2 
 
 CO CO 
 
 1 2 
 
 ind some algebra shows that ^(l - ^ h )(l - w , ) = p 
 
 2S0R 
 
 (ii) Since < kL < n < 1, ™ e can set u_ = sin a and \x = sin p, 
 
 re ^ a < p < - Then 
 
 CO 3 
 
 -b "' 1 + sin a sin P + cos a cos P " 1 + cos (p - a) 
 
 co = 
 
 b 1 - sin a sin P + cos OL cos P " 1 + cos (p + a) 
 
 and 
 
 co = 
 b 
 
 A 
 
 1 - \x c 1 + cos P 
 
 Sine C a < < £, it follows that cos (p + cc) ;_ cop p < cos (p - a), 
 
 ■ 
 
 and so 1 < w < co < co . 
 
 b - b ~ b 
 
 If u_ = 0, then 
 
 J 2S0R 
 
 i -Ni - H 2 
 
 ■N[ 
 
 ■ ? 
 
 p S0R 
 
 But C — p < for u. > 0, and therefore * < p for u > 0. 
 
 riu_ M iS0R U ^2S0R - p S0R L - 
 
25 
 (iii) This follows at once from the formulas which define the 
 
 quantities involved. 
 
 It is of interest to trace the movement in the complex plane of 
 
 the roots of <P (\) for fixed w end w satisfying 1 < w < u as n 
 
 decreases from + « to zero. It follows from (2.3.2) and (2.3.3) in the 
 
 proof of the preceding theorem that each such pair (w u) ) corresponds to 
 
 some pair (u_> u) satisfying < u^ < u. 
 
 (l-" 2 ) 
 
 \J(l-u) i )(l-w 2 ) 
 
 Figure 3 
 
 For u = + c» 5 the roots of cp (k) are zero and + <»„ As j_ decreases., the 
 roots approach each other until for u = u,, they are equal to 
 Yl - u n )(l - ^p^ ^ E ^ decreases from u to u_, the roots move around the 
 circle of radius m«, 1 - w , )(l - °0 in the complex plane, becoming equal 
 to - y'l - w)(l - w ) when u = ji. As u decreased from u_ to zero, the 
 roots move apart on the real axis so that for u ■--. ; one root is (l - u ) 
 and the other is (l - w )„ 
 
26 
 Theorem 2,3°^ Using the notation of lemmas 2„3ol and 2„3.>3> let 
 
 tL = 
 
 u_ if p = r - r 
 1 ' o e 
 
 if p > r - r 
 
 ' o e 
 
 > 
 
 Then 
 
 0^ ,u> } =P 2 S0R = mln pU co , W } • 
 
 Proof : If p = | r - r |, then by theorem 2<,3°1 the characteristic 
 
 polynomial of £ w is 
 1' 2 
 
 P ^ 
 $(\) = [(1 - w ) - M n cp n (M o 
 
 ^ i=l •*! 
 
 By theorem 2„3o3? the maximum of the moduli of the roots of II cp (\) 
 
 i=l **I 
 
 is minimal for either of the two choices for (u) u ) specified in the 
 
 theorem, and in fact the roots lie on a circle of radius P oono ° But the 
 
 2oUn 
 
 remaining roots of 0(\) lie inside this circle because 
 
 |1 " ^ul < Y-l " ->)(-'- " ^v) a P n nR j an ^ the conclusion of the theorem 
 
 follows „ 
 
 If p > | r - r I , then both (l - w ) and (l - w ) are roots of 
 1 o e 1 ' -b b 
 
 $(X), But these are the roots of cp (.\) if [i = 0, and so it follows that 
 
 o'£, w ) is minimized by taking ^ = 0„ 
 1-" 2 
 
27 
 2 , h Gen e ral VSOR : Extr apo lation Matrices 
 
 Let A be partitioned so that each diagonal submatrix A . is 
 
 11 
 
 square of order r. , All of the iterations considered in this thesis for 
 the solution of (2.1.l) are of the form 
 
 r m+1 
 
 = ft.[B. . Z 
 
 ~ m+l 
 
 + B 
 
 m 
 
 i i,i-l i-1 i.i+1 i+-l 
 
 ] 4 (I - fi.)Z 
 
 rn 
 
 + fl A K , 1 < i < q 
 l li l 
 
 (2.U.1) 
 
 where each ft. is a non-singular square matrix of order r. For example, 
 
 in the case of SOR or of the 2 factor VSOR of Section 2.3, each ft. is 
 
 simply a scalar multiple of the identity matrix. The iterations to be 
 
 investigated in Chapter 3 a re also special cases of (2„U,l), although 
 
 their implementation does not involve explicit multiplication by the ft.'s 
 
 i 
 
 If 
 
 o" 
 
 o 
 
 then the iteration matrix for (2.U„l) is 
 
 r, » [i - nL]" [nu + (i - n) 
 
 where L + U - M, the block. Jacobi matrix derived from A, If each ft. is 
 
 ' l 
 
 non-singular and if Z - Z. the true solution of (2.1.l) 5 then Z 
 also, 
 
 rm+i a - 
 
28 
 In case each fi. is a full matrix,, the implementation of (2„4 l) 
 evidently entails considerably more work than extrapolation by a scalar „ 
 Suppose, for example, that r. ~r, l<i<q Z. can be computed most 
 economically (even if fi. is a scalar multiple of the identity) by first 
 computing 
 
 1 
 
 I™'' B. . ^4 B. 2* +A ;!if. (2.U.2) 
 
 l i.i-l l-l i.i-ii-l hi 
 
 and then extrapolating: 
 
 1 
 
 2* 1 - a. [t* 2 - z m j 4 z m . (2.U. 3 ) 
 
 11111 ' 
 
 2 ""* 
 
 ; 2,4,3) involves r multiplications and r(r+l) additions, and since Z. 
 
 contains r unknowns, it follows that the extrapolation procedure for 
 
 2„^ol) requires r multiplications and r+1 additions per unknown of the 
 
 linear system (2.1„l)o Extrapolation by a scalar, however,, requires only 
 
 one multiplication and two additions per unknown. The direct use, 
 
 tnerefore„ of full extrapolation matrices is computationally advantageous 
 
 only if a large increase in the rate of convergence can be obtained over 
 
 iterations using scalar extrapolation factors. 
 
 Such a large increase is indeed possible. In fact,, for a 
 
 properly chosen set [fi.] of extrapolation matrices, £, is nilpotent (i„e 0! 
 
 p'Ju ' ='• Oj so that the rate of convergence of (2.4,l) is infinite, The 
 
 following theorem gives criteria for choosing such a set (fi„ }. 
 
29 
 
 Theorem 2 ; U, 1 , jV is nilpotent if ft , ft , ft satisfy one of 
 
 (i) fi. - I and ft. =. [I - B. . 0. B . J" 1 . 2 < i < q 
 1 i i,i-l i-I i-l,i J - — ^ 
 
 (ii i fi = I and ft. =. [I - B. . _fi. _B. . . ] , 1 < i < q-1 
 q i i,i+l l+l i+l,i } - - ^ 
 
 (iii) ft = I, ft =» I, and for some m such that 1 < m < q; 
 
 ft. = [I - B. . .ft. n B. . . ]" ls 2 < i < m, 
 l i, i-I i-I i-I, i - ' 
 
 ft. = [I - B. . fl, B. . .J" 1 , m+1 < i < q- 
 1 l,14l 1+1 i+l,i ' - - ^ 
 
 ft -• [I - B ,ft _B 
 
 m 
 
 B ft B 
 
 -1 
 
 m,m-l m-1 m-l,m m,m+l m+1 nH-l,m 
 
 Proof ; Suppose that 
 
 x -»il x I|--| x ^ 
 
 is an eigenvector of £,„ belonging to the eigenvalue K It follows from 
 
 (2.4,1) that for 1< i < q, KK. ..• ft, UB. n X. _ + B. . n X. , ] 4 I-ft >X. 
 
 - y l i 1,1-1 i-I i,i+l i+1 i i 
 
 and hence JX --■ 0, where G is 
 
 -u )-\ ft B 
 
 'll' 1 J 1 12 
 
 v 
 
 2 21 
 
 \ft B 
 
 q q,q-i 
 
 ft ,B . 
 
 q-i q-l,q 
 
 q q' 
 
 
 Here I. is the identity matrix of order r , 
 i i 
 
30 
 
 There exists a non-trivial solution of GX - if and only if 
 det(G) = 0, and since detiG) is a polynomial of degree n in \, it follows 
 that its roots are precisely the eigenvalues of £ „ We now show that 
 if the ft.'s satisfy (i), ii ; , or (iii), then det(G) ~ +\ 
 
 If (i) is satisfied,, then the block in the (l, l) position of G 
 
 is -VI „ We now perform the block elementary row operation of adding 
 
 Q 13 times the first r n rows of G to the next r rows , That is , we add 
 2 21 1 2 ' 
 
 n 2 B 21 (-M 1 ) to Mi 2 B 21 and ^ 21 ^ 12 to [(i - I^) - Xl^]. This operation 
 leaves det' G ; unaltered since it corresponds to premultipli cation of G 
 by 
 
 I 
 
 21 
 
 L 
 
 o 
 
 whose determinant is obviously 1„ 
 
 block in the (2,l) position is now null,, while the block 
 in the (2,2) position is [ft B ft B + (i-fl )] - \1 „ The quantity in 
 square brackets vanishes because ft satisfies the recursion relation 
 and so the diagonal block in the .2,2) position is -XI Hence., we 
 repeat the procedure, adding ^B times the second r rows to the next 
 r rows, The sub -diagonal block in the (3,2) position vanishes^ while 
 the diagonal block in the (3,3) position becomes -M since ft satisfies 
 
 t ■ . P} continuing in this waj. , all sub-diagonal blocks of G are made to 
 vanish (i.e., G is made block upper triangular}, and the i-t'h diagonal 
 block becomes -\I.„ This procedure leaves det(G) unaltered,, and 
 det (j) t \ „ Hence £, is nilpctent. 
 
31 
 
 case the ft.'s satisfy (ii), we proceed by block column 
 
 operations starting from the right to reduce 1 to the same upper 
 
 triangular form. 
 
 If (iii) is satisfied, we proceed by block row operations 
 
 starting from above to eliminate the sub-diagonal blocks in positions 
 
 (2,lV (3,2), (m,m-l) and by block column operations starting from 
 
 the right to eliminate the remaining sub-diagonal blocks, Again the i-th 
 
 diagonal block becomes -XI. . Hence, if one of (1 Li) or ^iii) is 
 
 satisfied, det(G) - + A. , and so JlL is nilpotent in all cases, 
 
 QoE D„ 
 
 If Jo is nilpotent, then the error in the successive iterates 
 
 eventually becomes zero in the absence of rounding error, and so the 
 
 exact solution to (2 a l„l) is obtained after a finite number of 
 
 iterations. By the Cay ley -Hamilton theorem at most n iterations are 
 
 required, and in fact it can be shown tnat q iterations suffice, Hence, 
 
 the iteration !2,U,l) with the fi. 's chosen according to the criteria of 
 
 of theorem 2„U,1 is more properly classed as a direct method than as an 
 
 iterative one, and direct methods are net the principal concern of this 
 
 tnesis, In Chapter 3- however, we will use tne criteria of theorem 2, 4,1 
 
 not to obtain the ft. 's themselves, but to obtain certain scalar 
 
 acceleration parameters used in the actual iteration, 
 
 2,5 Extra polat ion Matrices_When the_B lock Matrices Have a Common Basis 
 of Eigenv ectors 
 
 Consider now the linear system (2,1,1, arising from the 
 discretization of the Dirichlet problem on a rectangular domain R for a 
 self -aaj.omt , elliptic partial differential equation of tne special form 
 
32 
 
 ^ f x) ^^zl) _ d (f (y) ^lls 
 
 ox 1 Ox Oy 2 oy 
 
 - (c, ,x: - a_(y))u(x 5 y) = g^x,y), 
 
 (2.5.1) 
 
 where f x f ;. , a (x), o y are continuous in the closure of R 
 
 and satisfy f [x) > 0, f (y) > ; o (x) > 0, o (y) > 0. Let a uniform 
 
 mesh be imposed on R, and suppose that the mesh points are numbered 
 
 either by rows or by columns starting from the upper left. Thus, the 
 
 mesh points .and hence the unknowns in (2„l.l)) are grouped into q blocks 
 
 of r elements each. It follows that for I < i < q, the blocks B. , n and 
 
 - — i,i-l 
 
 B. . _ of the block Jacobi matrix derived from A are square of order r and 
 
 i i+l 
 
 all have a common basis of eigenvectors u„e„, a set of r linearly 
 
 independent eigenvectors.,, In this case the matrices [ft } can be chosen 
 
 in various ways so that they share this common basis of eigenvectors, and 
 
 consequently a useful factorization of the characteristic polynomial of 
 
 £ is possible o In Chapter 3 we will investigate some iterations which 
 
 are equivalent to the use of extrapolation matrices having the same 
 
 eigenvectors as the matrices B. . . and B . . . . The matrices {ft.} will 
 
 i.i-l i,i+l i 
 
 then be determined by specifying certain of their eigenvalues. 
 
 Let t . ., , 6 . . , , and ^ denote the eigenvalues of B. . _, 
 
 i . i - 1 i.i-l i i,i-l' 
 
 B. . , and ft. respectively to which the I r-dimensional) eigenvector i 
 1 5 1+ 1 l 
 
 belongs, and consider the scalar iteration 
 
 1 1,1-11-1 l,lrl i i' 
 
 1 < i < q„ 
 
33 
 Let Ju\ note the iteration matrix for (2.5.2), and let $ (\) denote 
 
 the characteristic polynomial of £r , We then have the following 
 
 decomposition theorem, to be used in Chapter 3° 
 
 Theorem 2„5,,1; If the matrices B. . .. , B. , _. and ft. have a common 
 
 <L L _^ 1,1+1 1 
 
 basis of eigenvectors for 1 < i < q, then the characteristic polynomial 
 
 r 
 of £ n is <|>(\) = n <B J (X). 
 
 Proof • If x .: [x , x , — , x ] is an eigenvector of iV\ belonging to 
 the eigenvalue \, then it follows from (2,5 » 2) that 
 
 \x. = w?(\£. x. _ + p. -X. _) 4 L-w?)x., 1 < i < q, 
 
 1 11, 1-1 1-1 L , 1+1 1+1 % 11 - - 
 
 Hence G J X --■■ 0, where 
 
 G J = 
 
 [(1-^ - . 
 
 2 21 
 
 q q-q-1 
 
 q-1 q-l,q 
 
 L-oo J )-\l 
 
 q 
 
 follows that A. is an eigenvalue of JlTj if and only if \ is a root of 
 0, and therefore $ J (\) = del 
 Let H be the square matrix of order r whose j-th column is % , 
 
 q 
 
 and let H= E©H, where © denotes the direct sum. The linear 
 independence of E , £ , — - £ implies that H is non-singular, and so is 
 letl j) } wn-re G was defined in the proof of theorem 2.^.1, and 
 it follows that 4>(\) = det (H G is block tridiagonal , and each 
 non-null block of G is square of order r Hence, the effect of 
 
3h 
 
 a of G by H is to transform each non-null block of G by H, 
 and diagonal! ze it. Thus, 
 
 H~Vb, 
 
 1, l-l 
 
 H = 
 
 and 
 
 H fi.B. n H 
 
 i i ; 1- 1 
 
 [(l-w})-\] 
 
 H* [(i-fi.)-X" 
 
 [(l-<)-\] 
 
 Hence, H is an shown in Figure k, . 
 
 o be the permutation which selects the first element of the 
 first block first, the first element of the second block second, and so on 
 through the first elements of the various blocks^ then likewise through 
 the second elements of the various blocks, etc. If P is the permutation 
 
 _T .Tl 
 
 matrix corresponding to a, then cj>' ,} ■■-- det(P H GHP )„ But 
 E H = £ © G J , and hence d>(\) = IT det(G J i = TT J (\) o 
 
 
 j=l 
 
 j-i 
 
 Q,E D, 
 
35 
 
 0) 
 
 M 
 
 •H 
 
36 
 
 3 c SEQUENTIAL EXTRAPOLATED IMPLICIT METHODS 
 
 3,1 A Model Problem 
 
 The ususl discretization of the Dirichlet problem for Poisson's 
 equation 
 
 dx dy 
 
 (3.1.1) 
 
 on a rectangle R using a five point star leads to a linear system whose 
 coefficient matrix A is block tri-diagonal with tri-diagonal diagonal 
 blocks. Let a uniform r x q mesh be imposed on the interior of R, and 
 suppose that the mesh points are numbered by columns starting from the 
 upper left „ Then 
 
 D -I 
 
 ■I D 
 
 where I is the identity matrix of order r, D is of order v, and 
 
 h -1 
 
 D ^ 
 
 1 k 
 
 The resulting linear system (2.1.1) will be referred to in what follows 
 as the model problem,, 
 
37 
 For the model problem, each non-null block of the block Jacobi 
 matrix derived from A is simply^ D , and the block Gauss-Seidel iteration 
 for the solution of (2.1 l) is 
 
 m+1 ^-IrTT ro+1 ^ m 
 
 Z 
 
 ^ x = D" X [Z u+ " , Z " + K.], 1< i < q, (3.1.2) 
 
 i i-I i+1 i ' - - * w ' 
 
 or 
 
 _J7 m+ 1 _ m+ 1 „ m v , / _ .. _ \ 
 
 DZ. =sZ._+Z._+K.. 1 < l < q. (3.1.3) 
 
 i i-1 l+l i' - - * 
 
 The block SOR iteration obtained by extrapolating (3.1.2) is frequently 
 called SLOR (successive line over-relaxation) since it consists of updating 
 simultaneously the unknowns corresponding to an entire vertical line of 
 mesh poin~. 
 
 3.2 The Sequential Extrapolated Implicit Method (SEl) 
 
 ~* m+1 ~* m ~~ 
 If Z - Z = Zj the exact solution of the model problem, 
 
 then for any scalar s, 
 
 / • l~ m+1 ^" m+ 1 ^m Zm rT . / \ 
 
 (D + sI;Z = Z . . + Z . _ + sIZ. + K. . 1 < l < q, (3.2.1) 
 
 or 
 
 Z m+1 = CD + si)" 1 [Z ^ + Z m , + sIZ m + K. ]. (3.2.2) 
 i i-I i+I i i ' ' 
 
 1 < i < q. 
 The iteration (3.2.2) can be written as 
 
 z m+1 = qd" 1 [z mtl + z m n h it. ] + (i - n)z m 2.3) 
 
 i i-l i+l i i J ~" 
 
 i < i < 0. ; 
 
38 
 
 where ft = (D + si)" D. That is, (3=2 2) is the VSOR iteration (2.1+.1) 
 
 i ■ "I 
 
 with ft. = (.D + si) D„ 1 < i < q„ We denote the iteration matrix for 
 
 f3»2.2) by S. 
 
 UJ 
 
 The most efficient way to perform SLOR is to use (3.1.2) to compute 
 
 1 
 — m + p 
 
 Z . " and then to compute 
 
 1 
 
 2 
 Z 1 = <4Z. "Zj+Z.^l^i^q, (3=2,4) 
 
 — . ni -t- — _ 
 m+1 ,, r T 2 "mi "m 
 
 Hence, the additional work required to perform SLOR over that required 
 for Gauss-Seidel iteration is one multiplication and two additions per 
 unknown, 
 
 •ation (3°2„2) requires that (D + si) " be computed instead 
 of D , but both of these matrix inversions require the same amount of 
 work since the diagonal of D is non-zero Therefore the only additional 
 work required to perform SET over that for Gauss-Seidel iteration is one 
 multiplication and one addition per unknown in order to compute the vector 
 on the right hand side of (. 3o2„l)„ Hence SEI requires one less addition 
 per unknown per iteration than SLOR„ 
 
 Any eigenvector of D is also an eigenvector of ft,, and since the 
 eigenvectors of D form a basis, the converse also holds. Hence, the 
 decomposition theorem 2 5»1 applies. 
 
 Let d be an eigenvalue of D and % an eigenvector belonging to 
 d, If (3 denotes the corresponding eigenvalue of D , then (3 = —, and the 
 eigenvalue of ft = (D +• si) D to which % belongs is - — r~ , which we 
 denote by w(g ;S ) 
 
 
39 
 
 Let JlT denote the matrix for the (scalar) iteration 
 
 x^ +1 = <o(p)p(x^ + x* ) + (l-u(p))x» L < i < q, (3«2 o5 ) 
 
 and let ty (\) be the characteristic polynomial of £T. It follows from 
 theorem 2.5.1 that if w(p) = w(p,s) for all 3, then the characteristic 
 polynomial of S^ is ty(\) = TT\|r (XJ, where the product is taken over all 
 
 eigenvalues of D 
 
 ft 
 Let P-Vo'i denote the maximum of the moduli of the zeroes of 
 
 ft ft 
 
 \|/ (\) and let w (ft) denote the value of w (f3) which minimizes p fa \. It 
 
 follows from the theory of SOR that for w(p) > w h (P); the zeroes of 
 
 \|r (\) are complex of modulus (w(ft)-i),, 
 
 For fixed s, let p(S, ,) denote the spectral radius of S, lC If 
 
 'CO c CO 
 
 §_ and ft are the smallest and largest eigenvalues respectively of D } 
 
 then p(S ) = max __ p , », and we let p = min p(S ). 
 R_ < ft < ft p ' by ox,x s 
 
 The eigenvalues of D are 
 
 and so 
 
 P, - 
 
 U - 2cos( 
 
 rJIL^ 
 
 r+l' 
 
 \ < & < Pj < P < | 
 
 1 < J < r > 
 
 Let 
 
 A = 
 
 1 
 
be square of order q The eigenvalues of A are CC. = 2cos —r } 1 < i < q_ f 
 1 since the block Jacob! matrix derived from A is M = A @ D where 
 @ denotes the tensor product^ it follows that the eigenvalues of M are 
 \ a a.B.j 1 < i < q_ s 1 < j < r Moreover,, the theory of SOR applied to 
 the iteration (3o2„5) shows that 
 
 V- p) ■ — 7— . (3.2.6) 
 
 1 + \[i-a a p 2 
 
 where OL =a max {CH.} „ Since the largest positive eigenvalue of M is ji = a Pj 
 
 we have 
 
 b b VH; 
 
 + \pV 
 
 r: VSOR iteration (2.U.1) becomes SOR lif fl. a w I 1 < i < q, 
 
 Ls case,, theorem 2 5 o l shows that the characteristic polynomial of 
 
 8 • B 
 
 £ is <£>(a) 3 TI O^'.a. 1 j, where cj> (A.) is the characteristic polynomial 
 
 a b B o 
 
 of JlT with 1 to The maximum of the moduli of the zeroes of $ (\) 
 
 B 
 
 J 
 
 t 
 
 is p w , 
 
 Let s^tO =: - b~7pT - x ? and denote s b (p)by s bD Then for 
 
 Lemma 
 
 j_b 
 B < B- u b (p) < w(B,s b ) < w fe (3) a u^, 
 
 oof s s (6) < since ^.(B) > 1, and so 1 + s J < 1 + s B if < B < P, 
 
 Hence, w(B ; s ) < w b (B) Substituting (3„2„6) into the formula for s (p) 
 
 and differentiating shows that -rr s. (b) < 0. and hence s, < s, (b) for 
 
 dB b ' b b 
 
 B < B . Then 
 
 b h s 
 
 an 
 
 d so u) (B) < w (B 5 s, ) 
 
 Q.E.D, 
 
Ul 
 
 Theorem 3.2.1 : If s = s, , then: 
 ' b 
 
 R R R 
 
 (i) r (X.) = $ (A.), and hence p^ ^ = p g0R 
 
 b ' 
 
 (it) For P <e, p^ (p) <ol {&)%) <P S0R 
 
 (iiij lim p , \ = for any s. 
 R.*0 ^> S) 
 
 Proof : 
 
 (i) By definition, s is that value of s such that 
 
 — — — 6* R 
 
 <*>(0,s) = w b (P)» But w b (P) = °* h > and so V (M = $ (^). Since the moduli 
 
 o 
 
 of the roots of <$> (\) are equal to p C( ~ D for all B < (3, it follows that 
 
 B 
 
 P w (0) - P SOR : 
 b 
 
 (ii) It follows from lemma 3«2.1 and the theory of SOR that the 
 
 roots of ty (\) are complex of modulus (a)(f3,s )-l) for < R, and 
 
 (w (B)-l)) < (w(0,s )-l) < ( w ^P)-l). The assertion (ii) then follows 
 
 from the definitions of the quantities involved, 
 
 (iii) lim w(B,s) = 1, and since p , \ = (w(0,s)-l), we have 
 fi->0 \P> S . 
 
 lim p / „\ = 0, 
 
 B^O W( ^ S) Q.E.D. 
 
 Corrollary 3,2,1 : P SEI * P S(}R . 
 
 Proof ; The theory of SOR applied to J^ shows that P^/q \ is minimal if 
 
 and only if w(r s ) = w (r) which is the case if and only if s a s. , It 
 
 b b 
 
 then follows from part (ii) of theorem 3.2.1 that max Pm(a 1 ^ s 
 
 R_ < R < R lP ' Sj 
 minimal if and only if s = s , and part (i) of the theorem shows that 
 
 P SEI ''"' p S0R 
 
 Q.E,D, 
 
1*2 
 
 t | ' normalized to have (Euclidean) length 1 be the 
 
 -1 ~"* I r 
 
 eigenvector of D belonging to P.j 1 < J £ r ° ^he vectors { | } form 
 
 J - 1 - 
 
 an ortho-normal basis for r-space since D is symmetric., and indeed 
 
 'r+1 ' r+1 ' 3 r+1 
 
 t e " denote the q-dimensional column vector having 1 in the i-th 
 position and 1 s elsewhere,, The vectors [e } _r form an ortho-normal basis 
 for q-space,, and therefore the vectors [e ' @ £ '\)^-S.^-S. (: lj^- < S.^'S. x 's 
 form an ortho-normal basis for qr-space,, where © denotes the tensor 
 product o 
 
 eigenvector of D " corresponding to f3 is | , and it 
 follows from part (i) of theorem 3.2,1 that if the initial error E lies 
 in the subspace spanned by {e @ % } then the error E in the m-th 
 iterate for SE. r is identical to the error in the m-th iterate for SLOR,, Vm, 
 
 Because of part of the same theorem, however,, it is to be expected 
 
 -* 
 that if E does not lie entirely in this subspace,, then a fixed number of 
 
 iterations with SEI will produce a greater reduction in the I norm of the 
 
 error than the same number of iterations with SLOR,, Numerical results 
 
 which support this expectation are presented in Section 3-^. 
 
 3 . 3 Cv clic Chebyshev SEI 
 
 The cyclic C hebyshev semi -ite rative method [6 j 11, Chap „ 5 ] for 
 the solution of the model problem is 
 
 z. = u) d [z._+z,_+K.J + (l-w ;z . , (3,3,1a; 
 1 1-1 1+1 1 1 
 
 1 - : t, 3 ,, 2 ,, ' J 
 
><3 
 
 Z m+1 = c 2m+2 D 4 [Z m+ ! + Z "*} + K.] + (l-^ m+2 )Z m , (3.3.1b) 
 1 l-l l+l i l' 
 
 i = 2, h, 6, --- , 
 
 r m i 
 
 where the sequence [u J satisfies 
 
 1 , 2 1 
 
 W = 1. W = 
 
 ) 
 
 -2 
 
 1- L 
 
 1-g-eJ" 
 
 (3.3.2) 
 
 and u = p(M) , 
 
 In what follows (3»3.l) with the sequence {w } satisfying 
 
 (3.3.2) will be called cyclic Chebyshev SLOR , while (3.3.l) with 
 
 m 
 
 10 = OJ 
 
 b' 
 
 Vm, will be called SLOR with the odd-even ordering . Iteration (2,2.3) 
 
 with w = to will be called SLOR with the natural ordering,, 
 b — " 
 
 If the sequence fw } satisfies (3.3.2), then lim w = w and 
 ^ l ' \ -><■-> I? m->°° b 
 
 it follows that the asymptotic rate of convergence of cyclic Chebyshev SLOR 
 is equal to that of SL0R< However, if E denotes the error in the m-th 
 iterate for cyclic Chebyshev SLOR and if E denotes the error in the m-th 
 iterate for SLOR (with any consistent ordering), then Golub and Varga [6] 
 have shown that 
 
 max | |E m | | < max| |E m | | , m < 2, 
 
 ;o - o 
 z z 
 
 where denotes the I norm. Moreover, the sequence (||E ||) is 
 
 mi 
 
 strictly decreasing., whereas the sequence (j|E j} may increase 
 initially. 
 
ts m } by 
 
 Let the sequence {w ) satisfy (3=3=2), and define a sequence 
 
 S m = I [-i - 1], m> 1 . (3,3,3) 
 
 For < & < f>. define the sequence {w (f3,s )} by 
 
 m,„ m N 1 
 
 > (p, s m ) - x m , m> l. (3.3.U) 
 
 1 + si 
 
 Then 
 
 m /o „m\ m v_/ 
 
 The iteration 
 
 m- 1 / 2m-l M -l r m _ m 2m+l _ m v -, 
 
 Z I D + s I J Z , n + Z , , + s Z . + K , J , 
 
 i l-l l+l i i ' 
 
 1 = 13 5— (3 < 3 ' 5a) 
 
 ~Z m- 1 /_ 2m+2 7 >-l r r" m+1 ^ nul 2m+2 r" m IT 1 
 
 Z DtS 1) [Z-t+Z -,+s Z . + K, J 
 
 i v l-l l+l l i 
 
 i . 2, k, 6, - (3.3.5b) 
 
 for the solution of the model problem will be called cyclic Chebyshev SEI , 
 Iteration ( 3 - 3 - 5 ) with s replaced by s,, Vm., will be called SEI with the 
 od d -ev en ordering; while (3,2,2) with s - s will be called SEI with the 
 natural ordering. 
 
 The remarks in Section 3 2 concerning the convergence properties 
 of SEI relative to those of SLOR apply also to the convergence properties 
 of cyclic Chebysnev SEI relative to those of cyclic Chebyshev SLOR D 
 
^5 
 
 3.^4 Leal Results 
 
 Numerical experiments were performed to compare the actual 
 rates of convergence of SLOR and SEI with both the natural and I 
 odd-even ordering, cyclic Chebyshev SLOR, and cyclic Chebyshev SEI. The 
 test problem selected was the model problem with f(x,y) = (i.e., 
 Laplace's equation) and homogeneous boundary conditions. The rectangle 
 R was taken to be a square containing N mesh points on a side. For 
 N = 10, 20, 30, and Uo the number of iterations with each method needed 
 to reduce the I norm of the error below 1 was determined. 
 
 Starting vectors were constructed as follows. Let Q, be the 
 
 N-dimensional vector whose j-th component is (-l) , and let J be 
 
 _ N _ , 
 the N-dimensional vector whose j-th component is j. Let r\ - ^ £ , 
 
 ^1^2 -* N -1 ^ =1 
 
 where | , £ s , £ are the eigenvectors of D . Experiments were 
 
 performed using each of the three starting vectors Z - J © Q , 
 
 Z = 10 Q © Q„ and Z = J © r\. 
 
 The graphs shown in Figures 5 , 6, and 7 were constructed by 
 
 interpolating linearly between the data points, It can be seen from 
 
 these graphs that SEI converges much more rapidly than SLOR for certain 
 
 starting vectors but that the relative advantage of SEI over SLOR is 
 
 strongly dependent on the starting vector chosen. 
 
ber of Iterations Needed to Reduce 
 ti of Error below 1 for N x N Test, 
 nlem, 
 
 j ®5 
 
 60 
 
 to 
 
 a 
 o 
 
 ■H 
 
 p 
 
 CD 
 -P 
 
 O 
 
 X2 
 
 50 - 
 
 SLOR, Natural Ordering 
 
 SEI„ Natural Ordering 
 
 SLOR, Odd-Even Ordering 
 
 SEI, Odd-Even Ordering 
 
 Chebyshev SLOR 
 o o 
 
 Chebysnev SEI 
 
 -> N 
 
 10 
 
 20 
 
 30 
 
 ko 
 
 Figure 5 
 
hi 
 
 nber of Iterations Needed to Reduce ft 
 Norm of Error Below 1 for N x'N Test 
 blem, 
 
 Z° -- 10 Q ® Q 
 
 2 
 
 CO 
 
 C 
 
 o 
 
 -I 
 -p 
 
 CD 
 S- 
 CD 
 -P 
 
 Cm 
 
 - 
 
 SLOR^ Natural Ordering 
 
 SEI, Natural Ordering 
 
 SLOR, Odd-Even Ordering 
 
 o o 
 
 SEI, Odd-Even Ordering 
 
 Chebyshev SLOR 
 o o o 
 
 Chebyshev SET 
 O- - 
 
 o 
 
 ..^ N 
 
 Figure 6 
 
48 
 
 Number of Iterations Needed to Reduce 
 Norm of Error Below 1 for N x N Test 
 
 Z - J <g n 
 
 w 
 
 G 
 O 
 
 CD 
 U 
 
 cu 
 -P 
 
 ■ 
 
 CD 
 
 SLOR„ Natural Ordering 
 
 SET, Nat-oral Ordering 
 
 SLOR, Odd -Even Ordering 
 
 D D O 
 
 SEJ , Odd-Even Ordering 
 
 Chebyshev SIOR 
 o o o 
 
 Chebyshev SE ' 
 
 o— — — -o — -o 
 
 ->N 
 
 Figure 7 
 
k 9 
 
 3 . 5 The Variable Sequential Extrapolated Implicit Method (VSE l) 
 
 Consider the following generalization of (3.2 2) for the 
 solution of the model problem. 
 
 _ m+1 , , ,_ T \-l r _ m+1 _ m T . „ ran 
 Z. u) (D + s.I ) [Z,_+Z J _+K + s.Z . J 
 
 1 i l l-l i+l i li 
 
 + (l - u.) Z . , 1 < i < q ; 
 
 (3.5.1) 
 
 where ( w .), and {s.}, are any two sequences of real numbers, The special 
 case of (3.5.l) in which w_ = ]_,, Vi, will be called 1-parameter VSEI 
 whereas the general case will be called 2-parameter VSEI „ 
 Iteration (3.5=l) can be written as 
 
 I* m+1 n -l r r*m+l I'm -, /_ \^ m 
 Z . u.u IZ._+Z._J+(I-fi.;Z. J 
 
 1 1 l-l i+l 1 1 
 
 (3.5.2) 
 1 < i < <L, 
 
 where SI = w.(D + s.l) D„ Hence ( 3.5.1*) is a special case of the VSOR 
 111 
 
 iteration (2,U.l) The iteration matrix for (3«5°l) will be denoted by 
 
 The work required to perform 1-parameter VSEI is the same as 
 for SET, The most efficient way to perform 2-parameter VSEI is to 
 compute 
 
 1 
 
 p» 2 fT . _ \ - 1 r _ m+ 1 _ m _ m „-, 
 
 Z = ' D + s. I) [Z n +Z.^+sZ.4-K.J, 
 
 1 i-1 i+l 11 1 
 
 (3.5.3) 
 
 1 < i < q , 
 and then to compute 
 
SO 
 
 1 
 
 z '" = ". [Z . m+2 - Z m ] + z m , l <: i < q , (3.5.U) 
 
 ii i i' — — 
 
 e the additional work required to perform 2-parameter VSEI over that 
 
 required for SLOR is one multiplication and one addition per unknown per 
 
 iteration. 
 
 The eigenvectors of W .(D + s.l) D are the same as those of 
 ° ii 
 
 -1 - 
 
 D j, and so the decomposition theorem 2„5ol applies. If £ is an eigen- 
 vector of D " belonging to the eigenvalue B, then the eigenvalue of 
 D + s.l) D to which f belongs is 
 
 i 
 
 1 -* s.B 
 
 which we denote by w . (B,s.) 
 
 Let £L denote the matrix for the (scalar) iteration 
 ft 
 
 x m+1 = uWB,s.)e(X m ^ f x m J 4 (l-U).(p,s.))X m , 
 
 1 11 1-1 1+1' 11 1 
 
 i < . i < q. , 
 
 (3.5.5) 
 
 R R 
 
 and let a- (\) denote the characteristic polynomial of £ „ It follows 
 
 from theorem 2,5=1 that the characteristic polynomial of S is 
 
 R -1 
 
 • IT ^ ( A -) ? where the product is taken over all eigenvalues of D 
 
 Let p denote the spectral radius of Jl~ „ and let P vqT? - denote 
 
 the spectral radius of S n „ Then p. .__... = max p , where 8 and B 
 
 £ < 3 < F 
 denote the smallest and largest eigenvalues of D " respectively , 
 
 w (B,s ) is specified for two distinct values of 8, 8 n and 
 l l 1 
 
 8„ say, then s„ and w are determined by the formulas 
 2 ii 
 
51 
 
 w. : (l 4 s i P 1 ) u i i 'S 1 ,s.) . (3,5.6b) 
 
 P 6 
 3,6 VSEI with d^ 1 - p Q 2 -- 
 
 The criteria of theorem 2 c U c i with fi = w/ft s ) and 
 
 1 L ' 1 
 
 q 
 
 B. . n = B. , = P can be used to determine sequences f w .(P n3 s )}_, and 
 i,i-l i,i+l i x 1 J i 1 
 
 - q P l P 2 
 
 & ,s.,')., such that £ and £, are nilpotent; i c e c1 such that 
 
 ^ J- -L 06 Ob 
 
 to be used in [3.5. l)« It then follows from (.3.5 »6b) and theorem 2 C 4,1 
 
 o = p --- o Then (3»5=6) determines the sequences {s.} and l.^." 1 
 
 that I -parameter VSEI is the special case (3. = of 2-parameter VSEI 
 
 It can be shown that p is the same function of p if the 
 
 criteria of either part (i) or part 11,, of theorem 2,4,1 are used to 
 
 determine ( w .(P, s s.)l and f w .(P ,s.)} . in which case (3c5„l) will be 
 i ' 1" £ '1 i 2 i 1 
 
 called uni -directional VSEI . If the criteria of part (iii) of the 
 tneorem are used, (3.5.l) will be called b i-directional VS E[ 
 
 /.'as found experimentally that if P : P then max p 
 
 increases as P increases while max p decreases, Moreover 
 
 P L < P < P Q 
 
 P P 
 max p„ -: pE-. it follows that if p.. is the values of p. such that 
 
 P < 6 " " lb X 
 
 _ _ _ ^ 
 
 P P P 
 
 max p = pE;, then max p^ is minimal for P n - p ni , 
 
 P x < P < P ° ° & < P < P X lb 
 
For the model problem, p = p^ /-v, where p^/ R \ was defined 
 
 b 
 
 3 3 
 
 nor 
 
 in Section 3=2 The theory of SOR shows that p ,-> ,— \ 
 
 -dp 
 
 < P < Pj and — p^ ,r\ = + oo for p = P + „ The approximate shape of 
 
 the graph of p as a function of P is shown in Figure 8, and the graph of 
 
 bh (3 = P o is shown for compari! 
 b 
 
 Fj gure 8 
 
 For the model problem. § > ?- and (3 < —j Vr„ Moreover 
 
 lim £ = g- and lim p - Hence P YSEJ 
 
 max 
 
 r 
 
 r -■>«. 
 
 ft. < P < & 
 
 u 
 
 P P 
 
 p < max p . 
 
 ~ 1 P - 1 
 
 6 - p 
 
 iteration I s) with u.(3,s. ) = ; ■ — , 1 < l < q. was performed to 
 
 i ' l 1 + s.P ; ~ - 
 
 ! 
 
 Q 
 
 determine p as a function of P for q = 10, 19, and 33, Formula (3.5.6) 
 
 with P = ~ W5£- used to determine ( s . } _r and f w .} , For the ^-parameter 
 c. 1 J i J 
 
 case p was determined by numerical search. The results of these 
 
 !■■■ r ' m< tit J: • x e ih esented Ln Dab] < 1-U 
 
53 
 
 Table 1. 1-Parameter, Uni-Direction VSEI 
 
 » 
 
 P VSEI 
 
 P S0R 
 
 ^SEI 
 
 R S0R 
 
 10 
 
 .55 
 
 .56 
 
 .60 
 
 .58 
 
 19 
 
 .7»» 
 
 .73 
 
 ,30 
 
 ,31 
 
 33 
 
 -85 
 
 = 83 
 
 .16 
 
 .19 
 
 Table 2, 1-Parameter, Bi-Directional VSEI 
 
 q 
 
 P VSEI 
 
 p S0R 
 
 R SEI 
 
 R S0R 
 
 10 
 
 M 
 
 ,56 
 
 .82 
 
 = 58 
 
 19 
 
 ,65 
 
 = 73 
 
 M 
 
 = 31 
 
 33 
 
 o78 
 
 ,83 
 
 = 25 
 
 = 19 
 
 Table 3, 2-Psrameter, Uni-Directional VSEI 
 
 q 
 
 \b 
 
 P VSEI 
 
 P S0R 
 
 R VSEI 
 
 R S0R 
 
 10 
 
 .Ul 
 
 M 
 
 .,56 
 
 .81< 
 
 -58 
 
 19 
 
 A7 
 
 ,6k 
 
 .73 
 
 M 
 
 .31 
 
 33 
 
 .*9 
 
 = 78 
 
 = 83 
 
 = 25 
 
 ,19 
 
54 
 
 Table h c 2«Parameter n Bi-Directional VSEI 
 
 q 
 
 p ib 
 
 
 
 P VSEI 
 
 p S0R 
 
 ^SEI 
 
 r i 
 
 SOR i 
 
 10 
 
 .3^ 
 
 o32 
 
 .56 
 
 1.1 
 
 ,a 
 
 19 
 
 .hi 
 
 s . .... , 
 
 ,73 
 
 o62 
 
 •31 
 
 1 r > 
 
 M 
 
 • ,69 .83 
 
 1 .i .... 
 
 .37 
 
 = 19 
 
 For bi-directional VSEI,, best results were obtained with m in 
 part 1 i of theorem 2 4ol equal to the greatest integer not exceeding 
 '^r~i and the results presented are for this case. In all cases except 
 that of 1-parameter, uni-directional VSEI, the rate of convergence of 
 VSE] for the model problem exceeds that of SLOR, and in the case of 
 2-parameter VSE; r , the improvement is substantial. 
 
. CONCLUSION 
 
 Summary 
 
 Several schemes have been investigated for impn the rate 
 of convergence of extrapolated Gauss-Seidel iteration by the introduction 
 of a multiplicity of extrapolation parameters to replace the single 
 scalar parameter used by SOR. In Chapter 2 it was shown that the use of 
 two optimally chosen extrapolation factors results in an improved rate 
 of convergence for linear systems whose coefficient matrices are 
 consistently ordered S-matrices In the case of linear systems arising 
 in the numerical solution of boundary value problems for elliptic- 
 partial differential equations the improvement is small because u_ for 
 these systems is small. It is evident, however, that there exist linear 
 systems for which u_ and u are more nearly equal, and for such systems 
 the use of two optimally chosen factors offers a substantial advantage. 
 
 It was also shown in Chapter 2 that by use of matrices 
 rather than scalars as extrapolation parameters, an extrapolated 
 Gauss-Seidel iteration having an infinite rate of convergence can be 
 constructed, although its implementation requires more work per 
 iteration than extrapolation by a scalar. 
 
 The SEI and VSEI methods of Chapter 3 were analyzed for a 
 more limited class of linear systems, namely those arising from the 
 discretization of the Dirichlet problem for (3.1.].) on a rectangu] 
 domain , In this case SEI was shown to be equivalent to several 
 simultaneous SOR iterations, each on a different sub 
 
 a different scalar extrapolation factor. Moreover, SEI was shown to 
 have convergence properties superior to those of SOR for the problems 
 considered and to require less work. 
 
56 
 .- investigation of VSEI was largely experimental, but the 
 results indicate that for the class of problems considered it is superior 
 to SOR not only with respect to average rate of convergence but also with 
 respect to asymptotic rate of convergence. 
 
 k c 2 Extensions 
 
 It seems probable that the results of Chapter 3 can be 
 extended to linear systems associated with the more general self-adjoint,, 
 elliptic partial differential equation (2„5.l) since the non-null blocks 
 of the Jacobi matrix in this case still have a common basis of 
 eigenvectors though different eigenvalues. Extensions to boundary value 
 problems for non-rectangular regions appear more difficult, but such 
 extensions would be very desirable. 
 
LIST OF REFERENCES 
 
 Arms. R. Zondek, B Method of Block Iteration ' 
 
 Journal of the Society for Industrial and Applied Mathematlr 
 Vol, U, (1956). pp c 220-229. 
 
 Forsythe, G ; E., Wasow, W R, Finite Difference Methods for Partial 
 Differential Equations , John Wiley and Sons, Inc. , New York, 1960, 
 
 Forsythe, G. E e , "Gauss to Gerling on Relaxation," Mathematic a l Ta bles 
 an d Other Aids to Computatio n-, Vol. 5 (1951), pp. 255-258, 
 
 Frankel. S. P., "Convergence Rates of Iterative Treatments of Part 1 - 
 Differential Equations/' Mathematical Tables and O ther Aids to 
 Computation , Vol. h (1950), pp, 65-75, 
 
 Geiringer, H„, "On the Solution of Linear Equations by Certain Iterative 
 Methods/ 1 Reiss ne r Anniversary Volume , J = W. Edwards. Ann Arbor, 
 Michigan, (l9^9), PP = 3^393- 
 
 Golub- G, H., Varga, R : S = , "Chebyshev Semi-Iterative Methods, 
 Successive Over-Relaxation Iterative Methods, and Second Order 
 Richardson Iterative Methods I, II," Numerische Mathematik , Vol. 3 
 (1961), pp r 1U7-168, 
 
 Jacobi, C e G. J., "Uber ein Leichtes Verfahren die in der Theorie der 
 Sacularstorungen vorkommenden Gleichungen numerisch auf zulbsen," 
 Jo urnal fur die reine und angewandte Mathe matik. Band 30 (18^6 
 
 Kahan, W. M,, "Gauss-Seidel Methods of Solving Large Systems of Linear 
 Equations," Doctoral Dissertation, University of Toronto, Toronto - 
 Canada ; 1958. 
 
 9-. Ostrowski, A, M. . "On the Linear Iteration Procedures for Symmetric 
 
 Matrices." Rendiconti di Matema tica e delle sue Applicazioni , Vol. lh 
 pp.. 1UO-163, 
 
 Seidell L. - "Uber ein Verfahren. die Gleichungen, auf welche die 
 Methode der kleinsten Quadrate fuhrt, sowie lineare Gleichungen 
 uberhaupt. durch successive Annaherung aufzulosen, " Abhandlungen der 
 m athematisch-physikalischen Classe der Bayerischen Akadamie de r 
 Wissenschaf ten , Band 11 (187M, 3-te Abtheilung, pp ; 11-108, 
 
 11, Varga, R, S,, Matrix Iterative Analysis , Prentice -Hall-, Inc., 
 Englewood Cliffs, New Jersey, I962, 
 
 12, Young, D, M., "Iterative Methods for Solving Partial Differential 
 Equations of Elliptic Type," Transactions of the Amer ican Mathe matical 
 Society, Vol. j6 (195*0, pp. 92-111, 
 
 13- Wachspress, E, L c , I terative Solution of Ellipt i c Systems an d 
 
 Appli c ations to the Neu tron Diffusion Equati ons of Reactor Physics . 
 entice-Hallj Inc. Englewood Cliffs, New Jersey, 1966, 
 
58 
 
 VITA 
 
 Leland Kitchin McDowell was born on March 15, 19^0? in 
 Tarboro,, North Carolina, From 1958 "to 1962 he attended North Carolina 
 State University } where he received Bachelor of Science degrees in Applied 
 Mathematics and in Electrical Engineering. In 1963 he received a Master 
 of Science degree in Mathematics from the University of Illinois, where 
 he pursued graduate work from 1962 to 19^7= During this period he held 
 appointments as teaching assistant in the Department of Mathematics and 
 as research assistant in the Department of Computer Science „ 
 
\