LIBRARY OF THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAICN blo.%4 ,/no.42l-42G> // /*-z^ Report No. 423 NUMERICAL STUDIES OF STONE'S SYMMETRIC FACTORIZATION AND THE ITERATION PARAMETERS, a AND T by RICHARD JOSEPH TARANTO January, 1971 THE LIBRARY OFJIHE NOV 9 1972 UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN URBANA, ILLINOIS Digitized by the Internet Archive in 2013 http://archive.org/details/numericalstudies423tara numerical studies of stone's symmetric factorization' and the iteration parameters, a and t BY RICHARD JOSEPH TARANTO A.B., Rutgers University, 19&9 THESIS Submitted in partial fulfillment of the requirements for the degree of Master of Science in Computer Science in the Graduate College of the University of Illinois at Urbana- Champaign, 1971 Urbana, Illinois Ill ACKNOWLEDGEMENT I would like to express my appreciation to Professor Paul E. Saylor for suggesting my thesis topic and patiently guiding my efforts and to my fiancee for her kind understanding. This thesis is dedicated to the memory of my mother and father iv TABLE OF CONTENTS Page I. INTRODUCTION 1 II . MATHEMATICAL PRELIMINARIES 3 III. DESCRIPTION OF TEST PROBLEMS 8 IV. THE STONE PROCEDURE 11 V. THE ADI PROCEDURE 23 VI . RESULTS AND CONCLUSIONS 26 LIST OF REFERENCES 3 1 * I. INTRODUCTION Let Au = s be a linear system of algebraic equations resulting from the discretization of an elliptic boundary value problem. Direct methods of solution that factor A into a lower triangular matrix, L, and an upper tri- angular matrix, U, such that A = LU, are computationally inefficient since the factorization yields unique L and U which are not sparse matrices. For this reason, an iterative procedure is more generally desirable. Iterative procedures have recently been proposed by several people, notably H. L. Stone [3] and [k-~\, and Dupont, Kendall and Rachford [1], based on the idea of forming a matrix (A+B ) which can be factored into the product of sparse triangular matrices, L and U , i.e., (1.1) A+B = L U a a a The solution of Au = s is then computed as the limit of the sequence defined by (1.2) ( A+B a } Vl = ( A+ VV T(A V s) • where a and i are iteration parameters. The procedures of Stone, and Dupont, Kendall and Rachford differ only in the method used to determine (A+B ) . Another efficient iterative technique for solving the linear system Au = s is the alternating direction implicit (ADl) procedure. ADI was ini- tially proposed by Peaceman and Rachford [2]. The idea is to partition A as the sum of two easily invertible matrices H and V, i.e., (1.3) A = H + V . The following iteration scheme is then used to obtain the solution, u: (H+w ih l)u n + l/2 = S " (V - W ih l)u n (V+u. I)u = s - (H-w. I)u . /r> iv u+1 iv ' n+1/2 where to., and w. are iteration parameters, lh iv In his paper [3], Stone described a factorization procedure that proved to have several advantages over those now in use. His studies indi- cated that, in addition to reducing computational effort, the convergence rate did not depend strongly on the nature of the coefficient matrix, A, or the shape of the region being considered and suitable iteration parameters could easily be calculated. For these reasons, his method is the subject of much current research, but, so far, attempts at rigorous analysis have been unsuc- cessful. One of the difficulties is that his factorization replaces A with an asymmetric matrix, (A+B ), and, moreover, this matrix, being dependent on the variable parameter a, changes frequently during iteration. A new factor- ization recently devised by Stone [k] defines a symmetric, positive definite matrix, (A+B,.), to replace the asymmetric one in the procedure. Since sym- metry greatly simplifies matrix algebra, it is reasonable to expect progress in the mathematical analysis of this symmetric procedure. If so, these findings might also be applicable to the asymmetric procedure. However, the symmetric procedure is not just an approach to the asymmetric procedure. If it is competitive to ADI, it too is of practical interest. Therefore, the purpose of this thesis is to test the effectiveness of Stone's symmetric fac- torization by comparing the results obtained with those from both Stone's asymmetric procedure and the ADI procedure . In addition, since proper usage of parameters is important to convergence rates, pertinent data will be compiled to help determine their properties. II. MATHEMATICAL PRELIMINARIES Consider the Dirichlet problem <«•« - 4 (*!«%) - 4 ( a2(x %) = s(x) ' xeD u(x) = f(x), xeD where = (x ,x ), D is the interior of a compact region, and the differential operator, < 2 - ia > - 4 i a i< x ^) - 4 v a 2 (x) 4) ' is elliptic in D. In order to solve this problem numerically we will let D be the interior of the unit square, see Figure 1, i.e., (2.2) D={( X;L ,x x ) :0 and let D , the grid system over D with uniform spacing, be defined by (2.3) I> h b ((ih,jh) : 1 < i, j < n) ; i, j integral. The grid point (in, jh) will be denoted by (i,j) and the boundary dD of EL by (h, 0) (nh, 0) (1,0) Figure 1. Grid Point System '(i,0) , < i < n + 1 (1, j) , < j < n + 1 (2.4) ^\ = < (i,l) , < i < n + 1 (0,5) , < j < n + 1 s. We will replace the differential operator (2.1a) with the difference operator (2.7) and replace u with a vector whose components correspond to the grid points of the system. This yields a set of simultaneous linear equations having a one to one correspondence with the grid points. The grid points will be ordered in two ways and the ordering may alternate from one iteration to the next. The first, called normal ordering, arranges the grid points in a left to right, down to up pattern. Thus, grid point (i, j) precedes grid point (k, j) if i < k, and precedes grid point (l,m) if j < m. The second, called reverse ordering, arranges the grid points in a left to right, up to down pattern. In this way, grid point (i,j) precedes grid point (k, j) if i < k, and precedes grid point (l,m) if j > m. Thus, normal ordering is down to up and reverse ordering is up to down. Note that components of the iteration vector, u. ., must conform to the ordering of the i, J grid points at all times. If u is the iteration vector with components u. ., (i, j) a grid i, J point, let A defined by (2.5) (Au). . = H. .u. . n + F. .u, _ . + E. .u. . i,J i,0 i,J-l i,J i-l, J i,J i,J + F. n ,U. .. . + H. . ..u. . . l+l, J l+l, j i,J+l 1, J+l be a discretization of the differential operator (2.1a). The important properties of A are that (i) H ,F < for 1 < i, j for 1,0 1,0 1,0 1+1, G 1,0+1 1 < i, j < n, if F ^ and H. . J ; (v) E. . > -(H. . + F. . + F. , . + H. . .,) for 1,0 i,j i,o i + i,o 1,0+1 1 < i. j < n. if F. . = or H. . = 0, - ' - ' 1,0 1,0 Note: If i = n, F. ,. . = and if j = n, H. . , = ' 1+1,0 ' 1,0+1 The above properties are a result of many possible discretizations of (2.1a), for example, (2.7) -D_ Xi (a 1 (x + e 1 j&D^uW) - D^fa^x^ |)D^u(x)) where e. is the unit vector in the direction of the positive i coordinate i axis and (2.8) hD u(x) = ±(u(x±he.) - u(x)) ix . l l This discretization yields H, , = -a.(ih, (j- hh) , if j = 1, H, . = ; i,J i>J H i,j + 1= -a 2 (ih,(J + |)h) , if j = n, H.^. +1 = ; F. . 1 = -a ((i- ±)h,jh) , if i = 1, F . . 1 i+l,J = " a !((i + 2 )h ^ h) « if i = n > F i + l,d = ° ; E. . = -(H. . + F. . + F. _ . + H. . J . i>J 1*J i,J 1+1, J i, J+l III. DESCRIPTION OF TEST PROBLEMS Four distinct problems, identical to those described by Stone in [3], w iH be studied. Each problem will be solved using the square grid system described in Figure 1. The solution, u(x), of (1.2) and (l.^) was chosen to be (3.1) u. . = sin(ifrh) • sin( jTfh) for i, J = 1, ...,n 1, J where n is the grid size and h = -> — =-y. The data vectors, s, of (1.2) and (lA) were computed from (3-2) (Au E ). . = s. . . 1,0 1,0 The initial vector, u. ., was chosen equal to zero. 1,0 The first test problem is the model problem, defined by letting the coefficients, a (x) and a (x), be equal to -.5 over the entire region described in (2.3). For the second problem, the generalized model problem, the ratio of a (x) to a (x) is equal to 100, e.g. a (x) = -.5, a (x) = -.005. The third problem is defined by a heterogeneous region containing homogeneous subregions; see Figure 2. In region A, the ratio of a^(x) to a (x) is unity with a (x) = -.5. In region B, this ratio is ten with a- [ (x) = -.5; in C it is one-tenth with a (x) = -.05, and in D it is one with a-^x) = -.05. The last problem is identical to the third with one difference. In region A, the values of the a. (x) were randomly chosen from the set {z:0 < z < 1 subject to the condition that (0, 1) 30 (0,0) 25 _ 20 _ * 15 -a 10 _ 5 — Figure 2. Heterogeneous Test Problem Geometry (Region Defined on Unit Square) (3.3) If z < - 1 set a -( x ) = 10 ~ otherwise set a i (x) = z , i = 1,2 10 11 IV. THE STONE PROCEDURE Direct methods to solve the linear system, Au = s, involving the L-U decomposition of the coefficient matrix, A, are computationally inefficient when the grid size is large; because L and U will possess n non-zero diagonals. To make L-U decomposition efficient, A is replaced by a matrix of the form, A + B , which can then be factored as the product of sparse lower and upper triangular matrices, L and U , respectively. An approximate solution of Au = s is then computed by stopping the iteration, (4.1) (A+B )u . = (A+B )u - t(Au -s) , v cr n+1 an n ' provided, of course, that convergence is sufficiently rapid to be practical. In order to increase convergence rates it is plausible that the vector norm, C*.2) l|B a u||, should be minimized. Stone suggested constructing L and U , see Figures 3 and 4, to be sparse matrices with three non-zero diagonals, each of which coincide with the non-zero diagonals of A, i.e., (*.3) (L u) . . = b. .u. . _ + c. .u. _ . + d. .u. . a 1, j i, j l, j-1 i,j l-l, J i,j i,o (U u) . . = u. . + e. .u. , , + f . .u. . , a i, o i,0 i>J i+l, J i,0 i,0+l • Figure 5 and equations (4.3) show that the components of (L U u) are given by 12 Figure 3. Lower Triangular Matrix, L Figure h. Upper Triangular Matrix, U a a Figure 5. Product Matrix L^Uq, = A+^, 13 (L U )u . . = b. .u. , , + b . .e. . ,u. , . _ + c. .u. , . - ddf Ji,j 1,0 i,o-l i,j 1,0-1 1+1,0-1 i, i-l, n. \,\ + (d. .+b. .f. . _+c. ,e. , .)u. . + d. .e. .u. . . (4.4) i,J 1,0 1,0-1 1,0 1-1,0 1,0 i,J i,J 1+1,0 + c. .f. , ,vl. . ... + d. .f. ,u. . , 1,0 1-1,0 i-i, 0+1 1,0 1,0 1,0+1 where the subscripts i, j in (4.3) and (4.4) are the grid point indices, not row-column indices. (Throughout this section, lower case lettering will be used to represent components of the product matrices, L and U . Capitals are reserved for labeling the components of the coefficient matrix, A, and the auxiliary matrix, B .) Figure 9 displays the components of the iteration vector, u, that appear in (4.4) and the following equations are Taylor series expansions for these components in terms of the value of u at the grid point (i,0): (i) u , ,., = u. - hu „(i,o) + hu (i,o) + 0(h ) p (ii) u = u - hu (i,o) + 0(h ) i-i, x } J -*■ 2 (iii) u. .,. su. . + hu (i,o) + 0(h ) 1,0+1 i,0 y (4.5) (iv) u. . = u. . 1,0 1,0 2 (v) u = u + hu (i,j) - hu (i,o) + 0(h ) i + i,o-i i, o x y (vi) u = u + hu (i,o) + 0(h ) 2 (vii) u. . _ = u. . - hu (i,o) + 0(h ) 1,0-1 i,0 y (viii) u. . = u. . 1,0 i,0 The matrix, B , is then constructed such that: the elements of any row mul- tiply only those components of u which appear in (k.k) and (U.5); and the vector norm, (h .2) , is minimized. Figure 7 depicts an asymmetric matrix, B , Ik u whose components, (B u), are given by (B u). . = -qC. .u. . . + C. .u. _ . _ - aG. .u. , . a i,j 1,0 1,0-1 i, j i+i, j -l 1,0 1-1,0 (+.6) + a(C. .+G. .)u. . - aC. .u. , . + G. .u. _ . _ '1,0 i, i,J 1,0 1+1,0 1,0 1-1,0+1 - aG. .u. . _ . i,0 1,0+1 Note that when a = 1.0, the quantity, ±G. ., is a multiplier for the first 1 , J four equations in (^.5) and the quantity, ±C. ., is a multiplier for the last i, four equations in (4.5)- Consequently, the value of (^.2) is reduced to second 2 order, 0(h ), terms, when the u. . are given as values of a smooth function. i, The algorithm for determining the elements of L and U for this to to a a asymmetric factorization is derived by equating like terms of the matrices (L U ) and (A+B ), Figures 5 and 8, respectively. The resulting equations are b. , = H. , ,/(l+ae. . ) i,j i,j-l 1,0-1 c. . = F, - ,/(l+af. . .) i,J i-l, J 1-1,0 (U.7) d. , = E. . + a(C. . + G. .) - b. f. . - c. e i,j 1,0 1,0 1,0 1,0 1,0-1 i,0 i-i,, e. , = (F. . - aC. .)/d. . i,j 1,0 1,0 1,0 f . . = (H. , - aG. . )/d. . 1,0 i,0 1,0 i,0 15 ;'-QC» C...-OG' Cg£, -QC...G' -QG 1 Figure 6. Coefficient Matrix, A Figure 7« Auxiliary Matrix, B (Asymmetric) a .• H. . -QC» C...F -OG' E. .-fQC'+oG' F. .-QC...G 1 H. .-QG' 1_1 >J i*J i,J i,J Figure 8. Asymmetric Altered Matrix, A+B Q 16 3>l,j i»J-l ifl,d-l Figure 9. Points of Grid System Corresponding to the j-th Row of the Auxiliary Matrix, B a 17 where i,j = l,...,n and (U.8) C. , = b. ,e. G. = c. f , i,j i,j l-l, j A second method of minimizing (k.2) is obtained by constructing the matrix, B , as shovn in Figure 11. In this case, if q.= 1, (B u) is first order if the components, u. ., are given as values of a smooth and continuous i »0 function. This factorization has the improved property that B is symmetric. The following equations for determining the elements of L and U are obtained & ^ & a a by equating like terms of the matrices, (L U ) in Figure 5 and (A+B ) in Figure 12,: b. . = d. . _f. . . i,j i,J-l 1,0-1 c. . = F. , . - aG . . , 1,0 i-l,0 1,0-1 (^•9) d. . = E. . + 2aG. . . - b. .f. , _ - c. ,e, . . i,0 1,0 1,0-1 isO i,0-1 i,0 1-1,0 i. . = (F. . - aG. . . _)/d. i,0 1,0 1+1,0-1 i, f. . = (H. . - aG. .)/d. . 1,0 '1,0 i,0 1,0 where i, j = 1, . . .,n and G. . . = b. . .e 1,3-1 1-1,0 1-1,0-1 (+.10) G. , . _ = b. ,e. l+l, j -1 i,0 1,0-1 G. . = c. .f. . . 1,0 1,0 1-1,0 18 '-0G. G„...-OG. +2QG n -0G-. 12 1 12 ■ G -QG,, Figure 10. Coefficient Matrix, A Figure 11. Auxiliary Matrix, B (Symmetric) Q H. . -0G. G^...F. , .-QG n E. .+2QG n F .-QG_...G n H. .-QG n i,j-l 1 2 i-l, J 1 i,J 1 i,j 2 3 i,j 3 Figure 12. Symmetric Altered Matrix, A+B a 19 If u is not smooth, cancellation is not effective and Stone suggests using a variable a . If < a < 1, then whenever the matrix modification introduces the value u. , . ., into the individual grid point equation, it is to be partially cancelled by adding the value a(u. . - u. - u .). Stone suggests using a set of m parameters, a., i = 1, ...,m, < a. < 1 computed by letting: (4.11) a = m a max 1 - min 2Ax n a 2 (x) A x£ 1 + a 1 (x)Ax 2 , 2Ax a x (x) Ax 2 1 + a 2 (x) Ax x and o: = a . =0; the intermediate parameters a., i = 2, . ..,m - 1 are computed from (4.12) i-1 a. = 1 - (1 - a ) m_1 , i = 2, . . .,m - 1 i max ' ' ' which insures that < a. < a. . < a 1, . . . , m - 1 i "i+1 — "max , For the test problems in this thesis involving Stone's procedure with parameters, a set of nine parameters was computed and each was used on two successive iterations. For all odd iterations beginning with the first, the solution was obtained using the normal ordering of the grid points as explained in Chapter II. For all even iterations beginning with the second, the solution was obtained using the reverse ordering of the grid points, also explained in Chapter II. Altering the ordering for each iteration requires that the coefficients of L and U , given by equations (^.7) and (4.9), be recomputed each time. Since each parameter was used on two successive itera- tions, the parameter cycle length is 18. Before the parameters were applied to the iteration scheme, they were numbered in order of increasing magnitude, 20 i.e., cl < ql < . . . < a , and then rearranged in the following manner: °ty %> Oy Vq> ^ ( \> a 1 - Stone [3] remarked that if the values of the coefficients E, F, and H in the matrix ,A, see Figure 10, represent sharp, point to point discontinuities, and if the parameter a is varied during the iteration scheme, then any values of x other than unity do not generally improve the method. Since no numerical or analytical justification for this statement appears in [3], tests were performed using other values of x. These tests show an improvement and the results are discussed in more detail in Chapter VI. Numerical studies to determine the properties of the iteration parameters, a and x, were performed using the symmetric factorization .pro- cedure. For each test, a and x remained fixed for a maximum of 32 iterations. For each change of a or x, the initial vector, u. ., was chosen as zero. A systematic search for the optimum value of x, i.e., the value that provided optimal convergence for 32 iterations, was made for a = 1.0, in the range C 2. For the model problem on a 30 x 30 grid, divergence resulted after a few iterations whenever a was less than unity. Since no value of x could effect convergence, a was fixed at 1.0 for the remainder of the tests. Figures 13 and Ik show the results of tests performed on the model problem and the generalized model problem for various grid sizes ranging from 5x5 to 30 x 30. Similar results do not appear for problems 3 and k because of prohibitive work and time requirements. Only the optimum values for x are shown. Although only grid sizes which are multiples of 5 were used, a con- tinuous curve was drawn to indicate a trend. Figures 13 and Ik indicate that (i) the optimum value of x decreases as the grid size increases; 21 1.0 T A U 10 15 20 25 30 GRID SIZE Figure 13 • Model Problem - Tau vs. Grid Size 1.0, GRID SIZE Figure ik. Generalized Model Problem - Tau vs. Grid Size 22 (ii) the optimum value of t is less affected by the grid size when the ratio of a (x) to a (x) is large. This ratio was equal to 100 in Figure lk. For all problems the calculated optimum value of x was less than unity, regardless of grid size. Other tests were made to determine the effect of reordering the grid points on alternate iterations. Results were conclusive in showing that more effective cancellation and, therefore, faster convergence is achieved. More details are contained in Chapter VI. The actual programming of the procedure described in (U.l) was done as follows : Define and Solve R = t(s - Au ) n n 6 , , = u - u . n+ 1 n+ 1 n L V = R a n for V and U 6 4.1 = v a n+1 for 5 J _, . Then set n+1 n+1 n n+1 . 23 V. THE ADI PROCEDURE A well known, efficient iterative method for solving large systems of algebraic equations of the form Au = s on rectangular grid regions, such as described in (2.2), is the alternating direction implicit (ADI) procedure. The coefficient matrix, A, is partitioned into the sum of two, real and symmetric matrices, H and V, where H and V are discretizations of (5-D -^V 1 '^ and (5.2) -3x7 (a 2 U) ^ ) •> respectively. The solution, u, of (H+V)u = s is computed using the following 2-step iteration scheme: (H+ V l) Vl/2 = - (V -V )U n + S (V+oj. I)u .. = -(H-oj. I)u .- /0 + s jv n+1 jv n+1/2 where u^ is an initial estimate of the exact solution, u, u„ and gj . are jh jv iteration parameters, and j = n(modulus t), where t is the parameter cycle length. In (5»3), 2*t iteration parameters are used in a cyclic fashion and changed at each step of the iteration procedure. The problem we are concerned with is how to choose parameters to minimize the computer time required to obtain the unknown vector, u, to some preassigned accuracy. The solution to this problem, see Wachspress [6], 21* requires that the matrices, H and V, commute and that A "be positive definite. If these two conditions are satisfied, optimum parameters that minimize the spectral radius, u(T , ), of the t-step iteration matrix, * -1 -1 T x = n (V+w. I) (H-o). l)(H+a) T) (V-oo.,1) ; t , =1 jv jv jh jh (5.M (u t _ u E ) = T t (u Q -u E ) , can be easily and efficiently computed by a well-known algorithm [5], which cal- culates 2*t optimum parameters by minimizing the following continuous function: t X-w. Y-w., I n jv jh y = max n u — • u — t ' . .X+w... Y + w. J=l Jh jv (5.5) a mm where a, b, c, and d are eigenvalue bounds for H and V, respectively. Since every eigenvalue of T is of the form t X-u. Y-w.i (5.6) n — ^ & j=l Jh jv it follows that after t iterations, the spectral radius, y(T ), is such that u, -u (5.7) y(T ) = Hi" < u. V ii E. | - P t * V U 1 2 25 No scheme has yet been devised for choosing optimum parameters when H and V do not commute. A theorem in [2] proves that, if the above calculated parameters are used, despite lack of commutation of H and V, any desired error reduction can be obtained if the cycle length, t, is sufficiently long. It is our experience, see also Wachspress [6], that repeated application of a shorter cycle has proved better when H and V do not commute. Convergence rates improved, dramatically, when t was reduced from 32 to 9, while the number of iterations remained 32. Since the calculation of iteration parameters requires the values of the eigenvalue bounds for H and V, a method must be available to compute them. The values of the upper eigenvalue bounds, b and d of H and V, respec- tively, are less critical and can be obtained by Gerschgoren' s Theorem. Values for the lower eigenvalue bounds, a and c of H and V, respectively, are very critical for good parameter calculation and must be determined more carefully. Wachspress [7] suggests an algorithm based on Newton's method applied to the characteristic polynomials, |H-\I | and |V-\l|, which can be evaluated at any point, X, by a simple recursion [J] . If A is positive defi- nite, it follows that H and V are positive definite ,and a reasonable starting value for Newton's method is zero. 26 1 1 E 1 i u ■ -u L 1 1 n 1 '2 1 1 E i l u • -u^ L 1 1 ' '2 VI RESULTS AND CONCLUSIONS Figures 15-l8 display results corresponding to the four problems described in Chapter III that were obtained by solving each problem on a 30 x 30 grid system. In each of the figures, the ordinate is the L relative error defined by (6.1) where u is the exact solution, u^ is the initial estimate, and u is the ' ' n estimate. For the AD I procedure, two half steps are equivalent to one Stone iteration. The abscissa is the number of iterations needed to obtain the given L norm ratio. Only straight line approximations to the curves are shown in order to facilitate visual comparisons. The actual curves con- tain oscillations caused by parameter cycles. Each of the four Figures, 15-l8, contains six graphs labeled A-F. Each letter designates one of the following methods of solution: A. ADI procedure with 9 parameter pairs used cyclically. B. Stone's asymmetric procedure with a parameters used cyclically as described in Chapter IV, with alternating normal and re- verse ordering. The parameter, t, is fixed at unity. C. Same as B, except t is fixed initially at the value that pro- vides optimal convergence for 32 iterations. In all figures, except Figure 16, t is 1.6. In Figure 16, x equals 1.2. D. Same as B, but with Stone's symmetric procedure. E. Stone's symmetric procedure with normal ordering on every iter- ation. The parameter, a, is fixed and equal to unity. The parameter, x, is fixed and chosen as described in C above. 27 For Figures 15, 16, IT and 18, the optimum values for x are .365, .866, .3^1 and .326, respectively. F. Same as E, except with alternate normal and reverse ordering. The optimum values for x are also different. They are .331, .78U, .319 and .32U, respectively. By studying the results shown in Figures 15, 16, 17 and 18, we are able to make the following observations and conclusions: (i) For all test problems, procedure B converged faster than pro- cedure D. Therefore, the symmetric factorization is not as effective as the asymmetric factorization of Stone. (ii) For the generalized model problem, the ACT procedure provided faster convergence than any of Stone's methods. This con- tradicts the results given by Stone in [3] • The reason for this discrepancy is that Stone did not use optimum iteration parameters. Improvement in the performance of the ADI pro- cedure is not as dramatic for the remaining problems as for problem 2 . (iii) In every case procedure C out-performed procedure B, disputing Stone's statement that x = 1.0 gives better results. It must be noted, however, that the optimal values for x used in procedure C were found by a tedious systematic search. To fix x at unity is as reasonable as an inefficient search for a better value . (iv) The results of procedure F are, for all cases, an improvement over those of procedure E. Since the only difference between these two procedures is the ordering of the grid points, the suggestion by Stone [3] to use alternate normal and reverse ordering to increase convergence rates is numerically justified. 28 -. r- - J _ c.co 7.00 14.00 21. NO. OF ITERATION'S 28.00 Figure 15- Model Problem - Dirichlet Boundary Conditions 29 10 •Q r IT r-j 10 Figure 16. 7.00 14.00 ?1.00 28.00 3 r J.C NO. Or" ITERATIONS Generalized Model Problem - Dirichlet Boundary Conditions 30 10 u m o -.r,-'4 :L 1 — - r ' — 1 i u S; F — > r~ ' f - ' ' : L ^. < — i — K ,-R r~ 1 O B= n-'d r~ -F 10' 10 -1 ^F i 1 ! io-^t 10 iu t io- 1S G.CO 28.00 7.00 14. OG c'I.Ul 1 NO. OF" ITERATION'S Figure 17 . Heterogeneous Problem with Homogeneous Subregions Boundary Conditions 3 r J.O0 Dirichlet 3i r j i L ES 2d 10" 10 10' 10' 10' 10' J -15 c.co 7.00 14.00 21.00 NO. OF ITERATIONS 28.00 ; r J.oo Figure 18. Heterogeneous Problem with Random and Homogeneous Subregions - Dirichlet Boundary Conditions 32 (v) Stone [3], in an attempt to show that one of the sequence of values of a should be less than unity, argued for large grid sizes, that if both components, F and H, of the coefficient 7rAx IT ^ V matrix, A, are large relative to — ^ — and — — — , respec- tively, then repeated application of a = 1.0 in the iteration scheme, (k.l), leads to divergence. Stone, however, fixes x to be unity. For a = 1.0, values of x other than unity can be selected to yield convergence, as procedure F in Figures 15, 17 > and 18 show, (vi) The results of procedure F, when compared with those of pro- cedure D, show that the iteration scheme defined by (^.l) converges faster when a sequence of values of a is used. The results from procedure F, for problems 1, 3> and K, are optimal using fixed a and x in (^+.1) because experiments indicate that the optimal values of x and a occur when a = 1.0. Values of a less than 1.0 cause divergence, independently of the value of x used for the run. Therefore, finding the optimal values of (x, a) reduces to finding the optimal value of (x, 1.0). Although the findings in observation (i), showing the symmetric procedure, D, to be less effective than the asymmetric procedure, B, are somewhat discour- aging, they do not conclusively show that further consideration of procedure D is useless. It is more important that the results of procedure D be similar to those of B as the following observations point out: (vii) For the complex problems, 3 and k- in Figures 17 and l8, proce- dure D out-performed the ADI procedure. The improvement be- comes more pronounced as the complexity of the problem increases, see Figure l8. Hence, Stone's symmetric procedure merits further analytical attention. 33 (viii) For the simpler problems, 1 and 2, in which there were no sharp, point to point discontinuities in the values of the coefficients, F and H, in the matrix, A, the symmetric procedure, D, was nearly as effective as the asymmetric routine, B, see Figures 15 and l6. Thus, for simple problems, any findings resulting from analytical studies of Stone's symmetric procedure may be applied to his asymmetric procedure. 3k LIST OF REFERENCES [l] T. Dupont, 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. D. W. Peaceman, and H. H. Rachford, Jr., "The Numerical Solution of Parabolic and Elliptic Differential Equations", SIAM Journal On Applied Mathematics , Vol. 3, 1955. [3] H. L. Stone, "Iterative Solution of Implicit Approximations of Multi- dimensional Partial Differential Equations", SIAM Journal On Numerical Analysis , Vol. 5, 1968. [k] , Private Communication, April, 19&9 • [5] E. L. Wachspress, "Extended Application of Alternating Direction Implicit Model Problem Theory", SIAM Journal On Applied Mathematics , Vol. 11, 1963. [6] , Iterative Solution of Elliptic Systems , Prentice-Hall, Inc., Englevood Cliffs, New Jersey, 1966. [7] , "Numerical Solution of Neutron Diffusion Problems", SIAM-AMS Proceedings On Numerical Solution of Field Problems in Continuum Physics , Vol. 2, 1970. NO^ 1 \#l