LIBRARY OF THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAICN 510.84- IJKor no. 403-4-08 cop. 2 The person charging this material is re- sponsible for its return to the library from which it was withdrawn on or before the Latest Date stamped below. Theft, mutilation, and underlining of books are reasons for disciplinary action and may result in dismissal from the University. To renew call Telephone Center, 333-8400 UNIVERSITY OF ILLINOIS LIBRARY AT URBANA-CHAMPAIGN one* 07 ^0|i| OK 3 RE I'D L161— O-1096 Digitized by the Internet Archive in 2013 http://archive.org/details/numericalstudies407beck Report No C &0 ' and let (9) x H = iA Xj x 2 . = jA x 1 < i < N Ax = l/N 1 < j < N Ax = l/N. Define the grid system D over the region D by h (10) D = {(ih, jh)}. Denote the grid point (i h, j h) by (i, j). The boundary dD of D, is defined by (ID 8D h = (i, 0) < i <_N+1, (1> j) < j < N+l, (i, 1) <_j < N+l, (0, j) < j < N+l. The grid system is shown in Figure 1. Order the grid points of D in a left-to-right, down-to-up fashion, and let u be a vector whose components u. . are ordered in the same sequence as the grid points 1 » J (0,1) (0,Nh) (0,h) (0,0) i (i,D r~ 1 i 1 i 1 i 1 1 ! S _ 1 I 1 1 1 1 1 1 I I 1 1 1 1 (h,0) (Nh.O) (1,0) # i Figure 1. Grid Point System (i,j). Thus, u. . precedes u., . if i < i', and u. . precedes u if 1, j i , J i» J ^» -t/ j < L Several discretizations of the expression are possible. For example, (13) h" 2 {[a^i.j) aj (i+1, j)] 1/2 [u(i, j) - u(i+l, j)] + [a^i, j) a x (i-1, j)] l/2 [u(i, j) - u(i-l, j)] + [a 2 (i, j) a 2 (i, j + l)] l/2 [u(i,j) - u(i,j+l)] + [a ,(i,j)a, (i,j-l)] l/2 [u(i,j) - u(i,j-l)]}, or h -2 (14) — [a^i.j) + a^i+l.j)] [u(i.j) - u(i+l,j)] + [a^i, j) + a^i-1, j)] [u(i,j) - u(i-l.j)] + [a 2 (i,j) + a 2 (i,j+l)] [u(i,j) - u(i, j+1)] + [a 2 (i,j) + a 2 (i,j-l)] [u(i,j) - u(i,j-l)]. No attempt has been made to consider the mathematical merits of one discretization over another. The criteria used result from storage problems, not mathematical considerations. The need 8 to minimize the storage requirements of the routines used in this thesis for testing the ADI and Stone algorithms force the computa- tion of the elements of A, at almost every step of execution, from the coefficients, a (x) and a (x). Computing the elements of A 2 this way allows the storage of only 2N + 4N values of the coeffi- 2 cients, a^x) and a (x), instead of the (at most) 5N - 2N non-zero elements of A. For these reasons, the discretization, (14), is the most practical for converting the boundary value problem to matrix form. The exact definition of these elements of A as obtained from (14) depends on the type of boundary value problem, Dirichlet or Neu- mann. Either problem may be written in the form Au = S, where A 2 2 2 2 is an N x N matrix if Dirichlet, and (N+2) x (N+2) if Neumann. The difference in the order of A, of course, is due to the fact that for the Dirichlet problem, components of u are given on the boundary, whereas for the Neumann problem, they are not. If (i, j) is a point of D for which (i+_l, j) or (i, j+_l) is not on the boundary, then the discretization (14) yields the following values for the elements of A for either the Neumann problem or the Dirichlet problem: B. = -a_(ih, jh) - a (ih, (j-l)h), i.J 2 2 D. . = -a.(ih, jh) - a ( (i-l)h, jh), i.J 1 l (15) E. . = -(B. . + D. . + F. . + H. .), F. = -a. (in, jh) - a ( (i+l)h, jh), i.J 1 ! H. . = -a_(ih, jh) - a_(ih, (j + l)h), i.J 2 2 and S. . = f(i, j). i» J The treatment of what happens at the boundary is straight- forward. Let (i, j) be a grid point for which an adjacent point is on the boundary of D, , that is, i = 1 or N or j = 1 or N. Suppose i = 1, h 1 < j < N. The other cases are the same. For the Dirichlet problem, (16) D . = i.J and (17) S. . = f(i, j) + [a (i, j) + a (i-l),j)] u(i-l, j). The other elements are defined as in (15). This completes the de- scription of the Dirichlet problem. For the Neumann problem the elements of A and the inhomogeneous part, S , are defined exactly i. J as in (15). Now, suppose i = 0, 1 < j <_N. Only the discretization cor- responding to the Neumann problem is defined. We have, (18) D. . = and (19) f(i,j) = -2 [a^i.j) + a^i+l.j)]. The other elements are defined as in (15). 10 11 3. STONE ALGORITHM 3. 1 The Factorization The factorization of A into the product of a lower triangular matrix, L, by an upper triangular matrix, U, such that A = LU, does not yield a sparse fadtorization. Stone [57] suggests constructing L and U to be sparse matrices with three non-zero diagonals, each of which coincide with appropriate diagonals of A. The product of L and U now defines an altered form of A, (20) A + B = LU, Of which makes the direct solution of (A + B ) u = s computationally easy. Let the elements of L and U be defined by (21) (Lu). . = b. . u. . + c. . u. , . + d. . u. ., i,j i,j i,j-l i,j i+l, j i, j i,j (Uu). . = u. . + e. . u. , . + f. . u. . ,. 1 >i !>J i.J i+l, J i, J i.j + 1 The product LU is defined by (22) [(LU)u]. . = b. . u. . , + b. . e . u. , i,j i,j i, j-1 l, j i,j i+l.j-l + c u . + (d. . + b. . f. . + c. . e. . .) u. . X »J 1_1 'J i,J i*J i» J-1 i,J i-l, J 1,3 + d. . e. . u. . + c. . f. , . u; , 1,3 1,3 1+1,3 i,J i-l, J i-l,3+! + d. . f. . u. . ,. 1,3 1,3 1,3 + 1 Stone computes the elements of L and U by the following al- gorithm: 12 b. . = D. . - ry G. ., i.J LJ i.J c. . = D. . - a C. ., l »J LJ i.J (23) d. . = E. .+ a (C. . + G. .) - b. . f. . . - c. . e. , ., LJ LJ i.J i.J i.J i.J-1 i.J i-l. J » J F. i .,)' •Of c. 1, j e . i, d. -.j . f. if J H. h j" a G. if j d. 9 I.J where < q, < 1. 3. 2 The Iterative Method The iteration defined in (4), (24) (A + N)u n+1 = (A + N)u n - T (Au n -s), is computed as follows (Stone [5]) for j = 1. Define (25) R n = s - Au 11 and / _ / . p n+l n+1 n (26) 6 = n - u . Solve (27) LV = R n for V and 13 (28) U 6 n+1 = V n+1 . , for 5 . Then set ._., n+1 n n+1 (29) u = u + 6 3. 3 Application Notes Stone [5] suggests a variation for alternate steps of the itera- tion which has been followed in this thesis. For all odd iterations be- ginning with the first, the solution is obtained at the (i, j) grid points by starting with the j = row and i taking on the values < i <_ N. Then j is incremented by 1 and i again assumes the above range of values, until j = N. For all even iterations beginning with the second, this process is reversed, that is, the solution is obtained at the points in the j = N row with i assuming the values < i <_ N. Then j is de- cremented by 1 and i again assumes the range of values < i < N. At the even iteration steps then, the computation of the solu- tion begins at the upper left-hand grid point, (1, Nh), and moves in a left-to-right, up-to-down order over the grid instead of the left-to-right, down-to-up order of odd iterations. This new ordering of the grid points necessitates recomputing the coefficients of the altered matrix, A + N, for the even steps. A detailed description of the Stone program is provided in Ap- pendix A. 14 4. API ALGORITHM Let H and V be discretizations of (30, -&:( a i (x) £:) and < 31 >-i^( a 2<*>iij)- respectively [8], so that (32) A =H + V. The solution u of Au = S is computed as the limit of the sequence de fined by (33) (H + co , I)u n+1 / 2 = S - (V - co I)u n , n, h n, v (V + co I)u n+1 = S - (H - co u I)u n+1/2 . n, v n, h Since H and V are real, tridiagonal, symmetric matrices, both (H + co I) and (V + co I) are positive definite and their inverses n, h n, v exist. The ADI algorithm is based on (33). A detailed description of the ADI program is provided in Appendix B. 15 5. CODING OBJECTIVES 5. 1 Storage Minimization Let A be the matrix formed by discretizing the differential operator (34) -&_/a ( X ) rM+ -*-fa,(x) -M. b x : \ i d Xl ; ax 2 \ z ax 2 ; Let (i, j) be a grid point, and let PS = -[a 2 (i, j-1) + a 2 (i,j)], PE = -[a (i-l,j) + a x (i, j)], (35) PW = -[a x (i,j) +a 2 (i+l,j)], PN = -[a 2 (i,j) + a 2 (i,j+l)], PO = -(PS + PE + PW + PN). If the discretization in (14) is followed, the row of A corresponding to the grid point, (i, j), 1 < i,j < N, is (36) PS . . . PE PO PW . . . PN, 2 2 where, if the order of A is N x N , PS and PW are N positions re- moved from PO on the main diagonal. All possible values of PS, PE, PO, PW, and PN form five diagonals. The maximum number of non-zero elements of these ,. . , 2 2 2 2 2 diagonals are respectively, N - N, N , N , N and N - N, totaling 16 5N - 2N non-zero elements of A. Actually, values of PE corres- ponding to i = 1 are zero. Accordingly, the number of non-zero ele- 2 ments of the diagonal formed by the possible values of PE is N - N. Since A is symmetric, the number of non-zero elements in the diagonal 2 2 formed by PW is also N - N. Therefore, the bound 5N - 2N may be replaced by 5N - 4N. As a result, storage of A is possible with only 5N - 4N loca- tions, although the treatment of the boundary points required to obtain this bound is a somewhat hazardous programming feat. By using the fact that the diagonals of A are computed from the coefficients, a (x) and a (x), it is possible to reduce the storage requirements even further, although at the expense of program ef- ficiency, by storage of the necessary values of a (x) and a (x). The necessary values of a (x) are values a (i, j) for which <_i < N + 1, 1 <_ j <_ N, and the necessary values of a (x) are all values for which 1 < i 0>y Qy Of g > a g > OL^, a 5 > Oi 2 > "/, ®7' °V a 4 - & 4 , ofj, cvj. The minimum parameter was always zero and the maxi- mum parameter was computed by (40) a = 1-min max 23 2 A x 1 2 A x. 1 + a (x) A x 2* L 1 + a (x) A x. a^x) A x ' *2 (XMX 1 The intermediate parameters, < I), n v, n h, n h, n v, n Accordingly, n <«> ll» 2n h <_ II TT t || 2 ||» n 2 . i=0 Assume that V and H commute and that they have the same eigenvectors 24 Let x be such an eigenvector with eigenvalues \ and ^ corresponding to H and V respectively. Then x is an eigenvector of n 77 t. 1=0 with corresponding eigenvalue n (45) TT a - w X - w , • ** v, 1 h, 1 1=0 ** v, 1 h, 1 Conversely, every eigenvalue of n IT T i 1=0 is of the form (45) where ^ and \ are eigenvalues of H and V. Ideally, the parameters ax . and go . should be chosen 7 r h, 1 v, 1 to minimize the absolute value of the expression in (45) when ^ an d X range through all the eigenvalues of H and V respectively. Such a choice, of course, minimizes n TT T ; 1=0 since this quantity is equal to the largest eigenvalue. Parameters, how- ever, are selected somewhat differently. Let a and b, c and d, be bounds on the eigenvalues of H and V, i. e. , < a < \ < b and blem. For the testing performed in this thesis, parameters a) . and co . were chosen to minimize the rational expression v, 1 25 max (46) *—t < o X Boundary Conditions ■ 2 • • 7. COMPUTATIONAL WORK (NO. OF ITERATIONS FOR STONE PROCEDURE) ogeneous Problem with Random Variations in Subregions 'hlit Boundary Conditions 10 . § 10 w w > I— I H $ w en I— I H w . , i 10 20 30 40 COMPUTATIONAL WORK (NO. OF ITERATIONS FOR STONE PROCEDURE) Figure 9. Generalized Model Problem - Dirichlet Boundary Conditions o « w w > l-H < w (NJ 32 Figure 10. COMPUTATIONAL WORK (NO. OF ITERATIONS FOR STONE PROCEDURE) Heterogeneous Problem with Homogeneous Subregions - Dirichlet Boundary Conditions O w w > H W f\i 1 1 . COMPUTATIONAL WORK (NO. OF ITERATIONS FOR STONE PROCEDURE) Heterogeneous Problem with Random Variations in Subregions " hlel Boundary Conditions < w W P w N »— i o z X 33 COMPUTATIONAL WORK (NO. OF ITERATIONS FOR STONE PROCEDURE) Figure 12. Model Problem - Neumann Boundary Conditions < 9 W Q W N i— i -1 < « o z ED >< <: COMPUTATIONAL WORK (NO. OF ITERATIONS FOR STONE PROCEDURE) Figure 13. Generalized Model Problem - Neumann Boundary Conditions < a w W Q W N i—i s OS o z X 2 10 10 -2 10" 3 L ID" 4 !. 10" 5 L 10 10 Stone ADI 34 10 20 30 40 Neu- C COMPUTATIONAL WORK (NO. OF ITERATIONS FOR STONE PROCEDURE) Figure 14. Heterogeneous Problem with Homogeneous Subregions mann Boundary Conditions < 9 W Q W N »— i < o z S X 2 COMPUTATIONAL WORK (NO. OF ITERATIONS FOR STONE PROCEDURE) ire 15. Heterogeneous Problem with Random Variations in Subregions Neumann Boundary Conditions on o on on w w > i— i H w on (M 35 20 30 40 COMPUTATIONAL WORK (NO. OF ITERATIONS FOR STONE PROCEDURE) Figure 16. Model Problem - Neumann Boundary Conditions O on on w w > 3 W on (NJ . 5 COMPUTATIONAL WORK (NO. OF ITERATIONS FOR STONE PROCEDURE) Figure 17. Generalized Model Problem - Neumann Boundary Conditions Pi o Pi w w > 5 w (M Figure 1 36 10 20 30 40 COMPUTATIONAL WORK (NO. OF ITERATIONS FOR STONE PROCEDURE) Heterogeneous Problem with Homogeneous Subregions maim Boundary Conditions Neu- O w w > H OS ^ i . 9 COMPUTATIONAL WORK (NO. OF ITERATIONS FOR STONE PROCEDURE) 19. Hf-'r-rogeneous Problem with Random Variations in Subregions imann Boundary Conditions 37 6. 3. 1 Dirichlet Problem Figures 4 and 8 present data for the model problem with Dirichlet boundary conditions. Note that ADI converges at a faster rate than the Stone method. Data for the generalized model problem with Dirichlet boundary conditions are presented in Figures 5 and 9. ADI maintains approximately the same rate of convergence as in the model problem; however, there is a significant increase in the rate of the Stone method. Figures 6 and 10 exhibit the results for the heterogeneous problem with homogeneous subregions with Dirichlet boundary conditions, The complexity of this and the next test case is comparable to the com- plexity of practical problems. The last case for the Dirichlet problem is the heterogen- eous region with random variations in subregions, and the data for this problem are presented in Figures 7 and 11. As in the previous case, the Stone method converges faster than ADI. Also, the ADI rate of con- vergence decreased as the test cases increased in complexity, whereas the rate of convergence of the Stone method has either increased (gen- eralized model problem) or remained nearly the same. Table 1 indicates the work required by either method to achieve a maximum normalized residual of 10"° for the four test cases. Where required, straight line extrapolation was used to obtain the data. Note that except for the generalized model problem, the Stone method was less sensitive to the nature of the problem itself. 38 method Stone ADI TABLE 1 Dirichlet Problem Work Requirements of Iterative Methods model 31 13 RATIO: ADI/STONE 0. 42 generalized model 1. 60 heterogeneous model homogen. random no. subregions subregions 27 34 1. 26 28 85 3. 12 6. 3. 2 Neumann Problem Figures 12 and 16 present the data for the model problem with Neumann boundary conditions. Again ADI converges at a faster rate than Stone. Note, however, that for the Neumann conditions, both rates have decreased slightly. For the generalized model problem with Neumann bound- ary conditions, Figures 13 and 17 show that the Stone algorithm is con- verging faster than ADI. For either method, convergence is slower than for Dirichlet boundary conditions. Figures 14 and 18 present the data for the heterogeneous problem with homogeneous subregions. In this case, with Neumann boundary conditions, the faster rate of convergence is again obtained by the Stone algorithm. I "inally, the heterogeneous model with random subregions esented in Figures 15 and 19 also shows the superiority of the Stone jorithm. Convergence rates have consistently decreased for both 39 algorithms with the increasing complexity of the test cases. Table 2 is a summary of the work required by both -5 algorithms to reach a maximum normalized residual of 10 for the test cases with Neumann boundary conditions. This table is patterned after a table of Stone [5], and straight line extrapolation was used where necessary. It is apparent that the Stone algorithm is less sensitive to the nature of the test case. TABLE 2 Neumann Problem Work Requirements of Iterative Methods generalized method model model Stone 57 33 AD I 27 46 RATIO. ADI/STONE 0.48 1.41 heterogeneous model homogen. random no. subregions subregions 41 38 45 66 1. 10 1. 74 6. 4 Summary (i) The Stone algorithm was more effective than ADI as the com- plexity of the problems being solved increased, although the parameters used for ADI were optimal and no method is known for determining optimal Stone parameters. ADI rates of convergence tended to decrease with increasing problem complexity while Stone rates were relatively stable. 40 (ii) Convergence rates for both methods were faster for problems with Dirichlet boundary conditions than for Neumann boundary- conditions . 41 7. CONCLUSION In view of the results obtained by this author, the Stone method is in fact superior to ADI methods. Most importantly, the re- sults obtained by Stone [5] have been verified. It is imperative to note, however, that since the problems considered in this paper were not ex- actly the same as Stone's, the numerical data did fluctuate somewhat. Nevertheless, the results were still verified; furthermore, the compari- sons for both Dirichlet and Neumann boundary conditions further solidifed them. 42 LIST OF REFERENCES [1] D. W. Peaceman and H. H. Rachford, Jr. , "The Numerical Solu- tion of Parabolic and Elliptic Differential Equations, " SIAM Journal on Applied Mathematics , " Vol. 3, (1955). [2] J. Douglas, Jr., and H. H. Rachford, Jr., "On the Numerical Solu- tion of the Heat Conduction Problems in Two or Three Space Variables, " Transactions of the American Mathematical Society , Vol. 82, (1956). [3] R. S. Varga, Matrix Iterative Analysis , Prentice -Hall, Inc. y Engle- wood Cliffs, New Jersey, (1962). [4] E. L. Wachspress, Iterative Solution of Elliptic Systems , Prentice -Hall, Inc. , Englewood Cliffs, New Jersey, (1966). [5] H. L. Stone, "Iterative Solution of Implicit Approximations of Multi- dimensional Partial Differential Equations, " SIAM Journal on Numerical Analysis , Vol. 5, (1968). 1 [6] T. Dupont, R. P. Kendall, and H. H. Rachford, Jr., "An Approximate Factorization Procedure for the Solution of Elliptic Difference Equa- tions, " SJAM_JJDuj-nalj)nJ^ Vol. 5, (1968). [7] T. Dupont, "A Factorization Procedure for the Solution of Elliptic Difference Equations, " SIAM Journal on Numerical Analysis , Vol. 5, (1968). [8] J. Douglas, Jr. , "A Survey of Numerical Methods for Parabolic Differ- ential Equations, u _A^vajT£e_s_Jn_Compjiters , Vol. 2, (1961). [9] C. W. Gear, "High Speed Compilation of Efficient Object Code, " Communications of the ACM , Vol. 8, (1965). t 10 ] Applied Numerical Methods, John Wiley and Sons, Inc. , New York, (1969). 43 APPENDIX A STONE PROGRAM (DIRICHLET) 44 C STONE PROCEDURE C FOR ITERATIVE SOLUTION OF IMPLICIT APPROXIMATIONS OF C TWO DIMENSIONAL PARTIAL DIFFERENTIAL EQUATIONS C C FOR THE LINEAR SYSTEM OF EOUATIONSt A*T»Q, THIS PROGRAM C COMPUTES THE MATRIX A+N WHICH IS FACTORED INTO THE PRODUCT C OF TWO SPARSE MATRICES, L AND U, SUCH THAT A+N=L*U. THE C SOLUTION OF (A+N)*T*Q IS COMPUTED AS THE LIMIT OF THE SEQUENCE C DEFINED BY C ( A+N ) *T ( N+ 1 ) « ( A+N ) *T ( H )-TAU* ( A*T < N )-Q ) , C WHERE TAU IS A PARAMETER TO ENHANCE CONVERGENCE. C IMPLICIT REAL*8(A-H,0-Z) DIMENSION AK960),A2(960),Q(900),RN(900),T<900),BL(870),CL(900>, XDL(900),EU(900),FU(900) ,PARAM(9) C INPUT M,A1,A2,Q,T,N,ITMAX C C NOTE: IT IS NECESSARY FOR THE USER TO INPUT THE FOLLOWING C VARIABLES INTO THIS PROGRAM: C M = ONE HALF THE NUMBER OF PARAMETERS PER CYCLE (THERE C ARE 2*M PARAMETERS PER CYCLE CALCULATED IN SUBROUTINE C PARAMS). C A1,A2 = THE VECTORS OF DIFFERENTIAL EQUATION COEFFICIENTS. C = THE VECTOR CONTAINING THE RIGHT HAND SIDES OF THE C LINEAR SYSTEM. C T = THE INITIAL VECTOR. C N = NUMBER OF GRIDPOINTS ON ONE SIDE OF THE SOLUTION C GRID. (THERE ARE N**2 GRIDPOINTS IN THE TOTAL SQUARE C GRID SYSTEM). C ITMAX = TERMINATION CRITERIA OR THE NUMBER OF ITERATIONS C TO BE PERFORMED. C C THE INPUT THAT FOLLOWS IS A SAMPLE OF THE NECESSARY DATA. ITMAX=32 N = 30 M = 9 DO 109 1=1,960 Al ( I )=-.5D00 109 A2( I )=-.5D00 DO 110 1=1,900 T( I )=0.0D00 110 0(1 )=0.0D00 C CALCULATE THE PARAMETERS FOR A CYCLE OF 2*M ITERATIONS FOR C CONSTANT OR VARIABLE Al AND A2. CALL PARAMS(N,A1,A2,M,RN,NSQP2N) C it IS UP TO THE USER TO ARRANGE THE PARAMETERS AS HE WANTS THEM C USED IN THE CYCLE IN THIS SECTION. THE PARAMETERS HERE ARE C ARRANGED IN THE FOLLOWING MANNER. IF THE PARAMETERS ARE r - RN( 1 ) ITER, BIG ITER=ITER+1 —SOLVE L*V=RN FOR V. CALL SOLVEL(N,BL,CL,DL,RN,NSQ,NSQMN) —SOLVE U*DELTA=V FOR DELTA. CALL SOLVEU(N,EU,FU,RN,NSQ,NSQMN) -COMPUTE T(N+1)=T(N)+DELTA. CALL NEWT(T,RN,NSQ> -TEST FOR END OF ITERATION SCHEME IF( ITER.GE.ITMAX) GO TO 20 -ODD ITERATION STEPS ARE PERFORMED HERE. THE GRIDPOINTS ARE —CONSIDERED IN A NEW ORDERING FOR THESE STEPS(SEE PAPER). -COMPUTE THE COEFFICIENTS OF THE ALTERED MATRIX, A+N, FROM A. CALL BALTRM(N,A1,A2,BL,CL,DL,EU,FU,ALPHA,NSQ,NSQP2N,NSQMN) —COMPUTE THE PRODUCT A*T. CALL PR0DMT(N,A1,A2,T,RN,NS0,NSQP2N) -COMPUTE THE RESIDUAL 0-A*T. CALL RESIDL(RN,Q,NSQ) -COMPUTE THE MAXIMUM NORMALIZED RESIDUAL. CALL MAXNR(NSO,RN,SUM,BIG) 46 C OUTPUT THE MAXIMUM NORMALIZED RESIDUAL AFTER ITER ITERATIONS. WRITE(6,905)ITER,BIG ITER«ITER+1 C SOLVE L*V*RN FOR V. CALL BSOLVL(N,BL»CLtDLtRN,NSQ,NSQMN) C SOLVE U*DELTA=V FOR DELTA, CALL BSOLVU(N,EUtFU,RNtNSOfNSOMN) C COMPUTE T(N+1)»T(N)+DELTA. CALL NEWT(T,RN,NSQ) C TEST FOR END OF ITERATION SCHEME IF( ITER.LT.ITMAX)GO TO 10 C COMPUTE THE PRODUCT A*T. 20 CALL PRODMT(N»Al t A2,T,RN r NSQ,NSQP2N) C COMPUTE THE RESIDUAL Q-A*T. CALL RESIDL(RN»Q,NSQ) C COMPUTE THE MAXIMUM NORMALIZED RESIDUAL FOR THE LAST ITERATION, CALL MAXNR(NSOtRN,SUM,BIG) WRITE(6,905)ITER,BIG C OUTPUT THE FINAL SOLUTION VECTOR T(I). WRITE(6,901) WRITE(6,902) ( I , T < I ) , I = 1 ,NSQ ) C FORMAT STATEMENTS 901 FORMAT! 'OTHE SOLUTION VECTOR T(I) IS:') 902 FORMAT (• ',/,4(' T(',I3,») = ', 016,8)) 903 FORMAT( '0'tl3X, 'MAXIMUM NORMALIZED RESIDUAL') 905 FORMAT(' I TERAT ION ' , I 3 1 6X, D16.8 ) STOP END SUBROUTINE PARAMS ( N, Al , A2t M,P ARM,NSQP2N ) IMPLICIT REAL*8( A-H,0-Z) DIMENSION Al(NS0P2N)tA2(NS0P2N) »PARM(M) C C CALCULATE THE PARAMETERS FOR THE USE IN ONE CYCLE OF 2*M C ITERATIONS FOR VARIABLE Al AND A2 BY TAKING THE ARITHMETIC C AVERAGE OF THE MAXIMUM PARAMETER FOR ALL THE GRID POINTS, C PTMAX=0.0000 NSQ=N*N DO 1 I=l t NSO IF( Al ( 1+1 > ,EQ.0.0D00.OR.A2( I +N ) . EO .O.ODOO ) GO TO 4 TEMP1=A1 ( I+l )/A2( I+N) TEMP2=A2 ( I+N)/A1( I+l) GO TO 3 A TEMP1=0.0D00 TEMP2=0.0D00 3 AMAX = 1.0D00-DMIN1( (2.0D00*( 1.0D00/(N+1 ) ) ) / ( 1 .0000+ ( TEMP2 ) ) « (2, X0D00*( 1.0000/(N+1 ) ) ) /( 1.0D00+( TEMPI ) ) ) PTMAX=PTMAX+AMAX I CONTINUE AMAX=PTMAX/DFLOAT(NSO) PARMI 1 )=0.0000 H1=M-1 47 DO 2 1=2, M 6X=DFLOAT< I-l)/DFLOAT(Ml) PARMU )*1.0DOO-( (1.0DOO-AMAX)**EX) 2 CONTINUE M2=2*M WRITE(6,100)M2,(PARM( I ),I«1,M) 100 FORMAT('0THE UNORDERED PARAMETERS FOR • XARE: • ,/,(D16.8) ) RETURN ENO SUBROUTINE ALT ERM( N ,A1 ,A2 ,BL ,CL,DL,EU, FU, ALPHA, NSQtNSQP2N ,NSQMN ) IMPLICIT REAL*8(A-H,0-Z) DIMENSION AKNSQP2N) , A2 ( NSQP2N ) ,BL ( NSOMN ) ,CL ( NSOMN ) ,DL < NSO ) ,EU ( NSQ XMN) ,FU(NSOMN) -CALCULATE THE ALTERED MATRIX M+N=L*U FROM THE GIVEN MATRIX M. L IS -A TRI-DIAGONAL MATRIX CONSISTING OF THE MAIN DIAGONAL, THE LEFT -ADJACENT DIAGONAL, AND A LOWER DIAGONAL AT A DISTANCE N FROM THE -MAIN ONE. THE ELEMENTS OF THESE DIAGONALS ARE STORED IN THE -VECTORS DL, CL, AND BL RESPECTIVELY. U IS A TRI-DIAGONAL MATRIX -CONSISTING OF THE MAIM DIAGONAL EVERY ELEMENT OF WHICH IS EQUAL -TO THE CONSTANT ONE, THE RIGHT ADJACENT DIAGONAL, AND AN UPPER -DIAGONAL AT A DISTANCE N FROM THE MAIN ONE. THE ELEMENTS OF THE -THE TWO OFF DIAGONALS ARE STORED IN THE VECTORS EU AND FU -RESPECTIVELY. VANISHING ELEMENTS OF THE DIAGONALS ARE NOT STORED -IN THESE AFORE MENTIONED VECTORS. COMPUTATION IS DONE IN A LEFT-TO-RIGHT AND DOWN-TO-UP ORDER -SINCE THE ELEMENTS OF L*U CALCULATED AT ANY GRID POINT ONLY -DEPEND ON PREVIOUSLY CALCULATED ELEMENTS AT PREVIOUS GRID POINTS. PS,PW,PO,PE, AND PN DEFINE THE LOWER DIAGONAL, LEFT ADJACENT -DIAGONAL, MAIN DIAGONAL, RIGHT ADJACENT DIAGONAL, AND UPPER -DIAGONAL OF THE MATRIX M. -DEFINE CONSTANTS N1=N+1 NN1=N1+N NN=N+N N21=NS0-N-N+1 NSQM1=NSQ-1 -GET M+N=L*U FROM M PE = A1 (2I+AK3) PN=A2(N1 )+A2(NNl) PW=A1 (1)+A1<2) PS=A2(1)+A2(N1) PO=-(PE+PN+PW+PS) -THE COMPUTATION IS PERFORMED AT THE FIRST POINT, THEN THE -INTERMEDIATE POINTS, AND THEN THE LAST POINT OF THE GIVEN -HORIZONTAL LINE IN ORDER TO TAKE ADVANTAGE OF THE VANISHING OF -COMPONENTS OF M ALONG THE BOTTOM AND TOP HORIZONTAL LINES, AND -THE LEFT AND RIGHT VERTICAL LINES OF THE GRID. THE ELEMENTS OF L*U -FOR THE FIRST GRID POINT ARE COMPUTED HERE. 0L( 1 ) = P0 48 EU(1)=PE/DL(1) FU(1)=PN/DL(1) DO 15 1=2, N C 1 FOLLOWS THE REMAINING GRID POINTS ON THE BOTTOM HORIZONTAL C LINE. THE ELEMENTS OF L*U FOR THOSE POINTS ARE COMPUTED HERE. I IS C THE BASIS FOR THE INDICES OF Al, A2, AND THE ELEMENTS OF L*U. PE=AHI + l>+AlU + 2) PN=A2< I*N)+A2< I+NN) PWsAl(I)+AHI + l) PS=A2(I )+A2U+N) PO=-(PE+PN+PW+PS) CLU-1 )«PW/Q.ODOO+ALPHA*FU(I-l>) DL(I )=PO*CL(I-l)*(ALPHA*FU(I-l)-EU(I-l>) FUU )=(PN-ALPHA*CL( I-1)*FU< 1-11 )/DLU ) C FOR THE LAST GRID POINT ON THE BOTTOM HORiZONTAL LINE THE RIGHT C ADJACENT DIAGONAL ELEMENT OF L*U VANISHES SO THE COMPUTATION OF C EU(N) IS INCORRECT AND WILL BE CORRECTED IN STATEMENT 16. EU(I ) = PE/DL( I) 15 CONTINUE 16 EU(N)=O.ODOO C L AND M ARE COUNTER VARIABLES THAT AID IN CALCULATING THE CORRECT C SUBSCRIPTS FOR THE COMPONENTS OF L*U. M2 INCREMENTS THE Al INDEX C AT EACH HORIZONTAL LINE. N1=N+1 ,N2 1=N**2-2*N+1 . M=0 M2=0 L = 2 DO 30 J=N1,N21,N C j MOVES THE COMPUTATION UP TO THE NEXT HORIZONTAL LINE OF THE C GRID. FOR EACH J CALCULATIONS ARE DONE AT THE FIRST POINT, THE C INTERMEDIATE POINTS, AND THE LAST POINT OF THE GIVEN HORIZONTAL C LINE. J IS THE BASIS FOR THE INDICES OF Al, A2 , AND THE ELEMENTS C OF L*U. THE CALCULATIONS FOR THE FIRST POINT ON EACH HORIZONTAL C LINE ARE DONE HERE. M2=M2+2 PE=A1(J+M2+1)+A1( J+M2+2) PN=A2( J+N)+A2< J+NN) PW=A1(J+M2)+A1( J+M2+1) PS=A2( J)+A2< J+N) P0=-( PE+PN+PW+PS) BL(J-N)=PS/(1.0D00+ALPHA*EU(J-N-M) ) DL( J )=PO+BL( J-N)*( ALPHA*EU( J-N-M)-FU( J-N) ) EU(J-1-M)=(PE-ALPHA*BL( J-N)*EU( J-N-M) )/DL(J ) FU( J )=PN/DL( J) J1=J+1 J2=J+N-2 DO 25 I=J1,J2 C 1 MOVES THE COMPUTATION ALONG THE REMAINING GRID POINTS ON THE C GIVEN HORIZONTAL LINE AND THOSE COMPUTATIONS OF THE ELEMENTS OF C L*U APE DONE HERE. I IS THE BASIS FOR THE INDICES OF Al, A2, AND C THE ELEMENTS OF L*U. PE=A1(I+M2+1)+A1( I+M2+2) PN=A2( I+N)+A2( I+NN) PW = A1 ( I+M2 )+Al( 1+M2-H) K0=A2( I )+A2( I+N) PO=-(PE+PN+PW+PS) BL(I-N)-PS/( 1.0D00+ALPHA*EU( I-N-M) ) CL ( I-L )=PW/< 1.0nOO+ALPHA*FU( 1-1 ) ) 49 DL(I )*PO+BL(I-N>*+CLU-L)* C C CALCULATE THE PRODUCT M*T WHERE M IS A MATRIX AND T IS A VECTOR. C THE COMPUTATION IS DONE IN A LEFT-TO-RIGHT AND DOWN-TO-UP ORDER C OVER THE GRID IN ORDER TO TAKE ADVANTAGE OF THE VANISHING OF C COMPONENTS OF M; I.E. THE COMPUTATION OF RN»M*T IS PERFORMED C AT THE FIRST GRID POINT, THEN THE INTERMEDIATE GRID POINTS, AND C -THEN THE LAST GRID POINT OF EACH HORIZONTAL LINE IN ORDER TO TAKE C ADVANTAGE OF THE VANISHING OF COMPONENTS OF M ALONG THE BOTTOM C AND TOP HORIZONTAL LINES AND THE LEFT AND RIGHT VERTICAL LINES OF C THE GRID. C PS,PW,PO,PE, AND PN DEFINE THE LOWER DIAGONAL, LEFT ADJACENT C DIAGONAL, MAIN DIAGONAL, RIGHT ADJACENT DIAGONAL, AND THE UPPER C DIAGONAL OF THE MATRIX M. C C DEFINE CONSTANTS N1=N+1 NN1=N1+N NN=N+N N21=NSQ-N-N+1 C BEGIN PRODUCT CALCULATION PE = A1 (2)+AK3) PN=A2(N1 )+A2(NNl) PW=A1(1)+A1(2) PS = A2U )+A2(Nl) PO=-(PE+PN+PW+PS) C THE CALCULATION OF RN=M*T AT THE FIRST GRID POINT IS PERFORMED C HERE. NOTE THE COMPONENT PS OF M VANISHES ALONG THE BOTTOM C HORIZONTAL LINE. RN(1 )=PO*T< 1)+PE*T(2)+PN*T(N1) DO 15 1=2, N c i MOVES THE COMPUTATION ALONG THE REMAINING GRID POINTS ON THE C BOTTOM HORIZONTAL LINE. RN=M*T AT THESE POINTS IS COMPUTED HERE. I c IS THE BASIS FOR THE INDICES OF Al, A2 , AND RN. PE = A1 ( 1+1 )+AH 1+2) PN=A2( I+N)+A2( I+NN) PW=A1 ( I )+Al( 1+1 ) PS=A2( I )+A2( I+N) PO=-(PE+PN+PW+PS) C FOR I=N THE VALUE OF RN(N) IS INCORRECT AND IS CORRECTED BELOW. RN(I )=PW*T( 1-1 )+PO*T( I )+PE*T( I+1)+PN*T( I+N) 15 CONTINUE C RN=M*T AT THE LAST GRID POINT ON THE BOTTOM HORIZONTAL LINE IS C CALCULATED HERE SEPARATELY FROM THE INTERMEDIATE POINTS ON THE C LINE IN ORDER TO TAKE ADVANTAGE OF THE VANISHING COMPONENT PE OF C *. THIS IS DOME BY SUBTRACTING PE*T(N+1) FROM THE RN(N) JUST C CALCULATED. KN(N)=RN(N)-PE*T(N+1 ) C N1=N+1, N21=N**2-2*N+1 M1=N1 M2 = N 30 J=N1,N21,N 51 c j MOVES THE COMPUTATION UP TO THE NEXT HORIZONTAL LINE OF THE c GRID. FOR EACH J RN=M*T IS CALCULATED AT THE FIRST GRID POINT, THE c INTERMEDIATE POINTS, AND THE LAST GRID POINT OF THE GIVEN c HORIZONTAL LINE. Ml, M2 AND J ARE THE BASES FOR THE INDICES OF Al, c A2 , AND RN. RN=M*T FOR THE FIRST GRID POINT OF THE HORIZONTAL LINE c CALCULATED HERE IN ORDER TO TAKE ADVANTAGE OF THE VANISHING C COMPONENT PW OF M. MlxMl+3 M2=M2+1 PE=A1(M1)+A1(M1+1) PN=A2(M2)+A2(M2+N) PW=A1(M1-1)+A1(M1) PS=A2(M2-N)+A2(M2) PO=-(PE+PN+PW+PS) RN(J)=PS*T( J-N)+PO*T( J)+PE*T(J+l)«-PN*T(J + N> J1=J+1 J2=J+N-1 DO 25 I=J1,J2 c 1 MOVES THE COMPUTATION ALONG THE REMAINING GRID POINTS ON THE C GIVEN HORIZONTAL LINE AND RN=M*T FOR THOSE POINTS IS COMPUTED C HERE. Ml, M2 AND I ARE THE BASES FOR THE INDICES OF Al, A2 AND RN. M1=M1+1 M2=M2+1 PE = A1 (Ml )+Al (Ml + 1 ) PN=A2(M2)+A2(M2+N) PW = A1 (M1-1)+A1(M1) PS=A2(M2-N)+A2(M2) PO=-(PE+PN+PW+PS) RN(I )=PS*T( I-N)+PW*T( I-l)+PO*T( I )+PE*T( 1+1 )+PN*T( I+N) 25 CONTINUE C T he COMPUTATION OF RN=M*T AT THE LAST GRID POINT ON THE GIVEN c HORIZONTAL LINE IS DONE HERE SEPARATELY FROM THE OTHER GRID POINTS c ON THE LINE IN ORDER TO TAKE ADVANTAGE OF THE VANISHING COMPONENT C PE OF M. THIS IS DONE BY SUBTRACTING PE*T(J2+1) FROM THE RN(J2) C JUST CALCULATED. RN(J2)=RN( J2)-PE*T( J2+1) 30 CONTINUE L=J2+1 C l MOVES THE COMPUTATION UP TO THE FIRST GRID POINT ON THE TOP C HORIZONTAL LINE AND RN = M*T AT THAT POINT IS COMPUTED HERE. Ml, M2 C AND L ARE THE BASES FOR THE INDICES OF Al, A2 AND RN. NOTE THE C COMPONENT PN OF M VANISHES ALONG THIS LINE AND PW VANISHES AT C THIS POINT. Ml=Ml+3 M2=M2+1 PE=A1 (Ml )+Al(Ml+l) PN=A2(M2)+A2(M2+N) PW=A1 (Ml-1 )+Al(Ml) PS=A2(M2-N)+A2(M2) P0=-( PE+PN+PS+PW ) RN(L )=PS*T(L-N)+P0*T(L)+PE*T(L+1 ) L1=L+1 DO 45 K=L1,NS0 C K MOVES THE COMPUTATION ALONG THE REMAINING GRID POINTS ON THE C TOP HORIZONTAL LINE AND RN=M*T AT THOSE POINTS IS COMPUTED HERE. C Ml, M2 AND K ARE THE BASES FOR THE INDICES OF Al, A2 AND RN. M1=M1+1 52 M2=M2+1 PE=A1(M1)+A1(M1+1) PN=A2-CL(I-1)*RN>/DL(I) 10 CONTINUE M = 2 L = N DO 30 1=2, N RN(L+1)=(RN(L+1)-BL(L+1-N)*RN(L+1-N>>/DL(L+1) RN(J°L)=(RN(J + L)-BL(J + L-N)*RN(J + L-N,-CL(J + L-M)*RN(J + L-1))/DL*DELTA< 1+1 )+FU< I ) *DELTA ( I+N ) =V ( I ) 54 C C NOTE: IN ORDER TO SAVE STORAGE DELTA IS RESTORED IN THE VECTOR V C WHICH IN TURN IS RESTORED IN RN. ALSO f REMEMBER THAT EU AND FU C ONLY CONTAIN NON-ZERO VALUES OF THE DIAGONALS SO THEIR SUBSCRIPTS C ARE NOT EXACTLY THOSE REPRESENTED IN THE ABOVE RECURSION FORMULA. C NM1=N-1 V(NSQ)=V(NSQ) DO 10 I=1»NM1 N1=NSQ-I N3=NSQMN-I+1 V(NI)=V(N1)-EU(N3)*V ,DELTA(NSQ) C C CALCULATE THE NEW VALUE OF THE VECTOR T FROM THE OLD VALUE AND C THE NEW INCREMENT, DELTA .( REMEMBER DELTA IS STORED IN V WHICH C IN TURN IS STORED IN RN) I.E. CALCULATE: C c T(N+1)=T(N)+DELTA(N+1)=T(N)+RN(N+1) C DO 1 1=1, NSO T( I )=T( I )+DELTA( I ) 1 CONTINUE RETURN END SUBROUTINE B ALTRM ( N , A 1 , A2 , B ,C , D, E , F , ALPHA , NS0,NSQP2N ,NSQMN ) IMPLICIT RFAL*8(A-H,0-2 ) DIMENSION AKNS0P2N) , A2 ( NS0P2N ) , B < NSQMN ) ,C(NSO) ,D(NSO) ,E(NSQ),F(NS /OMN) 55 C CALCULATE THE ALTERED MATRIX M+N»L*U FROM THE GIVEN MATRIX M. L IS C A TRI-DIAGONAL MATRIX CONSISTING OF THE MAIN DIAGONAL, THE LEFT C ADJACENT DIAGONAL, AND A LOWER DIAGONAL AT A DISTANCE N FROM THE C MAIN ONE. THE ELEMENTS OF THESE DIAGONALS ARE STORED IN THE C VECTORS D, C, AND B RESPECTFULLY. U IS A TRI-DIAGONAL MATRIX C CONSISTING OF THE MAIN DIAGONAL EVERY ELEMENT OF WHICH IS EQUAL C TO THE CONSTANT ONE, THE RIGHT ADJACENT DIAGONAL, AND AN UPPER C DIAGONAL AT A DISTANCE N FROM THE MAIN ONE. THE ELEMENTS OF THE C the TWO OFF DIAGONALS ARE STORED IN THE VECTORS E AND F C RESPECTFULLY. VANISHING ELEMENTS OF C AND E ARE STORED AS ZEROES C IN THESE AFORE MENTIONED VECTORS. c THE COMPUTATION IS PERFORMED IN REVERSE ORDER FROM SUBROUTINE C ALTERM. THAT IS, THE SOLUTION IS COMPUTED AT THE GRID POINTS c IN THE J=N ROW WITH I ASSUMING THE VALUES 1 UP TO N. THEN C J IS DECREMENTED BY 1 AND I AGAIN ASSUMES THE RANGE OF C VALUES 1 THROUGH N. THIS LEFT-TO-RIGHT AND UP-TO-DOWN C ORDERING OF THE GRID POINTS FOR EVEN ITERATIONS BEGINNING WITH C THE SECOND MAY BE CONSIDERED A REORDERING OF THE ORIGINAL POINTS. C PS,PW,PO,PE, AND PN DEFINE THE LOWER DIAGONAL, LEFT ADJACENT C DIAGONAL, MAIN DIAGONAL, RIGHT ADJACENT DIAGONAL, AND UPPER C DIAGONAL OF THE MATRIX M. C C COMPUTE CONSTANTS NM2=N-2 NP1=N+1 NSQP1=NSQ+1 NSQPN=NSQ+N NM1S0=NSQ-N-N+1 NSQMNl=NSOMN+l C THE SUBSCRIPTS OF Al AND A2 FOLLOW THE NORMAL ORDERING. PE=A1 (NSOPN)+Al (NSOPN+1) PN=A2(NSQP1)+A2(NSQP1+N) PW = A1 (NSQPN-1 )+AKNSQPN) PS=A2(NSOPl-N)+A2(NSOPl) PO=-(PE+PN+PW+PS) D(l )=PO E(1)=PE/D(1) F(l )=PS/D(1) C(l )=O.ODOO C OBSERVE THAT THE MATRIX M IS REORDERED WITH THE RESULT THAT FOR C A GIVEN ROW I, PN AND PS ARE THE COEFFICIENTS OF X(I-N) AND XU+N) C RESPECTFULLY. THIS REVERSES THEIR ROLE IN THE OLD ORDERING WHERE C PN AND PS ARE THE COEFFICIENTS OF XU+N) AND X(I-N). C THE ROLE OF PE AND PW IS UNCHANGED. C THE FOLLOWING LOOP COMPUTES C,D,E, AND F ALONG THE TOP C HORIZONTAL LINE OF THE GRID. DO 10 1=2, N C LOOP PARAMETERS FOLLOW GRID POINTS IN THE NEW ORDERING. THESE C PARAMETERS AVOID SUPERFLUOUS ADDITIONS OF I. IM1=I-1 NS0I=NSQ+IM1 NSQP1I=NSQP1+IM1 NS0PNI=NSQPN+IM1 PE=A1 (NSQPNI )+Al(NSQPNI+l) PN=A2(NSQP1I )+A2(NSQPlI+N) PW=A1 (NSOPNI-1 )+Al (NSQPNI ) PS=A2 (NSQP1I-N)+A2(NSQP1I ) 56 PO=-(PE+PN+PW+PS) C(I )=PW/(1.0D00+ALPHA*F( 1-1 ) ) D(I )=P0+C(I)*(ALPHA*F(I-1>-E(I-1)> E(I )=PE/0(I ) C FOR I=N THIS VALUE OF E IS INCORRECT AND IS CORRECTED IN C STATEMENT 11. FII ) = 10 CONTINUE 11 E(N)=O.ODOO C THE NEXT SEGMENT COMPUTES B,CtD,E, AND F FROM THE SECOND C HORIZONTAL LINE (IN UP-TO-DOWN ORDER) TO AND INCLUDING THE N-l C HORIZONTAL LINE. THE SEGMENT IS A NESTED LOOP. ONE EXECUTION OF C THE OUTER LOOP COMPUTES THESE COMPONENTS ALONG THE FIXED C HORIZONTAL LINE. CERTAIN PARAMETERS AND COMPONENTS ARE ZERO AT THE C ENDPOINTS OF HORIZONTAL LINES. ACCORDINGLYt THE COMPUTATIONS C FOR A FIXED HORIZONTAL LINE PROCEED AS FOLLOWS: c (D c= o AND PW =o AT THE LEFT ENDPOINT. THIS FACT ABOUT PW C REQUIRES THE COMPUTATION AT THIS POINT TO BE A SPECIAL CASE C NOTE THAT SINCE C=0 AT THE ENDPOINTS* IT IS POSSIBLE TO C AVOID RESERVING CORRESPONDING STORAGE IN THE C ARRAY. THIS C WAS NOT DONE IN ORDER TO SIMPLIFY THE CALCULATION OF C SUBSCRIPTS OF C. C (2) POINTS 2 THROUGH N-l ON THE GIVEN HORIZONTAL LINE ARE C COMPUTED ITERATIVELY. C (3) E=0 AND PE=0 AT THE RIGHT ENDPOINTS. PE=0 MAKES THIS A C SPECIAL CASE. NOTE ALSO, THAT E=0 AT THESE ENDPOINTS MAKES C IT POSSIBLE TO AVOID RESERVING CORRESPONDING STORAGE IN THE C E ARRAY. THIS WAS NOT DONE IN ORDER TO SIMPLIFY THE C CALCULATION OF THE SUBSCRIPTS OF E. C (4) NOTE THAT THE VALUES OF C AND E THAT ARE ZERO AT CERTAIN C ENDPOINTS ARE STORED AS ZERO SO THAT IF ONE DESIRES A C PRINTOUT OF THE ARRAYS FOR DEBUGGING PURPOSES THE CORRECT C VALUES OF C AND E WILL BE STORED FOR THE CORRESPONDING C ENDPOINTS. C THE OUTER LOOP MOVES THE COMPUTATION FROM ONE HORIZONTAL LINE TO C THE NEXT LOWER. Kl=-N-2 C Kl STEPS Al DOWN K2 = -N C K2 STEPS A2 DOWN. THE LOOP VARIABLE FOLLOWS THE INDEX OF POINTS C OF THE LEFT VERTICAL LINE IN THE NEW ORDERING. DO 30 I=NP1,NM1S0tN Kl=Kl+N+2 K2=K2+N Ll=NS0-2-Kl C LI IS THE INDEX OF Al AT THE LEFT GRID POINT OF THE CURRENT C HORIZONTAL LINE. L2=NS0-N+1-K2 C L2 IS THE INDEX OF A2 AT THE LEFT GRID POINT OF THE CURRENT r HORIZONTAL LINE. PE=A1 (LI )+Al (Ll+1 ) PN=A2(L2)+A2(L2+N) PW = A1 (Ll-1 )+AKLl ) PS = A2(L2-N)+A2(L2 ) PO=-(PE+PN+PW+PS) [MN«I-N h( lMN)=PN/( 1 ,ODOO+ALPHA*E ( IMN) ) 57 C(I )=0.OD0O D( I )=PO+B( IMN)*( ALPHA*E( IMN)-F(IMN)) E( I )= C E = AND PE = AT THIS ENDPOINT. ZERO IS STORED IN E(JJ). JJ=JJ+1 JMN=JJ-N B( JMN)=PN C SINCE E = C( J J )=PW/(1.0D00 + ALPHA*F( JJ-1 ) ) D(JJ)=PO+B(JMN)*(-F( JMN) ) +C ( J J ) * ( ALPHA*F ( J J-l )-E ( J J-l ) ) E( JJ)=O.ODOO C SINCE PW DOES NOT APPEAR. F(JJ)=(PS-ALPHA*C< JJ)*F( JJ-1 ) )/D( JJ) 30 CONTINUE C THE COMPUTATION IS NOW AT THE FINAL HORIZONTAL LINE. PE = A1 (2)+AK3) PN=A2(NP1 )+A2(NPl+N) PW = A1(1 )+AK2) PS=A2(NP1-N)+A2(NP1) PO=-(PE+PN+PW+PS) C NM1SQ=(N**2)-(2*N) + 1 B(NMlSQ)=PN/( 1.0D00+ALPHA*E(NM1S0) ) C NSQMN1 = (N**2)-N+1 D(NS0MN1 )=P0+B(NM1S0)*(ALPHA*E(NM1S0)-F(NM1S0) ) EINSQMN1 )=(PE-ALPHA*B(NM1S0)*E(NM1S0) )/D(NSOMNl) C(NS0MN1 )=0.0D00 C THE COMPUTATION AT THE INTERIOR POINTS OF THE FINAL HORIZONTAL C LINE IS D0NE HERE. THE LOOP VARIABLE INDEXES THE GRID POINTS 58 IN THE NEW ORDERING. K ADJUSTS THE INDICES Kl AND K2 OF Al AND A2. K=-l L1=NSQMN1+1 L2=NSQ-1 DO 60 J=LltL2 K = K + 1 Kl=3+K K2=N+2+K PE=A1(K1)+A1(K1+1> PN=A2(K2)+A2(K2+N) PW=A1 IK1-1)+A1(K1) PS=A2(K2-N)+A2(K2) PO=-(PE+PN+PW+PS) JMN=J-N . B(JMN)=PN/(1.0D00+ALPHA*E( JMN) ) C(J)=PW D(J)=PO+B( JMN)*( ALPHA*E(JMN)-F( JMN) ) +C ( J ) * ( -E ( J-l ) ) EU) = (PE-ALPHA*B< JMN)*E( JMN) )/D(J) 60 CONTINUE WE ARE NOW AT THE RIGHT ENDPOINT. NSOMN= ( N**2 )-N. K1=N+1 K2=N+N PE = A1 (Kl )+AKKl + l ) PN=A2(K2)+A2(K2+N) PW=A1(K1-1)+A1(K1) PS=A2(K2-N)+A2(K2) PO=-(PE+PN+PW+PS> 8(NSQMN)=PN C(NSO)=PW D(NSQ)=PO+B(NSQMN)*(-F(NSQMN) ) +C ( NSO ) * (-E ( NSO-1 ) ) E(NSQ)=O.ODOO RETURN END SUBROUTINE BSOLVL ( N f B t C » D ,RN ,NSO t NSOMN ) IMPLICIT REAL*8(A-H,0-Z) DIMENSION B( NSOMN) ,C(NSQ) ,D(NSO) tRN(NSQ) C C THIS SUBROUTINE SOLVES THE EQUATION L*V=RN FOR V BY DIRECTLY C SOLVING THE SYSTEM OF EQUATIONS GENERATED BY THE ABOVE EQUATION, C 1_ is THE TRI-DIAGONAL MATRIX CALCULATED IN SUBROUTINE BALTRM AND C RN IS THE VECTOR CALCULATED IN SUBROUTINE RESIDL. I.E. C THE FOLLOWING RECURSION FORMULA IS SOLVED FOR V STARTING AT THE C FIRST GRID POINT AND MOVING IN A LEFT-TO-RIGHT AND DOWN-TO-UP C ORDER ACROSS THE GRID: C C BL(I)*V(I-N)+CL(I)*VU-l)+DLm*V(I)=RN(I) C C NOTE: IN ORDER TO SAVE STORAGE V IS RESTORED IN THE VECTOR RN. C fN(l) =RN( 1 )/D( 1 ) DO 10 I=?,N Ml I ) = (kN( 1 )-C( I )*RN( 1-1 ) )/D( I ) 59 10 CONTINUE N1=N+1 00 20 I=N1,NSQ RN( I )=*V(N1+N) 20 CONTINUE RETURN END APPENDIX B 60 ADI PROGRAM (DIRICHLET; 61 C ADI PROGRAM c T HIS PROGRAM COMPUTES THE SOLUTION X OF THE LINEAR SYSTEM OF c EQUATIONS A*X = S AS THE LIMIT OF THE SEQUENCE DEFINED BY: c (H + W*I )*X(K+1/2)=S-(V-W*I )*X(K) c (V+W*I )*X(K+1)=S-(H-W*I )*X(K+l/2) C IMPLICIT REAL*8 ,X(900),S<900),H<900),ABND<6) ,WH(33)» XWV(33) ,THETA(6) C INPUT N,S»X» Alt A2tNIT,ABND»BBND,CBND f DBND C c NOTE: IT IS NECESSARY FOR THE USER TO INPUT THE FOLLOWING C VARIABLES INTO THIS PROGRAM: C N = NUMBER OF GRID POINTS ON ONE SIDE OF THE SOLUTION GRID C (THERE ARE N**2 GRID POINTS IN THE TOTAL GRID SYSTEM). C S = THE VECTOR CONTAINING THE RIGHT HAND SIDES OF THE LINEAR C SYSTEM. C X = THE INITIAL VECTOR. c A1,A2 = THE VECTORS OF THE DIFFERENTIAL EQUATION COEFFICIENTS c NIT = THE POWER OF 2 ITERATIONS TO BE PERFORMED BY THIS C PROGRAM. C ABND(l) = UPPER BOUND ON THE EIGENVALUES OF H. C BBND = LOWER BOUND ON THE EIGENVALUES OF H. C CBND = UPPER BOUND ON THE EIGENVALUES OF V. C DBND = LOWER BOUND ON THE EIGENVALUES OF V. C c THE INPUT THAT FOLLOWS IS A SAMPLE OF THE NECESSARY DATA. N = 30 NIT = 5 NSQ=N*N NSQP2N=NSQ+N+N DN=DFLOAT(N) DNl=DFLOAT(N+l) PI=3. 1415926535 89793 DO 200 I=1,NSQP2N Al ( I )=-.5D00 200 A2( I )=-.5D00 DO 201 I=lt900 X( I )=0.0D00 200 S( I )=1.0D00 ABND(1 )=4.0D00*DSIN(DN*PI/(2.0D00*DN1 ) )**2 BBND=4.0D00*DSIN(PI/(2.0D00*DN1 ) )**2 CBND=ABND(1 ) DBND=BBND C COMPUTE THE SUM OF THE POSITIVE RIGHT HAND SIDES FOR USAGE IN C THE NORMALIZING OF THE MAXIMUM RESIDUAL. SUM=O.ODOO DO 3 1=1, NSQ IF(S( I ) .LE.O.ODOOGO TO 3 SUM=SUM+S( I ) 3 CONTINUE C SET PARAMETERS MCNT=0 NSQ=N*N NSQMN=N*N-N NSQP2N=NSQ+N+N C THE OPTIMAL ITERATION PARAMETERS FOR MIT = 2**NIT ITERATIONS ARE 62 c COMPUTED HERE BY SUBROUTINE OPTW. THEY ARE RETURNED IN THE C VECTORS WH AND WV RESPECTIVELY, IN DECREASING ORDER AND ARE C APPLIED TO THE ITERATION STEPS IN THAT SAME ORDER, MIT=2**NIT MIT1=MIT+1 NIT1=NIT+1 CALL 0PTW(NIT1,MIT1,WV,WH,ABND,BBND,CBND,DBND,THETA,PMU) WRITE(6,102) c the ACTUAL ADI ITERATION BEGINS HERE. MCNT KEEPS TRACK OF THE C ITERATION PARAMETERS. FOR A DESCRIPTION OF THE FUNCTIONS OF EACH C SUBROUTINE AND ITS RELATION TO THE ITERATIVE PROCEDURE, SEE C the COMMENTS AT THE BEGINNING OF THE INDIVIDUAL SUBROUTINES. 4 MCNT=MCNT+1 c THE FIRST HALF STEP OF THE ITERATIVE SEQUENCE IS PERFORMED HERE. W=WH(MCNT) C COMPUTE THE PRODUCT S-( V-W*I ) *X ( K ) . CALL PROD V(N,S, NSQ, A2,W,X,NSQP2N) C COMPUTE THE INVERSE OF (H+W*I) BY TR I ANGULAR I ZAT ION . CALL HH1 (N,A1,NSQMN,W,H,NSQ,NSQP2N) C COMPUTE X(K+l/2) BY BACKWARD SUBSTITUTION. CALL PPKN,A1,NSQMN,W,H,NSQ,X,NSQP2N) C c THE SECOND HALF STEP OF THE ITERATIVE SEQUENCE IS PERFORMED HERE. W=WV(MCNT) C COMPUTE THE PRODUCT S- (H-W*I ) *X ( K+l /2 ) . CALL PRODH(N,S,NSQ,Al,NSQMN,W,X,NSQP2N) c COMPUTE THE INVERSE OF *< > 101 FORMATCOTHE OPTIMAL PARAMETERS WH IN DECREASING ORDER ARE:'/,(D16 X 8 ) ) 102 FORMATCOTHE OPTIMAL PARAMETERS WV IN DECREASING ORDER AREC/,(D16 X.8) ) 103 FORMATCOTHE SPECTRAL RADIUS =«,D16.8) RETURN END SUBROUTINE DECR ( NN ,WH , MIT1 > IMPLICIT REAL*8 IMPLICIT REAL*8(A-H,0-Z) DIMENSION S(NSO) , X ( NSO ) t A2 ( NS0P2N ) C c CALCULATE THE PRODUCT ( S- < V-W*I ) *X ) . COMPUTATION OF THE c COMPONENTS IS DONE IN DOWN-TO-UP AND LEFT-TO-RIGHT ORDER TO C TAKE ADVANTAGE OF THE VANISHING OF COMPONENTS OF V ALONG THE C BOTTOM OR TOP LINES OF THE GRID. C PS,PN,AND PO DEFINE THE LOWER DIAGONAL, UPPER DIAGONAL* AND c HA i N DIAGONAL ( SOUTHfNORTHt ORIGIN) OF V. C DO 2 1=1. N r i FOLLOWS THE LOWER HORIZONTAL LINE. FOR EACH I, THE COMPONENTS rn^ PRODUCT ARE COMPUTED ALONG THE CORRESPONDING VERTICAL LINE. 65 c THE FIRST COMPUTATION FOR EACH It ANO THE LAST, ARE SPECIAL C CASES OCCURRING WHEN PS ANO PN ARE ZERO. c I MOVES THE COMPUTATIONS FROM LEFT TO RIGHT ALONG THE C (BOTTOM) HORIZONTAL LINE. IPN=I+N PS=A2(I )+A2(IPN) PN=A2(IPN)+A2( IPN+N) PO=-(PS+PN) C SAVE X(I) TO USE AGAIN X1 = X(I ) X(I )«S(I )-< (PO-W)*X(I)+PN*X( IPN)) c THE VARIABLE OF THIS LOOP STEPS UP EACH VERTICAL LINE IN C ORDER FROM LEFT TO RIGHT. THE LOOP VARIABLE IS THE INOEX OF X. C HENCE IPN IS THE INDEX OF X AT THE LOWEST POINT AND LI IS THE C INDEX OF X AT THE HIGHEST POINT. L1=NSQ-N-N+I DO 1 J=IPN,LltN PS = PN C JPN IS THE INDEX OF A2 JPN=J+N PN=A2(JPN)+A2( JPN+N) PO=-(PS+PN) X2=X(J) X(J)=S( J)-(PS*Xl+(PO-W)*X( J)+PN*X< JPN) ) X1 = X2 1 CONTINUE C THIS SEGMENT COMPUTES THE COMPONENTS ALONG THE UPPER C HORIZONTAL LINE. PS = PN C K2 IS THE X INDEX C |_1 = N**2-2*N+K K2=L1+N C K3 IS THE A INDEX K3=K2+N PN=A2(K3)+A2(K3+N) PO=-(PS+PN) X(K2)=S(K2)-(PS*Xl+(PO-W)*X(K2) ) 2 CONTINUE RETURN END SUBROUTINE HH1 ( N, Al , NSQMN, W ,H ,NSQ , NSQP2N ) IMPLICIT REAL*8( A-H t O-Z) DIMENSION H(NSQ)»A1(NSQP2N) C C THIS SUBROUTINE COMPUTES THE SOLUTION OF PW = PE 4- CONTINUE I 1 = 1 I+N INM1=INM1+N M=M+N1 5 CONTINUE RETURN END MJBROUT INE PP1 ( N , A 1 , N SOMN , W ,H ,NSO , X ,NSQP2N ) IMPLICIT REAL*8 (A-H t O-Z» 67 DIMENSION H(NSO) , X ( NSO > , Al ( NSOP2N ) C c this SUBROUTINE CONTINUES THE GAUSSIAN ELIMINATION PROCESS c B Y CALCULATING THE RECURSION FORMULA FOR THE RIGHT HAND SIDES; c i.e. it CALCULATES P(M)=X(M) FOR M«1,...,N BY THE FOLLOWING C FORMULA: C C P(M)*X(M)«(X(M)-PW*X(M-l))/(PO*W+PW*H(M-n ) C WHERE C P(l)»X(l)=X(l)/(PO^W) c C AND WHERE PW, PE, AND PO ARE THE LEFT DIAGONAL, RIGHT DIAGONAL, C AND MAIN DIAGONAL ELEMENTS OF H AT POINT M. C COMPUTATION IS DONE IN A LEFT-TO-RIGHT AND DOWN-TO-UP ORDER C vER THE GRID. FOR EACH HORIZONTAL LINE P(M)=X(M) IS CALCULATED C AT THE FIRST POINT AND THEN AT THE REMAINING POINTS IN ORDER TO c TAKE ADVANTAGE OF THE FACT THAT P(M) DOES NOT DEPEND ON A PREVIOUS C H AT THE FIRST POINT AND DOES DEPEND ON A PREVIOUS H AT THE C REMAINING POINTS. THE I LOOP MOVES THE COMPUTATION UP TO C THE NEXT HORIZONTAL LINE OF THE GRID, C INN=1 11=2 IN = N M = l DO 2 1=1, N C COMPUTATION FOR THE FIRST POINT ON THE GIVEN HORIZONTAL LINE C is PERFORMED HERE. INN IS THE X INDEX AND IM IS THE Al INDEX. IM=INN+M PW = A1( IM-1 )+AK IM) PE=A1(IM)+A1( IM+1) PO=-(PE+PW) X( INN)=X( INN)/(PO+W) DO 1 K=II,IN C THIS LOOP PERFORMS THE COMPUTATION AT THE REMAINING GRID POINTS C ALONG THE GIVEN HORIZONTAL LINE. K IS THE X INDEX AND KM IS THE C Al INDEX. KM=K+M PW = PE PE=A1(KM)+A1(KM+1) PO=-(PE+PW) KK=K-1 X(K)=(X(K)-PW*X(KK) )/{PO+W + PW*H(KK) ) 1 CONTINUE INN=INN+N II=II+N IN=IN+N M=M+2 2 CONTINUE C THIS SECTION COMPUTES THE BACKWARD SUBSTITUTION FOR THE GAUSSIAN C ELIMINATION PROCEDURE; I.E. IT COMPUTES: C C X(N)=P(N)=X(N) C X(M)=P PO=-(PE+PW) X(L)=S(L)-(PW*XH-(PO-W)*X(L) ) M = M + 2 2 CONTINUE RETURN END SUBROUTINE HH2 ( N , A2 ,NSQMN , W ,H ,NSQ,NSQP2N ) IMPLICIT REAL*8(A-H,0-Z) DIMENSION H(NSQMN) ,A2(NSQP2N) C C THIS SUBROUTINE SOLVES (V+W*I)X=S BY GAUSSIAN C ELIMINATION. SINCE V IS A TRI-DIAGONAL MATRIX CONSISTING OF A C MAIN DIAGONAL AND TWO OFF DIAGONALS AT A DISTANCE N FROM THE C MAIN DIAGONAL, GAUSSIAN ELIMINATION CAN BE PERFORMED DIRECTLY C ON THE MATRIX. C FIRST THE FORWARD ELIMINATION PROCESS IS PERFORMED C TRANSFORMING THE TRI-DIAGONAL MATRIX INTO AN UPPER DIAGONAL MATRIX C (SUBROUTINE HH2 ) THEN (SUBROUTINE PP2) THE RIGHT HAND SIDE VECTOR C IS TRANSFORMED IN ACCORDANCE WITH THE UPPER DIAGONAL MATRIX; I.E. C THE P(I)=X(I) ARE CALCULATED. THEN THE BACKWARD SUBSTITUTION C (SUBROUTINE PP1) IS PERFORMED TO GET THE RESULT. C ALONG THE FIRST HORIZONTAL LINE OF THE GRID THE VALUES OF C THE UPPER DIAGONAL COMPONENT, H, DO NOT DEPEND ON ANY PREVIOUS C VALUES OF H. PS,PN AND PO DEFINE THE LOWER DIAGONAL, UPPER C DIAGONAL AND MAIN DIAGONAL ELEMENTS OF V. C THE COMPUTATION IS DONE IN A DOWN-TO-UP AND LEFT-TO-RIGHT C ORDER TO TAKE ADVANTAGE OF THE VANISHING OF COMPONENTS OF V ALONG C THE BOTTOM AND TOP LINES OF THE GRID. C DO 2 I=1,N C i FOLLOWS THE LOWER HORIZONTAL LINE. FOR EACH I THE VALUES OF H C FOR THE TRANSFORMED UPPER DIAGONAL MATRIX ARE COMPUTED ALONG THE C CORRESPONDING VERTICAL LINE. I MOVES THE COMPUTATIONS FROM LEFT TO C RIGHT ALONG THE (BOTTOM) HORIZONTAL LINE. I IS THE H INDEX AND C ip N is THE A2 INDEX. IPN=I+N PS=A2( I )+A2( IPN) PN=A2( IPN)+A2( IPN+N) PO=-(PS+PN) H(I )=-PN/(PO+W) C J STEPS UP EACH VERTICAL LINE IN ORDER FROM LEFT TO RIGHT. J IS C THE H i NDE x AND JPN IS THE A2 INDEX. Ll=NSO-N-N+I 70 00 1 J«IPN f LltN PSxPN JPN*J+N PN=A2(JPN)+A2(JPN+N) P0=-(PS*PN) H(J)=-PN/(PO+W+PS*H( J-N) ) 1 CONTINUE 2 CONTINUE C THE VALUE OF H FOR THE TRANSFORMED UPPER DIAGONAL MATRIX C ALONG THE TOP HORIZONTAL LINE OF THE GRID SO IT NEED NOT C COMPUTED AND STORED, RETURN END IS BE ZERO SUBROUTINE PP2 ( N, A2 ,NSQMN, W,H,NSQ , X ,NSQP2N ) IMPLICIT REAL*8( A-H,0-Z) DIMENSION H(NSQMN),X(NSQ),A2(NSQP2N) C C THIS SUBROUTINE CONTINUES THE GAUSSIAN ELIMINATION BY CALCULATING C THE RECURSION FORMULA FOR THE RIGHT HAND SIDES; I.E. IT CALCULATES C THE RIGHT HAND SIDES P(M)=X(M) BY THE FOLLOWING FORMULA: C P(M)=X(M) FOR M=1,...,N C P X(I )=X( I )/(PO+W) C J STEPS UP EACH VERTICAL LINE IN ORDER FROM LEFT TO RIGHT. J IS c the X INDEX AND JPN IS THE INDEX OF A2. L1=NSQ-N+I DO 1 J=IPN,L1,N PS = PN JPN=J+N PN=A2 ( JPN)+A2( JPN+N) PO=-(PS+PN) X( J ) = (X( J)-PS*X( J-N) )/(PO+W+PS*H(J-N) ) 1 CONTINUE ? CONTINUE r. THIS SECTION OE THE SUBROUTINE PERFORMS THE BACKWARD SUBSTITUTION 71 C FOR THE GAUSSIAN ELIMINATION PROCESS; I.E. IT COMPUTES: C C X(M)»P(M)=X(M) FOR M*N,...,N**2-N+1 C X(M)=P(M)+H+A2(NNl ) PW=A1 ( I )+Al (2) PS = A2( 1 )+A2(Nl ) PO=-(PE+PN+PW+PS) C THE CALCULATION OF RN=M*T AT THE FIRST GRID POINT IS PERFORMED C HERE. NOTE THE COMPONENT PS OF M VANISHES ALONG THE BOTTOM C HORIZONTAL LINE. RN(1 )=PO*T(l )+PE*T(2)+PN*T(Nl) DO 15 1=2, N C 1 MOVES THE COMPUTATION ALONG THE REMAINING GRID POINTS ON THE C BOTTOM HORIZONTAL LINE. RN=M*T AT THESE POINTS IS COMPUTED HERE. I C IS THE BASIS FOR THE INDICES OF Al, A2, AND RN. PE = A1 ( I + D+AK 1 + 2) PN=A2( I+N)+A2( I+NN) PW = A1 ( I )+AK 1 + 1) 72 PS=A2(I )+A2< I+N) PO=-(PE+PN+PW+PS) C F0R I=N THE VALUE OF RN(N) IS INCORRECT AND IS CORRECTED BELOW, RNU )=PW*T( I-l)+PO*T( I )+PE*T(I+l)+PN*T(I+N) 15 CONTINUE RN=M*T AT THE LAST GRID POINT ON THE BOTTOM HORIZONTAL LINE IS CALCULATED HERE SEPARATELY FROM THE INTERMEDIATE POINTS ON THE LINE IN ORDER TO TAKE ADVANTAGE OF THE VANISHING COMPONENT PE OF M. THIS IS DONE BY SUBTRACTING PE*T FROM THE RN(N) JUST CALCULATED. RN(N)=RN(N)-PE*T(N+1 ) N1=N+1, N21=N**2-2*N+1 M1=N1 M2 = N DO 30 J=N1,N21,N J MOVES THE COMPUTATION UP TO THE NEXT HORIZONTAL LINE OF THE GRID. FOR EACH J RN=M*T IS CALCULATED AT THE FIRST GRID POINT, THE INTERMEDIATE POINTS, AND THE LAST GRID POINT OF THE GIVEN HORIZONTAL LINE. Ml, M2 AND J ARE THE BASES FOR THE INDICES OF Alt A2, AND RN. RN=M*T FOR THE FIRST GRID POINT OF THE HORIZONTAL LINE CALCULATED HERE IN ORDER TO TAKE ADVANTAGE OF THE VANISHING COMPONENT PW OF M. Ml=Ml+3 M2=M2+1 PE = A1 (Ml )+Al (Ml + 1 ) PN=A2(M2 )+A2(M2+N) PW=A1 (Ml-1 )+Al (Ml) PS=A2(M2-N)+A2(M2) PO=-(PE+PN+PW+PS) RN( J)=PS*T( J-N)+PO*T( J)+PE*T( J + l )+PN*T( J+N) J1=J+1 J2=J+N-1 DO 25 I=J1,J2 1 MOVES THE COMPUTATION ALONG THE REMAINING GRID POINTS ON THE GIVEN HORIZONTAL LINE AND RN=M*T FOR THOSE POINTS IS COMPUTED HERE. Ml, M2 AND I ARE THE BASES FOR THE INDICES OF Al, A2 AND RN. M1=M1+1 M2=M2+1 PE = A1 (Ml )+AKMl + l) PN=A2(M2)+A2(M2+N) PW=A1 (Ml-1 )+Al (Ml) PS=A2(M2-N)+A2(M2) PO=-(PE+PN+PW+PS) RN( I )=PS*T( I-N)+PW*T( I-1)+P0*T( I )+PE*T( 1+1 )+PN*T( I+N) 25 CONTINUE THE COMPUTATION OF RN=M*T AT THE LAST GRID POINT ON THE GIVEN HORIZONTAL LINE IS DONE HERE SEPARATELY FROM THE OTHER GRID POINTS ON THE LINE IN ORDER TO TAKE ADVANTAGE OF THE VANISHING COMPONENT PE OF M. THIS IS DONE BY SUBTRACTING PE*T(J2+l) FROM THE RN(J2) JUST CALCULATED. PN(J2)=RN(J2)-PE*T(J2+1) 30 CONTINUE L=J2+1 L MOVES THE COMPUTATION UP TO THE FIRST GRID POINT ON THE TOP HORIZONTAL LINE AND RN=M*T AT THAT POINT IS COMPUTED HERE. Ml, M2 AND L ARE THE BASES FOR THE INDICES OF Al, A2 AND RN. NOTE THE COMPONENT PN OF M VANISHES ALONG THIS LINE AND PW VANISHES AT 73 THIS POINT. Ml=Ml+3 M2=M2+1 PE = A1 (Ml l+AKMl + 1) PN=A2(M2)+A2(M2+N) PW=A1 (Ml-1 I+A1IM1) PS=A2(M2-N)+A2(M2) PO=-(PE+PN+PS*PW) RN(L)=PS*T(L-N)+PO*T(L)+PE*T BL(J-N)=PS/( 1.0D00+ALPHA*EU< J-N-M) ) DL(J)=PO+BL( J-N)*(ALPHA*EU( J-N-M ) -FU ( J-N ) ) C ****»THIS CHANGE IS REOUIRED FOR THE NEUMANN PROBLEM. EU(J-1-M)=(PE+PE-ALPHA*BL( J-N )*EU< J-N-M) )/DL(J) FU(J)=PN/DL( J) J1=J+1 J2=J+N-2 DO 25 I=J1,J2 c i MOVES THE COMPUTATION ALONG THE REMAINING GRID POINTS ON THE C GIVEN HORIZONTAL LINE AND THOSE COMPUTATIONS OF THE ELEMENTS OF c L *U are DONE HERE. I IS THE BASIS FOR THE INDICES OF Al, A2, AND c the ELEMENTS OF L*U. PE = A1 ( I+M2 + D+AK I+M2+2) PN=A2( I+N)+A2( I+NN) PW=A1 (I+M2)+A1( I+M2+1) PS = A2( I )+A2( I+N) PO=-(PE+PN+PW+PS) BL( I-N)=PS/( 1.0D00+ALPHA*EU( I-N-M) ) CL( I-L)=PW/( 1.0D00+ALPHA*FU( 1-1) ) DL(I )=P0+BL(I-N)*(ALPHA*EU(I-N-M)-FU(I-N))+CL(I-L)*(ALPHA*FU(I-1) XEU( I-2-M) ) EU(I-1-M)=(PE-ALPHA*BL( I-N)*EU( I-N-M) )/DL(I ) FU( I )=(PN-ALPHA*CL( I-L)*FU( 1-1) )/DL( I ) 25 CONTINUE I=.J2 + 1 C THE COMPUTATION OF THE ELEMENTS OF L*U FOR THE LAST POINT ON THE c GIVEN HORIZONTAL LINE ARE DONE HERE SEPARATELY FROM THE OTHER GRID C POINTS ON THE LINE IN OROER TO TAKE ADVANTAGE OF THE VANISHING C OF THE COMPONENT EU OF L*U AT THIS POINT. 77 BL(J-1)=PS C *«#**THIS CHANGE IS REQUIRED FOR THE NEUMANN PROBLEM. CL(I-L) = (PW + PW)/(1.0DOO-»-ALPHA*FU(I-l)) DLU >=PO+BL< I-N)*(-FU< I-N) ) +CL ( I-L >*( ALPHA*FU( 1-1 ) -EU ( I-2-M ) ) FU( I )= D(I )«P0+C(I)*(ALPHA*FU-1)-EU-1>) E(I )»PE/D( I ) C FOR I»N THIS VALUE OF E IS INCORRECT AND IS CORRECTED IN C STATEMENT 11. C»****THIS CHANGE IS REQUIRED FOR THE NEUMANN PROBLEM. F(I ) = (PS + PS-ALPHA*C( I)*F< I-1))/D(I ) 10 CONTINUE 11 E(N)=O.ODOO C*****THIS CHANGE IS REQUIRED FOR THE NEUMANN PROBLEM. C(N)=(PW+PW)/<1.0D00+ALPHA*F(N-1) ) C THE NEXT SEGMENT COMPUTES B,C,D,E, AND F FROM THE SECOND C HORIZONTAL LINE (IN UP-TO-DOWN ORDER) TO AND INCLUDING THE N-l C HORIZONTAL LINE. THE SEGMENT IS A NESTED LOOP. ONE EXECUTION OF C THE OUTER LOOP COMPUTES THESE COMPONENTS ALONG THE FIXED C HORIZONTAL LINE. CERTAIN PARAMETERS AND COMPONENTS ARE ZERO AT THE C ENDPOINTS OF HORIZONTAL LINES. ACCORDINGLY, THE COMPUTATIONS C FOR A FIXED HORIZONTAL LINE PROCEED AS FOLLOWS: C (1) C = AND PW=0 AT THE LEFT ENDPOINT. THIS FACT ABOUT PW C REQUIRES THE COMPUTATION AT THIS POINT TO BE A SPECIAL CASE C N0 TE THAT SINCE C*0 AT THE ENDPOINTS, IT IS POSSIBLE TO C AVOID RESERVING CORRESPONDING STORAGE IN THE C ARRAY. THIS C WAS NOT DONE IN ORDER TO SIMPLIFY THE CALCULATION OF C SUBSCRIPTS OF C. C (2) POINTS 2 THROUGH N-l ON THE GIVEN HORIZONTAL LINE ARE C COMPUTED ITERATIVELY. c ( 3) E = AND PE = AT THE right ENDPOINTS. PE = MAKES THIS A C SPECIAL CASE. NOTE ALSO, THAT E = AT THESE ENDPOINTS MAKES C IT POSSIBLE TO AVOID RESERVING CORRESPONDING STORAGE IN THE C E ARRAY. THIS WAS NOT DONE IN ORDER TO SIMPLIFY THE C CALCULATION OF THE SUBSCRIPTS OF E. C (4) NOTE THAT THE VALUES OF C AND E THAT ARE ZERO AT CERTAIN C ENDPOINTS ARE STORED AS ZERO SO THAT IF ONE DESIRES A C PRINTOUT OF THE ARRAYS FOR DEBUGGING PURPOSES THE CORRECT C VALUES OF C AND E WILL BE STORED FOR THE CORRESPONDING C ENDPOINTS. C THE OUTER LOOP MOVES THE COMPUTATION FROM ONE HORIZONTAL LINE TO C THE NEXT LOWER. Kl=-N-2 C Kl STEPS Al DOWN K2=-N C K2 STEPS A2 DOWN. THE LOOP VARIABLE FOLLOWS THE INDEX OF POINTS C OF THE LEFT VERTICAL LINE IN THE NEW ORDERING. DO 30 I=NP1,NM1SQ,N Kl=KH-N+2 K2=K2+N Ll=NSQ-2-Kl C L1 is THE INDEX OF Al AT THE LEFT GRID POINT OF THE CURRENT C HORIZONTAL LINE. L2=NSQ-N+1-K2 80 •L2 IS THE INDEX OF A2 AT THE LEFT GRID POINT OF THE CURRENT -HORIZONTAL LINE. PE»A1 PW«Al PS»A2(L2-N)+A2(L2) PO=-(PE*PN+PW«-PS) IMN»I-N B(IMN)«PN/<1.0DOO+ALPHA*EUMN>) C(I)«0.0D00 D(I )«PO*B( IMN)*(ALPHA*E(IMN)-F(IMN>) C ***** T HIS CHANGE IS REQUIRED FOR THE NEUMANN PROBLEM. E(I)=/D+A1(L1J> PS=A2(L2J-N)+A2(L2J) PO=-(PE+PN+PW+PS) JMN=JJ-N BUMN)=PN/(1.0+ALPHA*E< JMN) ) C(JJ)=PW/(1.0D00+ALPHA*F( J J-l > ) D(JJ) = PO+B*H-C(JJ)*ULPHA*F ) E(JJ ) = 0.0D00 C SINCE PW DOES NOT APPEAR. FUJ )=(PS-ALPHA*C< JJ)*F< JJ-1) )/D< J J ) 30 CONTINUE C THE COMPUTATION IS NOW AT THE FINAL HORIZONTAL LINE. PF = A1 (2 H-Al (3) PN=A2(NPl H-A2(NP1+N) 81 PW«A1U)«-A1<2) PS»A2(NP1-N)«-A2 C(NSQMN1)=0.0D00 C THE COMPUTATION AT THE INTERIOR POINTS OF THE FINAL HORIZONTAL C LINE IS DONE HERE. THE LOOP VARIABLE INDEXES THE GRID POINTS c im T he NEW ORDERING. K ADJUSTS THE INDICES Kl AND K2 OF Al AND A2 K = -l L1=NSQMN1+1 L2=NSQ-1 DO 60 J=LltL2 K=K + 1 Kl=3+K K2=N+2+K PE=A1(K1)+A1(K1+1> PN=A2*(-E( J-l ) ) E (J )=(PE-ALPHA*B( JMN)*E( JMN) )/D( J) 60 CONTINUE C WE ARE NOW AT THE RIGHT ENDPOINT. NSQMN= ( N**2 )-N. K1=N + 1 K2=N+N PE = A1 (Kl )+AL(Kl+l) PN=A2(K2)+A2(K2+N) PW = A1 (K1-1I+AHK1) PS=A2(K2-N)+A2(K2) PO=-(PE+PN+PW+PS) C***«*THIS CHANGE IS REQUIRED FOR THE NEUMANN PROBLEM. B(NSQMN)=PN+PN C1) PN=A2(M2)*A2(M2*N) PW«AKMl-l)+Ai(Ml) PS*A2(M2-N)+A2(M2) PO«-(PE+PN+PW+PS) C FOR K=NSO THE VALUE OF RN(NSO) IS COMPUTED BELOW. IF GO TO 40 C*****JHIS CHANGE IS REQUIRED FOR THE NEUMANN PROBLEM. RN(K)=(PS+PS)*T(K-N)+PW*T(K-1)+P0*T(K)*PE*T(K+1) GO TO 45 c the COMPUTATION RN=M*T AT THE LAST GRID POINT ON THE TOP C HORIZONTAL LINE IS PERFORMED HERE SEPARATELY FROM THE OTHER GRID C POINTS IN ORDER TO TAKE ADVANTAGE OF THE VANISHING OF THE C COMPONENTS PE AND PN OF M. C *#**# T HIS CHANGE IS REQUIRED FOR THE NEUMANN PROBLEM. 40 RN(K)=(PS+PS)*T(K-N)+(PW+PW)*T DIMENSION H(NSQ),A1(NSQP2N) C C THIS SUBROUTINE COMPUTES THE SOLUTION OF (H+W*I)X«S BY GAUSSIAN C ELIMINATION, THE MATRIX (H+W*I) IS TRIDIAGONAL WITH THE MAIN C DIAGONAL AND THE TWO ADJACENT DIAGONALS ON EITHER SIDE OF C THE MAIN DIAGONAL. C COMPUTATION IS DONE IN THE FOLLOWING MANNER: FIRST FORWARD C ELIMINATION IS PERFORMED TRANSFORMING THE TRIDIAGONAL MATRIX C INTO AN UPPER DIAGONAL MATRIX (SUBROUTINE HH1). THEN (SUBROUTINE C PP1) THE RIGHT HAND SIDE VECTOR IS TRANSFORMEDt THEN THE C BACKWARD SUBSTITUTION IS PERFORMED TO SOLVE THE RESULTING C UPPER TRIANGULAR SYSTEM. H IS SET EQUAL TO ZERO ALONG THE C RIGHTMOST VERTICAL LINE OF THE GRID, WHICH RESULTS FROM THE C FACT THAT THE RIGHT ADJACENT DIAGONAL VANISHES AT THESE POINTS. C DO 1 K=N,NSQ,N H(K)=0.0 1 CONTINUE C NEXT THE COMPUTATION IS PERFORMED ON THE FIRST VERTICAL LINE OF C THE GRID SINCE H AT THESE POINTS DEPENDS ONLY ON THE DIAGONAL C ELEMENT AND THE RIGHT ADJACENT DIAGONAL ELEMENT, AND DOES NOT C DEPEND ON ANY PREVIOUS H. C PW,PE,AND PO DEFINE THE LEFT DIAGONAL, RIGHT DIAGONAL, AND C MAIN DIAGONAL ELEMENTS OF THE MATRIX H. NTOP=NSQMN+l M=l DO 2 I=l,NTOP,N C 1 IS THE H INDEX AND MOVES THE COMPUTATION UP THE FIRST VERTICAL C LINE. IM IS THE Al INDEX. IM=I+M PW=A1 ( IM-1 )+Al( IM) PE=A1 ( IM)+A1( IM+1) PO=-(PE+PW) C*****THIS CHANGE IS REQUIRED FOR THE NEUMANN PROBLEM. H(I )=-(PE+PE)/(PO+W) M = M + 2 2 CONTINUE C THE COMPUTATION OF H I S NOW PERFORMED IN LEFT-TO-RIGHT AND DOWN- C TO-UP ORDER; I.E. H IS CALCULATED ALONG THE REMAINING GRID POINTS C ON THE FIRST HORIZONTAL LINE STARTING WITH THE FIRST AND I MOVES C THE COMPUTATION UP TO THE NEXT HORIZONTAL LINE. N1=N+1 M=l I 1=2 INM1=N-1 DO 5 1=1, N PW = A1 ( I+M)+A1 ( I+M+l ) DO 4 J=l I t INM1 C J IS THE H INDEX AND Jl IS THE Al INDEX J1=J+1 PE = A1 (Jl )+Al ( Jl + 1 ) f>0»-( pf-f pw) H(J )=-PE/ ( PQ + W+PW*H( J-l ) ) ^--^ <* CONTINUE 87 II=II+N INM1=INM1+N M=M+N1 CONTINUE RETURN END 88 SUBROUTINE PP1 ( N, Al ,NSQMN, W,H ,NSQ,X,NSQP2N) IMPLICIT REAL*8(A-H,0-Z) DIMENSION H(NS0),X(NSQ),A1(NSQP2N) C C this SUBROUTINE CONTINUES THE GAUSSIAN ELIMINATION PROCESS C BY CALCULATING THE RECURSION FORMULA FOR THE RIGHT HAND SIDES; C I.E. IT CALCULATES P(M)=X(M) FOR M«1,...,N BY THE FOLLOWING C FORMULA: C C P(M)=X(M)=(X(M) -PW*X ( M-l ) ) / ( PO+W+PW*H ( M-l ) ) C WHERE C P(l)=X DIMENSION S(NSQ),X(NSQ>,A1(NSQP2N) C C CALCULATE THE PRODUCT ( S- (H-W*I ) *X ) • COMPUTATION OF THE COMPONENT C IS DONE IN LEFT-TO-RIGHT AND DOWN-TO-UP ORDER TO TAKE ADVANTAGE C OF THE VANISHING OF COMPONENTS OF H ALONG THE LEFTMOST AND C RIGHTMOST VERTICAL LINES OF THE GRID. C PW,PE»AND PO DEFINE THE LEFT DIAGONAL, RIGHT DIAGONAL* AND MAIN C DIAGONAL < WEST , EAST, OR IGEN ) OF H. C M = l NTOP=NSQMN+l DO 2 l=l,NTOP,N C 1 FOLLOWS THE LEFTMOST VERTICAL LINE OF THE GRID. FOR EACH I THE C COMPONENTS OF THE PRODUCT ARE COMPUTED ALONG THE CORRESPONDING C HORIZONTAL LINE. C THE FIRST COMPUTATION FOR EACH I, AND THE LAST ARE SPECIAL C CASES OCCURRING WHEN PW AND PE ARE ZERO. THIS FIRST SECTION C PERFORMS THE COMPUTATION AT THE FIRST POINT ON THE GIVEN C HORIZONTAL LINE. c i i S THE x INDEX AND IM IS THE Al INDEX. XI SAVES X(I) TO USE C AGAIN IM=I+M X1=X(I ) PW = A1 (IM-D+AH IM) PE=A1( IM)+A1( IM+1) PO=-(PE+PW) C*****THIS CHANGE IS REQUIRED FOR THE NEUMANN PROBLEM. X(I )=S(I )-( (PO-W)*X( I )+(PE+PE)*X(I+l ) ) C THIS SECTION PERFORMS THE COMPUTATION ALONG THE INTERMEDIATE C POINTS ON THE GIVEN HORIZONTAL LINE. 11=1+1 I2=I+N-2 DO 1 J=I1,I2 C J IS THE X INDEX AND JM IS THE Al INDEX. X2 SAVES X(J) TO USE C AGAIN. JM=J+M X2 = X(J ) PW = PE PE=A1( JM)+A1( JM+1 ) PO=-(PE+PW) X(J )=S(J)-(PW*X1+(P0-W)*X( J)+PE*X(J+1 ) ) X1 = X2 1 CONTINUE C THIS SECTION PERFORMS THE COMPUTATION AT THE LAST POINT OF THE C GIVEN HORIZONTAL LINE. L= I 2 + 1 LM=L+M C L IS THE X INDEX AND LM IS THE Al INDEX. PW = PE PE=A1 (LMJ+A1 (LM+1 ) PO=-(PE+PW) ^♦♦♦♦♦this change is required for the neumann problem. /'l)=s(l)-(( pw+pw)*x1+(p0-w)*x(l) ) m = m + ;> 2 CHNTINUF- 91 RETURN END 92 SUBROUTINE HH2 < N, A2 , NSOMN , W ,H , NSQ » NSQP2N ) IMPLICIT REAL*8( A-H»0-Z> DIMENSION H(NSQMN)»A2(NSQP2N) C C THIS SUBROUTINE SOLVES (V+W*I)X«S BY GAUSSIAN C ELIMINATION. SINCE V IS A TRI-DUGONAL MATRIX CONSISTING OF A C MAIN DIAGONAL AND TWO OFF DIAGONALS AT A DISTANCE N FROM THE C MAIN DIAGONALt GAUSSIAN ELIMINATION CAN BE PERFORMED DIRECTLY C ON THE MATRIX, C FIRST THE FORWARD ELIMINATION PROCESS IS PERFORMED C TRANSFORMING THE TRI-DIAGONAL MATRIX INTO AN UPPER DIAGONAL MATRIX C (SUBROUTINE HH2 ) THEN (SUBROUTINE PP2) THE RIGHT HAND SIDE VECTOR c IS TRANSFORMED IN ACCORDANCE WITH THE UPPER DIAGONAL MATRIX; I.E. C THE P(I)=X(I) ARE CALCULATED. THEN THE BACKWARD SUBSTITUTION C (SUBROUTINE PP1 ) IS PERFORMED TO GET THE RESULT. C ALONG THE FIRST HORIZONTAL LINE OF THE GRID THE VALUES OF C THE UPPER DIAGONAL COMPONENTf H, DO NOT DEPEND ON ANY PREVIOUS C VALUES OF H. PS»PN AND PO DEFINE THE LOWER DIAGONALt UPPER C DIAGONAL AND MAIN DIAGONAL ELEMENTS OF V. C THE COMPUTATION IS DONE IN A DOWN-TO-UP AND LEFT-TO-RIGHT C ORDER TO TAKE ADVANTAGE OF THE VANISHING OF COMPONENTS OF V ALONG C THE BOTTOM AND TOP LINES OF THE GRID. C DO 2 1=1, N C i FOLLOWS THE LOWER HORIZONTAL LINE. FOR EACH I THE VALUES OF H C FOR THE TRANSFORMED UPPER DIAGONAL MATRIX ARE COMPUTED ALONG THE C CORRESPONDING VERTICAL LINE. I MOVES THE COMPUTATIONS FROM LEFT TO C RIGHT ALONG THE (BOTTOM) HORIZONTAL LINE. I IS THE H INDEX AND C IPN IS THE A2 INDEX. IPN=I+N PS=A2(I )+A2( IPN) PN=A2(IPN)+A2( IPN+N) PO=-(PS+PN) C*****THIS CHANGE IS REOUIRED FOR THE NEUMANN PROBLEM. H(I )=-(PN+PN)/(PO+W) C J STEPS UP EACH VERTICAL LINE IN ORDER FROM LEFT TO RIGHT. J IS C THE H INDEX AND JPN IS THE A2 INDEX. L1=NSQ-N-N+I DO 1 J=IPN,LltN PS = PN JPN=J+N PN=A2( JPN)+A2( JPN+N) PO=-(PS+PN) H( J )=-PN/(PO+W+PS*H( J-N) ) 1 CONTINUE 2 CONTINUE C THE VALUE OF H FOR THE TRANSFORMED UPPER DIAGONAL MATRIX IS ZERO C ALONG THE TOP HORIZONTAL LINE OF THE GRID SO IT NEED NOT BE C COMPUTED AND STORED. RETURN END 93 SUBROUTINE PP2 ( N, A2 , NSQMN, W»H ,NSQ » X , NSQP2N ) IMPLICIT REAL*8( A-H»0-Z) DIMENSION H( NSQMN) , X ( NSQ > * A2 < NSQP2N ) C c THIS SUBROUTINE CONTINUES THE GAUSSIAN ELIMINATION BY CALCULATING c THE RECURSION FORMULA FOR THE RIGHT HAND SIDES; I.E. IT CALCULATE c THE RIGHT HAND SIDES P(M)*X(M) BY THE FOLLOWING FORMULA: c PIM)xX(M) FOR M*lt...»N C P(M)»X(M)= IS CALCULATED AT c THE FIRST POINT AND THEN THE REMAINING POINTS IN ORDER TO TAKE C ADVANTAGE OF THE FACT THAT P(M) DOES NOT DEPEND ON A PREVIOUS H C AT THE FIRST POINT OF ANY VERTICAL LINE. C DO 2 1=1, N C 1 FOLLOWS THE LOWER HORIZONTAL LINE. FOR EACH I THE VALUES OF c p-x ARE COMPUTED ALONG THE CORRESPONDING VERTICAL LINE. I MOVES C THE COMPUTATIONS FROM LEFT TO RIGHT ALONG THE (BOTTOM) HORIZONTAL C LINE. I IS THE INDEX OF P (P=X) AND IPN IS THE A2 INDEX. IPN=I+N PS=A2( I )+A2( IPN) PN=A2( IPN)+A2( IPN+N) PO=-(PS+PN) X(I )=X( I )/(PO+W) C J STEPS UP EACH VERTICAL LINE IN ORDER FROM LEFT TO RIGHT. J IS C THE X INDEX AND JPN IS THE INDEX OF A2. L1=NSQ-N+I DO 1 J=IPN»L1»N PS = PN JPN=J+N PN=A2(JPN)+A2( JPN+N) PO=-(PS+PN) C*****THIS CHANGE IS REQUIRED FOR THE NEUMANN PROBLEM. IFIJ.E0.L1 )GO TO 4 X(J)=(X( J)-PS*X( J-N) )/(PO+W+PS*H( J-N) ) GO TO 1 4 X 0.101195940- 0.722432800-02 JD.6C799540D-03 0.12 82 71540-03 0.498937000-04 0.936444 720-05 0.332202750-05 ■08 08 09 Q9_ 10 0.2183950 3D-11 0.12C555250 JLa_£_Q l.C 8 3 63 ■ 0.102915060- 0.264337650 0.355649 80 0- 0.117 20 26 10- 0.717648800- -0 _.3_4 43 2_3 3 D- 0.913163020- .389 883400- 0.107633010- ^.JJ^_9287_8_0r 101907850- _1116A06 2D- 1096 36810- _1Q4 341040- _0 _0_ 0.807247560- JL._S1A5_9_431I>: 11 12 12 13_ 14 i_4_ 15 15 16 16_ 16 I6_ 0.52734296D-06 JLl2_58151_310-06 0. 805263900-07 -P. • JL1 03 3 5 2 6 - 7_ 0.504483430-08 0.1718 19690 -08 0. 302390000-09" _0. 132522740-09 0.254808720-10 _0_.„7619 5 591D-11 0.27822733D-11 0.136527380-11 0.212719510-12 _0_. 95897955 0-13 0.355086330-13 0.21699 89 30-13 0.180295930-13 _0. 178457760-13 0. 112499390- _C_._78_720_9 5 00- 16 l_6 16 16 17 17_ 16 17 0.174681090-13 _0_._i7193l310-13 0.171073140-13 0_. 16446736 D- 13 0.164683530-13 0.164271250-13 0.16 250 3730-13 0_._158 757030-1 3 Results of Stone Algorithm on Generalized Model Problem with Dirichlet Boundary Conditions 98 PAX IMiJM NHRMAI T/FD RESIDUAL LI VFCT UR NO RM L2 NO RM RATIO I 0.797638270-03 0.132027410 02 0.79557158D 00 — 2 0*21 9 C 98 8 8D_r.O 3 0, 16331 76 3D_0 1 Q. 6 76 7 7 38 D 00_ 3 0.140652440-03 0.564218360 01 0.578310120 00 — 4 _X2.06.661 50-0 3__* 0.45Q3108 90 01 _ 0. 4922 7654 0_QlO_ 5 0.10169786D-03 0.359253310 01 0.405728460 00 — 6 0. 84191 32 40 - OA Q-__e.8ji___.L3 2 01 0.330909970 pp 7 0.677108250-04 0.230288240 01 0.26614257D 00 8 _U3133 79 890-04 0. 181339620„ Ol_„ 0*2 10604830 Q0_ 9 0.33077995D-04 0. 11 1 14725D 01 0.130136180 00 _ 1__ 0_. 201 132A20j^0A Q._66786606_0_ 00 0^786700.9 7 D-iOl 11 0.119375860-04 0.391223090 00 0.463211190-01 - 12 -L^__-_i_66 320-Q5 -_-£____A_6840 0. 2 6447227D-0 ] 13 0.381611070-05 0.121712130 00 0.145383 110-01 -14 _U______£l^l_5J_rJ_5 0.6 3930894LH0 1 _0. 766663370-02 15 0.103131300-05 0.320013560-01 0.385200880-02. ^ 6 •CljlA.S33.52_62 0ji.06. 0,15159.0400-01 0. 1831241 90-02_ 17 0.12128716D-06 0.353128240-02 0.436659050-03 -^ Q.?75 .9_33A.6J0_r__J 0. 789096120-03 0. 9691 5161 D-04 19 0. 576003310-08 0.159898130-03 0.197599490-04 -2& 0.108395390-08 Q_. 294 3 893.1 0-04 „ 0. 365741 11 D-05 21 0.182199790-09 0.481213610-05 0.600641480-06 -22 Cl*_26_69_422 3^-L0 . C__691 i_3433D--g6_ _0 x 866l53860-07_ 23 0.335122460-11 0.853194880-07 0.107299830-07 -2A Q_-.3517 5 573JU ZL L2 0_ y _3_82681 48L)-08 0.111347660-08 25 0.28986685D-13 0.72 3561600-09 0.908414150-10 -2& 0.22C26265D-14 O.A973 36520-10 _0. 578305550-1 1 27 0.149838090-14 0.162947840-10 6.2 74531 710- I 2~ 28 0^1.15762730- .14. .JA3JI 2 3 5 0- 1 0_ 0. 120402 32D-l_3_ 29 0.103671200-14 0. 1 53913900-10 0.135199540- 13 -^ Q.1338^650Jl_J-4 0.168981130-10 0. 18434904D- 1 3 31 0.153786330-14 0.169494620-10 0.171655260-13 32 0.142934340-14 __0. 160742210- 10 0.133213100-1 3__ Results of ADI Algorithm on Generalized Model Problem with Dirichlet Boundary Conditions 99 ITERATION MAXIMUM NORMALIZED RESIDUAL 0.675675680-01 L2 NORM RATIO 0.0 ITERATION ITERATION ITERATION I TERATION I ? 3 4 5 6 0. 943201420-02 0.21729C41D-02 0.43241693D-03 0.88390354D-C4 0.58071917D 00 0.384511640 00 0.33300406D 00 0.28885845D OC ITERATION ITERATION C.81952496D-04 0.76944366D-04 0.272475510 00 0.25716231D 00 ITERATION ITERATION 7 8 9 10 11 12 0.65054589D-04 0.57914335D-04 0.195799160 00 0.149466790 00 ITERATION ITERATION 0.51634272D-C4 - C. 408757680-04 0.13478101D 00 0. 12170359D 00 ITERATION ITERATION C. 416565630-04 0.37906C480-04 0.116314710 00 0.111237850 OC ITERATION ITERATION 13 14 15 16 17 18 C. 327552990-04 0.276181640-04 0.915899870-01 0.75465564D-01 ITERATION ITERATION 0.261786300-04 C.22081712D-04 0.698832840-01 0. 647777970-01 ITERATION ITERATION 0.226677680-04 0.211568260-04 0.624851140-01 0.60303519D-01 ITERATION ITERATION 19 20 C. 154642440-04 C.14249C730-C4 0.413547370-01 0.2860C5790-01 ITERATION ITERATION 21 22 2 3 24 0.114200C6D-04 0.871719260-05 0.24833732D-01 0.216134940-01 I TERATION ITFRATION 0. 791038830-05 0.674848950-05 0.2038 402 6D-0 I 0.19245750D-01 ITERATION ITERAT ION ITERATION ITERATION 25 26 27 23 29 30 0.552276400-05 0.479452C30-05 0.41538816D-C5 0.325878950-C5 0.328332930-05 0.296451740-05 0.147042170-01 . 0.112552510-01 0.101484940-01 0.916555210-02 ITERA T ION ITERATION 0.87539842D-02 0.837677770-02 ITERATION ITERATION 31 32 0. 252968720-05 0.212497900-05 0.689941040-02 0.563688060-02 Results of Stone Algorithm on Heterogeneous Model with Homogeneous Subregions with Dirichlet Boundary Conditions 100 : MAXIMUM N CRMALIZEQ RESIDUAL I 0.494798920-02 -2 Q_, 2843356 1 0-02 3 0.195144450-02 .4 0_ a 16 62 2 63 6_Dr_0_2_j 5 0.14742072D-02 -A 0.133137820-07 LI VECTOR NORM 12 NORM RATIO 9 JLQ_ 11 12_ 13 _L4_ 15 _L6_ 17 JLfl. 0. 120948990-02 ._£)__. .1 098 544 2 0- 2 0.995361680-03 _ JO ,.905 090820 -03 0.111484520-02 C.139C733Q0-Q2 C.20C499770-02 0.2 81 804 6 10-0 2 0.389338490-02 XL. 5 l27 2606 50-0 2 0.7017 05990-02 XL^_ 9 13630 1 5 0-02 19 _2.0_ 21 22 0.113329830-01 Q._15C3 72_90O-0l 0. 188353290-01 ,0.2 32626260-01 23 _2_4_ 0.28382900D-01 0.342741270-01 25 26 27 23 29 0.410170270-01 0.486 7463 00_-_0_l_ 0. 572637160-01 0. 66 75266 70-31 0.769726740-01 0.876125280-01 0.758910990 _0_. 57833 31.60 0.471624960 0.3 9746038D 0.34C83210D 0.295145260 0.258891210 02 .0.22 7617060 02 0.202 593350 02 _0_._18_400 377 02 0.1727116 20 02 0. 171604630 02 02 0.95094974D 00 .02 0_. 920531390 00 02 0.89465124D 00 02 0.8 700 3 394 00 02 0. 845064 C70 00 02 0. 818631020 00 0. 790042 760 00 0. 758386110 00 0.72295505D 00 0.682953190 00_ 0.63 7549090 00 0.535354460 00 0.18139187D 02 JL.202918170 02 0.236487040 02 _0. 2 81 2 568 10 02 0.339072860 02 0.413396430 02 0.504574450 02 _0.61240J340 02 0.735960170 02 0_.JT75_672_16O 02 6.103362 010 3 0. 121357690 03 0.527115090 00 0.461026670 00 0.383303870 00 0.3U32799D 00 0.234642720 00 0.165285670 00 0.11338167D 00 0.906478230-01 0. 95895432 D-01 0.113267370 00 0. 134655730 00 0. 15953446D 00 0.142056750 03 JD. 165944500 03 0.193189530 03 _0. 223804360 03 0.2 569 86350 3 0.290985500 03 0. 188666050 0.222440490 0.260939320 0. 303912070 31 32 0.98 13 13470-0 1 0.107725950 00 0.322851940 03 0.348697680 03 00 00 00 00 00 00 0.445377730 00 0.485431830 00 0.350490380 0.398764720 suits of ADI Algorithm on Heterogeneous Model with Homogeneous Subregions th Dirichlet Boundary Conditions 101 MAXIMUM NORMALIZED RESIDUAL L2 NCRN RATIO ITERATION . 67 1 76877Q-C 1 Q.Q ITERATION 1 C.e92877130-C2 0.580917300 00~ -IHRATjCN 2 C.15A21A 550-C2 0.38 75189<5D 00 ITERATION 3 C .25728335D-C3 ~07T359T5TlD~ 6d" _1_TERATI0N A 0. 83 0A70 15D-0A 0.2917921 CD OC ITERATION 5 .820 273700-CA " T72T53l5'35D~0C" ITERATION 6 C .7776 10A Q D-CA 0.26 C0AQ700 00 ITERATION 7 0.659789600^ 0TT9 827108C 00 ITERATION 8 0.574627 970-C4 0.15166685D CO ITERATION 9 C. A9 6C6 8 6 6D-CA ~07T369Tr2C'D~00" JLI E R AJUO N_ jm C.A1A92 CJB6 D-CA 0,1 21725. 8 6D 00 ITERATION 11 C.A161857?D-CA " T 0riTS3C5C90~00" ITERATION 12 0_.378A627CD-CA 0.113183130 OC IT ERA TICK n . 3 32836580-CA 0.9330e597~f>oT~ ITE RATION 1A Q . 2 7 9C956 C D^C A 0.77C02A71D-01 ITERATION 15 C.26C13916D-CA ~CT71 36 3~1 520-0 V JJER ATI CN 16 C 226C 8 A5 6D-C A g .661897580-01 ITERATION 17 C . 2288AC520- CA ~0~ 63 37 18 C 6 HJT~ ITERATION 18 C . 2 10C2 1 79D-CA .6 1660336D-Q1 ITERATION 19 C. 1 58C325CD-CA 0. A 2356229D-C I JTERATICN 20 C . 1A070 1 78D -CA . 293A5259D-C1 ITERATION 21 . 1082306 1D-0A 07255Tsl9T6^0T" ITERATICN 22 .8A6672860-C5 0.2223 1700D-Q1 ITERATION 23 C . 78 1 50AC30-C5 ~T720V8 _ 1387C^0T~ ITERATION 2A C .685 182 79C-C5 0. 198 18A96D-C 1 ITERATION 25 C . 5668C A 360-C5 0. 1 516AC86C-C1 JJER ATJT N_ 2b . A839 8599D-C5_ _ 0.116277A50-01 ITERATION 27 C .A05201 37C-C5 " ^Tl CVsTC SW- CT~ U ERA_11C_N_2_8 C .33255 9^ 2D- C 5 .9A8666060-C? ITERATION 29 .33381 8A2D-C5 ~^0T9C7b~A6ClC-02~ ITERATION 30 . 298893 33D-C5 0.867 776750-02 ITERATION 31 0.2628559C0^C5 0. 7156A31AD-02 ITERATION 32 C .2 18A8 1 7 1D-C5 .59C76955D-02 Results of Stone Algorithm on Heterogeneous Model with Random No Sub- regions with Dirichlet Boundary Conditions 102 5 6_ —MA.X1H UM NOR MALIZE: J RESIDUAL 1 0.564897290-02 _2 0_* 29 5 4 .4 3_6 5 0- .2_ 3 0.225026750-02 , 1_9£JL_Q 5000- C2_ ; 0. 180U133D-02 0.16 5742 7JLQr£L2 LI VECTOR NORM L2 NORM RAT| n 0. 693666390 Q , 5 327 79 7a D 0.434 A 8336 _0j 3669 700 ID 0. 3 1 58 31 78 D 0.274832070 9 .JO. 11 _L2_ 0. 1526 30330-02 _C-._i3 t 970702D-02_ 0.126531500-C2 0.115776280-02 0. 12 33 77210-02 -Q...13123 082 0-02 0.240801760 _0_,_212504 370 0.189578620 0.17539019D 0.170659230 13 15 _16_ 17 _L8_ .179791520 0.185037900-02 _ _Q_. 25 77 3939 0-0 2_ 0.350730460-02 _J). . 4 6 3 V. 18 6 0-0 2 0.6197 30880-02 0.8 08013300- 19 0.104750090-01 _ 2 _0_. 1319 930 70-0 1 21 0.163947000-01 JL2^„ _0.210 166770-01 23 0.258303940-01 _2it OJ. 14 1133 50-01 25 0. 37334343 J-OL 26 0.451510100-01 27 C. 533607850-01 2d 0^.623 735920-0 1 29 0.719957750-01 _1Q O.Bia345750-0 1 0.204563040 02 0. 244791330 02 0.302629430 02 _0.3„802 063 50_02 0.477167450 02 0.593006450 0? 0.728576570 0.381267470 0.104851150 .0.122891130 0.142281030 0.163263320 0.186038400 0.211098400 0.238659410 0.268506630 0.299893300 0.331496220 02 0.95195615D 00 ?___ D j. 9 21972520 00 02 6.896417590 00 02 0.87208958D_00_ 02 0.84 74 11140 00 .Qi .821 35 4290 00 02 0.79311204D 00 02_0,761960 39D 00 02 0.727186630 00 0_2 _0_. 6 88046970^00. 02 0.643 74633D 00 02 0. 593458430 00 0.53643826D 00 0._4_L23Q4710_00 0.40152653 00 0.326007650 00 0.249570580 00 0. 17828852D 00 02 0.121090990 00 02 0.9014 5 140 0-01 03 0.90285282D-0r 03 0.106_98389jD J)0_ 03 0.128591640 00 03 0.1 531 3 3350 00 03 0. 181304910 00 03 0.213611390 00_ 03 0.2 5010 8560 00 03_ 0^2903 9578 00 03 0.333408030 00 03 0.377115530 00 0.913102260-01 0,996320100-01 0.361179630 03 0.386320870 03 0.418308430 00 0.4527218 70 00 Results of ADI Algorithm on Heterogeneous Model with Random No. Subregions with Dirichlet Boundary Conditions 103 MAXIMUM NORMALIZED RESIDUAL L2 NORM RATIO _TTERATION 0.900000000 00 0.0 I TERATION I 0.36565756D CO 0.9996 71210 0~0~ ITERATION 2 0. 8489971 7D-CU 0.99887916D 00 ITERATION 3 0.165282650-01 I)7997Tr6~3 8 D~0"0" ITERATION 4 0. 1 77350 75D-01 0.99634342D 00 ITERATION 5 .37 3 1 72 I A I'D- 02 0T9957 C3~5~2D~ 00" ITERATION 6 .309260560-02 0. 9950782 9D 00 ITERATION 7 0. 809555210-02 0.992611170 00 ITERATION 8 . 933506220-C2 0.990061UD 00 ITERATION 9 6.242126280-02 "0Y9890 375lD~00" ITEPATION 1 CU1 18377 890- C 2 .9880 249 7D OC ITERATION 11 0. 824298 83D-C3 ~0."W7 579 70 7J~00 ITERATION 1 ? C. 781409880-03 0.98713668D 00 ITERATION 13 0. 1 31 07080D-02 0.98525454D 00 ITERATION 14 0.210 179240-02 0.983371950 00 ITERATIONS . 71 0662300- 03 - ^0. 98262 9 5 8 D~0 6~ ITERATION 16 0^646004690-03 0. 9818S087D 00 I TF RAT I ON 17 0.521390310-03 ~0~. 98T543 8 9 D~~00 ~ ITERATION 18 .5 1 1375 1 8D-03 0.981 19 8010 00 ITERATION 19 0.253U6880-C2 0.977439130 00 .ITERATION 20 0.39461 250D-02 0.97367703D 00 ITEPATION 21 0.13340077D-C2 ~0^9T232ib%D~QW J_T1. R . A JI0N. 22 0.59908996D-03 0.970971050 00 ITERATION 23 . 36484634D-C3 79T04T8 5 3 D~OlT I TERATION 24 0. 343600 26 0-03 0.96986 810D 00 ITERATION 25 0.95862 7690-0 3 0.967208210 00~ ITERATION 26 0. 15J35P86D-C2 0.96456122D 00 ITERATION 27 0.454951 340-03 " ~^796T5T377D~0C~ ITERATION 23 0.348 261640^0 3 0.96260936D 00 ITERATION 29 0.259921700-03 OT96Tl"84860 _ 00 _ ITERATION 30 0.255745800-03 0.961761360 00 ITERATION 31 0.394828340-03 0.959895600 00 ITERATION 3J 0.686508220-C3 0.958035180 00 Results of Stone Algorithm on Model Problem with Neumann Boundary Conditions V 104 —MAX I MUM-NO RMAUZ ED -RESXDUAL. 1 0.910332430-01 -2 n.flflr>Q4RARp-oi 3 5 -6 7 _8_ 0.851564820-01 --{L.82218099D-0L 0.792797160-01 _a.76 34133 3DrJ)l_ 0.734029500-01 n .7n464 S 67n-ni 9 11 _12 13 14 15 _16 17 -18 19 _2D_ 21 -22 23 -24 25 _26_ 27 -2B 29 30 31 -32- 0.675261840-01 J).6A5 873Q1D-QX 0.616494180-01 _0.58 7_11035D-rO.L_ 0.55772652D-01 0.57 fi34. AQn-m 0.498958860-01 ..£. . 4 6 9 5 75 3 D-XU 0. 440191200-01 -0 . 41 8 073 70-01 0.381423540-01 -0*352-0 39 7 1D_-Q_1 0.32265588D-01 _ 0.29327205D-01 _ 0.26388822D-01 0.23450439D-01_. 0.205120560-01 _0 *X7_5 736-7 aO^QX_ V FC T DR-NORM L 2_NORil_RiLrJil_ 0.98290694D 00 n.97sr>^nt;7n qq 0.819299230 0. 797853780 02 _02_ 0.766408340 02 -0..739 96 2 89D 0? 0.713517440 JL.f> 870.72000. 0.66062655D O.ft^i fli inn 0.96661858D Jl. 35.85 6.51 20 02 00 _QQ_ 0.950962010 00 0.607735660 JL._58129_.2LD 0. 554844 76D J). 5 2 83 99 3 ID 0.501953870 0_475'.034?n 0.44906297D JL. 4226 17530 -02 Q__S4_245 OJLLD _Q0_ 02 0.93385958D 00 J12 O.q. «_474.sn on 02 0.917734690 00 02 a. 90_9155_7_&Q; 00_ 02 0.901230150 00 02 . 5.9 34 2 93DB_„Q0__ 02 0.884943900 00 -02 0.877164030 00 02 .02. 02 0.39617208D J>.36_972663D_J)2 0.343281190 02 0.31 6R^S74 n 07 0.146352900-01 0. 1169690 7D-01 0.875852400-02 0.582014100-02 0.28817580D-02 J) ._ 5 6 62 5 0.0-9B. _04_ 0.290390290 02 _JL.263944_4D 02__ 0.23749940D 02 -Q_2L1Q5395D_Q2 0.184608500 02 JU15ai6306n 07 0.13171761D -Q. 105272160 0.78826716D 0.86887236D 00 J),. .86 0X48 8 10 00 0.852514520 00 JL. 8443961 4.0— 0_0_ 0.83658839D 00 Q.8?B0?633n 00 -0.523812690 0.259358220 02 02 01 01 01 -Q .5096?50RO-OJ 0.81971369D 00 JK81173325D J)Q__ 0.80353219D 00 JL__79 542 79__6D.___>Q__ 0.78761281D 00 0.7797364RO 00 0.771831610 00 0.76381693D 00- 0.756081540 00 0.74754817D 00_. 0. 739924650 00 0.73 i60f.77n nn Results of ADI Algorithm on Model Problem with Neumann Boundary Conditions 105 ITFRATH3N (. MAXIMUM NORMALIZED RESIOUAL ITERATION IT ERATI UN ITERATION ITERA TION ITERATIUN ITERATION 3 4 ITERATION 7 JTER AT. I ON 8. ITERATION 9 J.TEJLAHQN__LO- ITERATION 11 ITHRAXIiiLL_12_ ITER AT ION XXERAXION. 13 J_4_ ITERATION 15 XIERAXX_]N_16_ ITERATION 1/ XTFRAT IG N 13 ITERATIUN 19 XTJTRAXION. _20_ ITERATION 21 XX£RAXJUN_22_ ITERATION 23 ITFRAT inN__2it_ I TERATION _JJJ_RAX10N_ I TERATIUN _lXERA.TJ.Ufi.. ITERATION J_Xi_RATIO___ 25 2.6. 27 2_8_ 29 _____ ITERATION £TERAI_1l_N. 31 _3_L 0.900000000 00 L2 NORM 0.0 RATIO 0.107392520 00 C. 204319370 00_ 0.107602810 00 C.Z7 02MJL_5 0- 1_ 0.424807320-01 0.222449520-02 0. 997978260 00 0_. 9844078 70 00 0.964534080 00 J_L_93772 33_>D__00 0.93452406D 00~ 0.93169797D 00 0.161115760-01 C. L9684_9130-0l 0. 319523930-01 _-_-_-2612_J6_5_00- 2 0.442965420-02 0.2 91172130- :)3 0.197460260-03 _a_-lJ_5____ajU.:i^_22 0. 125540070-02 _Q _XS_<__i_3J_L3_C______3_3 0.62 910435 0-04 ___ . 281 142720-04 0. 835252460-02 J-_.J9_04O_6i3_9 00^.02. 0. 135616490-01 -0_^3_3l5376 7.50-02_ 0.525450290-02 Q.3d3480S2n-0^ C.16 7118890-03 _Q.-J3AXL613_L0-_O3_ 0.835947700-03 JJl-..3J3_a_5_5_854ir0_3. 0.922035060-04 Q. 75240 30 70-04 0.300983690-03 _0._22J.8L 20 .50- CJ3_ 0.376993090 00 0_^8_2 74 365 30__00_ 0.821257310 00 J>i-81A1_3 865 0_00_ 0.814826800 00 0.813532650 00 0.795536830 00 JO .7 7 744V1 *> qo 0.774758630 00 .0.772 1 5 15 0__0 0_ 0.77l099 7o0 00 0.770093940 00 0.706515680 00 -__ JO- 643 83 2.530 00 0.635508810 00 ._6 2 8_0 2 5 2 3 Q_ 00_ 0.626689300 00 -_--JL__L2^3_52980_00_ 0.597814320 00 0.5 72143 9 50 __0 0_ 0.568614610 00 _0_ .__56 50_ 83 3 1 0_ 6 0.564206370 00 Jl__JL6JL3J_68J3_0_fi0_ 0.549965840 00 .0.53726 3190 00 Results of Stone Algorithm on Generalized Model Problem with Neumann Boundary Conditions 106 --HA.XIHUM-NORMAL1 ZED-RESIDUAL I 1 VFf.TnR NDRM 1 0.817678490 00 . -2 0^zai3-a535a_oxD 5 9 -4 11 _12 13 14 15 -16 17 -18^ 19 20 0.76509221D 00 -XK73&7.9-9-07D -QO 0.712505930 00 -0.686212790-00. 0.65991965D 00 -0.6336265in 00 0.607333370 00 -0.581O4023D-QQ- 0.554747090 00 -0.52845395D-QQ 0.502160810 00 n. 47s q < s 7^7n on 0.449574530 00 ,0. 42 328 1390 -00 _ 0.396988250 00 -0.3.706951 ID. jO.O_ 0.344401970 00 _a^3-L8-L0 8830 00— 21 -22 23 -24 25 -2A- 27 -28^ 29 -30 31 -32- 0.29181569D 00 — Q.2 65522 5 5D-0O- 0.23922941D 00 0.212936270 00 0.186643130 00 .16034 939Q. 00 0.134056850 00 0.107763710 00 0.814705710-01 0.551774310-01 0.28884291D-01 0^-259-L L5 100-02 0.735910640 03 0.71??46ft?n 03 _12_N0RMJlAIJn. 0.989806940 00 n.QflflR?ns7n nn 0.68858299D 03 -0^64919-160-03 0.641255340 03 -0 .i> 175 9151 D_Q3_ 0.593927690 03 0. 570763860 03 0.98731858D 00 __ 0..9H61651 20 _QQ__ 0. 985462010 00 _ 0, 9fl.335jQ3_LD_QQ_ 0.982159580 00 o.QR n f>74?sn nn 0.54660003D 03 -0. 5229 36 210-03- 0. 499272380 03 -0L..4 756.08 5 6D_03.. 0.451944730 03 n.4?R7ftOQnn ni 0.979834690 00 _Q . 97B155_7£D- Q Q_ 0.977130150 00 _ .976229300 _Q0__ 0.974643900 00 0.973764030 00 0.404617080 03 _0. 380 95 325D__03- 0.357289430 03 0. 333.625600- 0.309961770 03 0.28629 795 03 0.972372360 00 -0 J .97_1148 8_10__0O_ 0.969814520 00 ..9.6.B 5.9.6 14D_Q0_ 0.967688390 00 n.Qf>fen?^ft3n nn 0.262634120 _Q.2389_7Q3.QD 0.215306470 .Q_.1_9164264D 0.167978820 0.14431499CL 03 0.964613690 00 3 . 9635 33 25D_ -QO. 03 0.962232190 00 03 -0.961027960 0Q- 03 0.960112810 00 .03, 0.9591364 8O-O0_ 0.120651170 0. 9698734 OD 0.733235140 J}. 4965 968 80 0.259958620 0.23 3203 59O. 03 0.958131610 00 02 _ 0.957016930^00. 02 0.95618154D 00 02 J). 9545481 7D 00 02 0.953824650 00 -OJ n.Q5740677n no Results of ADI Algorithm on Generalized Model Problem with Neumann Boundary Conditions 107 MAXIMUM NORMALIZEO RESIOUAL L2 NORM RATIO ITERATION o 0.90C0OOOOD 00 0.0 ITERATION ITERATION I 2 0.369512080 00 0.101449410 00 0.999659990 0.998884560 0.997601570 0.996144170 00 00 ITERATION ITERATION 3 4 5 6 0.211599460-01 0*250774230-01 00 00 ITERATION ITERATION 0.434708430-02 0.333185050-02 0.99540619D 0.994735040 00 00 ITERATION ITERATION 7 8 0.108194810-01 0.120142660-01 0.991891070 0.938906630 00 00 ITERATION ITERATION 9 10 0.294619460-02 0.134584380-02 0.987748060 0.986590460 00 00 ITERATION ITERATION 11 12 C.90835755D-03 0.826936230-03 0.986097510 0.9856 00660 00 00 I TERATIUN ITERATIUN 13 14 0.171005180-02 0.270927790-02 0.983456980 0.98127579D 00 00 ITERATION ITERATION 15 16 0.773537480-03 0.689625840-03 0.980450240 0.979615680 00 00 . ITERATION ITERATION 17 18 0.546359090-03 0.54744063C-03 0.979236300 0.978852270 0.974537240 0.97015980D 00 00 ITERATION ITERATION 19 20 0.321482210-02 0.473847870-02 00 00 I TERATIUN ITERATION 21 22 0.153938920-02 0.777612370-03 0.96864095D 0.967109430 00 00 ITERATION ITERATION 23 24 0.385158940-03 0.365654 740-03 0.966503410 0.96589334D 00 00 ITERATION ITERATION 25 26 27 23 0.117428000-02 0.18269337D-02 0.514232440-03 0.373301670-03 0.962871210 0.959833770 00 00 ITERATION ITERATION 0.95874951D 0.957660960 00 00 ITERATION ITERATION.. 29 30 0.269154350-03 0.268811860-03 0.957198530 0.956734630 00 . 00. ITERATIUN ITERATION 31 32 C.493722280-C3 0.654752350-03 0.95463500D 0.95252532D 00 00 Results of Stone Algorithm on Heterogeneous Model with Homogeneous Sub- regions with Neumann Boundary Conditions 108 MXIUJM NCRMALIZEO RESIOUAL Li VECTOR Ni 3RM 03 03 03 03 03 03 L2 NORM RA 0.996260690 0.995443060 0.994571860 0.993736510 0.992946200 0.992065030 no 1 2 0.846730210 C. 819488790 CO CO 0.76205719U 0.737539910 00 00 3 4 G. 79224737^ C.765CC5950 CO 00 ■ 0.713022630 0.6885053oO 00 00 5 0*737764530 0.71C523110 00 00 0.66398808J 0.639470800 00 00 7 0.6832S169D 0.65604C27O 00 00 0.014953520 0.59043624D 03 03 0.991175960 0.99030743J 00 00 9 10 C. 626758850 C. 601557430 C. 574316010 G.547074;>9D 00 CO CO 00 0.565918970 0.541401690 03 03 0.989503470 0.988615580 00 00 11 12 0.51688441D 0.492367130 03 C3 03 03 03 03 03 03 03 03 0.987793020 0.986^82930 0.98610439U 0.985296400 0.984437240 .983594880 0.98274145D 0.981899610 0.981086840 0.980202680 00 00 13 14 0,519833170 0.492591750 00 00 0.46784985J 0.443332580 00 00 15 16 0.465350330 .0.43S1C691Q CO cc 0.418815300 0.394298020 uo 00 1/ 13 C.41C66749D 0.38362607U CO CO 0.369780 740 0.345263460 00 oc 19 2Q C. 356384650 0.3291432 30 C. 301901810 0.274660390 00 00 0.320746190 C. 296228910 00 00 21 22 00 00 0.271711630 0.247194350 03 03 03 03 0.97934137D 0.978513330 0.977663220 0.976822300 00 00 23 ?4 C. 247418970 C.22C177550 00 00 0.222677070 0.19dl59800 00 00 25 PL C. 192936130 C. 165o94710 00 CO 0.173642520 0.149125240 03 C3 0.97601128D 0.97519365D 00 00 21 2fl C. 138453290 C. 111211870 00 00 0.1246C7960 C.IC0090680 03 03 0.974373160 0.973541690 00 OC 29 30 0.8397C4510- 0.567290^10- 0.294376110- C. 224619100- 01 01 0.755734C60 0.510561200 02 02 0.972738150 0.971854820 00 00 31 32 CI C2 0.265368500 0.202157190 02 01 0.9710624&0 0.97020067J 00 00 Results of ADI Algorithm on Heterogeneous Model with Homogeneous Sub- regions with Neumann Boundary Conditions 109 ITERATION MAXIMUM N0RMALI7ED RESIDUAL 0.9OC0OO00D 00 L2 NORM RATIO 0.0 ITERATION ITERATION 1 2 0.369929740 00 0.753381780-C1 0.99965111D 0.998801280 00 00 ITERATION ITERATION 3 A 0.284090740-Cl 0.247755C20-C1 0.997489930 C. 995950180 00 00 ITERATION ITERATION 5 6 0.386844190-C2 0.31873126D-C2 0.99522704C 0.994535150 00 oc ITERATION ITERATION 7 8 0.96959258D-C2 0. 109539940-01 0.9916P776D 0.988647160 00 00 ITERATION ITERATION 9 10 0.247845270-02 • 0.12013744D-C? 0.987490540 0.98631513D 00. 00 ITERATION ITERATION 11 1? 0. 833465960-03 0.793707690-C3 0.98582263D 0.98531846D 00 00 ITERATION ITERATION 13 14 C.15899757D-C2 0.25396574C-C2 0.983161690 0.99094702D 00 00 ITERATION ITERATION 15 16 0.725349670-03 0.657191000-03 0.980115810 0.979269300 00 00 ITERATION ITERATION 17 19 0.526397610-03 0.574259650-03 0.978887410 0. 978497990 00 00 ITERATION ITERATION 19 20 0.30667442D-02 0.459923400-02 0.974112090 0.969663180 00 00 ITERATION ITERATION 21 22 0.140454230-C2 0.61596664D-03 0.968122260 0.9665690 50 00 00 ITERATION ITERATION 23 24 0. 369477230-03 0.35R09762D-03 0.965953670 0.965334590 oc 00 ITERATION ITERATION 25 26 0. 111*71 960-C? 0.175541510-02 0.962259790 0.9591763*0 00 00 00 0-0 ITERATION ITERATION 21 28 0.465407510-C^ 0.368780590-C3 0.958074620 0.95696957D ITERATION ITERATION 29 30 C. 259936450-03 0.264754P«0-C3 0.95649999D 0.956028920 00 00 ITERATION ITERATION 31 32 0.471754570-C3 0.82321863D-03 O.953992 , 530 0.95175142D oc 00 Results for Stone Algorithm on Heterogeneous Model with Random No. Subregions with Neumann Boundary Conditions 110 MXIFUN r*CR^ALlZfcC RESIOUAL LI VECTOR NORM 12 NORM RATIO 1 ? 0.76G17056G C.73577C01D CC CO 0.684153520 0.662193010 03 03 0.995768300 0.995069690 00 00 3 4 0. 711369*40 C. 666968870 00 CO ■ 0.640232490 0.613271980 03 03 0.994462480 0.993825610 00 00 5 6 0.662568300 C. 636167730 00 CO 0.596311470 0.57*350950 03 C3 0.993129220 0.992503d5D 00 00 7 a C. 613767160 0.53936659D CO 00 C. 552390440 0.530429930 C3 03 0.991826770 0.991165720 00 00 9 10 .56*966020 0.5*3565450 oc CO 0. 508469420 0.466508900 03 03 0.990495620 0.96985629Q 00 00 ' 11 12 C.51t>164S8C C. ^91 764310 CC CO 0.464548390 0.442587d80 03 03 0.989229920 0.98859931D 00 00 13 14 C. 467363740 C. 442963170 CO CO 0.4206273o0 0.396666850 03 03 0.98790922U 0.987280900 00 00 15 16 C.4165626C0 • C. 394 162030 00 00 0.376706340 0.35*745b2U 03 03 G.9866121ti0 0.985915560 00 00 17 29 -3JL 31 JlL. 0.3697614t>0 00 0.34536C890 CO C. 76<;;>4617l,-C1 C.52554047U-C1 C. 281534770-01 C .37529C690-02 0. 332/65310 03 0.31C824800 03 0.985313680 00 0.984652t>80 00 19 20 G.32C960320 C.;9655975D CO CO 0.288864290 0.2O6903770 03 03 0.983964820 0.983302^4D 00 00 21 22 C.2721591dD C. 247758610 CC 00 0.24*9*3260 0.222982750 03 C3 0.9ti26*23lD 0.981959650 00 00 23 24 C.22335d0*0 0.198957470 CO 00 0.2C102223D 0.179C61720 C3 03 0.981319300 0.980657590 00 00 25 26 0.17*5569CD 0. 15015633U 30 CO 0.157101210 0. 135140690 03 03 0.960015350 0.979320860 00 00 27 2d 0.125755760 0.1013^5190 CC CG 0.113180180 0.912196680 03 02 0.978681690 0.978003170 00 00 0.692591550 02 0.472986*20 02 0.253381290 02 0.337761620 01 0.977360100 00 0.97b71029D 00 0.976054640 00 0.97537623D 00 Results of ADI Algorithm on Heterogeneous Model with Random No. Sub- regions with Neumann Boundary Conditions n 0l s$ W?g MAY 1 7 197A