# solution techniques for large eigenvalue problems in structural

## Comments

## Transcription

solution techniques for large eigenvalue problems in structural

'±~ctA ..... Lt/'. 'J.- ~ UILU-ENG-79-2006 ~.3 , CIVIL ENGINEERING STUDIES STRUCTURAL RESEARCH SERIES NO. 462 SOLUTION TECHNIQUES FOR LARGE EIGENVALUE PROBLEMS IN STRUCTURAL DYNAMICS II Metz Reference Room Civil Engineering Department BI06 C.E. Building University of Illinois Urbana" Illinois. 6180li By I.-W. LEE A. R. ROBINSON A Technical Report of Research Sponsored by THE OFFICE OF NAVAL RESEARCH DEPARTMENT OF THE NAVY Contract No. N00014-75-C-0164 Project No. NR 064-183 Reproduction in whole or in part is permitted for any purpose of the United States Government. Approved for Public Release: Distribution Unlimited UNIVERSITY OF ILLINOIS at URBANA-CHAMPAIGN URBANA, ILLINOIS JUNE 1979 50272-101 REPORT DOCUMENTATION PAGE 3. Recipient'. Acce••lon No. 2. 1_1./;REPORT NO. TTTT,TT ENG-79-2006 5. Repo·rt Date 4. Title and Subtitle June 1979 Solution Techniques for Large Eigenvalue Problems in Structural Dynamics - - - - - - - - - - - - - - - - ----------_._--------..------7. Author(s) 8. Performing Organization Rept. No. In-Won Lee and Arthur R. Robinson 9_ Department of University of Newmark Civil 208 N. Romine TT."..1-..... 12. SRS No. 462 Performing Organization Name and Address 10. Project/Task/Work Unit No. NR 064-183 Civil Engineering Illinois at Urbana-Champaign Engineering Laboratory Street T1...,; ...... ,....;c 11. Contract(C) or Grant(G) No. (C) N00014-75-C-0164 (G) h1Q(\1 SPo~~;ri'j;g' Organiz;ti~~ "N';me ~';d Address 13. Type of Report & Period Covered Material Science Division Structural Mechanics Program (code 474) Office of Naval Research (800 Quincy Street) Arlington, Virginia 22217 Technical Report r---------------------------~ 14. 15. Supplementary Notes --------------------_. --- --------------------1 ------------ 16. Abstract (limit: 200 words) This study treats the determination of eigenvalues and eigenvectors of large algebraic systems. The methods developed are applicable to finding the natural frequencies and modes of vibration of _large structural systems. For distinct eigenvalues the method is an application of the modified NewtonRaphson method that turns out to be more efficient than the standard competing schemes. For close or multiple eigenvalues, the modified Newton-Raphson method is generalizec to form a new process. The entire set of close eigenvalues and their eigenvectors are found at the same time in a two-step procedure. The subspace of the approximate eigenvectors is first projected onto the subspace of the true eigenvectors. If the eigenvalues are multiple, the results of the first stage indicate this fact and the process terminates. If they are merely close, a single rotation in the newly found space solves a small eigenvalue problem and provides the final results for the close set. The procedure for subspace projection can be expressed as a simple extremum proble~ that generalizes the known extremum property of eigenvectors. Computational effort and convergence are studied in three example problems. method turns Out to be more efficient than subspace iteration. The -----------------------------------------------------------------~ 17. Document Analysi5 a. Descriptof"l Eigenvalues Eigenvectors Numerical Analysis Dynamic Structural Analysis b. Identifiers/Open-Ended Terms c. COSAT! Fie!d/Group 19. Security Class (This Report) 18. Avaiiability Statement Approved for public release: Distribution unlimited Unclassified 21. No. of Pages 118 ~--------------~---------------20. Security Class (This Page) 22. Price Unclassified See ANSI-Z39.1S, See Instructions on Reverse OPTIONAL FORM 272 (4-77) (Formerly NTIS-35) Department of Commerce SOLUTION TECHNIQUES FOR LARGE EIGENVALUE PROBLEMS IN STRUCTURAL DYNAMICS by I.-W. Lee A. R. Robinson A Technical Report of Research Sponsored by THE OFFICE OF NAVAL RESEARCH DEPARTMENT OF THE NAVY Contract No. N00014-75-C-0164 Project No. NR 064-183 Reproduction in while or in part is permitted for any purpose of the United States Government. Approved for Public Release: Distribution Unlimited University of Illinois at Urbana-Champaign Urbana, Illinois June 1979 iii ACKNOWLEDGMENT This report was prepared as a doctoral dissertation by Mr. In-Won Lee and was submitted to the Graduate College of the University of Illinois at Urbana-Champaign in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Civil Engineering. The work was done under the supervision of Dr. Arthur R. Robinson, Professor of Civil Engineering. The investigation was conducted as part of a research program supported by the Office of Naval Research under Contract N00014-75-C0164, "Numerical and Approximate Methods of Stress Analysis. II The authors wish to thank Dr. Leonard Lopez, Professor of Civil Engineering, for his assistance. The numerical results were obtained with the use of the CYBER-75 computer system of the Office of Computer Services of the University of Illinois at Urbana-Champaign. iv TABLE OF CONTENTS Page 1. INTRODUCTION...... 1 . 1 Genera 1 . . . . . . 1.2 Object and Scope . . . . . 1.3 Review of Solution Methods. 1 . 4 No ta t ion . . . . . . . . . . . . . . . 2. 4. 6. 10 10 14 17 General . . . . . . . The Iterative Scheme . . . . . . . . . Convergence Rate and Operation Count . Errors in Approximate Eigenso1utions . Treatment of Missed Eigenso1utions 19 CLOSE OR MULTIPLE ROOTS . . . . . . . . . 21 3.1 3.2 3. 3 3.4 3.5 21 22 25 30 31 General . . . . . . . . . . . . . . . . . . . . . . . . . . Theoretical Background. . . . . . . . . . . . . . . . . Th e I t era t i ve Schern e . . . . . . . . . Treatment of Close Roots . . . . . . . Convergence Rate and Operation Count 34 APPROXIMATE STARTING EIGENSOLUTION . 4. 1 Genera 1 . . . . . . . . . . . . . 4.2 Subspace Iteration Method . . . . . . 4.2.1 The Iterative Scheme. 4.2.2 Starting Vectors. . . . 4.2.3 Convergence Rate, Operation Count, Estimation of Errors. . . . . . . 4.3 Starting Solution for the Proposed Method 5. 10 DISTINCT ROOTS . . 2. 1 2.2 2.3 2.4 2.5 3. 1 2 3 6 . . . . . . . . . and . 34 35 35 37 38 41 NUHER I CAL RESULTS AND COMPARISONS . . . 43 5. 1 5.2 5.3 5.4 5.5 43 44 45 46 Genera 1 . . . . . . . . . • . . . • . • . . . . . . . . •. Plane Frame. . . . . . . . . . . . . . . . . . . . . . . . Arch . . . . . . . . . . . . . . . . . . . . . . . . . . . . Plate Bending. . . . . . . . . . . . . . . . . . . . . . . Comparison between the Theoretical Convergence Rates and Numerical Results . . . . . . . . . . . . . . . . . . . . . SUM~ARY 6. 1 6.2 6.3 AND CONCLUSIONS . . . . . . . Summary of the Proposed Method Conclusions . . . ... Recommendations for Further Study . . . . . . . . . . . . . 47 49 49 50 51 v Page LIST OF REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . 52 APPENDIX A. NONSINGULARITY OF THE COEFFICIENT MATRICES OF THE BASIC EQUATIONS . . . . . . . . . . . . . B. CONVERGENCE ANALYSIS . . . . . . . . . . . . . . . . . . . . 61 B.1 B.2 C. 57 Case of a Di·stinct RDot Case of a Multiple Root 61 67 THE BASIC THEOREMS ON THE CONSTRAINED STATIONARY-VALUE PROBLEM . . . . . . . . . . . . . . . . . . 77 vi LIST OF TABLES Page Table 1. NUMBER OF OPERATIONS FOR EIGENSOLUTIONS . . . . . . . . 2. EIGENVALUES OF THE PLANE FRAME PROBLEM (DISTINCT ROOTS) . . . . . 90 3. EIGENVALUES OF THE CIRCULAR ARCH PROBLEM (DISTINCT ROOTS) . . . . 92 4. EIGENVALUES OF THE SQUARE PLATE PROBLEM (DOUBLE ROOTS) . . . 94 5. EIGENVALUES OF THE RECTANGULAR PLATE PROBLEM (CLOSE ROOTS) . 95 6. COMPARISON OF THE TOTAL NUMBER OF OPERATIONS 96 7. COMPARISON BETWEEN THE THEORETICAL CONVERGENCE RATES FOR EIGENVECTORS AND THE NUMERICAL RESULTS - FRAME PROBLEM (DISTINCT ROOTS) . . . . . . . . . . . . . . . . . . . . . . . . 97 8. COMPARISON BETWEEN THE THEORETICAL CONVERGENCE RATES FOR EIGENVECTORS AND THE NUMERICAL RESULTS - SQUARE PLATE PROBLEM (DOUBLE ROOTS) . . . . . . . . . . . . . . . . . . . . . 98 9. NUMERICAL CONVERGENCE RATES FOR EIGENVECTORS RECTANGULAR PLATE PROBLEM (CLOSE ROOTS) . . . . . . . . . . . . . 99 . . 86 nee Room t neering Departmen C iv~l Engi C E BUl.l~·ng ell . Metz .t:' Re~ere Bl06 .~ . f Illinoia Univers~vY ~ 's 6180] Urbana" Illl.nol. - - vii LIST OF FIGURES Figure Page 1. ESTIMATION OF ERRORS IN APPROXIMATE EIGENVECTORS. 100 2. TEN-STORY, TEN-BAY PLANE FRAME. . . . . . . . . . 101 1 1. 1.1 INTRODUCTION General Various engineering problems can be reduced to the solution of matrix eigenvalue problems. Typical examples in the field of structural engineering are the problem of determination of natural frequencies and the corresponding normal modes in a dynamic analysis and the problem of finding buckling loads in a stability analysis of structures. Since the advent of the digital com- puter, the complexity of structures which can be treated and the order of the corresponding eigenvalue problems have been greatly increased. Hence, the development of solution techniques for such problems has attracted much attention. For the dynamic analysis of a linear discrete structural system by superposition of modes, we must first solve the problem of free vibration of the system. The free vibration analysis of the linear system without damping reduces to the solution of the linear eigenvalue problem Ax = A Sx (1.1) in which A and B are stiffness and mass matrices of order n, the number of degrees of freedom of the structural system. A column vector x is an eigenvector (or normal mode), and the scalar A the corresponding eigenvalue (or the square of a natural frequency). Th2 matrices A and B are real and symmetric, and are usually banded and sparse. If a consistent mass matrix is used, the matrices A and B have the 2 same bandwidth [4,5J. be diagonal. semidefinite. If a lumped mass model of the system is used, B will The matrix B is positive definite, but the matrix A may be There are n sets of solutions of Eq. (1.1), that is, n eigen- values and their corresponding eigenvectors. Frequently, in practical eigenvalue problems, the order of A and B is so high that it is impractical or very expensive to obtain the complete eigensolution. On the other hand, to carry out a reasonably accurate dynamic analysis of the structure, it is possible to consider only a partial eigensolution. The partial solution of interest may consist of only few lowest eigenvalues and their eigenvectors, or eigenvalues in the vicinity of a given frequency and the corresponding eigenvectors. The method described in this study is aimed at effective solution of this type of problem rather than at a complete eigensolution. 1.2 Object and Scope The object of this study is to present an iterative method which is efficient and numerically stable for the accurate computation of limited number of eigenvalues and the corresponding eigenvectors of linear eigenvalue problems of large order. The method developed remedies the major drawbacks of the inverse iteration method with spectral shifting [13]: numerical instability due to shifting and slow convergence when eigenvalues are equal or close in magnitude. The proposed method converges rapidly and is numerically stable for any number of multiple or close eigenvalues and the corresponding eigenvectors. 3 The procedure for distinct eigenvalues is treated in Chapter 2, and a modified procedure for multiple or close eigenvalues in Chapter 3. Selection of initial approximate eigenvalues and eigenvectors by the subspace iteration method is described in Chapter 4. To show the efficiency of the proposed method, three sample problems are solved: vibration of a plane frame, of a plate in bending, and of an arch. Comparisons are made in Chapter 5 with a method which is generally regarded as very efficient, the subspace iteration method. 1.3 Review of Solution Methods Numerous techniques for the solution of eigenvalue problems have been developed. These techniques can be divided into two classes - techniques for approximate solution and techniques for "exact solution. ll The approximate solution techniques include well-known static condensation [2,3,24,25,27,42J, dynamic condensation [34J, Rayleigh-Ritz analysis [9,13,31,48J, component mode analysis and related methods summarized by Uhrig [50J. These methods are essentially techniques for reducing the size of a system of equations. The reduction of a system of equations eventually leads to a loss in accuracy of a solution. However, the advantage of lessened computational effort for a solution sometimes may compensate for the loss in accuracy. Moreover, an approximate solution found by these methods may serve as the starting solution for the exact methods, which will be discussed next. The exact methods are designed for the accurate computation of some or all the eigenvalues and corresponding eigenvectors. These methods consist 4 of vector iteration methods, transformation methods, the method based on the Sturm-sequence property, polynomial iteration method, and minimization methods. These methods are well described in Ref. 51. The methods differ in the choice of which mathematical properties of an eigenvalue problem are used. The vector iteration methods such as the classical vector iteration (power method) and simultaneous vector iteration deal with the form of equations Ax = A Sx. The transformation methods (LR, QR, Jacobi, Givens, and Householder methods) are based on the mathematical property that the eigenvalues of a system are invariant under similarity transformations. In the polynomial iteration method, the roots of det (A - AS ) = 0 are found, and minimization methods are based on the stationary property of the Rayleigh quotient [43J. In vector iteration methods and minimization methods, both the eigenvalues and corresponding eigenvectors are found simultaneously, but in other exact methods, only eigenvalues are computed or the computed eigenvectors are, in general, not suitable for use in the final solutions. In such methods, another method such as the vector iteration method with a shift may be used for finding the eigenvector corresponding to a computed eigenvalue. For a limited number of eigenvalues and corresponding eigenvectors df an eigenvalue problem of large order which we are concerned with in this study, the above methods have been modified or combined to take advantage of the useful characteristics of several of the methods. First, the determinant search method [7,9,22,23J combines the methods based on the Sturmsequence property, polynomial iteration, and inverse iteration. In this 5 method, eigenvalues in a specified range are approximately isolated by using the bisection method and the Sturm-sequence property and then located accurately by the polynomial iteration method. The corresponding eigen- vectors are computed by inverse iteration with a shift. By this method, eigenvalues in any range and corresponding eigenvectors can be found. However, it has the disadvantage that the matrix is factorized in each iteration to locate the eigenvalues of interest. Another method for the solution of large eigenvalue problems is the so-called subspace iteration method [6,15,32,39,47J, which is a combination of the simultaneous iteration method and a Rayleigh-Ritz analysis. In this method, several independent vectors are improved by vector inverse iteration, and the best approximation to the eigenvectors are found in the subspace of the iteration vectors by a Rayleigh-Ritz analysis. In this method, eigen- values at the end of the spectrum and the corresponding eigenvectors converge very rapi dly. Thi s method wi 11 be d-i scussed further in Chapter 4. The inverse iteration method with a shift is known to be extremely efficient for improving approximate eigenvalues and eigenvectors. However, as mentioned in the previous section, when the shift is very close to a true eigenvalue, the method exhibits numerical instability, yielding unreliable answers [13J. In addition, when the eigenvalues of interest are close to- gether, their convergence is very slow. Robinson and Harris [44J developed an efficient method to overcome the above difficulty for distinct eigenvalues by augmenting the equations used in the inverse iteration method by a side equation. While this method extracts eigenvalues and eigenvectors simul- taneously with a very high convergence rate, it has the disadvantage that the 6 algorithm is inefficient for problems with multiple or close eigenvalues. This method and some improvements on it will be discussed further in the next chapter. 1.4 Notation All symbols are defined in the text when they first appear. With regard to matrices, vectors, elements of matrices or vectors, and iteration steps, the following conventions are generally used: (1) Matrices are denoted by uppercase letters, as A, Band X. A column vector is denoted by a lowercase letter with a (2) x.. superior bar and a subscript, as ~.- J" -b.J and -- --- - J Elements of a matrix or vector are denoted by a lowercase (3) 1 e t te r with ado ub1 e sub s c rip t s, a sa. ., b.. and x. . . lJ lJ lJ Iteration steps are denoted by a superscript, as x(k), (4) x~k) J and x (k) ... lJ Increments are denoted by the symbol ~, as ~x~k) and ~X~~). (5 ) J Some symbols are assigned more than one meaning. 1J However, in the context of their use there are no ambiguities. - - A, a·, a·· J lJ = stiffness matrix, jth column vector of A, element of A *(k) A " = projection of A onto the subspace spanned a = radius of circular arch B, 6., b .. J *(k) B lJ . = mass matrlx, J. th co 1umn vec t or 0f by vectors B,e 1emen t 0 f B = projection of B onto the subspace spanned by vectors in y(k), B*(k) = y(k)T B y(k) 7 = expansion matrix of X(k), jth column vector of C(k), element of C(k), X(k) o = XC(k) = diagonal matrix, see Section 2.2 = plate bending stiffness, De = EH 3/12(1-u 2 ) = matrix for finding close or multiple eigenvalues and 0, 0. J eigenvectors, jth column vector of 0, see Eq. (3.24) = iteration matrix for Dafter k iterations, jth column vector of D(k), see Eq. (3.23) = Young's modulus E E, e., e .. J JJ = diagonal matrix, jth column vector of E, element of E, see Eq. (A.7) -* e.· -* E* , e., = h = thickness of plate J JJ diagonal matrix, jth column vector of E* , element of E*, E* = _E- l = number indicating rate of convergence of eigenvector, see Eq. (2. 13) I = moment of inertia of cross-section = identity matrix of order s i, j = indices of matrix elements k = L = lower triangular matrix superscript indicating number of iterations = Lagrangian, see Eq. ma , mb = (3.6) average half bandwidth of A, "oF R VI U = total number of operations required for finding eigenpairs by the proposed method, by the RobinsonHarris method, by the subspace-iteration method 8 = = = n p q order of A and B number of eigenpairs sought number of iteration vectors by subspace iteration method, q -r.(k) = J = max(2p, p+8) residual vector of approximation to jth eigenpair after k iterations s = number of close and/or multiple eigenpairs sought = number of iterations needed to find eigenpairs by proposed method, by Robinson-Harris method, by subspace iteration method - , x, x. J x .. lJ = matrix of eigenvectors (modal matrix), jth eigenvector, element of X -(k) , x~~) X(k), x· lJ J = approximation to X after k iterations, jth column vector of X(k), element of X(k) y(k), y~k), y~ ~) lJ J = matrix of iteration vectors improved from X(k) by simultaneous iteration method, jth column vector of y(k), element of y(k) = rotation matrix, approximation to Z after k iterations = error in A~k) or ~~~) J JJ = increment operator o.. = Kronecker delta e~k) = error in x~k) or y~k) = multiple eigenvalue = diagonal matrix of eigenvalues, jth eigenvalue, lJ J * 11, A. J J 11 = diag(Al' J A2' . . . , AS) 9 A(k) A(k) j = approximation to A, to Aj' after k iterations = shift applied in vector iteration method l.l (k) D(k) l.l' • , 1-1- • = element of 0, of P = mass density w = natural circular frequency, lJ lJ A = w2 10 2. 2.1 DISTINCT ROOTS General In this chapter, a method for finding a simple corresponding eigenvector will be presented. eigenvalue and the The method developed by Robinson and Harris [44J is modified here to save overall computational effort for finding an eigensolution. The Robinson-Harris method is an application of the Newton-Raphson technique for improving the accuracy of an approximate eigenvalue and the corresponding approximate eigenvector. In the proposed method, a modified form of the Newton-Raphson technique is applied instead of the standard one used in the Robinson-Harris method. In Section 2.2, the Robinson-Harris method will be discussed first; then the proposed method will be presented. The convergence rate of the proposed method and the number of operations per iteration will be given in Section 2.3. The estimation of error in an approximate solution is found in Section 2.4. A technique for the examination of the converged solution to determine whether the eigenvalues and corresponding eigenvectors of interest have been missed and a method for finding a missed solution will be presented in Section 2.5. 2.2 The Iterative Scheme Let us consider the following linear eigenvalue problem Ax. J = A. J Bx.J (j = 1, 2, . . . , n) (2.1) 11 where A and B are assumed to be given symmetric matrices of order nand B is taken to be positive definite. The A. and J X.J are the jth eigenvalue and the corresponding eigenvector. Let us assume that an initial approximate solution of Eq. (2.1), A. (0) J and x.(O), is available. J Denote an approximate eigenvalue and the corre- sponding eigenvector after k iterations by A/ k ) and x/ k ) (k = 0, l, . . . ). Then, we have - (k) - (k) r. J Ax. J (2.2) rj(k) is a residual vector. where The object is to remove the residual vector in Eq. (2.2). Raphson technique is applied for this purpose. The Newton- Let the (k + l)th approxi- mation be defined by A • (k+ 1) = J - (k+l) - (k) - (k) = x. + !:J.x. x. J (2.3) J J where ~A.(k) and ~x.(k) are small unknown incremental changes of A.(k) and J J J 1 1 x (k) Sub s tit uti ng A. (k + ) and X. (k + ) 0 f Eq . (2 . 3) for A. and X. i n j . J J J J Eq. (2.1) and discarding a nonlinear term ~A.(kLB~x.(k) as very small J J compared with the other, linear, terms, We get - (k) r. J where r.J (k) is the residual vector defined in Eq. (2.2). (2.4) 12 Note that in Eq. (2.4), there are n+l scalar unknowns (IlA.(k) and n J compone~ts of llX.(k)), but only n equations. Hence, it is required for the J solution of Eq. (2.4) that either the number of unknowns be reduced or one equation added. . Derwidue [16J and Ra11 [41J reduced the number of unknowns by setting the nth component of the vector lli.(k) or i.(k+l) at a preassigned J value - zero or one. ,J In these methods, it may happen that an unfortunate choice of one component results in failure of the procedure. \ Instead of reducing the number of unknowns, Robinson and Harris [44J added an extra equation (side condition) to the system of Eq. (2.4), to arrive at a set of n+l equations in n+l unknowns. T x.(k) . BllX. ( k) J This side condition is =0 (2. 5 ) J Equation (2.5) means that the incremental value llX.(k) is orthogonal to J the current approximate eigenvector x.(k) with respect to the matrix B. J The side condition prevents unlimited change in the xj(k). The resulting set of simultaneous linear equations may be written in matrix form as· - (k) J - Bx. -r. (k) llx- . (k) J J (2.6) = _ (k) T o B x. J llA.(k) J where the residual vector r.(k) is given in Eq. (2.2). J 0 The coefficient matrix for the incremental values is of order n+l and symmetric. it is nonsingular if A. is not multiple [44J. J Moreover, Equation (2.6) may be solved for llA.(k) and llX.(k) by Gauss elimination, or by any other suitable J J 13 technique. Note that the submatrix in the coefficient matrix (A - A.(k)B) J is almost singular when A.(k) is close to A.. J However, this does not cause J any difficulty in solving Eq. (2.6), since in the elimination process the last pivot element, in general, becomes very small. only Thus, the inter- change of columns and rows does not increase significantly the column height of the factorized matrix. computed from Eq. (2.3). The improved values, A. (k+1) and J X. (k+1), J are The procedure employing Eqs. (2.3) and (2.6) is repeated until the errors in the A.(k) and x.(k) are within allowable tolerJ J ances. The method of estimating these errors will be discussed in Sec- tion 2.4. The convergence of the above process for an eigenvalue and the corresponding eigenvector has been shown to be better than second order; the order has been found to be 2.41 [44J. However, the algorithm using Eq. (2.6) requires a new triangularization in each iteration, since the values of the elements of the coefficient matrix are changed in each iteration as a result of changing from Aj(k) to Aj(k+l). The number of operations (multiplications and divisions) required in such a triangularization is very large. To avoid the complete elimination procedure in each iteration, the following equations instead of Eq. (2.6) are used in the proposed method. A - A. (O)B J - (k) -Bx. J - (k) 6X. J - -r. (k) J (2.7) l- x.(k) TB J 0 6A.(k) J 0 14 where the residual vector r.(k) is defined in Eq. (2.2). J Equation (2.7) was obtained by introducing Eq. (2.3) into Eq. (2.1) and discarding a small linear term (A.(k+l) - 1...(0)) BllX.(k). J J J Note that Eq. (2.7) differs from Eq. (2.6) in such a way that the coefficient matrix in Eq. (2.6) has the submatrix (A - Aj(k)B), while the coefficient matrix in Eq. (2.7) has (A - Aj(O)B). The coefficient matrix in Eq. (2.7) is also symmetric, and non sin g~ 1a r i f AJ '1 S no t mu 1t'1 p1e. matrix will be prove~ The nonsingularity of the coefficient in passing, in Appendix A. From the form of the coefficient matrix, it can be seen that once the matrix is decomposed into the form LDL T, where L is lower triangular and 0 _is diagonal, only a small number of additional operations is required for the solution of Eq. (2.7) in the succeeding iterations, since only the vector Bx.J (k) in the matrix is changed in each iteration. The proposed method therefore considerably reduces the number of operations required in each iteration. On the other hand, the method lowers the convergence rate because of the neglect of the small linear term (A.(k+l) - 1...(0)) (BllX.(k)), J J which in turn increases the number of iterations for a solution. the overall co~putational effort for a solution does decrease. J However, It will be seen in Chapter 5 that the proposed method is actually more efficient than the Robinson-Harris method. 2.3 Convergence Rate and Operation Count The efficiency of a numerical method such as the one proposed here can be estimated given the convergence rate and the number of operations per iteration required in the process. The convergence analysis, which is given 15 in Appendix B, will be summarized as follows. - (k) Xj Let an approximate eigenvector be expanded in terms of the true eigenvectors xi' i.e., n x.(k) = r J ._-' c .. (k) 1J X. (2.8) 1 i=l ( k) 1. sac 0 e ff'1 c 1. en t were c.. h 1J 0f th eve cor t x,.. If Y ' (k) J 1'S the error ,. n A.(k) and e .(k) the error in x.(k), they may be defined as J J J A. - A.(k) y.(k) J = J (2.9) A,. J J 1/2 n ,-. ( "- .., (2. 10) c.1 J. (k)2 ) i=l v/here e.(k) is a measure of the angle between the vectors c.(k) and _J(k)T vJhere c. = (c1.(k), c .(k), ... , C .(k) and c~ J = c., and J (0, .. ,0, c .. (k),O, ... ,O). 2J nJ J JJ The geometric interpretation of 8 j (k) is illustrated in Fig. 1. With the above definitions, the errors in A.(k+1) and x.(k+1) may be J J J J written as (see Appendix B) (2.11) 16 e . ( k+ 1) = he. ( k) J J (2.12) where h = max m~j A - m < 1..(0) (m = 1, 2, . . . n) (2.13) J Equatiops (2.11) and (2.12) show that the convergence character of both eigenvalues and eigenvectors is linear. However, the eigenvalues converge much more rapidly than the eigenvectors. Note also that the closer A. is to another eigenvalue, the larger J is, yielding slow convergence. a Hence, the method is not suitable for finding close eigenvalues and the corresponding eigenvectors. Another important consideration which should be taken into account in estimating the efficiency of numerical methods is the number of operations per iteration. One operation is defined as one multiplication or division, which almost always is followed by an addition or a subtraction. expression of this number, let rna For the and mb be the half band-widths of the matrices A and B, and let n be the order of A and B. Let Tp be the number of iterations needed to find p eigenpairs by the proposed method and Tr by the Robinson-Harris method. Then, the number of operations for p eigenpairs, Np , required by the proposed method is -fi p = , 2 ~n ,2 _ ~ \ m +.:Sm a ~ a + Lm +_ b ?)\ '- + T n (5m + 2m + 6) b P a (2.14) 17 and by the Robinson-Harris method, N , is r Nr = 2 IT 2 r n (ma + 13ma + 6mb + 12) (2.15) It can be seen that the number of operations per iteration required by the proposed method is much smaller than for the Robinson-Harris method. The development of the above expressions is given in Table 1. 2.4 Errors in Approximate Eigensolutions An important feature of an iterative method such as the proposed method is some means of estimating the error in a computed solution. This permits one to terminate the iteration process at the point where a sufficiently accurate result has been obtained. It is important to have estimates in terms of numbers available in the calculation, since it is impossible to compare with the exact values. The error in A.(k), y.(k), can be estimated as follows: J J from Eqs. (2.9) and (2.11) A. = J A.(k+l) J (2. 1 6) Substituting Eq. (2.16) for Aj in Eq. (2.9) gives y.(k) = J A • (k) 1 - ---:J"'---_ Aj = (2. 17 ) 18 Since 0 < h « 1 and 0 < - y.Ck) « 1, from Eq. (2.17) J A • (k) 1 - J A ( k+1) (2.18) The error in x.(k), e.(k), can be approximated by [e.(k) - e.(k+1)] J J since e.(k+1) «s.(k). J J J J Furthermore, from Fig. 1, n L(~c .. (k¥ 1/2 lJ i=l i!:rj- - -_ __ - s. (k+ 1) : : : _ J ~x. (k) = J _ (k)T x·J T B~x. (k) 1/2 J (2.19) - (k) Bx. J Therefore, s . (k) J T ~x. (k) B~x. (k) J J T - (k) - (k) x· Bx. - J L J 1/2 (2.20) 19 The number of operations for the estimation of e.(k) is only about J n (2m + 3), which is small compared with the number of operations per iterb ation (see Section 2.3). 2.5 Treatment of Missed Eigensolutions Some of the eigenvalues and corresponding eigenvectors of interest may be missed when the initial approximations are not suitable. In order to check whether this occurs, the Sturm-sequence property [9,31,39,48,51J may be applied. The Sturm-sequence property is expressed as follows: if for an approximate eigenvalue A.(O), (A - A.(O)B) is decomposed into LOL T, J J where L is a lower triangular matrix and 0 a diagonal one, then the number of negative elements in 0 equals the number of eigenvalues smaller than Aj(O). A computed eigenvalue can be checked using the above property with negligible extra computation, since the decomposition of t~e matrix (A - A.(O)B) has already been carried out during the procedure for the J solution of Eq. (2.7). If some of the eigenvalues of interest are detected to be missing, finding them consists of three steps: finding approximations to the missed eigenvalues~ finding approximate eigenvectors corresponding to the missed eigenvalues, and improving the approximate eigensolutions. The approximate eigenvalues can be found by the repeated applications of the Sturm-sequence calculation mentioned above and the method of bisection [9,31,38,51J, or by the polynomial iteration method [7,8,9,38,51J, in which the zeros of the characteristics polynomial p(A) using variants of Newton's method. = det(A - AB) are found 20 In the second step, the approximation to the eigenvectors corresponding to the missed eigenvalues is found. Frequently, finding the eigenvectors corresponding to the missed eigenvalues is much more difficult than finding the missed eigenvalues. However, subspace iterations with a shift [6,32J, which will be discussed in Chapter 4, or dynamic condensation [34,42,50J may be used for this purpose. Finally, the approximate eigenvalues and corresponding approximate eigenvectors can be improved by the method of Section 2.2 if the eigenvalues are not multiple or close, or if they are, by the method of Chapter 3. 21 3. 3.1 CLOSE OR MULTIPLE ROOTS General As mentioned earlier, the method presented in Chapter 2 fails or exhibits slow convergence if it is applied to the solution for multiple or close eigenvalues and for their corresponding eigenvectors. The failure or slow convergence of the method is caused by impending singularity of the coefficient matrix for the unknown incremental values as the successive approximations approach the true eigenvalue and eigenvector. The method presented in this chapter overcomes this shortcoming. To accomplish this, all eigenvectors corresponding to multiple or close eigenvalues are found together. As in the method of Chapter 2, this method yields the eigenvalues and corresponding eigenvectors at the same time. The essence of the method consists first in finding the subspace spanned by the eigenvectors corresponding to multiple or close eigenvalues. The subspace is found using the Newton-Raphson technique in a way suggested by the Robinson-Harris method [44J. If the eigenvalues of interest are multiple, any set of independent vectors spanning subspace are the true eigenvectors, but if the eigenvalues are merely close together, the vectors must be rotated in the subspace to find the true eigenvectors. The eigenvalues are obtained as a by-product of the process of finding the subspace and any subsequent rotation. In this method, any number of close eigenvalues or an eigenvalue of any multiplicity can be found together with the corresponding eigenvectors. 22 The theoretical background of the method is presented in S~ction 3.2 The iterative scheme for finding the subspace of the eigenvectors corresponding to multiple or close eigenvalues is given in Section 3.3. The additional treatment required for close eigenvalues and corresponding eigenvectors is the subject of Section 3.4. The convergence rate and the number of operations per iteration are given in Section 3.5. 3.2 Theoretical Background Let us consider the system treated in Chapter 2, i.e., Ax.1 = A.1 BX.1 (i = 1, 2, . . , , (3.1) n) where A and B are symmetric matrices of order n, and B is positive definite. The -x. are ei genvectors, and the Ai eigenvalues in the order A1 <- A2 < - , ... , -<A n . 1 Let a set S consist of s integers p. (j = 1, 2, J . , s), that is, S = [Pl ' P2,···,P s J where 1 ~ p. < n. The s-dimensional subspace spanned by J the eigenvectors x. (jES) where none of the corresponding eigenvalues J A. (jES) are close or equal to eigenvalues A. (itS) is denoted by R. J 1 Let us take s vectors y. (jES) which are orthonormal with respect to B and are J in the neighborhood of the subspace R. This means that if the vector y. is J expanded in a series of true eigenvectors x. (i = 1, 2, ... ,n) 1 n y. J = L i=1 c .. lJ - v A. 1 (jES) (3.2) 23 then, the following relations must be met: \"' 2 L.. (jsS) c .. . S lJ (3.3) ls Hence, a vector Yj(jsS) needs not be close to one of the xj(jsS). With the above definitions~ the subspace R of the eigenvectors x.(jsS) J is characterized by the following constrained stationary-value problem: find the stationary values of n ~ w= L_, - T - y. J jsS Ay. (3.4) J subject to -T By. y. 1 J = o.. (i, jsS) lJ (3.5) where Qij is the Kronecker delta, i.e., Qij = 1 for i = j, and Qij = 0 for i f j. The function w could be regarded as a sum of Rayleigh quotients of the vectors y., J since by Eq. (3.5) the denominators of the Rayleigh quotients are equal to unity. The important result that the stationary property characterizes the subspace R is proved as Theorem 1 of Appendix C. The stationary-value problem may be treated by the method of Lagrange multipliers. 11·· = 11·' lJ Jl Introducing the undetermined multipliers ~ij (i ,jsS) and letting (see Eq. (3.5)), we have the Lagrangian L = / •. _r isS -T y. AY· 1 1 \' L--- '\"' L isS jsS 11· . lJ (-y.T By. - Q.. ) 1 J lJ (3.6) 24 The problem of Eqs. (3.4) and (3.5) is equivalent to that of solying the unconstrained stationary-value problem for the Lagrangian L. The problem is solved setting the first partial derivatives of L with respect to the unknowns y.J and ~ .. lJ aL - equal to zero, i.e., L ~ ..J AY·J = =0 'Oy. 1 By.1 (j sS) (3.7) (i,jsS) (3.8) isS J -T y. By. = ~=O a~ .. lJ 1 J o.. lJ Introducing the following notation o = (d ,d , ... ,d ) PI P2 Ps (3.9) we can write Eq. (3.7) in matrix form as AY·J = Bya.J (3.10) AY = BYD (3.11) or collectively 25 In the same way, Eq. (3.8) can be written as (3.12) where Is is the unit matrix of order s. Hence, the subspace R of the desired eigenvectors can be found by solving Eqs. (3.11) and (3.12). Note that Eqs. (3.11) and (3.12) are nonlinear in 0 and Y and that there are s (s + 1)/2 scalar unknown elements in 0, since 0 is symmetric, and s (s + 1)/2 independent equations in Eq. (3.12). In the next section, the solution of Eqs. (3.11) and (3.12) in the special case that (jsS) are all multiple or close eigenvalues will be discussed. 3.3 The Iterative Scheme In this section, the application of the Newton-Raphson technique to the solution of Eqs. (3.11) and (3.12) for multiple or close eigenvalues and their corresponding eigenvectors will be presented. To simplify the notation in this discussion, we take the set S = [1,2, ... ,sJ, that is, the s lowest eigenvalues are close together, or the multiplicity of the lowest eigenvalue is s. It should be emphasized that this is not restrictive, and the procedure is perfectly applicable to multiple or close eigenvalues in any range. Assume that the initial values for 0 and Y, 0(0) and y(O) are available (the solution for the initial values will be discussed in Chapter 4). Furthermore, we assume that the initial vectors in y(O) are in the neighborhood of the subspace of the eigenvectors X = [xl ,x2 , ... ,x s J and that they 26 T have been orthonormalized with respect to the matrix B, i.e., y(O)BY(O) = I . s With the above assumptions, we now apply the Newton-Raphson technique to the solution of Eqs. (3.11) and (3.12). For the general kth iteration step, let (3.13) where ~a.(k) and ~y.(k) are unknown incremental values for d.(k) and y.(k). J J J J Introducing Eq. (3.13) into Eqs. (3.10) and (3.12) and neglecting the nonlinear terms, we obtain the linear simultaneous equations for ~a.(k) and J (3.14) (3.15) By Theorem 3 of Appendix C, if the A. (j = 1,2, ... ,s) are multiple or J close eigenvalues, the off-diagonal elements of 0 are zero or very small compared with its diagonal ones, thus the last term of Eq. (3.14) may be . t e d by approxlma ~jj B ~Yj - (k) ,Yle . ld·lng Let us take (3.17) ,...,.., L/ Then, Eq. (3.15) becomes (3.18) which is the condition that the incremental vectors be orthogonal to the current vectors with respect to B. y.1 (k) altered so that the latest If the computational scheme is slightly is used at all times, the orthogonality condition is satisfied automatically provided that the initial vectors y.(O) 1 are orthogonal. What this means is that we use y.1 (k) (i 1,2, ... ,j - 1) = for the computation of y.(k+1). J The final equations to solve for ~a.(k) and ~y.(k) are Eqs. (3.16) J J and (3.18) along with the orthonormality condition, Eq. (3.17). These equations can be written in matrix form as A- 11 •• _ By(k) (k)B -- - (k) ~y. J JJ --- .- - (k) - rj (3.19) - T _ y(k) B 0 . ~a (k) J 0 where - (k) r. J (3.20) The coefficient matrix for the unknowns, dj(k) and Yj(k), is symmetric. Furthermore, it is nonsingular, as is shown in Appendix A. Thus, Eq. (3.19) 28 can be solved for ~a.(k) and 6y.(k), yielding improved values, a.(k+1) and J J J - (k+ 1) yj from Eq. (3. 13) . The algorithm using Eq. (3.19) requires a new triangularization in each iteration, since the coefficient matrix is changed in each iteration. It therefore seems useful, as in Chapter 2, to substitute (A - ~ .. (0)8) for JJ (A - ~jj(k)B) in Eq. (3.19) in order to save computational effort in the solution. That is, the basic equations for the increments are taken as = -.- - _. - o (3.21) o where the residual vector r.(k) is defined . as in Eq. (3.20). The coefficient J matrix in Eq. (3.21) is also symmetric and nonsingular (Appendix A). The equation (3.21) was obtained discarding a small linear term (~ .. (k) ~.. JJ (0) ~ )8 y. (k) J JJ of Eq. (3.19). . The procedure using Eq. (3.21) requires only partial triangularizations in each iteration, since only the vectors in y(k) are changed, reducing the number of operations per iteration. The pro- cedure depends, for its convenience, on the decoupling of the ~Yj(k) for the s vectors y.(k) (i=1,2, ... ,s). J the small linear terms i=l ifj The decoupling was possible only because 29 (see Eq. (3.14)) could be dropped for A. (j=l, 2, . . . , s) all J (3.21) for A. (j=l, 2, . , s) J close together. vihich are not Experience with Eq. close together indicates that satisfactory results cannot be obtained. Note that if s = 1, Eqs. (3.19) and (3.21) are equivalent to the equations used for distinct eigenvalues and corresponding eigenvectors: Eq. (3.19) becomes Eq. (2.6), the equations used in the Robinson-Harris method, and Eq. (3.21) becomes Eq. (2.7), used in the proposed method. With sufficient large k, the incremental values ~a.(k) and ~y. (k) J will vanish. J Then, from Eq. (3.21) 1i m k-)-OO r. (k) J = 1 i m (AY. ( k ) k~ (3.22) J Letting lim d.(k) dj = k-tco - _limy.(k) J J (3.23) AY = BYD (3.24) Yj - k-tco we write Eqs. (3.22) and (3.17) as where Y = (Y1'Y2'··· 'Ys), and 0 = (d 1 ,d 2 ,···,d s )· By Theorem 3 of Appendix C, if the eigenvalues A. (j=1,2, ... ,s) are multiple, the values of the offJ diagonal elements of 0 are all zero, and its diagonal elements have an equal 30 value which is the desired multiple eigenvalue. Yare the corresponding eigenvectors. Moreover, the vectors in However, if the eigenvalues are close but not equal, additional operations are required to find the desired eigenvalues and eigenvectors. These additional operations are the subject of the next section. 3.4 Treatment of Close Roots Once the converged solution 0 and Y has been found by the algorithm described in the previous section, but the values of the off-diagonal elements of 0 are not zero, the vectors in Yare rotated in the subspace of Y to find the true eigenvectors. a small eigenvalue problem. A rotation matrix is found by solving Furthermore, the eigenvalues of the small eigenvalue problem are the desired eigenvalues. small eigenvalue problem is as follows. in X = [x l ,x2,··.,xs J and The derivation of the The system with the s eigenvectors corresponding eigenvalues in A = diag 1 2 ,··· (A ,A ,AS) may be written as AX = BXA where A and B are symmetric matrices of order n. x= YZ where Z is the unknown rotation matrix of order s. (3.26) Now, let (3.27) Introducing Eq. (3.27) into Eq. (3.26), we get AYZ = BYZA (3.28) 31 Postmultiplying Eq. (3.24) by the matrix Z yields AYZ = BYDZ (3.29) Premultiplying Eqs. (3.28) and (3.29) and using yTBY = Is of Eq. (3.25), we obtain the special eigenvalue problem of order s DZ = ZA (3.30) where D is the converged solution found by the algorithm of the previous section. The matrix D is symmetric (see Eq. (3.24)) and of order s, the number of close eigenvalues, which is usually small. The absolute values of the off-diagonal elements of D are small compared with those of its diagonal elements (see Appendix C). The eigenvalue problem, Eq. (3.30) can be easily solved by any suitable technique such as Jacobi's method [31,51J, yielding the desired eigenvalues in A (AI' A2' which in turn . , As) and the matrix Z, gives the eigenvectors X by Eq. (3.27). The number of oper- ations required for the solution of Eq. (3.30) is very small compared with that of Eq. (3.21), since s is small. 3.5 Convergence Rate and Operation Count In this section, the convergence rates of a multiple eigenvalue and the corresponding eigenvectors found in Appendix B will be summarized. For convenience, we assume that the lowest eigenvalues are multiple, i.e., A* = Al = A2 = ... = As' Le t th e approxlma . t - ( k) (.J -- 1 , 2 , ... ,s ) e· elgenvec t ors Yj be expanded in terms of the ei genvectors xi (i = 1, 2, . . . , n), i. e. , 32 n -.(k) YJ = \"' c .. (k) X. i .. lJ J j = 1, 2, ... , s 1 (3.31) ;=1 where cij(k) is a scalar representing the components of the eigenvector xi on y- .(k). If y.(k) denotes the error in (k) and e.(k) the error in y- .(k), J J lljj J J then they may be defi ned by y. J rn t , \"' L ' !Li=s+l (c.. ( k) ( k 1J = 17[ lJ ') (3.32) A* 1/2 s \" (c.. 1J i=l 1-I ( k 2l1 /2 ') (3.33) J As shown in Appendix B, the error in ~jj(k+l) and ;j(k+l) m~ be written as y. J (k+1 ) (3.34) = e.(k+1) = he.(k) J (3.35) J where h = m~x i A*- ll"(O) JJ I ()) A. 1 u.!.!\V . JJ I I I « 1 i=s+1,s+2, ... , n; j=l, 2, ... , S (3.36) 33 t can be seen from Eqs. (3.34) and (3.35) that the eigenvalues and the corresponding eigenvectors converge linearly. However, the eigenvalues converge much more rapidly than the eigenvectors. The number of operations Np required for finding multiple or close eigenvalues and the corresponding eigenvectors is calculated in Table 1. Th is number is Np = 12 pn (m a2 + 3m a + 2m b + 2) + Tp n [( s + 4) rna + 2m b + 12 (s 2 + 7s + 4) ] (3.37) where s is the multiplicity of an eigenvalue or the number of close eigenvalues, and Tp is the total number of iterations required for a solution. It can be seen that if s = 1, the number of operations is equal to the number of operations required for finding a simple eigenvalue and the corresponding eigenvector (see Eq. (2.14)). 34 4. 4.1 APPROXIMATE STARTING EIGENSOLUTION General The iterative methods described in the previous chapters begin with an approximate starting eigensolution. In this chapter, a procedure to find the starting solution is presented. The approximate starting solution of an eigenvalue problem is often available either as the final answer in some approximate methods or as an intermediate result in other iterative methods. Numerous methods for approximate solutions have been developed. These include static or dynamic condensation [2,3,25,28,34,42J, Rayleigh-Ritz analysis [48,51J, component mode analysis [9,51J, and related methods summarized by Uhrig [50J. In all these methods, the approximate solution is found in a single step, and not in an iterative process. Hence, automatic improvement of the solution is not built into the procedure. Moreover, the success of the methods depends, to a great extent, on the engineer's judgment, which is difficult to incorporate into an automatic computer program. Another possible way for finding the approximate solution is to take the intermediate results from other iterative methods such as a method combining the Gran-Schmidt orthogonalization process [51J with simultaneous iteration method or combining Rayleigh-Ritz analysis [6,9,11,29,32,49J with simultaneous iteration method. The latter combined method is sometimes called the "subspace iteration rrethod" [6,9J. The subspace iteration method is used here to find approximate starting solutions because it has a better convergence rate than most others. selecting starting vectors. The method itself turns out to require However, a scheme to find starting vectors for 35 the subspace iteration method has been well established and is fairly routine (see Section 4.2.2). In the next section, the subspace iteration method will be di scussed. 4.2 Subspace Iteration Method 4.2.1 The Iterative Scheme The subspace iteration method is a repeated application of the classical vector iteration method (power method) and Rayleigh-Ritz analysis. Suppose that the p smallest eigenvalues A.1 (i = 1,2, ... ,p) and corresponding eigenvectors x.1 are required and that we have p initial independent vectors - ( 0 ) (.1 -- 1 " 2 ..• ,p ) spannlng . a p- dlmenSlona . . 1 su bspace ln . th e nelg . hbor h00 d xi of the subspace of the desired eigenvectors. If the approximate eigenvectors and corresponding eigenvalues after k iterations are denoted by x.(k) and 1 Ao 1 (k), X(k) = [x1(k), x2(k) , ... ,x p(k)J, . (k) ' A2 (k) , ... ,A (k) ), the subspace iteration method for and 0(k) = dlag (AI p the kth iteration may be described as follows: (i) Find the improved eigenvectors y(k) = [Y1 (k), Y2(k) , ... ,Y (k)] p by the simultaneous inverse iteration method; (4.1) (ii) Compute the projections of the operators A and B onto the subspace spanned by the p vectors in y(k); A(k) = y(k)T Ay(k) g(k) = y(k)T By(k) (4.2) 36 where A(k) and S(k) are pxp symmetric matrices. (iii) Solve the eigenvalue problem of reduced order p for the = diag (A1(k), A2(k) , ... ,Ap(k) Z(k) = [Zl(k), Z2(k) , ... ,Zp(k)]; eigenvalues in D(k) eigenvectors in and the (4.3) (iv) Find an improved approximation to the eigenvectors; (4.4) Then, lim O(k) = k~ lim X(k) = (4.5) k~ Note that Eqs. (4.2) through (4.4) represent a Rayleigh-Ritz analysis with the vectors in y(k) as the Ritz basis vectors, which results in x(k), the best approximation to the true eigenvectors in the subspace of y(k). More rapid convergence can be obtained by taking more iteration vectors than the number of eigensolutions sought. However, the more starting vectors are taken, the more computational effort is required per iteration. As an optimal number of iteration vectors, q, q = min (2p, p + 8) has been suggested [6,9J. Metz Referenoe Room Civil Engineering Departman~ BI06 C.E. Building 37 T~4~a~~4+~ UU~YV.U.VJ n~ v. Tl14"~~_ ~~~_~w~~ Urbana,,111ino1a 6180ll To find eigenvalues within a given range a < ~ < b and the correspond- ing eigenvectors, we may use, instead of Eq. (4.1), the inverse iteration with a shift [32J: (A - ~B) y(k) = BX(k-1) where ~ is a shift and can be taken as (a + b)/2. (4.6) It is clear from Eq. (4.6) that the eigenvectors corresponding to the eigenvalues in the vicinity of a shift ~ will converge rapidly. However, the convergence of other eigen- vectors may be slower than when the shift is not applied, since as a result of the application of the shift, the absolute values of some shifted eigenvalues may become closer. 4.2.2 Starting Vectors The number of iterations required for convergence depends on how close the subspace spanned by the starting vectors is to the exact subspace. If approximations to the required eigenvectors are already available, e.g., from a previous solution to a similar problem, these may be used as a set of starting vectors. If not, we may use one of the schemes for generating starting vectors which have been proposed as effective [6,11,32,47J. The scheme for establishing the starting vectors proposed by Bathe and Wilson [6,9J is used here because of its simplicity and effectiveness. scheme may be described as follows. The The first column of BX(O) in Eq. (4.1) is formed simply from the diagonal elements of Bo That is, if BX(O) is denoted by C, = b .. 11 (i = 1,2, ... ,n) (4.7) 38 This assures that all mass degrees-of-freedom are excited in miss a mode [6,9J. or~er not to The next (q-1) columns in C may each have all zeros except for a certain coordinate where a one is placed. are found in the following way. These coordinates First, compute the ratios a.·/b .. 11 11 (i = 1,2, ... ,n) and take the (q-1) s"s (j = 1,2, ... ,q-1) such that the J absolute values of the ratios aiilb ii for i (i = sl' s2,···,Sq_1) are smallest over all i. Then, c·1 , j-1 = 1 for =0 for = s. J i ~ s. J (i = 1,2, ... ,n) (j = 1,2, ... ,q-1) (4.8) If the absolute values of the ra ti os are close or equal, then it was recommended [6,9J that the s .IS (j J = 1,2, ... ,q-l) be chosen so that they are well spaced. 4.2.3 Convergence Rate, Operation Count, and Estimation of Errors With an adequate choice of the starting vectors, the subspace iteration method gives good approximations to the exact eigenvalues and eigenvectors even after only a few iterations. However, the subsequent convergence is only linear with the rates of convergence equal to Ai/Aq+1 (i = 1,2, ... ,p) for the ith eigenvector and (A i /A q+1)2 for the corresponding eigenvalue. These ratios indicate that for the higher eigenvalue convergence is slower. Hence, the convergence of the pth mode controls the termination of the iteration process. 39 One of the most important indicators of the effectiveness of numerical methods is the total number of operations required for finding a solution, which depends on both the rate of convergence and the number of operations per i terati on. This number for the subspace iteration method,N , (see s Table 1) may be expressed by where rna and mb are the half band-widths of A and B, and Ts is the total number of iterations required for the solution. The total number of iterations T , depends on the rate of convergence s and tolerances of the errors in approximate eigenvalues and eigenvectors. Bathe and Wilson [6,9J suggested use of the following formula for the estimation of errors in the ith eigenpair at the kth iteration: - (k) r. 1 (4.10) - ( k) Ax.1 where r. (k) = 1 (A - A. (k)B) , x. (k). 1 The error estimated by Eq. (4.10) is a function of both the approximate eigenvalues and eigenvectors. However, it may be more reasonable to estimate the errors in approximate eigenvalues and eigenvectors using separate formulas as follows: let Yi (k) and i (k) be the errors in the ith approximate eigenThen y. (k) may be estimated by value and eigenvector. 8 1 A. (k+ 1) _ A. (k) y. ( k ) ~ _,_ _.,....---=---,__ 1 A.(k+1) 1 (i = 1,2, ... ,p) (4.11) 40 For the estimate of 8,. (k) , we find the incremental vectors ~x. (k) from the , , re 1a ti ons X. (k+l) = a .. (k)x. (k) , X., (k) ", T B ~x. (k) , + ~x. (k) , =0 ( 4. 12) Then, (4.13) If some of the approximate eigenvalues Ai (i = PI' P2'··· ,ps) are equal or very close, we may then compute ~Xi (k) from the relations T Xj (k) Bt;x; (k) 8. , (k) =0 ::: ; (j (- (k)T X• " = Pl'P2'··· ,ps) B (Ps La.2lJ.(k)- J (k)TB_J' (k))1/2 - (k))1/2 X· X• j=p X (4.14) 1 For the purpose of comparison of the proposed methods of Chapters 2 and 3 with the subspace iteration method, the errors were computed using Eqs. ( 4. 11) to (4. 14 ) . 41 4.3 Starting Solution for the Proposed Method The intermediate results from the subspace iteration are used as the starting solutions for the proposed method. During the subspace iterations, the errors in approximate eigenvalues and corresponding eigenvectors can be estimated by the scheme described in Section 4.2.3. Furthermore, these errors can be used for estimating the number of iterations or the number of operations required for the solution by both the subspace iteration method and the proposed method. Hence, it is possible to estimate the optimal number of iterations to be carried out by the subspace iteration method. This optimal number of iterations is usually one or two. Let A~ and x~ (i = 1,2, ... ,p) be the intermediate solutions from the subspace iteration method after the optimal number of iterations. Then, if x.1 can be taken as the starting solutions the A.*1 are well separated, A.*1 and -* for the method of Chapter 2, A'(O) and i.(O). 1 1 However, if some of them, "\ *. and x-*.1 are t a ken as e.g., Ai* (.1 = PI ,P2'··· 'Ps ) are equa 1 or very close, 1\1 the starting solution for the method of Chapter 3 as y. (O) = x-* (0) = A~ 1 ~ii ~. i 1 . (0) = 0 1J for i f j 42 It should be noted that from Eqs. (4.3) and (4.4), the ite~ation vectors in the subspace iteration method are always orthogonalized with respect to B. Therefore, orthogonal;zation is not required for the first iteration of the proposed method. 43 5. 5.1 NUMERICAL RESULTS AND COMPARISONS General The relative efficiency of the methods developed in this study is illustrated in this chapter by the numerical results of the free vibration analyses of the following example problems: (a) Ten-Story, Ten-Bay Plane Frame (b) Two-Hinged Circular Arch (c) Simply Supported Plate. The problems were formulated using a stiffness method for the plane frame problem, a finite difference method for the arch problem, and a finite element method for the plate problem. No attempt has been made to present the solutions of eigenvalue problems of very large order, although the proposed method is developed for them. However, some trends can be inferred from the example problems presented here. The first two problems, with distinct eigneva1ues, were solved by the method discussed in Chapter 2 and the third one, with multiple or close eigenvalues, by the method of Chapter 3. The above problems were also solved using the Robinson-Harris method [44J and the subspace iteration method discussed in Chapter 4. The results are summarized in Tables 2 through 5. The numerical results given here are shown to be consistent with the convergence estimates of Appendix B. For each method, the total number of operations required for finding the desired eigenvalues and eigenvectors to the same accuracy was found. These are presented and compared in Table 6. Although a tolerance of]O -4 on the eigenvalues and eigenvectors should be sufficient for normal requirements, it 44 was taken as 10- 6 for the purpose of comparisons of the convergence characteristics of the methods. The numerical computations of the above problems were performed on the CDC CYBER 175 system of the Digital Computer Laboratory of' the University of Illinois, Urbana, Illinois. 5.2 Plane Frame The ten-story, ten-bay plane shown in Fig. 2 was taken as an example problem in order to test the method of Chapter 2 for problems with distinct eigenvalues. The problem was formulated by a stiffness method in which the axial deformations of the members are considered, but the shear deformations neglected [40J. The frame with three displacements per joint has a total of 330 degrees of freedom. The mass matrix is the consjstent mass matrix [4,5J with a maximum half-bandwidth of 35, equal to that of the stiffness matrix. The four smallest eigenvalues and their corresponding eigenvectors were computed by the proposed method, by the Robinson-Harris method, and by the subspace iteration method. The results are given in Table 2. For the subspace iteration method, ten starting vectors were formed by the technique suggested by Bathe and Wilson (see Section 4.2.2). The starting approximate eigenvalues and eigenvectors for the proposed method and for the RobinsonHarris method were established by performing two cycles of subspace iteration. Table 2 shows that even the eigenvalues calculated by two subspace iterations are already to only one or two figures. However, the eigenvectors are accurate In addition, the convergence of eigenvectors by the subspace iteration method is so slow, as discussed in Section 4.2.3, that 12 iterations were required for the convergence of both eigenvalues arid eigen- 45 vectors to the indicated tolerance. The proposed method and the Robinson~ Harris method required only two iterations for the convergence of eigenpairs except for that of the fourth mode, which required four iterations by the proposed method and three iterations by the Robinson-Harris method. The total number of operations to solve for all the desired eigenpairs by the proposed method is 3.50xl0 6 ; by the Robinson-Harris method, 4.57xl0 6 , and by the subspace iteration method, 9.27xl0 6 . Therefore, the Robinson-Harris method required 1.31 times as many operations as the proposed method did, and the subspace iteration method required 2.78 times as many operations, as shown in Table 6. 5.3 Arch A uniform 90 degree circular arch simply supported at both ends was analyzed for in-plane vibration behavior. thickness h, and the ratio a/h = 20. The arch has the radius a and the Melin and Robinson [36J investigated the free vibration behavior of such an arch as a part of a study of vibrations of a simply supported cylindrical shell using a finite difference method. arch was divided into 12 uniform segments giving 22 degrees of freedom. The The maximum half-bandwidth of the stiffness matrix is four and the mass matrix is a unit diagonal matrix. The problem was analyzed for the three smallest eigenvalues and their eigenvectors by the proposed method, by the Robinson-Harris method, and by the subspace iteration method. The results are summarized in Table 3. Five radial displacements were taken as master displacements for the iteration vectors of the subspace iteration method. Starting approximate eigenpairs for the proposed method and the Robinson-Harris method were established by carrying out just one cycle of the subspace iteration. 46 The comparison of the total number of operations for each method is given in Table 6. The proposed method needed 8.87xl0 3 operations, the Robinson-Harris method 9.77x.0 3 operations, and the subspace iteration method 1.76xl0 4 operations. Hence, the ratio of the total number .of operations by the Robinson-Harris method to that by the proposed method is 1.10, and this ratio for the subspace iteration method is 1.98. 5.4 Plate Bending A plate simply supported on all edges was analyzed in order to test the method presented in Chapter 3, for the solution of eigenvalue problems with multiple or close eigenvalues. theckness h. The plate has the lengths a and b, and the Two special cases were considered; an aspect ratio bla of 1.00 and bla equal to 1.01. The first case gives multiple roots, while the second one gives close roots. The problem was formulated by a finite element method, in which the plate was divided into 16 elements. Each unrestarined node has a deflection and two rotational displacements, giving a total of 39 degrees of freedom. The mass matrix is the consistent mass matrix [4,5J with a maximum half-bandwidth of 16, equal to that of the stiffness matrix. The four smallest eigenvalues and corresponding eigenvectors were computed for both cases by the proposed method and by the subspace iteration method. The results are summarized in Tables 4 and 5. the deflection at each node was taken as the master degrees of freedom, giving nine iteration vectors for the subspace iteration method. for the proposed method. Only one cycle of subspace iteration was performed The multiple eigenvalues of the square plate close eigenvalues of the rectangular plate were isolated by the method discussed in Chapter 3. 47 The total number of operations by the proposed method for both cases is 5 1 .27xl0 and by the subspace iteration method, 2.20xl0 5 , as shown in Table 6. Hence, the subspace iteration method needed 1.73 times as many operations as the proposed method did. 5.5 Comparison between the Theoretical Convergence Rates and Numerical Results It was shown in the previous chapters that in the proposed method, the convergence of eigenvalues is much faster than that of eigenvectors. Hence, the convergence of the eigenvectors governs the termination of process, when the tolerances on the eigenvalues and eigenvectors are same. Comparison between the theoretical convergence rates and numerical results was, therefore, carried out only for the eigenvectors. Comparisons between the proposed method and subspace iteration method are given in Tables 7, 8, and 9. The numerical convergence rates were computed by e~k+l)/e~k), where e~k) is the error on the ith approximate eigenvector at the kth iteration. These errors are given in Tables 2 through 5, showing that the numerical convergence rates for the proposed method and the subspace iteration method increase monotonically to approach the theoretical convergence rates as the number of iterations increases. A typical example for this is the convergence rates of the fourth eigenvector of the frame problem, as shown in Table 7. The number of iterations for this mode is large enough to provide a good comparison between the theoretical and numerical convergence rates. Tables 7, 8, and 9 show that in the proposed method, eigenpairs much faster than in the subspace iteration method. Note also that in Table 9, the numerical convergence rates for the proposed method are almost same as 48 those rates for the problem with double roots. Hence~ the expressions for the theoretical convergence rates for multiple eigenvalues also seem applicable to the case of close eigenvalues. 49 6. 6.1 SUMMARY AND CONCLUSIONS Summary of the Proposed Method Two iterative procedures for the solution of linear eigenvalue problems for systems with a finite number of degrees of freedom were discussed in Chapters 2 and 3. Chapter 2 developed a procedure for finding distinct eigenvalues and the corresponding eigenvectors, and Chapter 3 dealt with multiple or close eigenvalues and the corresponding eigenvectors. For distinct eigenvalues and the corresponding eigenvectors, the RobinsonHarris method [44J was modified to save overall computational effort by the use of a II modified" form of the Newton-Raphson technique. The modified method reduces both the number of operations per iteration and the convergence rates. However, the reduction of the number of operations generally compensates for the disadvantage of the decrease of the convergence rate, reducing the total number of operations. The procedure in Chapter 2 for finding a distinct eigenvalue and the corresponding eigenvector fails if the eigenvalue is one of multiple or close eigenvalues, because the matrix involved in the computation become ill-conditioned. This difficulty has been overcome by the new method of Chapter 3. In this mehtod, all eigenvalues close to an eigenvalue or a multiple eigenvalue and the corresponding eigenvectors are found in a group. In other words, a subspace spanned by the approximate eigenvectors is projected by iterations onto the subspace of the exact eigenvectors. the vec~ors If the eigenvalues are multiple, spanning the subspace are exact eigenvectors. However, if the eigenvalues are close, the exact eigenvectors are found by a simple rotation of the vectors in the subspace. The rotation matrix is found from a special 50 eigenvalue problem of small order s, the number of the close eigenvalues. The eigenvalues of the small eigenvalue problem are exact eigenvalues of the original system. The above procedures of the successive approximations require initial approximations to the eigenvalues and eigenvectors. These are available either as the final solution in some approximate methods such as static or dynamic condensation or as an intermediate result in an iterative method as the subspace iteration method described in Chapter 4. 6.2 Conclusions The method presented in this study is very efficient for finding a limited number of soltutions of eigenvalue problems of large order arising from the linear dynamic analysis of structures. The features of the method are summarized as follows. (a) The method has very high convergence rates for eigenvalues and eigenvectors. The method is more economical than the subspace iteration method, the advantage being greater in larger problems. For comparable accuracy, a ten-story ten-bay frame required only 36% of the number of operations need in applying subspace iterations. (b) A transformation to the special eigenvalue problem is not required. Thus, the characteristics of the given matrices such as the sparseness, bandness, and symmetry are preserved, mi~imizing operations. the storage requirements and the number of 51 (c) Any number of multiple or close eigenvalues and their eigenvectors can be found. The existence of the multiple or close eigenvalues can be detected during the iterations by the method of Chapter 2. (d) The eigenvalues in any range of interest and their eigenvectors can be found, if approximations to the solution are known. (e) The solution can be checked to determine if some eigenvalues and corresponding eigenvectors of interest have been missed, without extra operations. 6.3 Recommendations for Further Study Several possible areas of further study to improve the proposed method may be suggested. (a) The convergence rate may be improved by other modifications of the successive approximation method used for the proposed method. (b) Further improvements may be possible for the method of finding an initial approximation to the eigensolution, and for isolating the eigenvalues and their eigenvectors which may be missed by the proposed method. (c) The proposed method may be applied to other practical problems of our interest such as a stability analysis of structures. (d) The proposed method could be easily extended to the continuous eigenvalue problems if there were better ways of direct estimation of their eigensolutions. 52 LIST OF REFERENCES 1. Aitken, A. C., liThe Evaluation of Latent Roots and Vectors of a Matrix," Proceedings of Royal Society, Edinburgh, Vol. 57, 1937, pp. 269-304. 2. Anderson, R. G., Irons, B. M., and Zienkiewicz, O. C., "Vibration and Stabi 1i ty of Pl ates Us i ng Fi ni te El errents, II Interna ti ona 1 Journal of Solids and Structures, Vol. 4, No. 10, 1968, pp. 1031-1055. 3. Appa, K." Smith, G.o C. C. and Hughes, J. T., "Rational Reduction of Large-Scale Eigenvalue Problems," Journal of the American Institute of Aeronautics and Astronautics," Vol. 10, No.7, 1972, pp. 964-965. 4. Archer, J. S., "Consistent Mass Matrix for Distributed Mass Systems," Journal of the Structural Division, Proceedings of the American Society of Civil Engineers, ~ol. 89, No. ST4, August 1963, pp. 161-173. 5. Archer, J. S., "Consistent Matrix Formulation for Structural Analysis Using Finite-Element Techniques," Journal of the American Institute of Aeronautics and Astronautics, Vol. 3, No. 10, 1965, pp. 1910-1918. 6. Bathe, K. J. and Wilson, E. L., "Large Eigenvalue Problems in Dynamic Analysis," Journal of the Engineering Mechanics Division, Proceedings of the American Society of Civil Engineers, Vol. 98, No. EM6, December 1972, pp. 1471-1485. 7. Bat he, K. J. and Wi 1son, E. L., II Ei gens 0 1uti 0 n 0 f La rg e Str uc t ura 1 Systems with Small Bandwidth," Journal of the Engineering Mechanics Division, Proceedings of the American Society of Civil Engineers, Vol. 99, No. EM3, June 1973, pp. 467-479. 8. Bathe, K. J. and Wilson, E. L., IIS ol ution Methods for Eigenvalue Problems in Structural Mechanics,1I International Journal for Numerical Methods in Engineering, Vol. 6, 1973, pp. 213-226. 9. Bathe, K. J. and Wilson, E. L., Numerical Methods in Finite Element Analysis, Prentice-Hall Inc., Englewood Cliffs, New Jersey, 1976. 10. Bradbury, w. w. and Fletcher, R., "New Iterative Methods for Solution of the Eigenproblem," Numerische Mathematik, Vol. 9, 1966, pp. 259-267. 11. Corr, R. B. and Jennings, A., "A Simultaneous Iteration Algorithm for Symmetric Eigenvalue Problems," International Journal for Numerical Methods in Engineering, Vol. 10, 1976, pp. 647-663. 12. Courant, R. and Hilbert, D., Methods of Mathematical Physics, Vol. 1, Interscience Publishers, New York, 1953; 53 13. Crandall, S. H., IIIterative Procedures Related to Relaxation Methods for Eigenvalue Problems," Proceedings of Royal Society, London, A207, 1951~ pp. 416-423. 14. Crandall, S. H., Engineering Analysis, McGraw-Hill Book Co., Inc., New York, 1956. 15. Dong~ S. B. and Wolf, J. A., Jr., "On a Direct-Iterative Eigensolution Technique," International Journal for Numerical Methods in Engineering, Vol. 4, 1972, pp. 155-161. 16. Faddeev, D. K. and Faddeeva, V. N., Computational Methods of Linear Algebra, W. H. Freeman and Co., San Francisco and London, 1963. 17. Felippa, C. A., ilRefined Finite Element Analysis of Linear and Nonlinear Two-Dimensional Structures," Ph.D. thesis, University of California, Berkeley, Department of Civil Engineering, 1966. 18. Fox, R. L. and Kapoor, M. P., IIA Minimization Method for the Solution of the Eigenproblem Arising in Structural Dynamics,1I Proceedings of the Second Conference on Matrix Methods in Structural Mechanics, Wright-Patterson Air Force Base, Ohio, AFFDL-TR-68-150, 1968, pp. 332- 345. 19. Francis, J. G. F., "The OR Transformation, Parts I and 11,11 Computer Journal, VoiD 4, 1961, pp. 265-271, Vol. 4, 1962, pp. 332-345. 20 . Fr i ed, I., II 0Pt i rna 1 Gr ad i en t Min i mi z a t ion Sc heme for Fin i teE 1erne nt Eigenproblems,I' Journal of Sound and Vibration, Vol. 20, 1972, pp. 333- 342. 21. Givens, J. W., "Numerical Computation of the Characteristic Value of a Real Symmetric Matrix,1I Oak Ridge National Laboratory Report No. 1574, distributed by Office of Technical Services, Washington, D. C., 1954. 22. Gupta, K. K., "Vibration of Frames and Other Structures with Banded Stiffness Matrix," International Journal for Numerical Methods in Engineering, Vol. 2, 1970, pp. 221-228. 23. Gupta, K. K., "Eigenproblem Solution by a Combined Sturm Sequence and Inverse Iteration Technique," International Journal for Numerical Methods in Engineering, Vol. 7, 1973, pp. 17-42. 24. Guyan, R. J., IIReduction of Stiffness and Mass Matrices,1I Journal of the Arrerican Institute of Aeronautics and Astronautics,n Vol. 3, No.2, 19"65, pp. 380. 54 25. Henshell, R. D. and Ong, J. H., "Automatic Masters for Eigenvalue Economisation," Journal of Earthquake Engineering and Structural Dynamics, Vol. 3, 1975, pp. 375-383. 26. Householder, A. S., The Principles of Numerical Analysis, McGraw-Hill Book Co., Inc., New York, 1953. Republished by Dover Publications, Inc., 1974. 27. Irons, B. M., "Eigenvalue Economisers in Vibration Problems," Journal of the Royal Aeronautical Society, Vol. 67, August 1963, pp. 526-528. 28. Irons, B. M., "Structural Eigenvalue Problems: Elimination of Unwanted Variables," Journal of the American Institute of Aeronautics and Astronautics, Vol. 3, No.5, 1965, pp. 961-962. 29. Jennings, A., "A Di rect Iteration Method of Obtaining the Latent Roots and Vectors of a Symmetric Matrix," Proceedings of Cambridge Philosophical Society, Vol. 63, 1967, pp. 755-765. 30. Jennings, A., "~~ass Condensation and Simultaneous Iteration for Vibration Problems," International Journal for Numerical Methods in Engineering, V~l. 6, 1973, pp. 543-552. 31. Jennings, A., Matrix Computation for Engineers and Scientists, John Wiley &Sons, 1977. 32. Jensen, P. S., liThe Solution of Large Symmetric Eigenproblems by Sectioning," Journal of the Society for Industrial and Applied Mathematics, Numerical Analysis, Vol. 9, No.4, December 1972, pp. 534-545. 33. Lanczos, C., "An Iteration Method for the Solution of the Eigenvalue Problem of Linear Differential and Integral Operators," Journal of Research of the National Bureau of Standards, Vol. 45, 1950, pp. 255-282. 34. Leung, A. Y., \IAn Accurate Method of Dynamic Condensation in Structural Analysis," International Journal for Numerical Methods in Engineering, Vol. 12, 1978, pp. 1705-1715. 35. Meirovitch, L., Analytical Method in Vibration, The Macmillan Co., London, 1967. 36. Melin, J. W. and Robinson, A. R., liThe Analysis of the Dynamic Response to an Aboveground Simply Supported Cylindrical Shell Subjected to Blast Loading,1I Air Force System Command, Kirtland, Air Force Base, New Mexico, AFSWC-TR-61-55, August 1961. 37. Myskis, A. D., Advanced Mathematics for Engineers, MIR Publishers, Moscow, 1975. 55 38. Peters, G. and Wilkinson, J. H., IIEigenvalues ofAx=ABx with Band Symmetric A and B,II Computer Journal, Vol. 12, 1969, pp. 398-404. 39. Peters, G. and Wilkinson, J. H., IIAx=ABx and the Generalized Eigenproblem,1I Journal of the Society for Industrial and Applied Mathematics, Numerical Analysis, Vol. 7, No.4, December 1970, pp. 479-492. 40. Przemieniecki, J. S., IIMatrix Structural Analysis of Substructures," Journal of the American Institute of Aeronautics and Astronautics, Vol. 1, No.1, 1963, pp. 138-147. 41. Rall, L. B., Newtons , Method for the Characteristic Value Problem AX=ABx," Journal of the Society for Industrial and Applied Mathematics, Vol. 9, No.2, June 1961. 42. Ramsden, J. N. and Stoker, J. R., "Mass Condensation: a Semi-Automatic Method for Reducing the Size of Vibration Problems," International Journal for Numerical Methods in Engineering, Vol. 1, 1969, pp. 333-349. 43. Rayleigh, Lord, The Theory of Sound, Second edition, Vol. 1, Macmillan Co., London, 1894. Republished by Dover Publications, Inc., 1945. 44. Robinson, A. R. and Harris, J. F., "Improving Approximate Eigenvalues and Eigenvectors," Journal of the Engineering Mechanics Division, Proceedings of the American Society of Civil Engineers, Vol. 97, No. EM2, April 1971, pp. 457-475. 45. Rosen, R. and Rubinstein, M. F., "Dynamic Analysis by Matrix Decomposition," Journal of the Engineering Mechanics Division, Proceedings of the American Society of Civil Engineers, No. EM2, April 1968, pp. 385-395. 46. Rutihauser, H. and Schwarz, H. R., "The LR Transformation Method for Symmetric Matrices," Numerische Mathematik, Vol. 5,1963, pp. 273-289. 47. Rutihauser, H., "Computational Aspects of F. L. Bauer's Simultaneous Iteration Method," Numerische Mathematik, Vol. 13, 1969, pp. 4-13. 48. Schwarz, H. R., Rutihauser, H. and Stiefel, E., Numerical Analysis of Symmetric Matrices, Prentice-Hall, Englewood Cliffs, New Jersey, 1973. 49. Schwarz, H. R., "The Eigenvalue Problem (A-AB)x=O for Symmetric Matrices of High Order," Computer Methods in Applied Mechanics and Engineering, Vol. 3, 1974, pp. 11-28. 50. Uhrig, R., "Reduction of the Number of Unknowns in the Displacement Method Applied to Kinetic Problems," Journal of Sound and Vibration, Vol. 4,1966, pp. 149-155. II 56 51. Wilkinson, J. H., The Algebraic Eigenvalue Problem, Oxford ,University Press, London, 1965. 52. Zienkiewicz, O. C., The Finite Element Method in Engineering Science, Third edition, McGraw-Hill Book Co., Ltd., London, 1977.' 57 APPENDIX A NONSINGULARITY OF THE COEFFICIENT MATRICES OF THE BASIC EQUATIONS Consider the basic equations (3.19) used for finding multiple or close eigenvalues and corresponding eigenvectors of the system Ax = A Bx (A.1 ) in which A and B are symmetric and both of order n, and B is positive defi ni te. Let the coefficient matrix of Eq. (3.19) be denoted by F, that is A - ~ .. (0)8 By(k) 11 F= i = m, m+ 1, ... ,m + s- 1<n (A.2) o where fl.· 11 (0) (i = m,m+ 1, ... ,m+ s-l) are initial approximate values of the multiple or close eigenvalues Ai (i = m,m+ 1, ... ,m+ s-l), and the s vectors (k) _ - (k) in Y - (k) - (k) - [Y m ' Ym+l , ... 'Ym+s-1 vectors in X = [x m ,x m+ l' ... 'x m+s- lJ. symmetric matrix. . J are approxlmate values of the eigenNote that F is an (n+s)x(n+s) The determinant of F ;s a continuous function of the approximate eigenvalue and eigenvectors. Hence, if F is nonsingular when the approximate 58 values in F become the exact ones, then it will be also nonsingular for close enough approximations. It is therefore sufficient for our purpose to use the exact eigenvalue and eigenvectors in Eq. (A.2) to prove nonsingularity. Let us take m = 1 for the convenience of the following presen- tation, then the resulting matrix G will be A - A.B 1 G = - BX _. - - _. -- - -. .- .- -- - XTB where X = (i = 1,2, ... ,s) (A.3) o [x 1,x 2 , .. ·,xs J. To find the determinant of G, we follow the idea that Robinson and Harris [44J used for showing the nonsingularity of the coefficient matrix of Eq. (2.10), that is, we consider the eigenvalues y'S and corresponding eigenvectors u's of the system *Gu = yBu (A.4) * GU = BUD (A. 5) or collectively where B* B a 0 I s = U= CUI ,U 2 '··· ,us 59 Is = unit matrix of order s I t may be veri fi ed by di rec t subs ti tu ti on tha t the (n + s) eigenvectors u·s and their corresponding eigenvalues y·s of Eq. (A.4) are: x. x. J J u -e.* e. J j = 1,2, ... ,s k = s+ 1,s+ 2, ... ,n o J (A.6) y x.1 where A.1 and are the eigenvalues and eigenvectors of the system Ax = ABx. - and e. -* form the diagonal matrices E and E* such that The vectors e. J E = J [e 1 ' e2' ... , es J = diag ( e 11 ' e 22' ... , e s s ) -* ,e-* ,··· ,esJ -* = diag (-.-1- , --1- , ... , -=-L ) E* = eel ess 2 ell e 22 e .. = L .. = L.. +-JL~. + 4 J1 J1 2 JJ J1 A. J A. 1 j = 1,2, ... ,s i,j = 1,2, ... ,s (A. 7) 60 From Eqs. (A.5), (A.6), and (A.7) det G = (det B* ) (det D) = (_I)s (det B) (A.B) In a similar way, the determinant of G for general m > 1 is (A.9) det G = (_l)s (det B) where the set S = [m,m+ 1 , ••• ,m+ s-lJ. which implies that det B > The matrix B is positive definite, O. 'Thus, if Ai (isS) is not equal or close to any of Ak (k = 1,2, ... ,n; kiS), the determinant of G is never equal to zero or close to zero, independently of whether A. (isS) are close, multiple, or 1 distinct. From Eq. (A.2), if s = 1, the matrix F becomes the coefficient matrix of Eq. (2.10), and by Eq. (A.9), the determinant of the matrix can be approxi rna ted by n F = (-1) (det B) IT k=1 (A k - Am) (A.I0) k#m Therefore, if Am ~ Am-I and Am 1 Am+l' the matrix F is also nonsingular. 61 APPENDIX B CONVERGENCE ANALYSIS The convergence analysis of the methods introduced in Chapters 2 and 3 will be presented. The eigenvalue problem we deal with here is Ax. 1 = A. 1 BX. (i 1 = 1, 2 , . . . , n ) (B. 1) in which A and B are symmetric matrices of order n, and B is positive definite. The eigenvectors x. (i 1 ~ 1,2, ... ,n) are assumed to be ortho- normalized with respect to B. B.1 Case of a Distinct Root Let us rewrite Eq. (2.7) used for improving approximate values of a distinct eigenvalue Aj and the corresponding eigenvector xjOf the system represented by Eq. (B.1): (B.2) T -x. (k) B/),x. - (k) = 0 J (B.3) J where ).. (k) and 'X.(k) are approximate values of ).. and J J J X.J after k iterations, and /),)..(k) and ~x.(k) are unknown incremental values of A.(k) and x.(k). J J J J L~t the approximate eigenvector xj(k) and the incremental vector ~Xj(k) be expanded in a series of the true eigenvectors, i.e., 62 n - (k) x·J c. . (k) -X. ,--, \ = L lJ 1 i=l b.c. . (k)- (B.4) x·1 lJ i=l in which cij(k) and the vicinity of - 6C Since the x.(k) is in ij (k) are scalar coefficients. J x~, J c .. (k) ~' max i i#j A. (k) J The errors in = E: c .. « (B.5) 1 JJ and X. (k) , y. J J ( k) and e.(k), may be defined by J A. - A.(k) J J A. « 1 J « 1 (B.6) where the values of y.(k) and e.(k) are very small compared with unity. J J the vectors c.(k) and c~(k) are defined by J J If 63 _ (k) T c. J = ( c1 . (k) ,c 2 . (k) , ... ,c. (k) ) = (0,0, c .. (k), 0, ... ,0) _ *( k) T c. J J J nJ (B.7) JJ The e. (k) then represents very closely the angle between the vectors cj(k) J -*( k) and c. J (s'ee Fig. 1). The task here is to estimate y. (k+l) and e. (k+l), 1 the e r ro r sin ~. (k - 1 ) and X. (k - 1 ) . 1 J J Let us substitute Eq. (B.4) into Eqs. (B.2) and premultiply by x.1 T to obtain =(i = 1, 2, ... ,n) c .. lJ (k) (B.8) Substitution of Eq. (B.4) into Eq. (B.3) and use of the orthonormality of the eigenvectors with respect to B results in n \' c .. (k) L i=l lJ The unknown quantities, b.A.(k) and J and (B.9). b.c •• (k) lJ b.C • • (k) 1J = ° (B.9) will be found from Eqs. (B.8) 64 Frorl Eq. (B.8) (i = 1, 2, ... ,n) (B.10) Now introduce Eq. (B .10) into Eq. (B.9) to obtain (B.11) where n a . (0) L = i =1 i tj c. . A. - A. J lJ J A. - A.(O) 1 J c.. JJ (k)!2 ; ( k) / I [ ) . (B.12) Using Eqs. (B.5) and (B.6), we get I I « 1 Is I « 1 a (B.13) 65 where h = max i i~j A. J A. 1 (0) - A. J (0) A. - « 1 :::: 1 J A. - A.(k) § = max i i~j 1 J A. - A.(O) (i (B.14) J 1 Therefore, from Eq. (B. 11), the I1A. (k) J may be approximated by I1A.(k) = A. - A.(k) + (A. J = 1, 2, ... ,n) J A.(O))S J J (B.15) J or A. (k+ 1) = A. (k) + J J I1A. (k) J (B.16) Substitute Eq. (B.15) into Eq. (B.I0) to obtain c .. (k) lJ A. - A.(O) + J J A. _ A.(O) 1 (1 + S) c .. (k) lJ J (i = 1, 2, ... ,n) (B.17) 66 or c. . lJ (k+1) = A. - A'(O) J = (1 + S) c .. (k) J A' - A.(O) 1 (i = 1, 2, ... ,n) lJ J (B.18) from which, it follows that c . . ( k+ 1) = (1 + S) c.. ( k) JJ (B.19) lJ The measure of the error in x.(k+1), e.(k+1), can now be found: J J \ n 8. (k+ 1) = J c .. (k+1) \ L ( i=l i]ij 1J (k+ 2l1/2 1)) ( Cjj ) ( ) n Aj - A. '\ 1/2 0 )2 (c .. (k) \ 2 .:...----A~-r-::(O~) = [ ~) ( ---=A i=l 1 JJ J i]ij < h 8.(k) (B.20) 1 where h is given in Eq. (B.14) and is very small compared with unity. find y.(k+1), the measure of the error of A.(k+1), we use Eqs. (B.6) and J (B.16), giving J To 67 Yj = (k+ 1) A. - A. (k+ 1 ) J J A. J 0 A. - A • ( ) = J A. J S J A· < A. ( 0 ) J J A. J 9 e~(k) J (B.21) by which Yj (k+2) < A. - J A.(O) J A. (B.22) J Substitution of Eq. (B.20) into Eq. (B.22) and use of Eq. (B.21) results in y. J (k+2) = h2 y. (k+1) J (B.23) Hence, it can be seen from Eqs. (B.20) and (B.23) that the jth eigenvector and eigenvalue converge linearly with errors multiplied by h h2 respectively in each iteration. B.2 (h «1) and Case of a Multiple Root The convergence analysis of the method for finding a multiple eigen- value and the corresponding eigenvectors of the system given in Eq. (B.1) will now be presented. 68 For convenience, but without loss of generality, the s low~st eigen- values are assumed to be equal, and the eigenvalue of multiplicity s is denoted by A*, i.e., A* = Al = A2 = . . . = AS' Let us rewrite Eqs. (3.21) and (3.17), which are the basic equations for improving approximate values of the multiple eigenvalue and corresponding eigenvectors, i.e., (~ \J _. 1 .l, I) ,...\ ' - , . . . , ;:, J (j = 1, 2, ... ,s) (B.25) and (B.26) where Is is the unit matrix of order s, and Y(k) -d.(k)T J = [y= 1 (k) ' y- (k) y- (k)] 2 , ... , s (k) (k) (k) (111j' fl2j , ... 'flsj ) (B.27) The A* 11 .. JJ = Al (k) (j = 1, 2, ... ,s) are approximations to the multiple eigenvalue = A2 = ••• = "s' and the Yj(k) (j = 1,2, ••• ,5) are approximations 69 to the true eigenvectors -x. ( j = 1, 2, ... ,s). The J - (k) and ~d. - (k) are ~y. J J unknown incremental vectors for y.(k) and a.(k). J J . Let the approximate eigenvectors y.(k) and the incremental vectors J AYj(k) be expanded in a series of the true eigenvectors xi (i = 1, 2, ... ,n), as in Eq. (B. 4), i. e. , n - (k) y. J = L c .. (k) X.1 lJ i=l n - (k) ~y . = \"' ~c .. (k) X. J Denoting the errors in G lJ (j = 1, 2, ... ,s) 1 (8.28) i=l lJ •. (k) an d y. (k) by J JJ J A. J (k) = y. - lJ·· (k) and e . (k) , we have J (k) JJ A. J J s (k) j e . (k) = J y. a. J (8.29) « 1 ( k) where ( s a.(k) =! J L c~.(k)~ 1 l i =1 r s/k) =) 11/2 n L ( i=s+l i 1J J 11/2 c~.(k) \. 1J I J (8.30) 70 The task here is to estimate the values of y.(k+1) and e.(k+1). J J . Now let the vectors c.(k) and 6c.(k), and the matrix C(k) be J J -c (k)T j II-c.(k) J U C(k) T = (c (k) = (IIuC . (k) , uc I · (k) , ... , uC II . (k)) 1J 2J nJ 1j' c 2j , ... , -- [-c (k) ,c- (k) , ... , -C (k)] 2 1 s Then, defining the matrix X = [xl' - (k) y. = X c.(k) y(k) = X C(k) J (k) x2 ' ... , J (B.31) X ], we may write Eq. (B.28) as n (j = 1, 2, ... ,s) (B.32) or (B.33) and - (k) 6y. J = X 6C. (k) J (j = 1, 2, ... ,s) (B.34) Substitution of Eqs. (B.32) - (B.34) into Eqs. (B.24) - (B.26), and premultiplication of Eq. (B.24) by XT results in (B.35) (B.36) (B.37) 71 where In = unit matrix of order n. Let us find lIC j (k) and lId j (k) using Eqs. (B.35) - (B.37). From Eq. (B.35) Substitution of Eq. (8.38) into Eq. (8.36) leads to (8.39) where F(k) = C(k) T (A _ ~ .. (0) I )-1 C(k) JJ n (8.40) Note that F(k) is a symmetric matrix of order s. Using Eq. (8.33), we can show that (8.41) where (B.42) 72 and the (l, m)th element of the symmetric matrix E(k) n e mm (k) = )...* \"' L-. i=s+1 n (k) = e,m C - ~ .. (0) (0) ).... - ~ .. 1 JJ )...* - i=s+1 (Co (k); 1m ) JJ elm(k), is (m= 1, 2, ... ,s) \ ( k) \ am Ci,t(k)) ).... 1 ).... .. (0) 1 JJ al( k) ~ ! c. (k) \ 1m (k) am \ ) (l, m = 1, 2, ... , S; l ~ m) The a.(k) (i = 1, 2, ... ,s) of Eq. (B.42) are defined in Eq. (B.30). 1 (B.43) By Eq. (B. 29), the absol ute val ues of elm (k) (l, m = 1, 2, ... ,s) are very small compared with unity, thus 1 (0) (;.* _ . 2(k) 2(k) 2(k) d1 ag (a 1 ' a2 ' ... ,as ) (8.44) Similarly, the values of the s components of the vector g.(k) can be found J fro m Eq. ( 8 . 40), i. e. , g .. (k) JJ g .. 1J ( k) a. ='* )... (k) a· J - (k) (0) lJ·· JJ lJ .. JJ (0) llij (k) (i = 1, 2, ... ,s, i ~ j) (8.45) 73 where n n· . JJ (k) = /., (C . (k) = - m- L \ a· J A* A (c .(k)\(c .(k)' m1 I~. ( k)} k i . ( 0) m=s+1 Am - 11jj Since Eq. a/ 2 (k) l~ (0) m=s+1 Am - 11·· JJ n n· . 1J . A m a. /Ia. 1 (8.46) i \ J J k 1 T + ) = (lll/k+l), 112/ k+1), ... , lls/k+l) by definition (see (8.2 )), Eqs. (8.39), (8.4 ), and (8.45) result in 11jj (k+1) ~ A* + (A* _ 11 • • t-<JJ (0)) lljj (k) (j = 1, 2, ... ,s) (i,j =1, 2, ... ,s; i f j) Substitution of Eq. (8.43) into Eq. (8.38) resul ts in s f::.C •. 1J ( k) .--. 1 = A. - 11 .. 1 JJ A* = (8.47) (0) 1m 11 . rTlJ (k+1 ) m=1 Ai ---(--0-'-) C i . Ai - 11jj (k) C. c..~-J J (k) + A. 1 (k) .. (0) c 1J A·1 - li·· JJ (0 ) Ai - lijj (i = 1, 2, ... ,n) (8.48) 74 or A* _ (0) 11jj (k) c .. (k+ 1) ~ (0) C ij 1J A. - 11 1 "'jj from which, since A* = Al = A2 = (i =A n ' c .. (k+1) ~ c .. (k) 1J (8.49) it follows that (i 1J = 1, 2, ... ,n) = 1, 2, ... ,s) (8.50) Thus, 1/2 n e . (k+ 1) = J [ (c. \ k+ 1) i 2 .l=S+ 1 1 J s ~ (c. ~ k+1 ))2 '- 1J i=l :5 h e.(k) (8.51) J where Y. J (k) 1· S def1· ned A* h = max i 1· n Eq . (8 . 30) , and - lljj (0) A. _ (0) 1 lljj (i = s+ 1, To find y.(k+1), the measure of the error of J and (8.47) are used, which results in s+2, ... , n ) 11 .. (k+1), JJ (8.52) Eqs. (8.29) 75 A* - y. (k+l) JJ = J (k+l) fl .. A* A* - fl .. (0) JJ ~ n· . JJ A* where n .. ( k) JJ I i s gi ven in Eq. (8.46). n .. ( k) JJ 12 ( k) (8.53) Its absolute value is max (i . 1 = s+l, s+2, ... ,n) (8.54) Therefore, from Eq. (8.53) 2(k) (8.55) z.: 8. J where A. z.: = max i 1 ( 0) . (i = s+l, s+2, ... ,n) (8.56) (j = 1, 2, ... ,s) (8.57) Ai - fljj From Eqs. (B.55) and (8.51), y.(k+2) = z.: J 8~(k+l) J = h2 y. J (k+ 1) 76 Hence, it can be seen from Eqs. (B.5n) and (8.51) that the multiple eigenvalue and the corresponding eigenvectors converge linearly with errors multiplied by h 2 (h« 1) and h respectively in each iteration. From Eqs. (B.47) and (B.49), (k) = A* for i = j = 0 for i "1 j 1i m c .. (k) = k-+oo lJ 0 lim k-+oo ll .. lJ (i,j = 1, 2, ... ,s) (B. 58) and (i = S+ 1, s+2, ... , n ; j = 1, 2, ... ,s) (8.59) which shows that as k-+oo, the vectors x.(k) (j = 1, 2, ... ,s) span the subspace J of Xj (j = 1, 2, ... , s) whose corresponding eigenvalue is multiple. Thus, the -x. (k) (J._- 1, 2, ... ,s) themselves are a set of true eigenvectors, orthoJ normalized with respect to 8. 77 APPENDIX C THE BASIC THEOREMS ON THE CONSTRAINED STATIONARY-VALUE PROBLEM Three theorems used for the development of the method for finding multiple or close eigenvalues and the corresponding eigenvectors will be presented. For convenience, two definitions will be given first. Definition 1 Let S be a set of positive integers p. (i 1 = 1,2, ... ,s) which are smaller than or equal to n, the order of the matrices A and B of the system Ax = ABx, i.e., S = (PI' P2"" ,ps). subspace spanned by the s eigenvectors x.1 Then, R is defined as the (isS). Definition 2 If no vector in the subspace R is orthogonal to all vectors - Yi (i = 1,2, ... ,s) (with respect to B), then the set of the vectors Yi (i = 1,2, ... ,s) is said to be an admissible frame with respect tb the subspace R. Theorem 1 With the above definitions, if none of the eigenvalues A.1 (isS) is equal to any A. (jiS), then among all admissible frames of vectors J y.1 (i = 1,2, ... ,s) a frame which renders w extremum in the following con- strained stationary-value problem spans the subspace R, and its stationary value is the sum of the eigenvalues A. (isS): 1 78 Find the stationary value of 5 w= L -T y.Ay. 1 (C.1) 1 i=l subject to -T - (i,j y .By. = 8 .. 1 J 1J = 1,2, ... ,5) (C.2) where 8 .. is the Kronecker delta. lJ Proof: S = x. 1 For convenience, but without loss of generality, the set (1,2, ... ,5 n) is taken, that ;5, R is the subspace spanned by (i = 1,2, ... ,5). Let the vectors series of the eigenvectors Xk (k = y.1 (i = 1,2, ... ,5) be expanded in a 1,2, ... ,n): n y.1 = L (k = 1,2, ... ,n) (C.3) k=l It will be shown that a solution of Eqs. (C.1) and (C.2) yields for k = 5+ 1, 5+2, •.. , n Sub s tit uti 0 n 0 f Eq . ( C. 3 ) into Eq 5 • w = 5 n .-- \' L L ;=1 ;=1 and 2 Ak c ki ( C. 1) (C.5) and (C. 2) re 5 u1ts ; n (C.6) 79 n L. k=l •• c ki c kj = cS lJ (i,j 1,2, ... ,s) = (C.7) Let us use the method of Lagrange multipliers to solve the stationary-value problem of Eqs. (C.6) and (C.7). 11·· lJ (i ,j = 1,2, ... ,s) and letting L= s n ~ \' L L i=l k=l s s L L 2 AkCki - \ i=l j=l Introducing the undetermined multipliers 11.· lJ = 11.·, Jl we have the Lagrangian n 11· . lJ (L C Ckj ki 8 .. ) k=l lJ Since the first derivatives of L with respect to the unknowns c · and k1 should vanish, s aL \"' - - = 2 ( Akc ki - L aC ki j=l n aL - - = L\"' ckic kj all·lJ. k=l 11 •• c ·) lJ kJ = a .. = a - alJ (C.8) 11 •. 1J (C.9) (C.10) We may write Eqs. (C.9) and (C.10) in matrix form as AC = CD (C.11) (C.12) where Is = unit matrix of order s 80 and the (k,i)th element of the nxs matrix C is c 1. (k = 1,2, ... -,n; k i = 1,2, ... ,s), and the (i ,j)th element of the sxs matrix D is ]1 •• lJ (i,j = 1,2, ... ,s). Let the matrix C be partitioned into·the two submatrices Cs and Cn-s where Cs is the sxs matrix having the elements of the first s rows of C, and C s is the remaining (n-s)xs matrix. I n- Then, Eqs. (C.11) and (C.12) may be written as AS Cs = CsO An-s Cn-s = Cn-s 0 CTC + cT C s s n-s n-s = Is (C.13) (C.14) (C.15) where Since ]1 •• lJ = ]1 •• Jl was taken (C.16) by which from Eq. (C.13) (C.17) Postmul tipl ication of Eq. (C.14) by C and use of Eq. (C.17) leads to A C CT n-s n-s s = C OCT n-s s (C.18) 81 Let u= C CT n-s s (C. 19) Then, Eq. (C.18) yields (i = 1,2, ... ,n-s; j = 1,2, ... ,s) (C.20) Since u = Cn-s CsT But the set of vectors y (i = = (C.21) 0 1,2, ... ,s) is an admissible frame with respect to B, from which it is not difficult to show that det CsT 1- 0 (C.22) From Eqs. (C.21) and (C.22), we obtain C n-s = 0 (C.23) or (i = 1,2, ... ,s; k = s+l,s+2, ... ,n) (C.24) This shows that the subspace spanned by the vectors Yi (i = 1,2, ... ,s) is the subspace of the eigenvectors xi (i = 1,2, ... ,s), which is to be proved here. Furthermore, from Eq. (C.6) we obtain ... Yetz Reference Room C1Y~1 Engineering Department BI06 C.E. BUilding University of Illinois Urbana, Illinois 61801 82 s w= ~ s .---. .:.._, L i=l ~. k=l s s .-- \' L- = L k=l 2 AkC ki ki=l c2 ki s = L (C.25) Ak k=l which implies that the stationary value is the sum of the eigenvalues Ak (k = 1,2, ... ,s). Theorem 2 If a frame of s vectors Yi (i = 1,2, ... ,s) which are mutually ortho- normal with respect to B spans the subspace R, then Eqs. (C.13), (C.14), and (C. 15) are satisfied, i.e., w is stationary. Proof: Since the vectors Yi (i = 1,2, ... ,s) are in the subspace R, and are orthonormal with respect to B, and or Hence, Eqs. (C.14) and (C.15) are satisfied. (C.26) Furthermore, we have 83 As Cs = (C s CsT) As Cs C = (C.27) (C T A C ) s s s s T We now define D by D = Cs AsC s ' then s Cs = Cs 0 (C.28) A which is equivalent to Eq. (C.13), i.e., Eq. (C.13) is also satisfied. Theorem 3 If the s vectors Yi (i = 1,2, ... ,s) span the subspace R of -xi (isS), then the Lagrange multipliers have the following properties: ~ .. lJ (i,j = 1,2, ... ,s) defined in Theorem 1 if the eigenvalues Ai (isS) are close together I ~i j I « I ~i i I for i ~ j (C.29) and if the eigenvalues are multiple, i.e., A* = A'1 (isS) ~ ~ Proof: .. = 0 1J for i f j .. = A* 11 For convenience, we take the set S = (1,2, ... ,s). (C.30) 84 Then, from Eq. (C.13) o= CT s II. C (C.31) 5 5 or 5 11· . lJ L = AkCkiCkj k=l s s L C C + [ =A ( Ak - Ai) Cki Ckj ki kj i k=l k=l (C.32) From Eq. (C.7) 5 r-"'"' L, C .C = k, KJ I • o.. (i,j IJ = 1,2, ... ,s) k=l Thus 5 11· • lJ = [ (A k - A. ) 1 C for i :f j C ki kj k=l 5 11 •• 11 = A. + 1 L (A 2 k (C.33) - A. ) C ki 1 k=l If the eigenvalues A. (i = 1,2, ... ,5) are close together, i.e., 1 I Ak - A; I «Ai (k:f i), then Eq. (C.33) implies that I 11; j I « I 11i i I for i :f j (C.34) 85 Furthermore, if all the eigenvalues A. (i 1 * = A A I = 1,2, ... ,5) are multiple, i.e., = A2 =, ... , = AS' then from Eq. (C.33) ll· . 1J =0 ll· . 11 = A.1 = A* (i = 1,2, ... ,5) (C.35) 86 NUMBER OF OPERATIONS FOR EIGENSOLUTIONS TABLE 1. Method Operation Proposed Method of Chapter 2 Multipl ication A-A.(O)B Factorization LOU Calculation Number of Operations n (m b + 1) J t A. (0) B =A- J nm a (ma + 3) Iterati on Multiplication Ax. (k) n (2m + 1) a J Bx . (k) n (2m b + 1) J r . ( k ) = Ax. ( k) _ II A. (k) Bx . (k) J J Factorization J J = F( k) LOU n (m + 1) a Solve Eq. (2.7) for llX.(k) and llA.(k) J Tota 1 2n (rna + 1) J 12 pn (m 2a Np n + 3m + 2m a b + 2) + T n (5m + 2m + 6) p a b Where A-L(O)B J - . : _Bx.(k) I .f 1 -----I [ _x.(k) B' J J r 0 i . Total number of iterations by the proposed method· Tp Total number of operations by the proposed method. > Tr . 87 TABLE 1. (Continued) Method Operation Proposed Method of Chapter 3 Multiplication A -fl .. (O)B Factorization LOU Calculation Number of Operations n (m + 1) b JJ =A - fl .. ( JJ 0) B Iterati on Multiplication Ay. (k) J Multiplication BYi(k) (i n (2m = 1,2, ..• ,s) a + 1) sn (2m + 1) b s r.(k) = Ay.(k) _ \' (k) B (k) J J L flij Yi Multiplication sn i=1 LOU Factorization = F( k) sn [rna + Solve Eq. (3.21) for t:,y. (k) and ~d. (k) n (2m J J \.Jhere y(k) [(k) y1 ,y 2 (k) ' ... ,y s (k) J a t (s + I)J + s + 1) 88 TABLE 1. (Continued) Method Operati on Robi nsonHa rri s Method Multiplication A-L(k)B Multiplication r. (k) = (A - >..(k)B) x.(k) Multiplication Bx. (k) Factorization LOU Calculation Number of Operations n (m a + 1 ) J J J J n (m a + 1) n (m + 1 ) b J i nma (ma + 5) + n = G(k) Solve Eq. (2.6) for lIx.(k) and lI>..(k) J J Total 2n (m a + 1) N r Hhere A-.\ . (k) B J : I I -Bx . (k) J G(k) -x· J Tr = Total (k) T B 0 number of iterations by the Robinson-Harris method. Tr Nr = Total number of operations by the Robinson-Harris method. < Tp. 89 TABLE 1. Method Opera ti on Subspace Iteration Method Factorization (Conti nued) Number of Operations Calculation LOU =A nma (m a + 3)/2 Iteration Multiplication BX(k-1) Solve for y(k) Ay(k) = BX(k-1) Multiplication A*(k) = y(k)T BX (k-1) Mu 1t i p1i cat ion By(k) Multiplication B*(k) Sol ve for Z(k) A*(k)Z(k) qn (2m b + 1) qn (2m a + 1) qn (q + 1)/2 qn (2m b + 1) = y(k)TBy(k) qn (q + 1)/2 = B*(k)Z(k)O(k) a (q 3 ) neglected and O(k) nq 2 Multiplication qn (2m + 4m a b Subtota 1 + 2q + 4) Sturm Seguence Check Multiplication A- A (k) B p Factorization LOU = A- /. (k) B Total Note: q max n (m + 1) p 2 Ts qn a (2m)m + 4 b + 2q + 4) + n (m + 3m + mb + 1) a a (2p. p+8) Ts Total number of iterations by subspace iteration method. NS Total number of ooerations by subspace iteration method. It is assumed that ma : m . b b 90 TABLE 2. EIGENVALUES OF THE PLANE FRAME PROBLEM (DISTINCT ROOTS) Method of Analysis Iteration Number Proposed Method 0 Eigenvalues (1) 2 (2) (3) 0.474744xlO O (0.36xlO- 2 )* 0.474744xlO O (0.44xlO -5 ) 0.443880xl0 1 (0.28xlO -1 ) 0.443876xl0 1 (0.45xlO -3 ) 0.132953xl0 2 (0. 77xlO- 1 ) 2 0.132921xI0 (0.38xlO -2 ) 0.474744xlO O (O.82xlO- 13 ) 0.443876xlQl -9 (O.27xlO ) 0.132921xl0 2 -6 (O.17xlO ) 3 0.284745xl0 0 (0.22xlO ) 0.284091xl0 -1 (0.27x10 ) 0 1 2 2 0.284091xlO 2 (0 .18x 10 -3 ) 0.284091xl0 2 (0.19xlO -5 ) O.474744x10 0 (0.36xlO -2 ) O.474744xlO O (0.44xlO -5 ) 0.443880xlO 1 (O.28xlO- 1) 0.443876xl0 1 (0.45xlO -3 ) 0.132953xl0 2 (0.77xlO- 1 ) 0.13292lxl0 2 (0.38xlO -2 ) O.284745xl0 2 (0.22xlO a 0.284091xl0 2 (0.27xlO -1 ) 0.474744xlO O (0.82xlO- 13 ) 0.443876xl0 1 (0.27xlO -9 ) O.132921xl0 2 (0.17x10 -6 ) 0.284091xl0 2 (0.18xlO -3 ) 0.284091xl0 2 3 ( 0.13xlO -8 ) * 2 O. 284091x 10 2 (0.23xlO -7 ) 4 Robinson Harris Method (4) Numbers in parentheses indicate errors in the approximate eigenvectors x.(k). J 91 TABLE 2. Method of Analysis Eigenvalues Iteration Number Subspace Iteration Method 1 2 3 4 5 6 7 8 9 10 11 12 ----* (Continued) (1) (2 ) (3) 0.476915x10 0 0.474744x10 0 (0.36xlO- 2 )* 0.465927x10 1 0.443880x10 1 (0.28x10 -1 ) 0.153636x10 0.474744x10 0 (O.44x10 -5 ) 0.474744x10 0 (0.llx10 -7 ) O.443876x10 1 (0.44x10 -3 ) 0.443876x10 1 (0.22x10 -4 ) 0.474744x10 0 (0.83x10- 1O ) 0.474744xlO O (0.71xlO- 12 ) 0.443876x10 1 (0.19x10 -5 ) 0.443876xlO l (0.16xlO -6 ) 0.474744xlO O (O.llxlO- 13 ) O.474744xlO O (O.13xlO- 13 ) 0.443876x10 1 (O.57xlO -8 ) 0.443876x10 1 (0.26xlO -9 ) 0.474744xlO O (0.16XlO- 13 ) 0.474744xlO O (0. 13xlO- 13 ) 0.443876xlO l (0.17x10- 10 ) 0.443876xlO l (0.93xIO- 12 ) 0.474744xlO O (0.13x10- 13 ) 0.474744xlO O (O.28xlO- 13 ) 0.443876xlO l (0.52xl0- 13 ) O.443876xlO l (O.50xl0- 13 ) (4 ) 2 0.132953x10 2 (0.77x10 -1 ) 2 0.132921x10 (0.35xlO -2 ) 0.132921x10 2 -3 (O.33x10 ) 0.132921x10 2 (0.90xlO -4 ) 0.132921xlO 2 (0.28xlO -4 ) 0.132921xl0 2 (0.28xlO -5 ) 0.132921xl0 2 -6 (0.28xlO ) 0.132921x10 2 (0.53xlO -7 ) 0.132921xl0 2 (0.91xlO -8 ) 0.132921xl0 2 (0.12xlO -8 ) O.132921xlO 2 (O.18xlO -9 ) :wnDers 1n parentheses i nd i ca te errors in the approximate floenvectors x.(k). - J 2 0.369761x10 2 0.284745x10 (0.22xlO 2) 0.284116x10 (0.19xlO -1 ) 0.284099x10 (0.48xlO -2 ) 0.284094xl0 (0.34xlO -2 ) 0.284092x10 (0.36xlO -2 ) 2 2 2 2 0.284091xl0 2 (O.llxlO -2 ) 2 0.284091xl0 -3 (0.14xlO ) 0.284091x10 2 (0.32xlO -4 ) 0.284091xl0 2 (0.99xlO -5 ) 0.284091xl0 -5 (0.33xlO ) 2 O.284091xlO 2 -5 (0.12xlO ) 92 TABLE 3. Method of Analysis Proposed Method EIGENVALUES OF THE CIRCULAR ARCH PROBLEM (DISTINCT ROOTS) Eigenvalues/A O* Iteration Number (1) O.102714xlO- 2 0.102640xlO- 2 (O.59xlO- 3 )** 0 1 O.102640xlO- 2 (0.26xlO -8 ) 2 3 (2 ) 0.9l95l6xlO- 2 0.909467xlO- 2 -2 (O.53xlO ) 0.909467xlO- 2 -5 (O.lOxlO ) 0.909467xlO- 2 (O:lOxlO- 8 ) O.363073xlO- l O.34l422xlO- l (O.35xlO -1 ) O.34l448xlO- l (O.2lxlO -3 ) O.34l448xlO- l (O.17xlO -5 ) O.34l448xlO- l -7 (O.33xlO ) 4 RobinsonHarris Method (3) O.102714xlO- 2 O.102640xlO- 2 -3 (O.59xlO ) 0.102640xlO- 2 (0.26xlO -8 ) 0 1 2 3 2 0.9l95l6xlO- 2 O.909467xlO- 2 -2 (0.53xlO ) 0.909467xlO- 2 (O.lOxlO -5 ) O.363073xlO- l O.34l422xlO- l -1 (0.35xlO ) O.34l448xlO- l (O.2lx10 -3 ). 0.909467x10- 2 (0. 45xlO- 12 ) O.341448x10 -1 -8 (0.20x10 ) 2 * . AO = E/pa (1-\1 ). ** Numbers in parentheses indicate errors in the approximate eigenvectors x.(k). J 93 TABLE 3. Method of Ana lysi s Ei genva 1ues/ AO * Itera ti on Number Subspace Iterat; on Method 2 3 4 5 6 7 (Continued) ( 1) (2) 0.102714x10- 2 0.102640x10- 2 (0.59x10- 3 )** 0.919516x10- 2 0.909468x10- 2 (0.52x10 -2 ) 0.102640x10- 2 (0.41x10 -6 ) 0.102640x10- 2 (0. 58xlO- 9) 0.909467x10- 2 (0.62x10 -4 ) 0.909467x10- 2 (0.12x10 -5 ) O.102640X10- 2 (0.12xlO- ll ) 0.102640xlO- 2 (0.23x10- 13 ) 0.909467x10- 2 (O.28xlO -7 ) 0.909467x1O- 2 (0.79xlO -9 ) 0.102640x10- 2 (0. 12x10- 13 ) 0.909467x10- 2 (0.25xlO- 10 ) (3) 0.363073x10- 1 0.341476x10- 1 -1 (0.32x10 ) 0.341448x10- 1 (0.26x10 -2 ) 0.341448x10- 1 -3 (0.23x10 ) 0.341448x10- 1 -4 (0.21x10 ) 0.341448x10- 1 -5 (0.19xlO ) 0.341448x10- 1 (O.17xlO 2 * . AO = E/pa 2 (I-v). ** Numbers in parentheses indicate errors in the approximate eigenvectors xj(k). -6 ) 94 TABLE 4. EIGENVALUES OF THE SQUARE PLATE PROBLEM (DOUBLE ROOTS) Method of Analysis Iteration Number Eigenvalues/a* Proposed Method 0 2 Subspace Iteration Method 2 3 4 5 .. :l = n (2 ) (I) (3) O.375840x1O I O.2302I2x10 2 O.375838xlO l O.230I2lx10 2 (O.80xlO- 3 )** (O.lOXlO- l ) 2 O.375838xlO l O.230I2lxlO ( O. 44 xl 0-11 ) (O.46xIO- 7) O.2302I2x10 2 2 O.230I2lxlO -1 (O.lOxlO ) O.230l2lxlO 2 (O.46xlO- 7 ) O.530627xlO 2 O.529932xlO 2 O.375840xlO l O.375838xlO l (O.80XIO- 3 ) O.2302l2xlO 2 O.230I2lxlO 2 (O.lOxIO- l ) O.2302l2xlO 2 O.230l2lxlO 2 (O.IOXIO- I ) O.530627xlO 2 O.529932xlO 2 (O.2IxlO -1 ) O.375838xlO l (O.65xlO -6 ) O.375838xlO l (O.53xlO- 9) O.230l21xlO 2 (O.l1XIO- 3 ) O.230l2lxlO 2 (O.13xIO -5 ) O.230l2lx10 2 (O.llxlO- 3 ) 2 O.230l21xlO -5 (O.13xIO ) O.529932x10 2 (O.43xlO -3 ) O.529932xlO 2 (O.90xlO -5 ) O.375838xlO l (O.44xIO- 12 ) O.230I2lxlO -7 (O.15xlO ) 2 O.230l21xlO 2 (O.15xlO -7 ) O.529932xlO 2 (O.19xlO -6 ) 40 /(a 4p). where 0 = Eh 3/12(1 - }). e e in parentheses indicate errors in the approximate e1genvectors Yj(k). ~umbers (4) (O.2lxlO -1 ) O.529932xlO 2 (O.59xlO -6 ) 95 TABLE 5. EIGENVALUES OF THE RECTANGULAR PLATE PROBLEM (CLOSE ROOTS) Method of Analysis Eigenvalues/o.* Iteration Number (1) (2 ) (3) 0.368470x10 1 0.268468x10 1 (0.80xlO -3 ) O.222957xl0 2 O.222868xl0 2 (O.lOxlO -1 ) 0.228454xl0 2 0.368468x10 1 (0.44xlO- ll ) 2 0 Proposed Method Subspace Iteration Method 3 4 5 (4) 2 0.222868x10 2 (0.47xlO -7 ) 0.228264x10 2 (O.IOxlO- I ) 0.228364xl0 2 (O.45xlO -7 ) O.5l9533x10 (0.59xlO -6 ) 0.368470x1O l 0.368468x1O l (0.80xlO -3 ) O.222957xl0 2 0.222868xl0 2 (O.lOxlO -1 ) 0.228454x10 2 2 O.228364x10 (O.lOxIO -1 ) 0.5202l5x10 2 2 O.5l9533xlO (0.2lxlO -1 ) 0.368468xlO I (0.65xI0 -6 ) 0.368468x1O l (0.53xlO -9 ) 0.222868x10 2 (O.llxlO -3 ) 0.222868xl0 2 (O.13xIO -5 ) 0.228364xl0 2 (O.llxIO -3 ) 0.228364xl0 2 (O.13xlO -5 ) O.5l9533xl0 (0.43xlO -3 ) 0.519533x10 (O.90xlO -5) O.368468x1O l (0.46xlO- 12 ) 0.222868xlO (0.15xlO -7 ) 2 O.228364x10 2 (0.14xlO -7 ) O.5l9533xlO -6 (O.19xI0 ) h3/ l2(l - ':r- e i i,a 4 " ) ,w~ere 0e -- E \! 2) • r,urtJers In parentheses indicate errors in the approximate elqenvectors y.(k). - J O.520215x10 2 2 O.519522x10 -1 (0.2lxlO ) 2 2 2 2 96 TABLE 6. COMPARISON OF THE TOTAL NUMBER OF OPERATIONS Input Data Problem Type Frame n rna mb p q 330 35 35 4 10 4 o3 NLmber of I terati ons Tp Tr T s 10 9 8 Arch 22 5 9 Plate 39 16 16 4 9 8 Number of Operations Np N r Ns Ratio Nr/Np N/N p 3.50xl0 6 4.57xl0 6 9.72xl0 6 1. 31 2.78 3 4 1.10 1. 98 7 8.87xl0 9.77xl0 3 1. 76xl0 5 5 5 1. 27xl0 1. 73 2.20xl0 12 Note n Order of stiffness and mass matrices rna Average half bandwidth of stiffness matrix mb Average half bandwidth of mass matrix p Number of eigenvalues and eigenvectors sought q Number of iteration vectors, q Tp Number of iterations by the proposed method Tr Number of iterations by the Robinson-Harris method Ts = max (2p, p+8) Number of iterations by the subspace iteration method N Total number of operations by the proposed method N Total number of operations by the Robinson-Harris method Ns Total number of operations by the subspace iteration method p r 97 TABLE 7. COMPARISON BETWEEN THE THEORETICAL CONVERGENCE RATES FOR EIGENVECTORS AND THE NUMERICAL RESULTS - FRAME PROBLEM (DISTINCT ROOTS) Method of Analysis Proposed Method Subspace Iteration Method Eigenvalue Number I tera ti on Number 1 2 3 Theory 2 3 4 5 6 7 8 9 10 11 Theory * 2 3 4 -7 4.4x10 -5 -5 3.6x10- 4 -3 6.7xlO 1.1xlO -2 1.2x10- 2 1. 2x10 -2 1. 6xlO -2 5.0xlO- 2 8.7xlO- 2 8.4xlO -2 3.6xl0-2 4.6xl0- 2 6.5xlO- 2 5.5xl0- 2 5.6xl0 -2 4.5xlO -2 -2 9.4xlO 2.7xl0 -1 3.1xlO- 1 1. Oxl0- 1 1. Ox10- 1 1 1.9xI0-.L 1. 7xl0- 1 1.3xl0- 1 1.5xlO- 1 1.6xlO- 1 8.6xlO -2 2.5xlO -1 7.1xlO- 1 -0 1. Ix 10 -1 3.1x10 -1 1.3xl0 -1 2.3xl0 3.1xl0 -1 3.3xl0 -1 3.6xlO -1 3.3xl0 -1 1. 9x10 -8 6.0x10 1. 3xlO- 8 1. Ox10 1.2x10 -3 2.5xlO- 3 7.5xlO- 3 8.6xl0- 3 1. 5xl0- 3 * * * * * 5.6xl0 -3 * 5.2xl0 -2 Errors too small for comparison because of round-off error. 98 TABLE 8. COMPARISON BETWEEN THE THEORETICAL CONVERGENCE RATES FOR EIGENVECTORS AND THE NUMERICAL RESULTS - SQUARE PLATE PROBLEM (DOUBLE ROOTS) Method of Analysis Iteration Number Proposed Method 1 Theory Subspace Iteration Method 2 3 4 Theory Eigenvalue Number 1 5.5xlO -9 1. Ox 10-6 8.1x1O- 4 8.2x1O -4 8.3x1O -4 -2 1.1xlO 2 3 4.6xlO -6 4.7x1O -4 1.lx1O -2 -2 1. 2x1O 1. 2x1O- 2 6.8xlO --2 4.6xlO -6 4.7x1O- 4 1.1x1O-2 1. 2x1O- 2 -2 1. 2x1O . -2 6.8x1O 4 2.8x1O -5 2.3x1O -3 2.0x1O -2 -2 2.1x1O 2.1x1O -2 -1 1. 6x1O 99 TABLE 9. Method of Ana 1ys i s NUMERICAL CONVERGENCE RATES FOR EIGENVECTORS RECTANGULAR PLATE PROBLEM (CLOSE ROOTS) Eigenvalue Number Itera ti on Number Proposed Method 'Subspace Iteration Method 2 3 4 Theory 2 3 4 5.5x1O-9 4.7x1O -6 4.5x1O -6 2.8x1O -5 8.1x1O -4 8.2x1O -4 8.7x1O -4 1.1x1O -2 1. 1xlO-2 1.2xlO-2 1. 2x1O -2 6.8x1O- 2 1.lxlO-2 1. 2x1O -2 1.lxlO-2 -2 7.0x1O -2 2.0x1O -2 2.1x1O -2 2.1x1O -1 1. 6x1O 100 ~=====-----------------------------~Xj ~~~~--------~+-----------------~Cj e~ k+1) J FIG. 1 ESTIMATION OF ERRORS IN APPROXIMATE EIGENVECTORS 101 E LO o r() II E LO Q r() /77T 7'''77 777r 77'?7 77 77.'17 10 at 6.1 m "~ = 61.0 m For All Beams and Columns Area of Cross-Section Moment of Inertia of Cross-Section Young's Modulus Mass Density 7777- )'7."7 A = 2.787xlO- l m2 I = 8.631xlO- 3 m4 E = 2.068xl0 10 Pa 3 4 p = 1.602xl0 kg/m FIG. 2 TEN-STORY, TEN-BAY PLANE FRAME 77. ~ )7.~ 474:NP:7l6:lab 78u474-6l9 ONR Code 474 May 1980 DISTRIBUTION LIST for UNCLASSIFIED TECHNICAL REPORTS The ONR Structural Mechanics Contract Research Program This list consists of: Part 1 - Government Activities Part 2 - Contractors and Other Technical Collaborators Notes: Except as otherwise indicated, forward one CODY of all Unclassified Technical Reports to each of the addressees listed herein. Where more than one attention addressee indicated, the individual copies of the report should be mailed separately. ~s 474:NP:7l6:lab 78u474-6l9 Part 1 - Government Administrative and Liaison Activities Office of Naval Research Department of the Navy Arlington, Virginia 22217 Attn: Code 474 (2) Code 471 Code 200 Director Office of Naval Research Eastern/Central Regional Office 666 Summer Street Boston, Massachusetts 02210 Director Office of Naval Research Branch. Office 536 South Clark Street Chicago; Illinois 60605 Director Office of Naval Research New York Area Office 715 Broadway - 5th Floor New York, New York 10003 Director Office of Western 1030 East Pasadena, Naval Research Regional Office Green Street California 91106 Naval Research Laboratory (6) Code 2627 Washington, D.C. 20375 Defense Technical Information Center (12) Cameron Stat ion Alexandria, Virginia 22314 Undersea Explosion Research Division Naval Ship Research and Development Center Norfolk Naval Shipyard Portsmouth, Virginia 23709 Attn: Dr. E. Palmer, Code 177 Navy (Con't.) Naval Research Laboratory Washington, D.C. 20375 Attn: Code 8400 8410 8430 8440 6300' 6390 6380 David W. Taylor Naval Ship Research and Development Center Annapolis, Maryland 21402 Attn: Code 2740 28 281 Naval Weapons Center China Lake, California 93555 Attn: Code 4062 4520 Commanding Officer Naval Civil Engineering Laboratory Code L3l Port Hueneme, California 93041 Naval Surface Weapons Center White Oak Silver Spring, Maryland 20910 Attn: Code R-10 G-402 K-82 Technical Director Naval Ocean Systems Center San Diego, California 92152 Supervisor of Shipbuilding U.S. Navy Newport News, Virginia 23607 Navy Underwater Sound Reference Division Naval Research Laboratory P.O. Box 8337 Orlando, Florida 32806 Chief of Naval Operations Department of the Navy Washington, D.C. 20350 Attn: Code OP-098 474:NP:7l6:lab 78u474-6l9 Navy (Con't.) Navy (Con' t.) Strategic Systems Project Office Department of the Navy Washington, D.C. 20376 Attn: NSP-200 Commander and Director David W. Taylor Naval Ship Research and Development Center Bethesda, Maryland 20084 At tn: Code 042 17 172 173 174 1800 1844 012.2 1900 1901 1945 1960 1962 Naval Air Systems Command Department of the Navy Washington, D.C. 20361 Attn: Code 5302 (Aerospace and Structures) 604 (Technical Library) 320B (Structures) Naval Air Development Center Warminster, Pennsylvania 18974 Attn: Aerospace Mechanics Code 606 u.S. Naval Academy Engineering Department Annapolis, Maryland 21402 Naval Facilities Engineering Command 200 Stovall Street Alexandria, Virginia 22332 Attn: Code 03 (Research and Development) 04B 045 14114 (Technical Library) Naval Sea Systems Command Department of the Navy Washington, D.C. 20362 Attn: Code 05H 312 322 323 05R 32R Naval Underwater Systems Center Newport, Rhode Island 02840 Attn: Bruce Sandman, Code 3634 Naval Surface Weapons Center Dahlgren Laboratory Dahlgren, Virginia 22448 Attn: Code G04 G20 Technical Director Mare Island Naval Shipyard Vallejo, California 94592 U.S. Naval Postgraduate School Library Code 0384 Monterey, California 93940 Webb Institute of Naval Architecture Attn: Librarian Crescent Beach Road, Glen Cove Long Island, New York 11542 Commanding Officer (2) U.S. Army Research Office P • o. Bo x 12211 Research Triangle Park, NC 27709 Attn: Mr. J. J. Murray, CRD-AA-IP 474:NP:7l6:lab 78u474-6l9 Army (Con't.) Air Force (Con't.) Watervliet Arsenal MAGGS Research Center Watervliet, New York 12189 Attn: Director of Research Chief, Civil Engineering Branch WLRC, Research Division Air Force Weapons Lab-oratory Kirtland Air Force Base Albuquerque, New Mexico 87117 u.s. Army Materials and Mechanics Research Center Watertown, Massachusetts 02172 Attn: Dr. R. Shea, DRXMR-T u.s. Army Missile Research and Development Center Redstone Scientific Information Center Chief, Document Section Redstone Arsenal, Alabama 35809 Air Force Office of Scientific Research Bolling Air Force Base Washington, D.C. 20332 Attn: Mechanics Division Department of the Air Force Air University Library Maxwell Air Force Base Montgomery, Alabama 36112 Other Government Activities Army Research and Development Center Fort Belvoir, Virginia 22060 NASA National Aeronautics and Space Administration Structures Research Division Langley Research Center Langley Station Hampton, Virginia 23365 National Aeronautics.and Space Administration Associate Administrator for Advanced Research and Technology Washington, D.C. 20546 Air Force Wright-Patterson Air Force Base Dayton, Ohio 45433 Attn: AFFDL (FB) ( FBR) ( FBE) ( FBS) AFML (HBM) Chief Applied Mechanics Group U.S. Air Force Institute of Technology Wright-Patterson Air Force Base Dayton, Ohio 45433 Commanda.nt Chief, Testing and Development Division U.S. Coast Guard 1300 E Street, NW. Washington, D.C. 20226 Technical Director Marine Corps Development and Education Command Quantico, Virginia 22134 Director Defense Research and Engineering Technical Library Room 3C128 The Pentagon Washington, D.C. 20301 Dr. M. Gaus National Science Foundation Environmental Research Division Washington, D.C. 20550 Library of Congress Science and Technology Division Washington, D.C. 20540 Director Defense Nuclear Agency Washington, D.C. 20305 Attn: SPSS 474:NP:7l6:lab 78u474-6l9 Other Government Activities (Con't) Universities (Con't) Mr. Jerome Persh Staff Specialist for Materials and Structures OUSDR&E, The Pentagon Room 3D1089 Washington, D.C. 20301 Dr. Harold Liebowitz, Dean School of Engineering and Applied Science George Washington University Washington, D.C. 20052 Chief, Airframe and Equipment Branch FS-120 Office of Flight Standards Federal Aviation Agency Washington, D.C. 20553 National Academy of Sciences National Research Council Ship Hull Research Committee 2101 Constitution Avenue Washington, D.C. 20418 Attn: Mr. A. R. Lytle National Science Foundation Engineering Mechanics Section Division of Engineering Washington, D.C. 20550 Picatinny Arsenal Plastics Technical Evaluation Center Attn: Technical Information Section Dover, New Jersey 07801 Maritime Administration Office of Maritime Technology 14th and Constitution Avenue, NW. Washington, D.C. 20230 PART 2 - Contractors and Other Technical Collaborators Universities Dr. J. Tinsley Oden University of Texas at Austin 345 Engineering Science Building Austin, Texas 78712 Professor Julius Miklowitz California Institute of Technology Division of Engineering and Applied Sciences Pasadena, California 91109 Professor Eli Sternberg California Institute of Technology Division of Engineering and Applied Sciences Pasadena, California 91109 Professor Paul M. Naghdi University of California Department of Mechanical Engineering Berkeley, California 94720 Professor A. J. Dure11i Oakland University School of Engineering Rochester, Missouri 48063 Professor F. L. DL~aggio Columbia University Department of Civil Engineering New York, New York 10027 Professor Norman Jones The University of Liverpool Department of Mechanical Engineering P. O. Box 147 Brownlow Hill Liverpool L69 3BX England Professor E. J. Skudrzyk Pennsylvania State University Applied Research Laboratory Department of Physics State College, Pennsylvania 16801 Professor J. Klosner Polytechnic Institute of New York Department of Mechanical and Aerospace Engineering 333 Jay Street Brooklyn, New York 11201 Professor R. A. Schapery Texas A&M University Department of Civil Engineering College Station, Texas 77843 474:NP:7l6:lab 78u474-6l9 Universities (Con't.) Universities (Con't) Professor Walter D. Pilkey University of Virginia Research Laboratories for the Engineering Sciences and Applied Sciences Charlottesville, Virginia 22901 Professor A. C. Eringen Princeton University Department of Aerospace and Mechanical Sciences Princeton, New Jersey 08540 Professor K. D. Willmert Clarkson College of Technology Department of Mechanical Engineering Potsdam, New York 13676 Dr. Walter E. Haisler Texas A&M University Aerospace Engineering Department College Station, Texas 77843 Dr. Hussein A. Kamel University of Arizona Department of Aerospace and Mechanical Engineering Tucson, Arizona 85721 Dr. S. J. Fenves Carnegie-Mellon University Department of Civil Engineering Schenley Park Pittsburgh, Pennsylvania 15213 Dr. Ronald L. Huston Department of Engineering Analysis University of Cincinnati Cincinnati, Ohio 45221 Professor G. C. M. Sih Lehigh University Institute of Fracture and Solid Mechanics Bethlehem, Pennsylvania 18015 Professor Albert S. Kobayashi University of Washington Department of Mechanical Engineering Seattle, Washington 98105 Professor Daniel Frederick Virginia Polytechnic Institute and Department of Engineering Mechanics Blacksburg, Virginia 24061 Professor E. H. Lee Stanford University Division of Engineering Mechanics Stanford, California 94305 Professor Albert I. King Wayne State University Biomechanics Research Center Detroit, Michigan 48202 Dr. V. R. Hodgson Wayne State University School of Medicine Detroit, Michigan 48202 Dean B. A. Boley Northwestern University Department of Civil Engineering Evanston, Illinois 60201 Professor P. G. Hodge, Jr. University of Minnesota Department of Aerospace Engineering and Mechanics Minneapolis, Minnesota 55455 Dr. D. C. Drucker University of Illinois Dean of Engineering Urbana, Illinois 61801 Professor N. M. Newmark University of Illinois Department of Civil Engineering Urbana, Illinois 61803 Professor E. Reissner University of California, San Diego Department of Applied Mechanics La Jolla, California 92037 Professor William A. Nash University of Massachusetts Department of Mechanics and Aerospace Engineering Amherst, Massachusetts 01002 474:NP:7l6:lab 78u474-6l9 Universities (Con't) Universities (Con't) Professor G. Herrmann Stanford University Department of Applied Mechanics Stanford, California 94305 Professor Joseph L. Rose Drexel University Department of Mechanical Engineering and Mechanics Philadelphia, Pennsylvania 19104 Professor J. D. Achenbach Northwest University Department of Civil Engineering Evanston, Illinois 60201 Professor S. B. Dong University of California Department of Mechanics Los Angeles, California 90024 Professor Burt Paul University of Pennsylvania Towne School of Civil and Mechanical Engineering Philadelphia, Pennsylvania 19104 Professor H. W. Liu Syracuse University Department of Chemical Engineering and Metallurgy Syracuse, New York 13210 Professor S. Bodner Technion R&D Foundation Haifa, Israel Professor Werner Goldsmith University of California Department of ~echanical Engineering Berkeley, California 94720 Professor R. S. Riv1in Lehigh University Center for the Application of Mathematics Bethlehem, Pennsylvania 18015 Professor F, A. Cozzarelli State University of New York at Buffalo Division of Interdisciplinary Studies Karr Parker Engineering Building Chemis~ry Road Buffalo, New York 14214 Professor B. K. Donaldson University of Maryland Aerospace Engineering Department College Park, Maryland 20742 Professor Joseph A. Clark Catholic University of America Department of Mechanical Engineering Washington, D.C. 20064 Dr. Samuel B. Batdorf University of California School of Engineering and Applied Science Los Angeles, California 90024 Professor Isaac Fried Boston University Department of Mathematics Boston, Massachusetts 02215 Professor E. Krempl Rensselaer Polytechnic Institute Division of Engineering Engineering Mechanics Troy, New York 12181 Dr. Jack R. Vinson University of Delaware Department of Mechanical and Aerospace Engineering and the Center for Composite Materials Newark, Delaware 19711 Dr. J. Duffy Brown University Division of Engineering Providence, Rhode Island 02912 Dr. J. L. Swedlow Carnegie-Mellon University Department of Mechanical Engineering Pittsburgh, Pennsylvania 15213 474:NP:7l6:1ab 78u474-6l9 Universities (Con't) Universities (Con't) Dr. V. K. Varadan Ohio State University Research Foundation Department of Engineering Mechanics Columbus, Ohio 43210 Professor V. H. ~eubert Pennsylvania State University Department of Engineering Science and Mechanics . University Park, Pennsylvania 16802 Dr. Z. Hashin University of Pennsylvania Department of Metallurgy and Materials Science College of Engineering and Applied Science Philadelphia, Pennsylvania 19104 Dr. Jackson C. S. Yang University of Maryland Department of Mechanical Engineering College Park, Maryland 20742 Professor T. Y. Chang University of Akron Department of Civil Engineering Akron, ohio 44325 Professor Charles W. Bert University of Oklahoma School of Aerospace, Mechanical, and Nuclear Engineering Norman, Oklahoma 73019 Professor Satya N. Atluri Georgia Institute of Technology School of Engineering and Mechanics Atlanta, Georgia 30332 Professor Graham F. Carey University of Texas at Austin Department of Aerospace Engineering and Engineering Mechanics Austin, Texas 78712 Dr. S. S. Wang University of Illinois De?artment of Theoretical and Applied Mechanics Urbana, Illinois 61801 Professor J. F. Abel Cornell University Department of Theoretical and Applied Mechanics Ithaca, New York 14853 Professor A. W. Leissa Ohio State University Department of Engineering Mechanics Columbus, Ohio 43212 Professor C. A. Brebbia University of California, Irvine Department of Civil Engineering School of Engineering Irvine, California 92717 Dr. George T. ~ahn Vanderbilt University Mechanical Engineering and Materials Science Nashville, Tennessee 37235 Dean Richard H. Gallagher University of Arizona College of Engineering Tucson, Arizona 85721 Professor E. F. Rybicki The University of Tulsa Department of Mechanical Engineering Tulsa, Oklahoma 74104 Dr. R. Haftka Illinois Institute of Technology Department of Mechanics and Mechanical and Aerospace Engineering Chicago, Illinois 60616 Professor J. G. de Oliveira Massachusetts Institute of Technology Department of Ocean Engineering 77 Massachusetts Avenue Cambridge, Massachusetts 02139 Dr. Bernard W. Shaffer Polytechnic Institute of New York Route 110 Farmingdale, New York 11735 474:NP:7l6:lab 78u474-6l9 Industry and Research Institutes Industry and Research Institutes (Con't) Dr. Norman Hobbs Kaman AviDyne Division of Kaman Sciences Corporation Burlington, Massachusetts 01803 Dr. T. L. Geers Lockheed Missiles and Space Company 3251 Hanover Street Palo Alto, California 94304 Argonne National Laboratory Library Services Department 9700 South Cass Avenue Argonne, Illinois 60440 Dr. M. C. Junger Cambridge Acoustical Associates 54 Rindge Avenue Extension Cambridge, Massachusetts 02140 Mr. J. H. Torrance General Dynamics Corporation Electric Boat Division Groton, Connecticut 06340 Mr. William Caywood Applied Physics Laboratory Johns Hopkins Road Laurel, Maryland 20810 Dr. Robert E. Dunham Pacifica Technology P.O. Box 148 Del Mar, California 92014 Dr. M. F. Kanninen Battelle Columbus Laboratories 505 King Avenue Columbus, Ohio 43201 Dr. J. E. Greenspon J. G. Engineering Research Associates 3831 Menlo Drive Baltimore, Maryland 21215 Dr. A. A. Hochrein Daedalean Associates, Inc. Springlake Research Road 15110 Frederick Road Woodbine, Maryland 21797 Newport News Shipbuilding and Dry Dock Company Library Newport News, Virginia 23607 Dr. James W. Jones Swanson Service Corporation P.O. Box 5415 Huntington Beach, California 92646 Dr. W. F. Bozich McDonnell Douglas Corporation 5301 Bolsa Avenue Huntington Beach, California 92647 Dr. Robert E. Nickell Applied Science and Technology 3344 North Torrey Pines Court Suite 220 La Jolla, California 92037 Dr. H. N. Abramson Southwest Research Institute 8500 Culebra Road San Antonio, Texas 78284 Dr. R. C. DeHart Southwest Research Institute 8500 Culebra Road San Antonio, Texas 78284 Dr. M. L. Baron Weidlinger Associates 110 East 59th Street New York, New York 10022 Dr. Kevin Thomas Westinghouse Electric Corp. Advanced Reactors Division P. O. Box 158 Madison, Pennsylvania 15663 Dr. H. D. Hibbitt Hibbitt & Karlsson, Inc. 132 George M. Cohan Boulevard Providence, Rhode Island 02903 Dr. R. D. Mindlin 89 Deer Hill Drive Ridgefield, Connecticut 06877 474:NP:7l6:lab 78u474-6l9 Industry and Research Institutes (Conrt) Dr. Richard E. Dame Mega Engineering 11961 Tech Road Silver Spring, Maryland 20904 Mr. G. M. Stanley Lockheed Palo Alto Research Laboratory 3251 Hanover Street Palo Alto, California 94304 Mr. R. L. Cloud Robert L. Cloud Associates, Inc. 2972 Adeline Street Berkeley, California 94703