HYBRID COMPUTER SOLUTION 0F LINEAR STATE MODELS Thesis for the Degree of Ph. D. MICHIGAN STATE UNIVERSITY WILLIAM C. ELLSWORTH 1969 This is to certify that the thesis entitled HY BRID COMPUTER SOLUTION OF LINEAR STATE MODELS Date 0-169 presented bg William C. Ellsworth has been accepted towards fulfillment of the requirements for M; degree in H;— i ';1' '. .1 " I { *“‘ A fl ) J n "' ’_.’1LhL“‘/‘L*~fix%; ’ Major professor I I January 22, 1969 '71 M, UV ' LIFT, 4R" Mich; .u “Scare . University 'V'o ‘ v l 33.71% ABSTRACT HYBRID COMPUTER SOLUTION OF LINEAR STATE MODELS BY William C. Ellsworth The solution of a set of differential equations of the form ngUc) = Ax(t), where X(t) is an n-dimensional vector function of time and A is an (n x n) real square matrix, is desired in many engineering problems. By solving the system of equations on a hybrid computer, the digital computer can be assigned the arithmetic and automated analog set-up, and the analog com— puter can simultaneously integrate all the state variables, thus utilizing the best features of both machines. The pr: cedure can be programmed to require not more than (2n—ll Potentiometers nor more than 2n operational amplifiers (n bipolar integrators) on the analog computer. This require— ment can be met by tridiagonalizing the A matrix on the digital computer as well as having the digital computer scale the transformed system in time, amplitude and initial condi- tions before integrating the transformed system on the analog computer. William C. Ellsworth In this thesis, an improved method for tridiagonalizing an arbitrary real square matrix is developed based on Lanczos' method.1 After transforming the given matrix A by unitary Householder transformations into a similar upper Hessenberg matrix AH, the latter is transformed into a tri— diagonal matrix T = RAHC by appropriate matrices R {with -1,. rows R.) and C = R \Wlth columns C.). The vectors R. 1 1 1+1 and Ci+1 are obtained recursively in the Lanczos algorithm starting with vectors R1 and C1 which are supposed to be arbitrary. But, if R1 and C1 lie in certain unknown sub- spaces, the algorithm breaks down. By choosing C1 = [1 0 ... 01T’ and obtaining Ci+1 directly from column (i+1) 3:1 become unit upper triangular in the regular case, and the 1 of the idempotent matrix U — Z CjR" the matrices R and C vector C.+ is never indeterminate. Whenever the computation 1 1 Of R.+ from the recurrence relation breaks down, this new 1 method describes a procedure for continuing the algorithm. The tridiagonalized system is scaled in time and ampli— tude by an automated procedure. Time scaling multiplies the tridiagonal matrix T by a constant fi. Amplitude scaling transforms ET by a diagonal matrix into a final tridiagonal matrix in which at least {n—l) of the off-diagonal elements adjacent to the main diagonal are scaled to O, 0.1, 1, 10, or 100 in magnitude. Thus, at most (Zn-1) potentiometers are needed on the analog computer. William C. Ellsworth The tridiagonalization, as well as the setting of the potentiometers, can be done automatically by the digital computer. The digital computer can also sample the analog solution and transforms it back to produce the solution X(t) for the given initial conditions X(O). 1R. L. Causey and R. T. Gregory, “On Lanczgs' 2:33:122m for Tridiagonalizing Matrices,“ Soc1ety O4 antober 1961), and Applied Mathematics Revigwg 111. NO- ( C ' 322-328. HYBRID COMPUTER SOLUTION OF LINEAR STATE MODELS BY William C3 Ellsworth A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Electrical Engineering 1969 ACKNOWLEDGMENTS The author acknowledges the advice and guidance fur— nished by Dr. J. S. Frame throughout the period of time devoted to this thesis. He also thanks Dr. J. B. Kreer for his assistance in selecting the thesis topic as well as his counsel on the analog computer aspects of this thesis. Furthermore, the author wishes to thank Dr. H. E. Koenig, Dr. R. C. Dubes, and Dr. Y. Tokad for their advice and guidance during the writing of this thesis. Finally, be- cause of her devotion and loyalty, the author is indebted to his wife Barbara for her assistance and understanding. ii II III IV VI TABLE OF CONTENTS INTRODUCTION . . . . . . . . . . . . . . . KNOWN TRIDIAGONAL TRANSFORMATION METHODS . . . Lanczos' Method . . . . . . . . . . . Elimination Method Q . . . . . . . . . . . Kublanovskaya's Method . . . . . . . . . . T Algorithm . . . . . . . . . . . . . . . La Budde's Method . . . . . . . . . . . . NNNNN U'IDhOOND-l A NEW SCHEME FOR TRIDIAGONALIZING AN ARBITRARY REAL SQUARE MATRIX . . . . . . . . . . . . . . 1 Analysis of the Problem . . . . . . . . . 2 Householder Transformation to Upper Hessenberg Form . . . . . . . . . . 3 A Modified Lanczos Transformation from Hessenberg to Tridiagonal Form . . . . 3. 3. 3. Recurrence Formulas . . . . . . . . . Krylov Factorization . . . . . . . . . Computational Technique . . . . . . SCALING . . . . . . . . . . . . . . . . . . . 4.1 Time Scaling . . . . . . . . . . . . . . 4.2 Amplitude Scaling . . . . . . . . . . 4. 3 Initial Condition Scaling . . . . . . . PROGRAM AND EXAMPLES . . . . . . . . . . . . 5.1 Program Flow . . . . . . . . . . . . . Tridiagonalization . . . . . . . . Scaling . . . . . . . . . . . . . 5.2 Examples . . . CONCLUSIONS . . REFERENCES . . . . . . iii Page 11 14 16 20 20 24 26 27 28 33 43 44 45 49 51 51 52 57 58 64 67 I INTRODUCTION In many engineering problems, it is necessary to solve linear state models of the form d a;x(t) = AX(t) ‘ (1.1) F(t) = CX(t) where X(t) is an n—dimensional vector function of time, A is an (n x n) constant matrix, F(t) is a p—dimensional vector function of time, and C is a (p x n) constant matrix. The homogeneous system of differential equations (1.1) includes the class of linear state models in which the drivers can be generated as the solution of the set of linear differential equations with constant coefficients shown below. 37;st arm +2.5.(t) . (1.2) F(t) = Eflit) +.2§(t) where §(t) is a particular solution of £1--1=«:(t)= grit) By combining (1.2) and (1.3), we obtain the system 2 am =32- [5558} [%§][§I:I} Axm ' 1.4) F(t) = [Q 2][§)E(]= CX(t) ( which has the form (1.1). The object of this thesis is to present a method by which the hybrid computer may be used to solve equations of the form (1.1) in a manner which is numerically stable and optimal in the following sense: a. It should use a minimum number of operational ampli- fiers. b. It should use a minimum number of potentiometers. To solve equations of the form (1.1) on the analog computer, one potentiometer is needed for every entry of A not 0, 1, 10, or 100. To simplify the computation and meet the criterion of optimality, the matrix A is transformed into BSAS’I, a scalar multiple of a matrix similar to A, having at most (3n—2) non-zero entries of which as many as 10, or 100 in magnitude. However, for possible are 1, reasons of numerical stability, we wish to avoid methods requiring the solution of the characteristic equation, such as the transformation to Jordan form. The proposed method will require only rational operations and the extraction of Square roots. The change of variables t = Br, 2(1) = SX(t) = sx(er) (1.5) transforms the differential equations in (1.1) into 3 32—132“): sans-12h) (1-6) After this transformation has been implemented on the digital computer, the resulting equations (1.6) are solved on the analog computer. The digital computer can then sample the continuous analog solution periodically, re—transform the solution back to the problem solution X(t), and print the entries of X(t) and F(t). The entire process, from reading the equations (1.1) to printing the solution, can be made fully automatic using the hybrid computer. In the solution procedure described in this thesis, the best properties of both the digital and analog computers are utilized; the digital computer does the arithmetic, bookkeeping, and analog set-up, and the analog computer integrates all the transformed state variables simultaneously. The means proposed for determining an optimal matrix, similar to A, without computing the eigenvalues is to trans- form the matrix A into a tridiagonal matrix T having at most (3n-2) non-zero entries on or adjacent to the main diagonal. Chapter 2 presents five ways for tridiagonalizing a matrix that are described in the recent literature. Each of these methods has significant shortcomings which, for algo— rithmic computation, involve undue complexities and possible breakdowns. One type of shortcoming is a numerical instability due to division by a small quantity. The other type is called a breakdown and results when the algorithm cannot be continued 4 without modifications which may or may not be possible. Chapter 3 presents an original tridiagonalization algo— rithm based on Lanczos'[1] method, but which is better adapted to automation. The new algorithm produces a unit upper triangular transforming matrix in the regular case, starting with a computer—determined initial vector which the computer itself modifies in the irregular cases. Two theorems are proved which show that the algorithm can always be implemented, even if a breakdown occurs. Once the tridiagonal matrix is obtained, it must be scaled both in time and amplitude. Furthermore, the initial conditions must also be scaled before the system (1.6) is ready for solution on the analog computer. As another original contribution in this thesis, Chapter 4 presents an automatic scaling procedure whereby as many as possible of the entries on the three diagonals are assigned the values 0, 1, 10, or 100 in magnitude. II KNOWN TRIDIAGONAL TRANSFORMATION METHODS Five recently published methods for transforming an (n x n) real but non—symmetric matrix into a similar tri- diagonal matrix serve as a background for the theoretical development in this thesis. Following a description of the techniques and shortcomings of each of these methods, a new and improved method will be presented in Chapter 3. 2.1 Lanczos' Method[1] Given the arbitrary, real, square matrix A of order n, Lanczos' method attempts to construct non—Singular trans— formation matrices R and C such that (2.1.1) (2.1.2) RAC = T H RC U where U is the unit matrix and T is tridiagonal of the form i—tl q2 0 0 0 i e2 t2 Q3 °°' 0 0 T = 0 e3 t3 0 0 (2.1.3) 0 0 0 tn-l qn _D O 0 en tn_ 6 th ‘ row of R and Ci the ith column Let Ri denote the i of C; then Ri and Ci are computed recursively startirg with somewhat arbitrary vectors R1 and C1 satisfying the relation Equations (2.1.1) and (2.1.3) require that the entries in T have the following forms: = = l = 0 o 0 \ tii ti RiACi' i 1, , n (2.1.4) ti+1,i = e1+1 = Ri+1ACi’ 1 = 1’ '°" n-l (2'1°5> ti,i+1 : qi+1 = RiACi+1' l = 1’ ...’ n-1. (2°1'6) Using (2.1.2), it can be seen that (2.1.1) may be written in two other forms RA = TR (2.1.7) (2.1.8) ll 0 *3 AC It then follows that 1 = 1, 11—1 (2.1.9) R q 1+1 Ri(A — tiU) — eiRi_1, 1+1 (A — timci — qui—1' i = 1, '~-, n-1, (2.1.10) e . . 1+1C1+1 where q1 = el = 0, R0 = 0, and C0 — 0. As an intermediate step for subsequent discussion, let the vectors X and Y at stage i be defined as X = Ri(A - tiU; - eiRi_1 (2.1.11) and Y=(A-qm—qgkl 9AA? At the jiil iteration, the quantities ti, ei+1, CET1’ - O ‘ I I L ' v d - . to be determined in the order indicate . ql+1, and Rl+1 are First, t. is obtained from (2.1.4). If the column vector Y 1 in (2.1.12) is not zero, set ei+1 = 1, Ci+1 = Y. (2.1.13) If, on the other hand Y = 0. set ei+1 = 0 (2.1.14) and choose some vector orthogonal to R1, ---, R. as C. i i+1° If the row vector X in (2.1.11) is not zero, then qi+1 is computed using (2.1.6) and Ri+1 is then determined from R- = X . (2.1.15) However, if X = 0, set q. = 0 (2.1.16) and choose some vector orthogonal to C1, °°-, Ci as Ri+1' normalized so that R. 1+1Ci+1 = 1. These choices of Ci+1 and Ri+1 are consistent with (2.1.6) and (2.1.5), since the orthogonality criteria require Ri+1(ACU = Ri+1(ei+1ci+1 + tiCi + q1C1-1) = ei+1(Rir1Ci+1' %(Ri+1ci) + q1(Ri+1ci- ’5 =e. +0 +0 (2.1.17) 1+1 \ (RiA)Ci+1 = (qi+1Ri+1 + tiRi + eiRi-I/Ciri . C “ = qi+1(Ri+1Ci+1) + t1(RiCi+1) + e 1(Ri— 1 1+1 = q, .+ 0 + 0, (2.1.18) 1+1 - d'rectl' in The vectors Ri+1 and Ci+1 can be expressed 1 / terms of R1, C1, and A with the help of LanCZOS‘ polynomials. a [3] Consider the characteristic matrix [AU - T]. Householder 8 calls the leading principal jfifl order minors of {AU - T] Lanczos polynomials. Letting P_1(A) = O and P0(A) = 1, the expansion of the minor Pi(k) by cofactors of the last column yields the recurrence relation Pi(1) = P i_1(1)(1-ti) — Pi_2(>\)eiqi (2.1.19) Assuming that qi+1 # 0 and ei+1 # 0 for i = 2, --~, n and using (2.1.19), we see by induction that (2.1.9) and (2.1.10) may be written as q2 --- qi+1Ri+1 = R1Pi(A) (2.1.20) e2 ... eific1+1 = Pi(A)C1. (2.1.21) . l ‘_1 If P.(k) = Z P..13 ' , then i ._ 13 3—1 1 j 1 1 1 . ‘ = ..F. 2.1.22 R1Pi(A) '2 pinlA _§ 913 3 ( ) 3—1 3 1 where F. = RlAj-1 is row j of the Krylov matrix F. F1 R1 F R A F = 2 = 1 . (2.1 23) —1 Fh RlAn _ L I— _. Hence, from (2.1.20) and (2.1.22), the matrix R is related to the Krylov matrix F by the simple formula _ 2.1.247 R - LTF, ( where L is a lower triangular matrix whose entries in T . . j - 3. ' ' . row (1+1) are the coeffic1ents of A in Pi(h) d1V1ded by the cumulative product QZ "° qi+1° 9 In Lanczos' method, the vectors R1 and C1 are chosen arbitrarily subject to the condition that R1C1 = 1. Some choices of R1 and C1 lead to a breakdown because one or both of the products in (2.1.20) and (2.1.21) vanish, but these choices cannot be predicted in advance because the coeffici- ents in the polynomial Pi(A) depend on R1 and C1. 2.2 Elimination Method In the Elimination Method, described by Strachey and Francis[7], the given square matrix A of order n is first transformed to the lower Hessenberg form [L11 H12 0 "' 0 _ h21 h22 h23 "' 0 H = h31 h32 h3;3 0 (2.2.1) USing pivotal condensation with row and column interchanges. This Hessenberg matrix H. similar to A, is then transformed into tridiagonal form, when possible, by a speCial case of Lanczos' Method. Definition; An elimination transformation of tne matrix A . —1 is an elementary similarity transformation EAE where the elementary matrix E consists of the unit matrix plus one Off-diagonal element eij' which is chosen so that EAE has one more zero entry than A. 10 For a 3 x 3 matrix, the reduction to Hessenberg form is accomplished by one elementary transformation as follows, assuming a12 # 0. F1 0 O iF-all a1? a13 F1 0 O - -1 _p 0 ll a31 a32 a33 O O 1 A '_ 2.2.2 Féll a12 0 i ( ) = a21 352 3&3 = A'- I a31 332 a33‘ The primed elements of Al are those that were altered by the transformation. Note that a zero is produced if e23 = a13/a12. However, if a12 = O, a transposition of rows and columns must precede this step. By allowing row and column interchanges during the transformation to lower < ' ‘1 Hessenberg form, we can ensure that ieij‘ _.1 for optimum numerical stability. In the elimination procedure that reduces the lower Hessenberg matrix to a tridiagonal matrix. the key step 18 ../h. .; times row (j+1) from row i(i > j+1) 13' 3+]r3' to produce a zero in the ij pOSition. to subtract (h Row and column interchanges cannot be used at this stage since they would Hence. alter the form of the upper triangular portion. some of the e..'s may have magnitudes greater than one. 1 This may lead to numerical instability. 11 Additional difficulties arise if some of the sub— diagonal entries are zero. Consider an example in which H has the form F1111 hlz 0 o- H = ' . (2.2.3) h41 h42 h43 h44_ Since h21 = 0, this method cannot be used to reduce h31 and h41 to zero. An analysis of the transformation from Hessenberg to tridiagonal form by this method reveals that the transforma- T tion is identical with Lanczos' Method in which R1 = C1 = [1 O ... O], and, consequently, R and C are both lower triangular. T4 2.3 Kublanovskaya‘s MethodL ] 21 . . An earlier method of Hessenberg[ J is essentially equivalent to the first half of the elimination method in Give: Producing an upper Hessenberg matrix H similar to A. an n'th order square matrix A, find a lower triangular, non-singular matrix C such that C-IAC = H, (2.3.1) where In 1 n22 hzs hzn H = 0 1 h... h... . (2.3.2) The columns Ci of C and the entries hij of H are computed as follows. Let (2.3.1) be written as AC 2 CH. (2.3.3) By investigating column i on both sides of (2.3.3), it can be seen that Aci = .i cjhji + Ci+1 (2.3.4) 321 . or 1—1. Ci+1 = (A - hiiU Ci - jEI thji . (2.3.5) In general. the first i entries of column i in (2.3.4) are used to find hki for k = 1, ...3 i. The last (n—i) entries of (2.3.5) are used to find Ci+1 since C is lower triangular. The process is then repeated with i replaced m I . by 1+1. The author uses [1 O ... 0]‘ for C1. ConSider the following example. 13 F511 312 313 Fa O 0 F1 0 0 [£11 h12 h13 _331 a32 33:3) 0 C32 C31) 0 C32 Casi) 0 1 h33_) 11-) {fill—i _} -9 h11 a31_) Cazfl i = 2:AC2 = C1h12 + C2h22 + C3 igizczz + a13C32 )- h12 .} ‘2 h12 a22022 + a23C32 = C22h22 } ‘9'h22 if C22 # 0 _?32C22 + 333C32J C32h22 + Cay } ‘9 C3 1 = 3:AC3 = C1h13 + C2h23 + C3h33 ($13c33 F. h13 } '-> h13 a23C33 = C22h23 } ‘“> h23 if C22 # 0 _?33C33_) C32h23 + C33h33) } ‘“> h33 1f C33 ; O In this method, C1 is an arbitrary vector, which may lead to a breakdown whenever cii — 0. Kublanovskaya describes a transformation to tridiagonal form that is equivalent to the Elimination Method and has the same difficulties with stability and breakdown. Given an upper Hessenberg matrix H of the form (2.3.1), find an uPper triangular, non-singular matrix R such that ‘1 , (2.3.6) 14 where T is tridiagonal of the form i—tll 1 O ... O t21 t22 1 O T = O t32 t33 ... 0 (2.3.7) 0 O _ O tnnJ The . ' rows R1 of R and the entries tii and ti+1,i of T are computed recursively as follows. Let (2.3.6) be written as RH = TR' (203.8) Upon investigating row i on both sides of (2.3.8), we have R.H = t R (2.3.9) .. +t...+. i i,i—1 llRl Ri+1 i—1 or (2.3.10) Ri+1 2 Ri(H ’ tiiU) ‘ ti,i~1Ri-1' In general, the (i-l)St and ith entries of (2.3.9) are used to find ti i-1 and tii respectively. The last (n-i) I entries of (2.3.10) are used to find Ri+1° The process is again repeated with i replaced by (i+1). Once again, note that R1 is an arbitrary vector, which may lead to a breakdown when rii = 0. 2-4 T Algorithm[8] When possible, the T Algorithm transforms a square com- Plex matrix to tridiagonal form by a sequence of similarity 15 transformations called quasi-rotation matrices. It does so by using the (i+1, i) entry and the (i, i+1) entry to annihi- late the (j, i) and (i, j) entries at the same stage, where j > (i+1). Definition: A matrix R with entries rij is called a quasi- rotation matrix if and only if det R = i 1 and rij éij for qq H ' ' ' f , , r , a d With the pOSSible exceptions o rpp rpq qp n any 9 and q (p # q)- The suitable quasi-rotation matrices used at the (i, j) stage are determined from ten categories listed in a table. The categories are listed according to whether the elements in the (i, j), (j, i), (i, i+1), and (i+1, i) positions are zero or non—zero. Therefore, a total of sixteen cases re- sult. The author admits that this algorithm might involve ex— tensive programming because of the table look—up. He also gives a sufficiency theorem for tridiagonalization of a matrix by the T Algorithm. For three of the sixteen possible cases, the algorithm can fail; however, the author presents a modified T Algo— rithm which yields a matrix that is almost lower triangular. Application of the T Algorithm to certain matrices leads to numerical instability. When the algorithm does work, it leads to a unique transformation matrix. 16 2.5 La Budde's Method[5] This method, like the others, seeks to find non—singular transformation matrices R and C such that for a square matrix A of order n, RAC = T (2.5.1) and II RC U (2.5.2) where T is tridiagonal. Furthermore, the matrices R and C are each the products of (n—2) matrices; i.e., R = R(n-2) coo R(1) (2.5.3) C = C<1>C< At stage j in this method, the matrix A has been trans— formed into a similar matrix A(3) of the form F I " : o T : _._- j rows A”) = I' wT (2.5-4) _—_—_-—T _____ L _____ 0 : V i B n—j rows _ ' I J where T is tridiagonal, V and W are (n-j) dimensional column vectors, and B is an (n—j) X (n-j) submatrix. Starting with two column vectors X and Y of dimension I .\ Kn‘j), we construct the elementary matrices R(J} and C(3) where 17 (2.5.5) T 0 Un-j +aXY 0 . \ . In order that R(J’C(3) = U, the vectors X and Y must be related to the scalars a and b by the equation T _ a + b) YX—- 1Tb~' (2.5.6) . n1 . /. a‘s . The matrix Ab+ > = R(J)A‘J)C(J) differs from A(3) by replacing V, WT, and B by V', W‘T, and B' where - = / T v \Un_j + aXY )v W'T = WT(Un_j + bXYT) (2.5.7) I _ T / T B — (Un_j + aXY )B(Un_j + bXY ). The scalars a and b and the vectors X and Y are chosen so that (2.5.6) is satisfied, all but the first entries of v' and W' are zero, and p = w‘Tv' = wTv. (2.5.8) The choice is not unique, but it must avoid values of a and b near zero. Scalars c and d are defined by T , (2.5.9) ll 2 ?< Dali-l ll |< < l. c From (2.5.7) and (2.5.9) we obtain d , x E c X = (V'-v)-5, Y = (w -w)g (2.5.10, and substituting (2.5.10) into (2.5.9), we have 18 1 t d I '.C .r ' g = (wlvl ~ 9);; '5 = (wlvl - Prg \2‘5-11) Then, (3— + p)(2— + p) = (w v')’w‘v ‘ = pa a .’2 5 12) Cd \Cd , 1 1 \ 1 1/ j+1lj j’j+1 ) ‘ ' Th 1 t' fo 3— 's e so u ion r cd i 1 {-p(a+b)b/p2(a+b)2 +4abp(aj’j+1aj+1tj-P)) (2.5.13) Ed ‘ 2ab In order that all but the first entries of W' and V' be zero, it is necessary that d {C for k = 3+2, °", n. Then Xj+1 and yj+1 are expressed as X = L + p-aj 'j+1aj+1j —g_____ j+l Cd 3 a. .+1 3'3 (2.5.15) y. = l__+ p-aj,j+1aj+i,j c 3+1 Cd b aj+1.j In order to solve (2.5.14) and (2.5.15) for the x‘s and y‘s, it is necessary that the scalars a. b, c, and d be finite and non—max» For some matrices, it can happen that the scalar product p is zero at some stage: then the algo— rithm breaks down because %5 = 0. Otherwise. c or d is arbitrary in (2.5.13), there is no loss of general- since either ity if we let d = 1. To solve for c, the product ab as well as p must be non-zero. The signs of a and b can be chosen so that the discriminant in (2.5.13) is positive, and the 19 sign of the radical is chosen so that c is finite. The dif- ficulty in choosing a and b is to avoid having p = 0 at the next stage. Wang and Gregory[91 point out that such a choice is not always possible. Parlett[6] has shown that when A is an unreduced lower Hessenberg matrix, this method is identical with the Elimina- tion Method. Lanczos' Method appears to be the most general since the Elimination Method and Kublanovskaya‘s Method are special cases of the former. Furthermore, Lanczos' Method appears to be better than the T Algorithm or the method of LaBudde. The T Algorithm involves a table look—up with sixteen cases. LaBudde‘s Method involves rational operations and the solu- tion of a quadratic equation at each stage, but with no assurance of avoiding a breakdown at the next stage. Becasue of the generality of Lanczos' Method, a modi- fication of his method is presented in Chapter III. This modified procedure is then used to tridiagonalize the matrix A. III A NEW SCHEME FOR TRIDIAGONALIZING AN ARBITRARY REAL SQUARE MATRIX In order to bring the new tridiagonalization procedure into perspective, it is necessary to examine the ways in which Lanczos' Method can break down. This examination, as well as an introduction to the new two-step procedure, is presented in Section 3.1. The first step of the procedure is discussed in Section 3.2. Section 3.3 describes the second step and includes a lemma and two theorems which show that it is always possible to transform an arbitrary real non-symmetric matrix into tridiagonal form without using the eigenvalues. 3.1 Analysis of the Problem In order to overcome the difficulties associated with Lanczos' Method, it is necessary to examine the ways in which his method can break down. At any stage, three cases may occur. The regular case, or Case 1, is the case in which the algorithm proceeds from one stage to the next without a breakdown. There are four ways in which the algorithm can break down; these are grouped into two irregular cases, Case 2 20 21 and Case 3. The reason for this grouping will become ap— parent when we consider how to proceed with the algorithm following a breakdown. A breakdown occurs when, at any stage of the algorithm, we have e. q = 0. In terms of the vectors X 1+1 i+1Ri+1Ci+1 and Y in equations (2.1.11) and(2.1.12), a breakdown occurs whenever XY = 0. We will say that a Case 2 breakdown occurs when either X = 0, or Y = 0, or both X = 0 and Y = 0. A Case 3 break— down occurs when X # 0, Y # 0, but XY = 0. To show that Case 3 can indeed occur, consider the following example. Let A, R1, and C1 be as shown. 1 -1 1 a 17 0 1 0 1 1 A = , R1 = [1 O O 0}, c1 = , 0 0 1 1 1 L0 0 0 1 L}- For i = 1, we would compute t1 = RIACL = 1. Then x = R1(A - tlU) — 0 -1 1 0} 2 0 and 731 1 Y = (A - t1U)C1 - 7! 0 1 L0-) but 22 Let m be the degrse cf the minimal polynomial. A =“t Case 2 breakdown will always occur at the (In-+1;D stage. The breakdown, in which x a 0 and Y = 0. is unavoidable. ry We will show that a simple procedure can e used at the stage in which a Case 2 breakdown occurs to permit continu- ation with Case 1. Furthermore, a Case 2 breakdown in which X # 0 and Y = 0 will never occur, because the new algorithm does not generate the columns of C in the same manner as the rows of R. As will be pornted out in a lemna, the columns Ci of C are computed recrusively as functions of rows R1, ~-~: Ri-l of R and coltmns C1, ...; C1-1 of C. As a first step in transforming the matrix A to tri- diagonal form I, we will transform it to upper Hessenberg f3, . - form AH by uSing Householder“ * transformations to be de— scribed in Section 3.2 such that A, = HAH* 3.1.1 . < ) HH* = U (H is unitaryV. (3.1.2) The upper Hessenberg matrix AH has the form 1—- a a a ~~ a hii h}; his “in c a a “ i“ o a- j A = 0 e. a _ “U“‘a _ - .1. - H w has rrf. (3 3/ 0 0 0 ... ar’, In. "‘ .... 'Ihe reason for this tran‘formatiox is to introduce the desired zeros below the subdiaqoral by stable computation. 23 A unitary transformation preserves the lengths of vectors and does not involve division by small quantities. In the second step, a matrix R is computed such that 1 RAHR- = RAHC = T is tridiagonal, and has the same sub— diagonal entries as AH. Whenever AH has a subdiagonal with nonvanishing entries, Theorem 1, in Section 3.3, shows how to find an initial vector R1 such that there will be no breakdown in the algorithm. The occurrence of a subdiagonal zero entry in AH invalidates the hypothesis of Theorem 1, but may or may not cause a breakdown. When a Case 3 breakdown does occur, the algorithm necessitates returning to the last occurrence of a Case 2 breakdown (whereby the first stage is treated like a Case 2 breakdown). At this point, the algorithm shows, by Theorem 2 in Section 3.3, how to progress one stage past the point at which the Case 3 breakdown occurred. The combined two steps of the procedure for tridiagonal- ' ' ' ' ' , re re— lZlng an arbitrary real non—symmetric matrix A can be p sented by the change of variables yo.) = _S_X(t) (3.1.4) where s = RH (3.1.5) The original state model g—X(t) = AX(t). F(t) = CX(t) (3.1.6) dt ” then becomes §__. 3 = /t., F t) = C'Y(t) (3.1.7) dt X(t, TY) ). ( -1 ' _ 5.1 Where T = §A§_ and C — C~ . 24 3.2 Householder Transformation to Upper Hessenberg Form[3] As mentioned in Section 3.1, our first aim is to find a unitary matrix H such that HAH* = AH is in upper Hessen— berg form. The matrix H is found as the product of (n—2) Householder matrices Hi; i.e., n_2 ... H1 . (3.2.1) Let A = A(1) be partitioned as (3.2.2) A51) A(}) 1 . where Afg) is an (n—l) dimensional row vector, A51) is an 1) . (n—l) dimensional column vector, and Aéz' 15 an (n-l) x (n—l) submatrix. Let a1 be defined by L‘— / -)l~ 01 = +:JA§:) A51), (3.2.3) and let X(1> be an n dimensional column vector defined by F07 F071.“ --—h d---‘ —--- X51) idl 1 (1) _ (1) - 0 (3.2.4) X( ) = X3 - A21 0 1-“.1 _ 1.3 (1) , . ~ - . The where the sign of a1 is chosen to maXimize [X2 I first Householder matrix H1 is defined by 25 H1 = U- ‘ 221 (3.2.5) where 21 is the following Hermitian idempotent of rank one; (1) (1)* X X X(l) X(l/ It is easily shown that Transforming A(1) by this H1, a partially transformed . /2 matrix A) ) is obtained: A(2) = H1A(1)H§ . (3.2.8) In general, for i = 1, "°, n—2, we have I— .. ' .— Afw : Al: A(i) = —-—]—-—-%-—-—-—~ , (3.2.9) OlAéPIAii) l l . . i where Aii) is an ith order upper Hessenberg matrix, AEZ) is an (i) x (n—i) dimensional submatrix, Agi) is an (n—i) dimensionalcxflumn vector, and Ag:) iS an (n—i) x (n-i) sub— matrix. Let / (-)*A(i) (3.2.10) and 26 {—0—} to TFO) . 6 (3 6 X(l) = —;(;) z _____ _ _;;-_ 1 (3.2.11) 1+1 _ i : Ali" : xii) ° _. _ 1-”.1 ( . , , , i where the Sign of ai is chosen to maXimize )xi+2 . As before, the Householder matrix Hi is defined by - -* X(l)x(l) H.=u-2( ... .), (3.2.12) 1 X(i)x(i) and A(l) is transformed into A(l+1) = H.A(1)Hl- (3.2.13) 1 1 When i = (n—2), the resulting matrix A(n-1) is in upper Hessenberg form; i.e., A — A(n-1> - H A(n_2)H = HAH* . (3.2.14‘ H " — n—2 n—2 _ ’ It is noteworthy at this time to point out that if A is symmetric, the well known method of Givens is precisely this Householder transformation applied to A, and the final matrix is tridiagonal. 3.3 A Modified Lanczos Transformation from Hessenberg to Tridiagonal Form Given the upper Hessenberg matrix A, our goal is to describe an algorithm to construct a unit upper 27 triangular matrix R (in the event that a Case 3 occurs, R might not be unit upper triangular) that will transform A into a similar tridiagonal matrix T; i.e., RAR-l RAC = T . (3.3.1) RA, this becomes In terms of the matrix G TRzRAzG . (3.3.2) The matrices A, G, T, and R have the forms: — G. p—- _1 all 312 a13 ”° ain 911 912 913 "' 91m e2 a22 a23 " azn 92 922 923 "' an A = 0 e3 a33 °°° asn . G = 0 e3 933 °°° gsn . . .0 o 6 O 0.. e g b-O O 0 en arm L n ng (3.3.3) _ _ T _ t1 Q2 0 °" 0 1 r12 r13 ° rin 82 t2 q3 " O O 1 r23 °° an T: 0 e3 t3°°’0.R= 0 0 1- ran ' ' °. 0 6 0 - 1 L__O O 0 en tQ. _ 4 Recurrence Formulas ~ We equate the ith rows of TR and G in (3.3.2). . = 5.. (3.3.4) \St " " ',it Comparing the ith, (i+1) (3 > 1+1) entries .th , and 3 can be seen that 28 i: eiri-1,i + ti = gii (3.3.5) .+ : ' ' . ' ' ' = l 1 elrl—1,l+1 + tlrl,l+1 + qi+1 gi,i+1 (3°3'6) j: eiri-1,j + tirij + qi+1ri+1,j = gij' (3.3.7) Since the subdiagonal entries of T are identical to those of A, then (3.3.5) and (3.3.6) may be used to compute the coefficients t. and q. in T as functions of R. , R., and 1 1+1 1-1 1 A; i.e., . = .. - . . . = a.. + e. r. . - e.r. . t1 gll elrl-1,l 11 1+1 1,i+1 1 1-1,i (3.3.8) . = . . - e.r. . — t.r. . 3.3.9) q1+1 g1,i+1 1 i—1,i+1 1 1,i+1 ( If qi+1 # 0, then (3.3.4) may be used to compute the vector Ri+1; 1 i+1 R . = — . . — t.R. . 3.3.10) 1+1 q (G' e1R1 ) ( The case where qi+1 = 0 will be treated later. Krylov Factorization Let F be the Krylov matrix whose rows are the iterates of R1 under A. R1A (3.3.11) In the regular case, F has a factorization as the triple 29 product of a unit lower triangular matrix L, a diagonal matrix D, and a unit upper triangular matrix V; i.e., F = LDV (3.3.12) where L, D, and V have the forms: F-l 0 --- 0 i I d1 0 --- 0 £21 1 ... 0 0 d2 ... O L = . D = . ° ° .. ° 6 6 ... d _}n1 gnz ° 1J .1 _ (3.3.13) G V12 Vln 0 1 V2n V = . O 0 ll 1'2""' 1—1’1 b the minor of F formed from LEt f(1’2’... i‘1Ij) e I rows 1,2,°--,i—1,i and columns 1,2,--°,i-1,j; further, let 1'2'°'°’i_1’i“ be the ith leading principal minor f’ = f(1,2,---,i-1,i/ Of F. It can then be shown that the entries in L and V are expressible as: ' ' 1 2 oo..i—1Ii ,1,2,”’,]-1;l f r I ..'. . \ fillzl...lj-1J-_> V = (1.2,. ,1-1,J «3.3.141, Eij = fj ’ ij f1 Since D and F have equal leading principal minors, it follows that f. d = f l where f0 = 1. (3.3.15) 30 It will now be Shown that the upper triangular matrix V in (3.3.12) is identical to the matrix R in (3-3.1). The vector Ri+1 is a linear combination of tie vectors R1, RIA, --, RlAi, in which RlAi has a non—vanishing coefficient. Hence, the matrix R is a left multiple of F by a lower tri- angular matrix, say (L'D')—1, where L' is a unit lower triangular matrix and D' is a diagonal matrix. Then, F = L'D'R = LDV. (3.3.16) By rewriting (3.3.16) as (L')'1LD = D'RV-l, (3.3.17) where the left side is lower triangular and the right side is upper triangular, it can be seen that both sides must be diagonal. Hence, (L')'1L ~ RV-1 = U and D = D'. There- fore, the factorization (3.3.12) is unique and V — R. -1 -1 , . Note that (L'D') = (LD) 15 the matrix LT referred to in equation (2.1.24). It can be shown that di = Q2q3 "' qi’ so the qi+1's are related to the entries in D and F by the formulas . f. = d1+1 = fi—l 1+1 (3.3.18) q1+1 di f2 1 and i i i+1‘j — - f — 1 (3 3 19) f. = II d. = II q. where d1 — f1 — 0 - . . . l j=2 J j:2 j . The t.'s are related to the subdiagonal entries in L by 1 - - . for i = 2,---,n-1, (3.3.20) t1 = £21’ ti — £i+4,j. £1,1-1 31 and tn must be computed using either (3.3.8) with en+1rn,n+1 = 0 017- (2.1.4). At stage (i+1), the computation of Ri+1 requires that qi+1 # 0. This condition is met if and only if fi+1 # 0. Hence, the non-vanishing of the leading principal minors of F is necessary and sufficient for the continuation of the algorithm in the regular case. We will now investigate the minors fi as functions of R1 and show, by means of two theorems, that it is always possible to transform a given Hessenberg matrix to tridiag- onal form without involving the eigenvalues. Furthermore, the existence proofs are constructive. Although they may not be computationally optimal, they do include methods for determining the vector R1 which can be written as R1 = [1 r2 r3 °°° rn]. (3.3.21) Let the notation fij(r2, 5k) imply that the entries fij of the Krylov matrix F depend linearly on r2, °-°, rk for 2 f.k < (i+j), and if (i+j) = (k+1) : (n+1), they involve rk explicitly in a term e2e3---ekrk. We emphasize this latter point by underlining rk in the nota- tion fij(r2, ---, Ek). Hence, the matrix F can be written in the following form: f21(£2) f31(r2:£3) fn1(r2,...,£fl) Theorem_1: r2 f22(r2,£3) f32(r2'r3'£A) f33(r2,...,r5).. fn2(r2,...,rn) fn3(r2,...,r 32 r3 f23(r2,r3,£4)...f2n(r2, .0 (3.3.22) .0. r n 7 .,rn) .f3n(r2,...,rn). .,r ) “J n)...fnn(r2,.. Given the Hessenberg matrix A in the form of (3.1.3) in which none of the ei's are zero (or close to zero), there exists an initial row vector R1 as in (3.3.21) leading to a unit upper triangular matrix R that transforms the matrix A to tridiagonal form. More specifically, all fi's can be made non-zero. Proof: Form the Krylov matrix F as in (3.3.22). It can be seen that f1 = 1 # 0. Using functional notation, it can also be seen that f2 = f2(r2,r3). The method proceeds as follows. For f2, assume temporarily that r3 = 0, and select r2 such that f2 is not zero. Once selected, consider r2 to be fixed from this point on; i.e., fij(r2,r3,...rk) be— comes fij(r3,...rk) in general. Following this line of thought, we see that f3 = f3(r3,r4,r5). Again assume that r4 = r5 =10, and select r3 such that f2(r3) # O and f3(r3) # 0. Note that at this stage there are four values Of r3 to be avoided since f2 is linear in r3, and f3 is cubic in r3. Continuing as before, we now assume that r3 33 In general, at stage i we are finding values of ri such that fj # O for j :.i where fj = fj(ri...). This is always possible since we are dealing with a matrix F of finite order; hence, there are only a finite number of values for ri to be avoided at each stage, and since for fi the term e2...eiri appears in which the ej's are non-zero. Therefore, all f.'s can be made non—zero. 1 Q.E.D. It is important to point out that the ei's may be completely arbitrary (including zero) and yet have the factorization work without breakdown. It is only necessary that one be able to define ri so that the fj's involving r. are not zero for j j,i. We shall examine the computa- i tional technique in order to motivate the second theorem. Computational Technique Assume that the algorithm for factoring the Krylov matrix F has computed the vectors R2,...,Rk, but breaks The action taken at down at stage (k+1) because fk+1 — 0. this stage depends on whether the breakdown occurs under Case 2 or Case 3. Whether or not a breakdown occurs, column Ck+1 is ' ..., accord- computed as a function of R1,..-,Rk and C1, Ck ing to the following lemma. Lemma- Let R be a unit upper triangular matrix of order n. —1 ' = can be computed Then the column Ck+1 of the matrix C R ' f R and recursively as a function of rows R1""’Rk o 34 columns C1""’Ck of C by the formula k c (3.3.23) iRi)Uk+1 Ck+1 = (U - 1:1 where Uk+1 is the (k+1)St column of the unit matrix U. Proof: Since C = R_1, then k n U = RC = CR = z c.R. + 2 c.R. . (3.3.24) i=1 1 J. i=k+1 l l The second summation in (3.3.24) defines an idempotent matrix W + where k 1 k n W = U - Z c.R. = z c.R. . (3.3.25) k+1 i=1 1 l i=k+1 l The last (n-i) entries of Ci are zero, and the first (i-l) entries of Ri are zero since R and C are both unit upper triangular. Hence, the matrix Wk+1 has the following form in which C appears as column (k+1): k+1 F6 no. 0 C1 ,k+1 * coo * . . . £ £ 0 . . O Ck,k+1 . — * * . 3.3.26) wk+1 — 0 ... 0 1 ( O .. 0 O 1 ... * 3 ... 3 O 0 .. H The (k+1)St column of Wk+1 iS Wk+1Uk+1 = Ck+1. Hence, (3.3.23) follows from (3.3.25). Q.E.D. The reason for introducing the above algorithm for computing the columns Ck+1 of the matrix C will become 35 apparent when we next discuss the actual technique of ap— plying the Krylov matrix factorization. Consider the Krylov matrix F to be initially partitioned as follows d1 FR1 F = F1 z , (3.3.27) FLl H1 where F is an (n—l) dimensional row vector, FL1 is an R1 (n-l) dimensional column vector, and H1 is an (n—l) X (n-I) submatrix. Since d1 is not zero (d1 = 1 since r11 = 1), then F1 may be factored into the product of three matrices F" -' — - r—l /d-) 1 0 d1 0 F 1 F1= R1 = L1D1V1 (3.3.28) F _Ll. 0 U d1 Un_1 0 Fa _ n_14 where lel - FLIFRI (3.3.29) — I F2 — d1 and Where R1 = [1 FRi/d1] is read from V1. Again, let F2 be partitioned as d2 FRZ F2: I FL2 H2 (3.3.30) If d2 = q2 # 0, then F2 may be factored as follows: — qr — 1 fl) d 1 0 d2 0 1 FRZ/ 2 = L2D2V2 (3.3.31) F2 = SEE. U 0 F3 0 Un_2 36" where R2 = [O 1 FRz/dz] from V; and where <3sz —- F, F 2 In general, at the (k+1)St stage, we have Fk+1 par- titioned as . T '— dk+1 FR,k+-1 19k+1 = (3.3.33) FL,k+1 Hk+i_ where we have assumed that di % O for i = 1, ---, k. Note that all the previous factorizations may be combined into a single factorization '—-— '- r ('1 : 0 d1 "ff: 0 d2 0 I--I u I . F=LDV= FLII I d1. u I"_F'_— r-"T-- IFLZI .n 1| 0 Ldk' o :‘7: :--:--- k I I l I—g—k—I Un-k 0 :0 ,Fk+1 TI _ 1 I FRI/d1 0'1 F d2 1......132 —————— f - (3.3.34) ._____..'__ :1 : FRk/d‘k r-i----~ o I U_ _ 0 ' I “k3 37; Now, assume that dk+1vu O. This is the first indica- tion of a breakdown. To determine whether it is a Case 2 or Case 3 breakdown, we must also examine the entries of FR,k+1 in (3.3.33). The first row of Fk+1 can be expressed in terms of the entries qi of the matrix T and the row Rk+1 of the matrix R as [0°'°Odk+1FR,k+1] = q2q3H°q1<+1Rk+1 ' (3'3°35) The product q2---qk is non—zero since we assumed H H 0 K1 di % O for i Hence, dk+1 = 0 implies that qk+1rk+1,k+1 = 0. In Case 2 where FR,k+1 is the zero vector, we choose qk+1 = O, and determine Rk+1 as the sum of the (k+1)St row of the idempotent Wk+1 and any desired linear combination of lower rows. This choice, like the initial choice of R1, is not unique. We then continue with the factorization as before, with the exception that for i > (k+1), [o---o di FRi] = qk+2°--qiRi (3.3.36) where di = qk+2---qi. It is noteworthy to point out at this time that if m is the degree of the minimal polynomial of A(m < n), then ‘st . there will always be a breakdown at the (m+1) stage. Since - ... n'1 A is upper Hessenberg, the Krylov matrix [C1,AC1, ,A C1] h diagonal entry e2e3°-°em. If = O. is upper triangular with m m < n, then (m+1) columns are dependent, and em+1 38 If however, Fle+1 has some non—zero entry, one cannot set qk+1 = O. This causes a Case 3 breakdown in which rk+1,k+1 = Rk+1ck+1 = 0. It is the occurrence of a Case 3 breakdown that motivates the second theorem. If a Case 3 breakdown occurs, we must revert to the last occurrence of a Case 2 breakdown, say at stage j (or to the start if no such breakdown occurred after j = 1), and attempt to select a new vector Rj with the usual con- straint that rji = O for i = 1, ---, j—l and rjj = 1. Under some circumstances, it is not possible to find a vector Rj for the given Hessenberg matrix A such that fk+1 % O, as the following example indicates. Example: Let the matrix A be as follows, and let R1 be represented as indicated in the first row of the Krylov matrix F. F0 0 1— F1 x y "l A = 0 o 1 , F = 0 o 1+x+y . (3.3.37) 1+ + _C O 1“J _P 0 x y Here it is obvious that f2 = O for every choice of R1. Under these circumstances, when an Rj cannot be found such that fk % O, a series of similarity transformations +1 are performed before the factorization is again used. These ' ' ' ' and transformations conSist of uSing the matrices Ra’ HR' H to be discussed below. The final matrix R may no longer C be unit upper triangular as will be shown later in equation (3 .3 .46). 39 Let Ra be a matrix formed from rows R1, --°, R of R k and the last (n—k) rows of U. Furthermore, let the matrix W be defined by W = RaARa (3.3.38) where R11 R12 Ra = - (3.3.39) 0 U F . l :--1-- lel W12 w: ...—£4}. = , (33%) l O W l _. .1 where T11 is tridiagonal, and where B = [bk+2°--bn] has at least one non—zero entry, which appeared in a different form in F in equation (3.3.33). Note that the reason R,k+1 W21 = 0 is that had there been an ek+1 in the upper right corner, then Theorem 1 indicates that there was an Rj which would not have led to this Case 3 breakdown; hence, the assumption that W21 = O is valid. The matrix HR is a Householder matrix of the form U 0 H = I (3.3.41) R 0 H ...R of B Which is used to collapse all the non-zero entries into q and to force bi = O for (k+2).: i.: n; see equation 40 (3.3.43). Referring to Section 3'2"§R is row-determined rather than being column-determined. Letting the matrix W' be defined by ' * 3342 w — HRWHR , ( . . ) then W' has the form: r ' ‘ T11 ' 0 I l-'-l_'-'— I l q l O Wil w12 w'= —.._,___-.. = . (3.3.43) I 0 W22 0 I W52 ' .1 Finally, the transformation defined by A' = ch'Hé , (3.3.44) in which HC is a column-determined Householder matrix, has the effect of transforming W52 into upper Hessenberg form. Hence, the form of A' is r. __ -1 I 0 , 0 [_ I 0 i l I ' Tll' I Til [ I_-I__-_ ' I I ' Iq IO ' A' = 0 I t '.E = ' E. . ___L._.J——-— ——--|--"--"-"(3..3,45) I I ' e ' I e I O I; IAéz '3' A22 I . I l 6 I _ ,0, y _ . . I Where the tridiagonal block Til is (k+1) x (k+1) and A52 is in upper Hessenberg form. 41 The result of theae three transformations is that the order of the leading tridiagonal block has been increased by one. If e = O in (3.3.45), then A' has the form of W in (3.3.40), and the above process would be repeated if.§ # 0. However, if §.= 0, then A' could be repartitioned to in- crease the size of the tridiagonal block by one. If e # O in (3.3.45), then Theorem 1 shows that there is an initial row vector R1 for which the Krylov factoriza— tion again applies. The following theorem has just been proved. Theorem 2: When a Case 3 breakdown occurs at stage (k+1) of the factorization defined by (3.3.27) through (3.3.33), it is possible to progress to at least the (k+2)nd stage in three steps: (1) transform A into RaAR;1 = W where Ra is constructed from R as described above, (2) transform W into HRWH; = W' by row—determined Householder transformations with product HR’ (3) transform W' into HCW'HE = A' by column-determined Householder transformations with product HC. The above transformations can be abbreviated by letting 3=HHR: (3.3.46) hence, A' =3Ag . (3.3.47) 42 The procedure described above transforms an arbitrary real n by n matrix into tridiagonal form without requiring a knowledge of the eigenvalues. The operations are linear except for the square root extractions in the Householder transformations. The rare exceptional cases are fully ac— counted for. IV SCALING In order to solve a set of linear differential equations on the analog computer, it is necessary that the variables lie within the operating range of the computer. Time, ampli- tude, and initial condition scaling transform the problem variables to solution variables amenable to the analog com— puter. The procedure for time and amplitude scaling of the tridiagonal system of equations described in this chapter is unique in that for the first time, an explicit method for scaling is presented which does not rely on the trial and error method used in the past. This new procedure is facilitated by the system of equations being in tridiagonal form. Let 51 and 52 (O < 51 < 52) be the low and high values of the operating range. 51 is determined by the tolerable drift of the operational amplifiers while 52 is determined by the capabilities of the read—out devices associated with the analog computer. As will be shown in Section 4.2, ampli- tude scaling does not affect the diagonal entries of T. Furthermore, all terms dleiqil are unaffected by amplitude scaling; only the individual entries ei and qi are affected. Therefore, the time scale factor B is used to transform all 43 44 diagonal entries and all.JTezq:TD's to values that are less than or equal to 82 in magnitude, and amplitude scaling places the individual off—diagonal entries of the time scaled matrix (BT) within the operating range. Finally, initial condition scaling ensures that, for stable systems of equations, the state variables will also lie within the operating range. 4.1 Time Scaling In terms of the elements ti’ ej, and qj of the tri— diagonal matrix T, let the scalars b1 and b2 be defined by II B H. :3 7—: H. + (D U. 3 LJ. b1 b2 = max(]ti],+-d|ejqjl) for i - 1,---,n and j = 2,---,n. The effect of the time (4.1.1) scale factor 5 is to place b2 at the high end of the Oper- ating range; i.e., 5 is computed by the formula _ E2 5 - “13-2— . (4.1.2) The resulting inequality then becomes 0 f. 61:31 1 (3b,, = 82 (4.1.3) where either 0 lqil c3 3 E dieiqil . It is now necessary to specify values for 51 and 52 which dePend on the analog computer. For the Applied Dynamics AD‘4, let 51 = 0.1, 52 = 100, and let aj be defined by '—1 a. = 103 51 (4.2.6) for j = 1,...,4; i.e., d1 = 61 and 04 — E2° 47 The following table lists the scale factor ratios in th f . . ’ e orm (kl/k1_1) for c1e and c1q(c2q and Cze) depending on the value of c3. ki/ki-1 if if C1 = C1e c1 = Ciq C3 :-01 quil/Gi 01/5191) (11 < C3 :N/Qlaz al/B‘ei‘. B‘qil/Ql [/02 az/Bjei] (12 < C3 :Vazaa Qz/Ble-l quil/Gz ~h12a3 < c3 i.a3 BIQ-l/Ga 03/5191) (13 < C3 :Va3a4 Q3/filei‘ B‘qil/Qa 'V C1304 < C3 ._<.. (14 quil/a4 a4/Blei‘ Although the foregoing scheme appears complicated, it iS quite easy to program and ensures that (n—l) of the xi's and 21's will have magnitudes equal to aj. Furthermore, each xi and zi equal to a2, o3, or a4 in magnitude reduces by one the number of potentiometers needed. If either e. or qi is zero, but not both are zero, then we will have c1 = O and c2 equal to 6 times the magni— tude of the non-zero entry. The ratio (kl/ki-l) will be defined by 48 qui . k. Q1 1f ‘qil # O l = (4.2.7) ki-i 1 B‘ell if lei] # 0 l = 1. (4.2.8) Now that we have (n—l) ratios of the form (ki/ki 1), we solve for the ki's by letting k1 = 1 and compute the re- ma' ' ' ining ki s by ki = (ki/ki_1)ki_1 (4.2.9) for i = 2,-°-,n. Once the ki's are known, we compute T' by (4.2.2) . Any parameter that is transformed to a value less than 0.1 in magnitude can be replaced by zero due to the limita- tions in accuracy of the analog computer. The reason for this is that the circuitry for this value would consist of a potentiometer, set for the value, coupled into a gain of one amplifier whose output would be too noisy to be signifi- cant or reliable. Consequently, there are some problems that cannot be solved on the analog computer. Once again, in terms of a change of variables, let 2(1) be defined by Z(T) = KY(T). (4.2.10) The system of equations (4.1.6) then becomes d -1 -Z(r) = K(BT)K 2(1). (4.2.11) dT 49 or in terms of the original matrix A. we have -5%z(w) = fiSAS—JZ(T) (4 2.12) where S = KRH. (4.2.13) 4.3 Initial Condition Scaling Initial condition scaling ensures that the dynamic range of the state variables, as solved on the analog com— pUter, lies within the operating range by imposing a rela— tionship between one volt and one unit of solution. At this point, there appears to be no well defined method for scaling the initial conditions. since the behavior of the state variables also depends on the system eigenvalues. One possible method for scaling the initial conditions is to use an iterative approach which is dependent on the length of the time interval over which the solution is de- sired. Let the scalar c4 be defined by c, = max )zi(0)l (4 3.1) . . . , y for i = 1,°‘~,n. Starting With an initial \alue for the scale factor a as (“I 2 : (4.3.2) Q 2 4 I 0 . attempt to solve the system of equations. If no saturation . _ _ _ . '.'+' Of the operational amplifiers occurs, then the lnlulal choice is sufficient. If however, saturation does occur decrease o by a small amount and try again. 50 The actual set of equations to be solved by the analog computer then becomes . —1 2% [oz m] = 8 SAS [82(7)], (4.3.3) or under the change of variables Z'(T) = QZ(T), we have %:g z' (T) = B 8218* z'(1). (4.3.4) V PROGRAM AND EXAMPLES A discussion of the digital computer program is pre- sented in Section 5.1 followed by three examples in Section 5.2. The material concerning the program is not intended to be comprehensive, but rather to show the program struc— ture in general, as related to the tridiagonalization and the scaling. The examples are intended to show numerically how the program implements the theory of this thesis. 5.1 Program Flow The logical flow of information in the digital computer PrOgram is indicated in the following two outlines where each level of indentation represents a decision or action taken as a result of the decision. These outlines show only the major operation involved and are not intended to be comprehensive. The first outline indicates the logic for the tridiagonalization, and the second outline indicates the logic for the scaling. Furthermore, each similarity transformation is combined with S and 8-1 when the transformations are generated; i.e._ S is replaced by RS. and S_1 is replaced by S-lc. These combinations are not indicated in the outlines. 51 52 Tridiagonalization At this time it is important to point out that an ex— haustive search for the proper initial row vector R1 is not practical, in time and prOgramming, unless necessary for an extremely rare occurrence of a specialized Case 3. There is at least one other method, although not foolproof, which has been demonstrated to be worthwhile for dealing with the Case 3 breakdown. The procedure for implementing this alternate method will now be discussed; this alternate method was used in the program that computed the examples in Section 5.2. Let j be the stage of factorization in which the most recent Case 2 occurred. After exhausting the limited choices for R. in attempting to progress past the stage (k+1) in which the Case 3 breakdown occurred, the follow— ing three transformations are performed: - ‘1 5.1.1' w - RaARa . ( ) w3 = PWP-l, (5.1.2) and A“ = H'W'(H')*. (5.1.3) The matrix R is formed from the first (j-l) rows of ' a R and the last (n-j+1) rows of U. It can be shown that W has the form (5.1.4) 53 where T11 is tridiagonal, and A28 is zero except for the entry ej in the upper right corner. Let E be the column in F of (3.3.33) in which the R,k+1 first non—zero entry occurred with the most recent Rj (now denoted Rj). The matrix P is the unit matrix except for the entries pfii’ i = k+1,°--,£-1; these entries are -1's. The matrix W' then has the form: w'= WP = . (5.1.5) 0 3 A21 W52 The effect of this transformation is to attempt to make fk+1 # 0 when the Krylov factorization, using_R_j to start with, again reaches the (k+1)St stage. Since W52 in (5.1.5) is not necessarily in upper Hessenberg form, then the Householder transformation H' in (5.1.3) is applied so that A52 in (5.1.6) is in upper Hessenberg form. (5.1.6) 0 H' A21 A52 ‘ The Krylov factorization is again attempted treating Ag? in (5.1.6) as A and 33 as R1. The resulting R matrix is now in the form 0 lw where the first row of 3 is Bj- H1 L‘J 54 For a given n and i, 1e: RT = [r. l l°--rn] be a row vector of length k = n-i+1 whose entries ri = 1 and rj = 0 or 1 for i < j i.n. Ri can be considered to be the binary representa— tion of Bi where 2k_1 j,£i j_2k»1. Example: If n = 5 and i = 2, the possible values of £1 lie in the range 23 = 8 j'fli‘: 24-1 = 15. The digital program for the tridiagonalization uses the above values for Ri’ when Ri is the result of a Case 2, starting at the low end of the range and ranging through all binary numbers to the high end of the range before it performs the partial transformation at a subsequent Case 3 stage. It is conceivable that one of these possible Ri's might not cause a Case 3 to occur, although a previous Ri did. In the outline for the tridiagonalization,the follow— ing symbols (combinations of letters; i.e., THS is a matrix name) are used to denote matrices: A, H, H*, B, THS, R, W, C, F, RA, CA, WT, FT, and Z. With the exception of B, THS. W, WT, FT, and Z, the remaining symbols parallel the ma rices referred to earlier in the thesis. B, THS, WT, FT, and Z are used for temporary matrix storage, and W 15 used to store the successive idempotents as they are computed. Since . T . A is assumed to be real, H* is actually H . Furthermore, R is row I of R F I is row I of the factored Krylov matrix I ’ R F (see equation(3.3.33)with k+1 replaced With I)' CI ls column I of C, and W is the idempotent. 55 As scalars, the symbols N, I, IT, and JT are used where N is the order of the matrix, I is a stage of factorization in which R1+1 is computed, IT + 1 is the most recent Case 2 row in the factorization, and JT (normally equal to zero) is set to a one when there is a Case 3 at a certain stage with all the limited choices for the prior Case 2 RIT+1'S being depleted. Outline (10) Given the matrix A and order N Perform Householder transformation (HAH* = B) Set THS = B (save Hessenberg matrix B) Compute R1 from binary vectors Compute C1 and W2 Compute Krylov matrix F Set I = 0, IT = 0, JT = 0 (I is the current stage of factorization, IT + 1 is the most recent Case 2 stage, JT is changed to one when a Case 3 has occured with all possible choices of RIT+1 depleted). Go to (c) Does I = N—l? Yes Go to (h) No Set I = 1+1 Factor Krylov matrix in place Test cases Case 1 56 (d) Set R : 1+1 FR,I+1/dI+1 Com ute P CI+1 and W1+2 Go to (b) Case 2 SetRA=R,CA=C,WT=W,andFT=F(save present R, C, W, and F) - Set IT = I and JT = 0 (e) Compute RI+1 from binary vectors Go to (d) Case 3 IT Set I Is JT 0? Yes I Is IT = 0? Yes ' Are all binary vectors for R1 used? Yes ' Set JT = 1 Go to (f) No Go to (a) No ' - o Are all binary vectors for RIT+1 used, Yes Set JT = 1 Compute (RA)(THS)(CA) = B Go to (g) 57 No SetR=RA,C=CA,W=WT,andF=FT Go to (e) No (f) Set B THS (restore upper Hessenberg or Partial tridiagonal matrix) (9) Compute PBP = Z Perform Householder transformation HZH* = B Set THS = B (save partial tridiagonal matrix) Compute Krylov matrix using most recent RIT+1 Go to (d) (h) Compute Rb(THS)Cb = T Scaling In the outline for scaling, T is used to denote the tridiagonal matrix. B, Q, C, and RK are used to denote vectors where E is the subdiagonal, Q is the superdiagonal, C is the geometric mean of the sub and super diagonals, and RK is initially the ratio of the amplitude scale factors and lastly the vector of amplitude scale factors. The scalars E1, E2, B2, BETA, and N correspond to 51, 92, b2, 6, and n in Chapter IV; Outine Given T, N, E1, E2 Set EI lt1,I—1l‘ 01 = ltI_1,1 : and C1 =:JE;6; for I = 2,°°°,N 58 Find 82 = max (‘tI,I" CI) Compute BETA = E2/B2 (time scale factor) Compute (BETA)T = T (time scaled) Compute (BETA)EI = EI, (BETA)QI = Q for I = 2,---,N I I, and (BETA)CI = c I ki/ki-i for l = I = 2'°",N Given k1 = 1, compute ki in (RK)I for i = I = 2,°°°,N Compute (RK) Compute (RK)('I‘)(RK)-1 = T (time and amplitude scaled) 5.2 Examples Three examples are given in this section to illustrate the tridiagonalization and the time and amplitude scaling procedures that were presented in Chapters III and IV. Example 5.2.1 illustrates the Householder transforma— tion which transforms an aribitrary real non—symmetric matrix A to upper Hessenberg form AH where, AH = HAHT. (5.2.1) Examples 5.2.2 and 5.2.3 illustrate the tridiagonali- zation and the scaling of an upper Hessenberg matrix AH where T' = fiSAHS-l. (5.2.2) Example 5.2.2 illustrates Case 2 operations, and Example 5.2.3 illustrates Case 3 operations. The IBM 1800 computed these examples using extended precision which is thirty-one bits of mantissa or roughly nine decimal places; however, printed answers are rounded to two or three decimal places. In these examples, 0_j is 59 used to denote a number whose characteristic is 10-j (0 < j); i.e., small order of magnitude. Zeros are shown where an exact zero was produced by the computer. Also, the symbol (—>) is used meaning "is factored into". The Krylov matrix F is factored in place, where F = LDV is stored as (L-U)+D+(v-U). \ LDV \\D\ \ '11 ll (5.2.3) Example 5.2.1: This example illustrates the Householder transformation of a real non-symmetric matrix A to upper Hessenberg form AH. A = 2 —2 2 5 _2 -1 —4 9 AH = HAHT (Householder transformation) f— -1 1.000 0 0 0 0 -0.333 —0.667 —0.667 H = o -o.133 -O.667 0.733 ' 0 -0.933 0.333 0.133 P 1 4.0 —1.0 3.0 1.0 —3.0 1.0 —4.4 -0.8 AH = 5.0 0 —2.0 _ 0.8 -1.0 4.9) 60 Example 5.2.2: This example illustrates the tridiagonali- zation procedure with a Case 1 followed by two Case 2 fac- torizations of the Krylov matrix F. It also illustrates the results of the time and amplitude scaling procedure. F- — 1 2 3 4 0 0 0 O A = H o 1 3 0 _P 0 0 i (Krylov factorization) F- _p F_' _[ A__ (_ 1 1 O O 1 1 0 0 1 1 O O 1 2 3 4 1' 1 3 4 1 1 3 4 F = T> : §> I—--——— , 1 5 12 16 1 , 4 12 16 1 4 I o o I l 1 14 39 52 1 :13 39 52 1 13 I 0 9 _ _ r 1 _ 1 1 1 0 O 1 1 0 0 1 1 0 0 1 1 3 4 1 1 3 4 1 1 3 4 _____ §> I 1 4 :1 o 1 4 1 0 1 4 1 o I l‘" _1 13 IO (1| _1 13 0 l 9.1 _1 13 O H T = RA C (Tridiagonalization) -H r- — F1 1 O 0 1 1 0 O 0 1 3 4 0 3 0 O R = I T = 0 0 1 0 0 1 0 -4 O 0 O 1 0 0 0 3 S = 33.33 61 T' = K(BT)K-1 (Scaling) F1.0 0 0 0 F33.3 0.1 0 0 0 333.3 0 0 0 100.0 0 0 K = , T' 0 0 1.0 0 0 0.1 0 -0.1 __0 0 0 1333.3) 0 o 0 100.@ S = KR (Transforming matrix) [1.0 1.0 0 0 0 333.3 104 1333.3 S: 0 0 1 0 __0 0 0 1333.3 Example 5.5.3: This example illustrates the tridiagonali- zation procedure with a Case 1 followed by a Case 3 factori— zation of the Krylov matrix F. Instead of changing the initial row vector R1, this vector is kept so as to illus— trate the use of the P matrix in (5.1.5) followed by a Householder transformation. A new Krylov matrix F' is computed using R1 and is factored with three Case 1's occurring. i 1 1 0 2 1 1 1 AH = 0 1 1 0 .9 0 1 U (Krylov factorization) 1.1 0 0 0 r1 0 0 0 [-1 0 0 | ——————— 1 1 1 0 1 l1 1 0 1 1 1 F = _1_> I _> r-— 3 3 3 1 3 I3 3 1 3 3 3 ,0 l l 9 9 10 I I _ ‘1) f .9 10 fl _9 9 I 1 —1 w' = PAHP F-l 0 0 0 [-1 1 1 01 0 1 0 0 3 2 2 1 P = t W' :: 0 0 1 0 0 1 1 O __-1 -1 —1 _1_J -3 —3 -2 0) A' = H'W'(H')T (Householder transformation) F1.00 0 0 0 1 [—1.00 -0.71 0_9 0 -0.71 0 0.71 -4.24 2.00 -0.82 I'Il = , A1 = 0 0.58 -0.58 0.58 0 1.22 0—9 0 0.41 0.82 0.41J 0 0__9 0 (Krylov factorization) [-1.00 0 0 0 1 1.00 0 0 r' _______ 1.00 -0071 O_9 1.22 1000:-Oo71 0_9 F = -> 4.00 -2.12 0.58 5.31 1 4.00{-2.12 0.58 I 13.00 -6.36 1.73 18.32} 13.00:-6.36 1.73 62 1.227 -4.04 -0.71 1.00 63 r - - 1.00 0 0 0 1.1.00 0 0 0 1.00 1.00 0_9 -1.73 1.00 1.00 0_9 -1.73 ————————— —> 4.00 3.00 :0.58 1.63 1 4.00 3.00 1.00 2.83 I _____ 13.00 9.00 :1.73 7.35J 13.00 9.00 3.00 {2.45 ... '" l . -1 _ , . . . . T — RbA Rb — RbA Cb (Tridiagonalization) . -r [— .— 1.00 0 0 0 1.00 -0.71 0-13 0 l 0 1.00 0_9 —1.73 —4.24 2.00 0.82 0_.18 Rb: I T: 0 0 1.00 2.83 0 1.22 0.9 4.24 __ 0 0 0 1.00 __ 0 0__9 0_18 1.00 8 = 50.0 T' = 1<(s'1*)1<'1 (Scaling) I—‘ " '— ‘— 1.00 0 0 0 50.0 75.0 0.16 0 0 0.47 0 0 -100.0 100.0 25.0 0_10 K = I T': 0 0 0.77 0 0 100.0 0_8 0.1 S = KRbH'P (Transforming matrix) F' 1.00 0 0 0 _ 0.9 —0.67 —0.57 0_,9 1.33 0 0.9 1.33 -666.67 0 666.67 666.67 These three examples have been presented in an attempt to illustrate the tridiagonalization and the time and amplitude scaling procedures of Chapters III and IV. VI CONCLUSIONS In this chapter, a summary of the original results is presented followed by an indication of some additional areas of study which have presented themselves through the develop- ment of this thesis. In many engineering problems, the need to solve a sys— tem of simultaneous homogeneous linear differential equa- tions with constant coefficients arises. This thesis showed how this system of equations, Sci-{X(t) = Ax(t), in which X(t) is an n—dimensional vector function of time and A is an (n x n) real square matrix, can be solved on the hybrid computer using no more than 2n operational amplifiers (n bipolar integrators) nor more than (Zn—1) potentiometers, as opposed to the more conventional means of solution using n2 potentiometers, on the analog computer. By tridiagonalizing the matrix A, the former n2 pos- sible number of non-zero parameters is reduced to at most (3n-2). Since one potentiometer is usually needed for each non-zero entry, by appropriately scaling the tridiagonal matrix, the number of potentiometers needed is further re— duced to (Zn—1). 64 65 In case automatic patching is desired, usually ac— complished using one reed relay for each non-zero entry in the coefficient matrix, a considerable savings in the number of relays needed is affected by tridiagonalizing the coef— ficient matrix. A new procedure for tridiagonalizing an arbitrary real non-symmetric matrix A, without using the eigenvalues, was presented in Chapter III. After transforming the matrix A by unitary Householder transformations into a similar upper Hessenberg matrix A the latter is transformed into a HI tridiagonal matrix T = RAHC by appropriate matrices R (with rows Ri) and C = R“1 (with columns Ci)' The rows R1 are computed recursively by factoring the Krylov matrix F whose rows are the interates of R1 under A The columns H. Ci are computed recursively from an idempotent matrix i-1 W. = U - Z C.R. by C. = W.U. where U. is the 1th column 1 j=1 3 3 1 1 1 1 of the unit matrix. In Section 3.3, Theorem 1 was proved which showed that the tridiagonalization is always possible without breaking down when AH has all non-zero subdiagonal entries. If AH has at least one zero subdiagonal entry, the hypothesis of the theorem is invalidated in which case the procedure may or may not break down. .When the procedure does break down, it does so under either Case 2 or Case 3. Case 2 is handled by obtaining Ri from the idempotent Wi_1- For Case 3. Theorem 2, proved in Section 3.3, shows how to advance at least one additional stage. 66 A lemma was also proved in Section 3.3; this lemma showed a new method for inverting a unit upper triangular matrix. For the first time, an explicit procedure for time and amplitude scaling was presented. This procedure, pre- sented in Chapter IV, is facilitated by the system of equations being in tridiagonal form, and ensures that no more than (Zn—1) potentiometers are needed on the analog computer. This thesis will be concluded by indicating some areas in which additional investigation could be done. It would be useful to extend this solution procedure to en- compass other forms of differential equations such as those that are non-linear. Also, the scaling of the initial conditions could perhaps be accomplished by means other than the method described in Section 4.3. For instance, the feasibility of using time—varying initial conditions over the desired interval of solution, in the form of a series of steps or as a continuous function. could be in- vestigated. REFERENCES Causey, R. L., and Gregory, R. T. "On Lanczos‘ Algorithm for Tridiagonalizing Matrices," Society of Industrial and Applied Mathematics Review, III, No. 4 (October, 1961), 322-328. Faddeyev, D. K., and Faddeyeva, V. N. Computational Methods of Linear Algebra. English translation. San Francisco: W.H. Freeman and Company, 1963, 241-247. Frame, J. S. Matrix Theory and Linear Algebra. Un— published book. Kublanovskaya, V. N. "Reduction of an Arbitrary Matrix to the Tridiagonal Form," U.S.S.R. Computational Mathematics and Mathematical Physics, IV, No. 3(1964), 202—203. . LaBudde, C. D. "The Reduction of an Arbitrary Real Square Matrix to Tri—Diagonal Form Using Similarity Transformations," Mathematics of Computation, VII (1963), 433-437. Parlett, Beresford. "A Note on LaBudde‘s Algorithm," Mathematics of Computation, XVIII (1964), 505—506. Strachey, C., and Francis, J. G. F. ”The Reduction of a Matrix to Codiagonal Form by Elimination,“ Computer Journal, IV(1961-1962), 168-176. Waltman, William L., and Lambert, Robert J. "T— Algorithm for Tridiagonalization," Journal of Society of Industrial and Applied Mathematics, XIII, No. 4 (December, 1965), 1069-1078. Wang, H. H., and Gregory, R. T. "On the Reduction of an Arbitrary Real Square Matrix to Tridiagonal Form,“ Mathematics of Computations, XVIII(1964), 501—505. 67 10. 11. 68 Wilkinson, J. H. "The Calculation of Eigenvectors by the Method of Lanczos," Computer Journal, I (1958), 148-152. Wilkinson, J. H. "Instability of the Elimination Method of Reducing a Matrix to Tridiagonal Form," Computer Journal, V (1962-1963), 61-70. “1"!leth 3