LIBRARY OF THE 
 
 UNIVERSITY OF ILLINOIS 
 
 AT URBANA-CHAMPAIGN 
 
 510 /c4 
 Cop. 2. 
 
Digitized by the Internet Archive 
 in 2013 
 
 http://archive.org/details/divergenceofston478diam 
 

it/* 
 
 JL Report No. 478 
 
 / 
 
 THE DIVERGENCE OF STONE'S FACTORIZATIONS 
 WHEN NO PARAMETERS ARE USED 
 
 by 
 
 MARTIN DIAMOND 
 
 September, 1971 
 
 1EUIBRARYOEIHI 
 
 NOV 9 1972 
 
 UNIVERSITY OF ILLINOIS 
 AT URBANA-CHAMPAIGN 
 
 DEPARTMENT OF COMPUTER SCIENCE 
 UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 
 
 URBANA, ILLINOIS 
 
Report No. 1+78 
 
 THE DIVERGENCE OF STONE'S FACTORIZATIONS 
 WHEN NO PARAMETERS ARE USED* 
 
 by 
 
 Martin Diamond 
 
 September 7, 1971 
 
 Department of Computer Science 
 University of Illinois at Urbana-Champaign 
 Urbana, Illinois 61801 
 
 *This work was supported in part by the National Science Foundation under 
 Grant No. GJ812. 
 
i. iwrroDucTioN 
 
 The solution of a system of equations AX = q by factorization 
 techniques can be separated into two segments. The first step is the 
 definition of an auxiliary matrix, B, which is chosen so that A + B = LU 
 where L and U are sparse lower and upper triangular matrices respectively. 
 In practical computer programs the elements of B are not computed; instead 
 the elements of L and U are computed. Then by using the factors L and U, 
 systems of equations whose matrix of coefficients is A+B can be easily 
 solved. 
 
 The second part of the method consists of an algorithm which 
 uses an iteration of the form 
 
 (A+B)X n+1 = (A+B)X n - (AX n -q), n = 0, 1, 2, ... 
 
 to generate a sequence of vectors {X } which converges to the solution AX = q. 
 
 The convergence of the iteration is dependent upon the definition of the 
 auxiliary matrix (through the definition of L and U) and a sequence of 
 parameters. 
 
 One of two types of parameters is employed in the iteration. The 
 first type multiplies the term (AX -q) in the iteration as follows, 
 
 (A+B)X ,_ = (A+B)X - t (AX -q). (l.l) 
 
 n+1 N ' n n n 
 
 an analysis of this parameter and an algorithm based on the iteration (l.l) 
 can be found in Diamond [ 1] . 
 
 The second set of parameters may be used to alter the auxiliary 
 : matrix for each n, producing the iteration 
 
 (A+B )X . = (A+B )X - (AX -q). (1.2) 
 
 n n+1 n n n 
 
 (Although there are algorithms which use a sequence of auxiliary matrices; 
 jthere has been no rigorous mathematical basis for the definition of an 
 
 jiuxiliary matrix, or a sequence of matrices, which will make (1.2) rapidly 
 
 i. 
 
 j convergent. 
 
Stone has developed in [3] an algorithm of this type which has 
 
 proved to be rapidly convergent on a number of test problems. The 
 
 fundamental idea of Stone's factorizations is to design the auxiliary 
 
 matrix so that the terms in each component of B X cancel. When B X and 
 
 n n n 
 
 B X are small it is intuitively plausible that one iteration of (1.2) 
 
 will yield a vector which is approximately equal to the solution of AX = q. 
 
 This cancellation is the basis of the factorization; however it 
 doesn't explain the convergence of the iteration. Indeed, experimental 
 results have shown that repeatedly using the auxiliary matrix B _, 
 
 (see 3*2 and 3*3 for precise definition) which produces the greatest 
 canceling yields a divergent iteration. 
 
 In this paper we explain the divergence of iteration (1.2) if 
 the auxiliary matrix B , is used on every iteration. We then present 
 
 some experimental results that indicate the eigenvalue of the iteration 
 matrix which has largest modulus grows without bound as the dimension of 
 A increases. 
 
3 
 
 2. ORIGIN OF THE PROBLEM 
 
 Consider the Dirchlet problem 
 
 - st: < a i» IB - i: ^ 1-) = .w, xed, 
 
 u(X) = i|r(X), for XeBD (2.1) 
 
 where x = (x,, x ), D is the interior of a compact region with "boundary 
 
 dD and where the differential operation in (2.1) is elliptic, that is 
 a (x) and a (x) are strictly positive on D. A common method of solving 
 
 (2.1) is to approximate the differential operator with a difference 
 operator. A system of linear equations results and the problem is then 
 to find a solution to the system of equations. 
 
 The system of equations is derived as follows. Assume D lies 
 in the first quadrant; cover D with two families of parallel lines, 
 
 x 1 = jh j = 1, . . . , n 1 
 
 and 
 
 x = kh k = 1, . . . , n . 
 
 2 
 
 Denote the points of intersection (jh, kh) by (j,k). The points (j,k) 
 are called the mesh points or gridpoints . Let dD, be the polygon formed 
 
 by joining those mesh points (i, j) for which one of (j"+l,k) or (j,k+l) 
 is not in D. Similarly let D, be the open and connected set bounded by 
 
 o\D, . At each gridpoint of D, the differential operator is approximated by 
 
 a difference operator and the equations are in one-to-one correspondence 
 
 with the gridpoints. 
 
 For simplicity of notation we shall consider the region D to be 
 
 the unit square buonded by the x, and x axes and the lines x, = 1 and 
 
 T 
 x = 1. Letting x = (x n _. x, _, ... x _, x n , .... x ) where 
 2 ° 1,1' 1,2' n, 1' 2,1' ' n,n 
 
 n, = n = n, we obtain the matrix equation 
 
AX = q (2.2) 
 
 defined by (AX) . , = q. ., where the subscript corresponds to the gridpoint 
 
 (j,k) and q. , is derived from q(jh, kh) and the boundary conditions. 
 
 The difference operator can be chosen so that the matrix A is 
 diagonally dominant and positive definite with a. . > and a. . < for 
 
 i = j. One method for deriving such an A is to approximate ^ — (a. (x) ^ — u( 
 
 1 1 
 
 by 
 
 -2 h h 
 
 h {a i (x+ -e i )[u(x) - u(x+he i )] + a i (x- g e i ) [u(x) - u(x-he i 
 
 where e-. and e p are unit vectors along the x, and x axes respectively. Thi; 
 defines A by 
 
 (AX) . . =B. .X. , - + D. . X. _ . +E. n x. n +D. n .x... 
 '3,k j,k J,k-1 j,k o-l, k j,k j,k j+l,k j+l,k 
 
 + B j,k+l X j,k-KL' 0,k = 1, ..., n, where 
 
 5 j k = " a 2^ h ' ( k -2^ 11 ^ 
 
 (2.3) 
 
 \k = '^(ti^M, 
 
 and 
 
 E j k = a 2 (jh, (k-|)h) + a 2 (jh,(k+|)h) + ai ((j-|)h,kh) 
 
 + a 1 ((j+ 2 -)h,kh). 
 
 When j = 1 or k = 1 we define D. , = or B . , =0 respectively. 
 
 J,k j,k 
 
 The following properties of A and the quantities defined in (2.3) will be 
 
 used below. 
 
and 
 
 (AX) . . = E . . X . . + B . . X . . ,+D. n X. nn + D . . . x . _ . 
 
 j,k+l j,k+l 
 
 B. , < , B. . = 
 3,^ ~ j,l 
 
 V'*° > v =o 
 
 (2A) 
 
 E. . > -2(B. . + D. . ) with equality holding if B. . ^ 
 
 and D. , / for j = 1, 2, ..., n and k = 1, 2, . .., n. 
 
 The subscript (j,k) refers to the gridpoint (jh, kh). 
 
3. THE FACTORIZATIONS OF STONE 
 
 Stone has defined two factorizations both of which depend on a 
 parameter a. One yields a symmetric auxiliary matrix and the other yields 
 a non symmetric auxiliary matrix. 
 
 For the symmetric case, [4] the factors L and U of A + B are 
 defined as 
 
 where 
 
 and 
 
 (LX) . . = b. . x. , -,+cx + d. X. , 
 3f* 3,k D,k-1 g,k o-l, k j,k j,k' 
 
 (UK) . . = x. , + e. ,x. _ . + f . .x. . _, 
 ; 0,k j,k j,k j+l,k o,k 3,3*4-1' 
 
 (3-D 
 
 0,k " j,k" ac j,k-l j-l,k-l' 
 
 c . = D. -ab . - , e. , (3-2) 
 
 j,k o,k a-i,k o-i,k-i' 
 
 0,k o,k o,k-l o,k j-l,k g,k j.,k-l j-l,k-l 
 
 + ab . i , e . , , -, 
 0-1, k o-l, k-1, 
 
 d . . e . , = D . , _ . - ab . , e . . , , 
 
 0,k j,k J+l,k j,k 3,k-l' 
 
 df = B. , ,, - ac ,f. t , i 
 0,k o,k o,k+l g,k o-l,k 
 
 In addition we define 
 
 C . . = Id . . e . , . and G . , = c . , f . . , . 
 3,k o,k j,k-l j,k oA j-l,k 
 
7 
 
 The nonsymmetric case [3] is more general in that it is not 
 necessary to have A defined as in (2.^). To avoid any confusion the 
 original matrix will "be called M and the auxiliary matrix N. M may be 
 nonsyrametric as is the matrix which naturally arises in the solution 
 of a problem with Neumann boundary conditions. Let M be defined as 
 follows : 
 
 (MX). , *= B. ,X. , . + D. ' X. _ . + E. .X. . 
 3,k j,k j,k _1 3,^ D-l,k j,k j,k 
 
 , j,k j+i,k j,k 3,k+l 
 
 The auxiliary matrix N is defined by factors L and U which are in the form 
 (3.1) just as in the symmetric factorization. Instead of (3.2) however, the 
 elements of L and U are defined by 
 
 j,k "" j,k " j,k' 
 
 c j,k " j,k " j,k' 
 
 d i k + b -i k f i k-1 + c -i k e i-l k = E i k + ^i k + ^-i k> 
 
 d. ,e. . = F. . - qC . . , 
 3,k 3,3s 3,k 3,k' 
 
 (3-3) 
 
 3,k 3,k . 3,k " 3,k> 
 
 where C . , = b . , e . „ - and G . n = c . , f . _ , . 
 3,k 3,k 3, k-1 3,k 3,k 3-l> k 
 
 Since there will usually be no confusion as to whether the auxiliary 
 matrix is defined by (3*2) or (3-3), B will be used to denote either of them. 
 With L and U of the form (3-1) A + B (or M + N) is given by 
 
((A + B)X) - C(LU)X) . >k - b Jf ^ jM + \ k e j, k -l* j+ l, k -l 
 
 c j,k x j-l ; k j,k j,k d,k-l C j,k e j-l,k' x j,k 
 
 + d. , e . . x. , . + c . . f . _ , x. _ . _ + d. . f . . x. . 
 
 j,k j,k j+l,k j,k j-l,k j-l,k+l j,k j,k j,k + 1 
 
k. WHY STONE'S FACTORIZATIONS ARE EXPECTED TO CONVERGE 
 
 In this section we explain why the factorizations defined in 
 (3*2) and (3*3) with a = 1 are expected to make the iteration 
 
 (A+B)X n+1 = (A+B)X n - (AX n -q) 
 
 rapidly convergent. 
 
 The matrix M which multiplies the error vector E = X - X on 
 
 n n 
 
 each iteration is called the iteration matrix. It is well known that an 
 
 iteration is convergent if the spectral radius of the iteration matrix is 
 
 less than one. Moreover as the spectral radius decreases, the rate of 
 
 convergence increases. Therefore we will show how the factorizations (3*2) 
 
 and (3*3) are designed to produce an iteration matrix with a small spectral 
 
 radius . 
 
 Throughout the following the eigenvalues of an n x n matrix M 
 
 will he denoted by \. (M) with 
 
 W M) = V M) ■> V M) ••• ^V M) = W M) - 
 
 The iteration matrix of 
 
 (A+B)X n+1 = (A+B)X n - (AX n -q) (k.l) 
 
 is M = I - (A+B) A. The spectral radius of this matrix is small if the 
 \. ((A+B) A) are close to 1. Thus the iteration is rapidly convergent if 
 
 (A+B) A is approximately equal to the identity I. That is if 
 
 lX-(A+B) AX || ig small fQr all x ^ Q> ([) _ >2) 
 
 PI 
 
10 
 
 The factorization (3*3), suggested by Stone, attempts to satisfy 
 (^.2) by making ||(A+B)X - AX|| small in the sense that 
 
 (BX). = o(h 2 ) (k.3) 
 
 where h is the meshsize of the grid used in the discretization of the 
 differential problem, i.e., deriving (2.2) from (2.1). 
 
 To see that BX satisfies the condition (k-3) we have the following. 
 
 Let x = (x_ .., .... x ) x where the x. . are values of a smooth function 
 1,1' ' n, n j,k 
 
 taken at the gridpoints (jh,kh). If B is the matrix given by (3«3) and C_ 
 
 3,k 
 
 definition of A in (2.1 
 
 to the definition of A+B in (3-3) and (3-^) yields 
 
 and G. are also as in (3.3), then comparing the definition of A in (2.k) 
 
 [(A+B)X] . . = (AX) . . + C. . [x. in , 1 + a(x. . -x.^- , -x. , ' 
 J 0,k y j,k j,k L j+l,k-l j,k j+l,k j,k-l' 
 
 + G. n [ x . _ . , _ + a(x . , -x . , . -x . .)]. 
 j,k j-l,k+l 3 ,k j-l,k o,k+l /J 
 
 V/hen a = 1, Taylor's Theorem implies the bracketed terms of (k.k) satisfy 
 
 2 
 
 x + (x. -x... , -x. _ ) = 0(h ) 
 
 j+l,k-l j,k j+l,k j,k-l 7 
 
 and 
 
 x. n . . + (x. , -x. _ , -x. , -, ) = 0(h ). 
 
 ,]-l,k+l j,k g-l,k j,k+l ; 
 
 Stone's symmetric factorization (3.2) satisfies (^.2) in the sense 
 BX) . = 0(h). This also follows from applying Taylor's Theorem to 
 
 (BX) = ((A+B)X-AX) . 
 
 Of*- d)^~ 
 
11 
 
 5. WHY STONE'S FACTORIZATIONS DIVERGE WHEN NO PARAMETERS ARE USED 
 
 The meaning of small used above is "weak since (^.2) and (k.3) 
 are not equivalent. In fact, empirical results show that both factorizations 
 (3.2) and (3.3) yield a B such that p(l-(A+B) A) > 1 when A is derived from 
 the Laplacian [a (x) = a (x) = 1 in (l.l)] and n = 30. Thus even though 
 
 (BX) . . = 0(h) or 0(h 2 ), ||X-(A+B) -1 AX|| is not small. It isn't even less than 
 
 1. In this section we will see how the above situation is possible. 
 
 We shall begin by deriving some properties of the eigenvalues of 
 
 (A+B)"^. 
 
 Lemma 5.1 : 
 
 If A is given by (2.k) and B is defined by either of Stone's 
 factorizations, (3*2) or (3-3)> then one is an eigenvalues of (A+B) A. 
 
 Proof : 
 
 If B is given by (3»2) or (3*3) then the first row and column are 
 null. Hence B has a zero eigenvalue. Let w be the corresponding eigenvector. 
 
 Then (A+B)w = Aw and consequently (A+B) Aw = w which completes the proof. 
 
 Theorem ^.1 : 
 
 The vector w is an eigenvector of (A+B) A iff Bw = or Aw = o"Bw 
 where <j is a scalar other than or -1. Furthermore the corresponding 
 
 eigenvalue is -2=r if Bw J and is 1 if Bw = 0. (5-1) 
 
 CT+-J- 
 
 Proof : 
 
 Suppose Aw = CT Bw a ^ -1, 0. Then 
 
 w + A Bw=w + -w= (- — )w. 
 a a ' 
 
12 
 
 Hence 
 
 (A+B) _1 Aw = (A+B)" 1 A(-^ T )(w+A" 1 Bw) 
 
 a+-L 
 
 = (- 2 T )(A+B)" 1 (Aw+Bw) = -2-w. 
 CT+-L a+J- 
 
 Now suppose (A+B) Aw = —rr v, a 4 ~lj 0. Then 
 
 cr+J- 
 
 Aw = -MA+B)™ 
 
 CT+ 1 
 
 and therefore Aw = a Bw. 
 
 The case Bw = was taken care of in Lemma 5.1 and the proof is 
 therefore complete. 
 
 For simplicity we shall now consider only the symmetric 
 factorization (3.2). The eigenvalues \. = A.. (l-(A+B)~ A) are real and 
 
 have magnitude greater than one iff 
 
 °i 1 
 
 A.. I = ll - =- I where -2 < a. < -—. 
 
 i ' cr. + 1 ' —i—2 
 
 If A is defined "by (2.k) and (A+B) is defined by (3 .2) then it is easily 
 shown that A.. ( (A+B)"^) > 0, see Diamond [1]. This and Theorem 5.1 imply 
 
 a i 1 
 T-T > 0. Thus if a- < - o then a. + 1 must also be negative and cr. < "!• 
 
 CT . + J. 1 <d 1 1 
 
 Hence the iteration (^.1) is divergent if there exists some w such that 
 
 Aw = crBw, -2 <C cr <; -1. Moreover as approaches -1, (l - — ~t) approaches 
 — — a + ± 
 
 infinity and the iteration diverges rapidly, i.e. the error grows rapidly. 
 
13 
 
 The iteration is rapidly convergent "when \\. | is small or, in 
 Dther words, when \a-\ is large. The fact that (Bw) . is small is therefore 
 Dnly significant if it is small relative to (Aw) . If for some w ||Aw|| is small, 
 
 say 0(h ), then selecting a B which satisfies (^.3) does not guarantee (^.2). 
 
 rhis difficulty arises on even the simplest problems as seen from the following 
 
 considerations. It is well known (see for example Frankel [2] that for A the 
 
 2 2 
 (n x n ) matrix derived from the Laplacian, 
 
 X . (A) = — 5 = 0(h 2 ) where h = 
 
 min 2(n+l) 2 ' ' ' " n + 1 
 
 Thus in this case the fact that B as defined by Stone satisfies (EX) = 0(h~ 
 or 0(h ) doesn't imply that ||x - (A+B)" 1 AX|| is small for all X ^ 0. 
 
Ik 
 
 6. EMPIRICAL RESULTS CONCERNING X^^ 
 
 In this section we shall present some empirical results which 
 
 indicate \ >r „, r ((A+B) A) grows without "bound as the dimension of A increases, 
 MAX 
 
 n 
 This is done by showing the existence of vectors w , where the n corresponds 
 
 to the dimension of A, such that 
 
 < Aw ,w > 
 lim l = co 
 
 n-*» < (A+B)w *,w 
 
 Then using the equality 
 
 we conclude lim ^ MAY ((A+B) A) = 0. 
 n-*» 
 
 n 
 The results of Section 5 suggest the vectors w will satisfy 
 
 Aw n = a Bw 11 (6.2) 
 
 n 
 
 where 
 
 lim a = -1. 
 
 n 
 n-*» 
 
 The vectors do not satisfy (6.2). However, as n increases 
 
 (Aw ) S -(Bw n ) , for j / l,n and k J- l,n: 
 
 and 
 
15 
 
 < Aw , w > = a < Bw ,w > 
 ' n ' 
 
 •where (6.3) 
 
 lim a = -1« 
 n 
 
 n-*» 
 
 From (6.3) we see that 
 
 . . n n ,. 
 < Aw ,w > 
 
 1 
 
 ff n 
 
 < (A+B)w n ,w n > 
 
 CT n 
 
 1 + *n 
 
 and this approaches infinity as a ^ -!• Furthermore the vectors w do 
 not maximize 
 
 < AX,X > 
 < (A+B)X,X > ' 
 
 Instead they give a lower bound for 
 
 $ < (A^')X,X > 
 
 and thus, from (6.1), they give a lower hound for A. . ((A+B) A). 
 
 The matrix A is derived from the Laplacian and Bis the auxiliary 
 matrix defined in (3*2) with a = 1. (The results also apply to the auxiliary 
 matrix defined in (3*3) since the elements of the two matrices approach the 
 same values as the order of A increases. This can be seen by noting that the 
 proof of Theorem 5-1 in [1] can be applied to the elements of the nonsymmetric 
 factorization as well as the symmetric factorization. ) 
 
 A number of different vectors were tried on a trial and error basis 
 
 in attempting to find a w such that < Aw , w > = - < Bw , w >. The process 
 of selecting vectors was started by using the algorithm of Diamond [1] to 
 
16 
 
 calculate the eigenvector w corresponding to the largest eigenvalue of 
 
 (A+B) A. Next a function f(x , x ) was found such that f(jh, kh) = w . 
 
 A vector F = (f , ..., f ), where f = f(jh,kh) and h = l/(l+n), 
 
 was defined and used as the first trial vector. A more or less random 
 
 procedure was employed to alter the trial vector F. and form the succeeding 
 
 n 
 vector F. _ • The final vector is denoted by w , where as before the n 
 
 corresponds to the dimension of A. 
 
 From the eigenvector, which was calculated to select the first 
 trial vector, the following information was inferred about the function 
 f(x x ,x 2 ). 
 
 The function should have a maximum at the center of the 
 
 rl 
 , 2 i 
 
 ,1 In 
 
 region, (-, -J. 
 
 2. The function is symmetric with respect to the lines x = -x , 
 
 1 i .1 
 
 X_ - r>) X p - 1 anC! - X "I _ O* 
 
 3. The function approaches zero as (x , x p ) approaches the 
 
 boundary. 
 k. The function decreases faster along the line x = x than 
 
 along x = 1 - x . 
 
 Functions with the above properties were tried with the additional 
 possibility that they were periodic in the x and x directions with period 
 
 Tr/m , m an integer. In this case the function exhibited the above propert: 
 
 in each period. 
 
 n n n n 
 
 The intention was to determine < Aw , w > and < Bw , w > when n 
 
 was very large. That is when B may be partitioned as 
 
 B = 
 
 B l,l B l, 
 
 '2, 1 \, 2/ 
 
IT 
 
 where B„ is an {W X HT) submatrix with almost constant diagonals. Using 
 the limiting values of the elements of B as derived in [ 1] we see that 
 
 B 2,2 " B 2,2 
 
 Instead of computing < Bw ,w > the quantity < B p , p w ,w > was 
 computed. We have < B , w , w > = < B , w , w > = < Bw , w > where 
 
 2 ' 2 
 
 2' 2 
 
 n 
 w = 
 
 w 
 
 and is an (n -IT) -component zero vector „ Similarly < Aw n ,w n > was computed from 
 
 . A n n ^ ^ . N N , 
 
 < Aw , w > = < A _w _, w > where 
 
 A 
 
 \l A l,2 
 
 k 2, 1 A 2, 2 
 
 Then instead of computing 
 
 
 a n : 
 
 / A n 
 
 < Aw , 
 
 < Bw , 
 
 n ^ 
 w > 
 
 n 
 w > 
 
 < A 2 g w , w > 
 
 < B 2 2 w , w > 
 
 the quantity 
 
 
 
 
 
18 
 
 ^ < A 2 2 w , w > 
 
 ° N = < w, ^ > 
 
 was computed. 
 
 N 
 The vector w , as described above, was chosen to be 
 
 v k N = exp {-(|jh-|| + |kh-|l + ±f [2 - cos( (jh-kh)^|)]} 
 
 cos 2 °((jh-kh)^|). 
 
 -x -x- 
 For N = 100 and N = 250 the values of a were a _ = -1.0918 and 
 
 -X -x- 
 
 a„.-~ = -1.066l. Values of N between 100 and 250 produced values of u_ T 
 2pU N 
 
 -X "X* 
 
 which varied monotonicly from cr-, m to Open • This indicates that as N 
 
 -x 
 increases cr gets close to -1 and 
 
 1 + a -, * 
 
 n 1 + a N 
 
 grows without bound. Thus X . ((A+B) A) grows without bound as the 
 dimension of A increases. 
 
LIST OF REFERENCES 
 
 [1] M. A. Diamond;, "An Economical Algorithm for the Solution of 
 Elliptic Difference Equations Independent of User-Supplied 
 Parameters/' Ph.D. Thesis, Department of Computer Science, 
 University of Illinois at Urb ana- Champaign, Urbana, Illinois, 
 to be printed. 
 
 [2] S. P. Frankel, "Convergence Rates of Iterative Treatments 
 of Partial Differential Equations," MAC, Vol. k, 1950. 
 
 [3] H. L. Stone , "iterative Solution of Implicit Approximations 
 of Multi- dimensional Partial Differential Equations," SIAM 
 Journal On Numerical Analysis , Vol. 5j> 19&8. 
 
 [If.] , Private Communication, April, 1969- 
 
Hov 
 
 Bq 
 
 19?$