# Perturbation Methods for the Eigencharacteristics of Symmetric and Asymmetric Systems

Authors: Philip D. Cha; Austin Shin
Journal: Shock and Vibration
 (2018)

---

## Abstract

Dynamic analysis for a vibratory system typically begins with an evaluation of its eigencharacteristics. However, when design changes are introduced, the eigensolutions of the system change and thus must be recomputed. In this paper, three different methods based on the eigenvalue perturbation theory are introduced to analyze the effects of modifications without performing a potentially time-consuming and costly reanalysis. They will be referred to as the straightforward perturbation method, the incremental perturbation method, and the triple product method. In the straightforward perturbation method, the eigenvalue perturbation theory is used to formulate a first-order and a second-order approximation of the eigensolutions of symmetric and asymmetric systems. In the incremental perturbation method, the straightforward approach is extended to analyze systems with large perturbations using an iterative scheme. Finally, in the triple product method, the accuracy of the approximate eigenvalues is significantly improved by exploiting the orthogonality conditions of the perturbed eigenvectors. All three methods require only the eigensolutions of the nominal or unperturbed system, and in application, they involve simple matrix multiplications. Numerical experiments show that the proposed methods achieve accurate results for systems with and without damping and for systems with symmetric and asymmetric system matrices.


---

## Body

## 1. Introduction

In order to analyze the vibration of any structural system, one usually first calculates the system’s eigensolutions. After the initial analysis has been performed, if any design modifications are then subsequently introduced, one would need to recalculate the perturbed system’s eigensolutions. However, when the degrees-of-freedom for the system is large, solving the eigenvalue problem multiple times for each set of perturbations can be computationally expensive and undesirable. In these cases, eigenvalue perturbation theory can be applied to obtain the approximate eigencharacteristics of the system without resolving an entirely new eigenvalue problem. The expressions for the approximate eigensolutions of the perturbed system consist solely of matrix multiplications, which can be efficiently performed by any computer.Over the years, there has been a large amount of interest in applying the perturbation theory and sensitivity analysis to study the effects of design changes on the eigensolutions of structures. Because of the large number of studies that have been published in these areas, only selected references will be addressed in this paper. Fox and Kapoor [1] derived expressions for the first derivatives of eigenvalues and eigenvectors with respect to a system’s parameters. The expressions were applicable only to self-adjoint systems, and the numerical examples used only contained real eigenvalues and were restricted to small perturbations in the system. Romstad et al. [2] formulated a power-series perturbation approach and applied the method to structural systems in order to find perturbed distinct or repeated eigenvalues and eigenvalue sensitivities. However, the systems they analyzed were symmetric and undamped. Plaut and Huseyin [3] derived expressions of the derivatives of eigenvalues and eigenvectors of non-self-adjoint systems. Rudisill [4] further explored sensitivity analysis by deriving simpler expressions of a system’s eigenvalues and eigenvectors using only one left-hand and right-hand eigenvector. These expressions could also be extended to any order of derivative of the eigenvalue and eigenvector, though it was not explicitly done so in the paper. Meirovitch [5, 6] obtained the complete response of gyroscopic systems using the perturbation method, but his method required an extra step of reducing the generalized eigenvalue problem to its standard form in order to take advantage of the fact that the gyroscopic matrix is skew symmetric. Rudisill and Chu [7] built upon the work of Rudisill and introduced two numerical methods of computing the derivatives of a system’s eigenvalues and eigenvectors. The first method consisted of an iterative method that could be used to compute the first partial derivatives of the eigenvalues and eigenvectors of a self-adjoint system and only the first partial derivative of the largest eigenvalue and its associated eigenvector of a non-self-adjoint system. The second method was an algebraic method that could be used to compute all orders of the eigensolutions’ derivatives without using the left-hand eigenvectors. However, unlike in Rudisill’s work, this method required the solution of n+1 equations as opposed to simply n equations. Nelson [8] presented a simplified procedure of calculating eigenvector derivatives of an eigensystem of any order using only the left and right eigenvectors and their associated eigenvalues. However, these expressions were derived based on only the standard eigenvalue problem, and the original eigensystem’s matrix of rank n−1 must be converted to rank n. Chen and Wada [9] formulated a first-order perturbation solution using a power-series expansion to conduct structural dynamic analysis. However, their analysis was restricted to systems with symmetric matrices. Meirovitch and Ryland [10] examined gyroscopic systems with small internal damping and utilized a second-order approach along with the Cholesky decomposition and linear transformations to reduce the generalized eigenvalue problem to a standard one. Because the damping considered was small, they used a perturbation approach to obtain the response of the system by considering the undamped gyroscopic system as the unperturbed system and treated the damping as perturbation. Belle [11] presented expressions of higher order sensitivities, allowing for more accurate approximations of a perturbed system’s eigensolutions than when using only the first partial derivatives. However, the paper did not quantify the improvement in accuracy with these higher-order sensitivities. Meirovitch and Ryland [12] extended their aforementioned perturbation approach to include external damping and used Rayleigh’s quotient to enhance the accuracy of the approximate eigenvalues used in the modal analysis. Baldwin and Hutton [13] surveyed various structural dynamics modification techniques, and in their analysis of perturbation techniques, they described a first-order approach restricted to systems without any damping. Chung and Lee [14] developed a second-order perturbation solution and analyzed heavily but weakly proportionally damped systems. Bickford [15] proposed a scheme that improved the accuracy of the approximate eigenvalues obtained using the usual first-order perturbation approach. It led to results whose accuracy approached those using the second-order perturbation approach. The method, however, was strictly valid for symmetric eigenvalue problems. Cronin [16] used the perturbation approach to analyze nonclassically damped dynamic systems. The proposed method involved a partial diagonalization of the homogeneous equations of motion. A perturbation quantity based on the off-diagonal terms of the partially diagonalized damping matrix was then defined. The eigenvalues and eigenvectors for the damped system were described in terms of power-series in the perturbation quantity. Wang [17] introduced two approximate methods for calculating a eigenvector’s derivatives when only a truncated set of the system’s mode shapes is available. The first was the explicit method, which approximated the contribution of higher modes. The second method was the implicit method, which assumed that the eigenvector derivatives were spanned by the truncated mode shapes and a residual static mode. Eldred et al. [18] developed an improved mass normalization method for expressing eigenvalue derivatives of structural systems. This improved method performed accurately only for slightly less restrictive parameter changes (up to 10%). In addition, Eldred et al. [19] explored the perturbation approach and developed a solution of the full nonlinear eigensolution perturbation equations that took into account all orders of the perturbation expansion, whereby the frequency and mode shape perturbation equations were coupled. The method presented in this paper will show that these coupled equations are not necessary in order to achieve accurate results with large perturbations. Kwak [20] proposed a power-series perturbation approach to express the perturbed eigensolutions of lightly damped systems. However, the analysis depended on the original system being undamped. Chen et al. [21] extended the perturbation approach to not only deal with distinct or repeated eigenvalues but also close eigenvalues or eigenvalue clusters. Their analysis, however, was constrained to deal with relatively small perturbations, and the system they analyzed contained only symmetric system matrices and no damping. Wang and Kirkhope [22] used a perturbation approach to analyze damped rotor systems. The solution scheme they developed, however, was restricted to gyroscopic systems with no cross-coupling other than gyroscopic effects. Tang and Wang [23] applied a power-series perturbation approach to analyze systems with distinct, repeated, and nearly equal frequencies. Their technique utilized complex operations such as subspace condensation and orthogonal decompositions, and they only considered undamped systems with small perturbations. Liu and Chan [24] developed an efficient matrix perturbation technique to approximate the changes in the eigensolutions of the system due to structural modifications. The solution scheme required one to perform a subspace condensation and an orthogonal decomposition in order to obtain the lower-order perturbed eigensolutions. The matrix singular value decomposition was then subsequently used to compute the higher-order perturbations of the eigensolutions. Trisovic [25] applied sensitivity analysis to two structural systems, a shaft of an electromotor and a cantilever beam, utilizing only the first partial derivatives of the eigenvalues and eigenvectors to alter the systems’ natural frequencies. Cha and Solberg [26] applied the eigenvalue perturbation theory to analyze various problems in structural dynamics. The method they used, however, was restricted to a first-order analysis and only valid for symmetric systems. ElBeheiry [27] analyzed weakly, nonproportionally damped systems, where this damping stemmed from gyroscopic moments or external damping factors. He utilized a second-order perturbation approach and represented the original system as a conservative gyroscopic system. However, the perturbations in the system were restricted to be one order of magnitude smaller than the values in the original system matrices.In this paper, simple methods are introduced that can be used to obtain the approximate eigenvalues and eigenvectors of various systems, including those with asymmetric mass, stiffness, or damping matrices and those with varying magnitudes of perturbations. A simple first-order perturbation analysis is sufficient if the modifications introduced into the system are small. By including perturbations up to the second-order, the accuracy achieved with the straightforward method becomes better than that obtained using a first-order perturbation approach. When the design changes are sufficiently large, the straightforward perturbation method, even with a second-order analysis, fails to approximate the structural vibration eigencharacteristics. In this case, the incremental perturbation method can be employed, which returns accurate eigensolutions because the system is minimally perturbed from one iteration to the next. Finally, a simple triple product method is proposed that further enhances the accuracy of the approximate eigenvalues of the modified systems. Numerical examples are presented to demonstrate the accuracy of the proposed perturbation methods for undamped and damped systems, and for structures with symmetric and asymmetric system matrices.


## 2. Theory

The use of eigenvalue perturbation theory in the design of mechanical systems has been well explored. However, most papers have only applied this method to systems with small damping parameters or to systems with symmetric mass and stiffness matrices. In the following, the first-order and second-order approximations of the perturbed eigencharacteristics are derived based on the perturbation theory. This formulation can be applied to systems with asymmetric matrices, such as gyroscopic systems, and to systems with large damping. An iterative scheme based on perturbing the system matrices incrementally is also proposed, which allows one to obtain accurate approximations of the eigensolutions of the perturbed system when large modifications are introduced. Finally, knowing the perturbed eigenvectors, a third approach is developed that yields approximate eigenvalues that are more accurate than the previous methods.


### 2.1. Straightforward Perturbation Method

Consider a system whose modes of vibration are governed by the solutions of the following generalized eigenvalue problem:(1)B0x0=λ0A0x0,where B0 and A0 denote the nominal system matrices, both of size p×p, λ0 represents a nominal eigenvalue, and x0 is the corresponding nominal right eigenvector. Assume all the unperturbed eigenvalues of the system are distinct, and B0 and A0 are asymmetric. The initial system’s left eigenvectors, denoted by y0, satisfy(2)B0Ty0=λ0A0Ty0.Assume thejth left eigenvector and the ith right eigenvector have been properly normalized such that they satisfy the following orthogonality conditions:(3)y0jTA0x0i=δij,y0jTB0x0i=λ0iδij,where δij denotes the Kronecker delta and i,j=1,…,p.After an initial analysis has been performed, design changes are introduced so that the system matrices become(4)A=A0+δA,B=B0+δB,where A0 and B0 denote the unperturbed matrices, δA and δB the perturbation matrices, and A and B the perturbed matrices. Moreover, elements of δA and δB are assumed to be an order of magnitude smaller than those of A0 and B0. Because the system parameters have been perturbed, the eigensolutions of the modified system will undergo changes as follows:(5)λj=λ0j+δλj+δ2λj,xj=x0j+δxj+δ2xj,yj=y0j+δyj+δ2yj,where δλj, δxj, and δyj denote the first-order eigensolution perturbations, and δ2λj, δ2xj, and δ2yj represent the second-order eigensolution perturbations. Resolving the eigenvalue problem associated with the perturbed system can be time-consuming and computationally taxing. Instead, the perturbation theory is used to find the approximate eigencharacteristics of the modified system. For completeness, the eigenvalue perturbation theory will be presented in detail.Consider thejth right and left eigenvectors of the perturbed system that satisfy(6)Bxj=λjAxj,yjTB=λjyjTA.Substituting Equations (4) and (5) in Equation (6), expanding, and keeping only the first-order and second-order terms, the following equations are obtained:(7)B0x0j+B0δxj+B0δ2xj+δBx0j+δBδxj=λ0jA0x0j+λ0jA0δxj+λ0jA0δ2xj+λ0jδAx0j+λ0jδAδxj+δλjA0x0j+δλjA0δxj+δλjδAx0j+δ2λjA0x0j,(8)y0jTB0+δyjTB0+δ2yjTB0+y0jTδB+δyjTδB=λ0jy0jTA0+λ0jδyjTA0+λ0jδ2yjTA0+λ0jy0jTδA+λ0jδyjTδA+δλjy0jTA0+δλjδyjTA0+δλjy0jTδA+δ2λjy0jTA0.Equating the first-order terms in Equations (7) and (8) gives(9)B0δxj+δBx0j=λ0jA0δxj+λ0jδAx0j+δλjA0x0j,(10)y0jTδB+δyjTB0=λ0jy0jTδA+λ0jδyjTA0+δλjy0jTA0.Similarly, equating the second-order terms in Equations (7) and (8) yields(11)B0δ2xj+δBδxj=λ0jA0δ2xj+λ0jδAδxj+δλjA0δxj+δλjδAx0j+δ2λjA0x0j,(12)δyjTδB+δ2yjTB0=λ0jδyjTδA+λ0jδ2yjTA0+δλjy0jTδA+δλjδyjTA0+δ2λjy0jTA0.Because the unperturbed eigenvalues are distinct, the unperturbed eigenvectors are linearly independent and form a basis. Thus,δxj and δ2xj can be expressed as a linear combination of x0r as follows:(13)δxj=∑r=1pεrjx0r,(14)δ2xj=∑r=1pηrjx0r,where εrj and ηrj are small coefficients that are to be determined.Consider Equation (9), which consists of the first-order terms. Substituting Equation (13) in Equation (9) gives(15)B0∑r=1pεrjx0r+δBx0j=λ0jA0∑r=1pεrjx0r+λ0jδAx0j+δλjA0x0j.Premultiplying Equation (15) by y0jT and noting the orthogonality conditions of Equation (3), one finds an expression for δλj, the first-order eigenvalue perturbation, as follows:(16)δλj=y0jTδB−λ0jδAx0j.Premultiplying Equation (9) by y0iT, where i≠j, one can readily solve for εij:(17)εij=y0iTδB−λ0jδAx0jλ0j−λ0i.Consider now Equation (11), which consists of the second-order terms. Substituting Equation (14) in Equation (11) gives(18)B0∑r=1pηrjx0r+δBδxj=λ0jA0∑r=1pηrjx0r+λ0jδAδxj+δλjA0δxj+δλjδAx0j+δ2λjA0x0j.Premultiplying Equation (18) by y0jT and exploiting the orthogonality conditions, one obtains the following expression for δ2λj, the second-order eigenvalue perturbation:(19)δ2λj=y0jTδB−λ0jδAδxj−δλjy0jTA0δxj+δAx0j.Premultiplying Equation (18) by y0iT, for i≠j, one obtains(20)ηij=y0iTδB−λ0jδAδxj−δλjy0iTA0δxj+δAx0jλ0j−λ0i.To determineεjj and ηjj, consider(21)yjTAxj=1.Substituting Equations (4) and (5) in Equation (21), expanding the resulting expression, and keeping only terms up to the second-order, one has(22)y0jTA0δxj+y0jTA0δ2xj+y0jTδAx0j+y0jTδAδxj+δyjTA0x0j+δyjTA0δxj+δyjTδAx0j+δ2yjTA0x0j=0.Equating the first-order and second-order terms in Equation (22) gives(23)y0jTA0δxj+y0jTδAx0j+δyjTA0x0j=0,(24)y0jTA0δ2xj+y0jTδAδxj+δyjTA0δxj+δyjTδAx0j+δ2yjTA0x0j=0.Expressingδyj and δ2yj as a linear combination of y0r, one obtains(25)δyj=∑r=1pε^rjy0r,(26)δ2yj=∑r=1pη^rjy0r.Substituting Equations (13) and (25) in Equation (23) gives(27)y0jTA0∑r=1pεrjx0r+y0jTδAx0j+∑r=1pε^rjy0rTA0x0j=0,and because of the orthogonality conditions, Equation (27) reduces to(28)εjj+y0jTδAx0j+ε^jj=0.Assumingεjj=ε^jj, then(29)εjj=ε^jj=−12y0jTδAx0j.Substituting Equations (14) and (26) in Equation (24) leads to(30)y0jTA0∑r=1pηrjx0r+y0jTδAδxj+δyjTA0δxj+δyjTδAx0j+∑r=1pη^rjy0rTA0x0j=0,which simplifies to(31)ηjj+y0jTδAδxj+δyjTA0δxj+δyjTδAx0j+η^jj=0,due to the orthogonality conditions. Assumingηjj=η^jj, then(32)ηjj=η^jj=−12y0jTδAδxj+δyjTA0δxj+δyjTδAx0j.In summary, the perturbed eigenvalues are given by(33)λj=λ0j+δλj+δ2λj,δλj=y0jTδB−λ0jδAx0j,δ2λj=y0jTδB−λ0jδAδxj−δλjy0jTA0δxj+δAx0j.The perturbed right eigenvectors are(34)xj=x0j+δxj+δ2xj,δxj=∑r=1pεrjx0r,δ2xj=∑r=1pηrjx0r,εij=y0iTδB−λ0jδAx0jλ0j−λ0ifori≠j,εjj=−12y0jTδAx0j,ηij=y0iTδB−λ0jδAδxj−δλjy0iTA0δxj+δAx0jλ0j−λ0ifori≠j,ηjj=−12y0jTδAδxj+δyjTA0δxj+δyjTδAx0j,and the perturbed left eigenvectors are(35)yj=y0j+δyj+δ2yj,δyj=∑r=1pε^rjy0r,δ2yj=∑r=1pη^rjy0r,ε^ij=y0jTδB−λ0jδAx0iλ0j−λ0ifori≠j,ε^jj=εjj,η^ij=y0jTδB−λ0jδAδxi−δλjy0jTA0δxi+δAx0iλ0j−λ0ifori≠j,η^jj=ηjj.Equations (33)–(35) can be used to obtain the approximate first-order or the second-order eigencharacteristics of an asymmetric system whose system parameters have undergone design changes or modifications. For a first-order analysis, the perturbed eigensolutions are given by(36)λj=λ0j+δλj,xj=x0j+δxj,yj=y0j+δyj,and for a second-order analysis, the perturbed eigensolutions are given by(37)λj=λ0j+δλj+δ2λj,xj=x0j+δxj+δ2xj,yj=y0j+δyj+δ2yj.While these equations appear to be complicated at first glance, they are easy to code and efficient to run because they only involve simple matrix multiplications. Finally, for a symmetric system whereA0=A0T and B0=B0T, the left and right eigenvectors are identical, i.e., yj=xj. Thus, in applying the perturbation equations, only the perturbed right eigenvectors need to be solved.


### 2.2. Incremental Perturbation Method

While the perturbation method is applicable only when elements ofδA and δB are assumed to be an order of magnitude smaller than those of A0 and B0, the method can be easily extended to analyze systems that are undergoing large changes. Specifically, assume the modifications introduced in the perturbation matrices δA and δB are large. A simple iterative procedure is proposed whereby the perturbations are accounted for incrementally. The proposed iteration scheme is applicable to either a first-order or a second-order analysis. The proposed iterative procedure is as follows:(1)
Solve for the eigensolutions,λ0j, x0j, and y0j, of the nominal system whose matrices are given by A0 and B0.(2)
Divide the large perturbation matrices intoN smaller increments δA/N and δB/N, where N is chosen such that δAi,j/N and δBi,j/N are less than 10% of A0i,j and B0i,j, respectively.(3)
Obtain the perturbed eigensolutions of the modified system whose matrices are given byA′=A0+δA/N and B′=B0+δB/N using Equation (36) or Equation (37), depending on if a first-order or a second-order analysis is desired. Denote the perturbed eigenvalues as λj′ and the perturbed right and left eigenvectors as xj′ and yj′, respectively.(4)
Assemble the perturbed right and left modal matricesXpert and Ypert, whose columns correspond to xj′ and yj′, for j=1,…,p.(5)
Normalize the perturbed modal matricesXpert and Ypert such that the diagonal elements of YpertTA′Xpert equal 1.(6)
For the next iteration, let the nominal system matrices beA0=A′ and B0=B′, and let the nominal eigensolutions be λ0j=λj′, x0j=xj′, and y0j=yj′.(7)
Repeat Steps 3 to 6 untilA′=A and B′=B.Using the proposed scheme, the error between the exact and perturbed eigensolutions at each iteration becomes negligible ifN is sufficiently large. The perturbed eigensolutions are obtained iteratively until the smaller perturbations have accumulated to equal to the total modifications. It is also important to note that after each iteration, the perturbed eigenvectors must be properly normalized such that the diagonal elements of YpertTA′Xpert are equal to 1. This step is necessary because the perturbed eigenvectors must satisfy Equation (3).Of great interest is also the accuracy of the perturbed eigenvectors, which can be determined by checking their self-compatibility using the orthogonality conditions. For the generalized eigenvalue problems, the following orthogonality check may be used:(38)I′=YpertTAXpert,where I′ denotes an orthogonal matrix, A is the system matrix, and Xpert and Ypert are normalized such that the diagonal elements of I′ are identical ones. For an ith-order analysis, the column vectors of the perturbed modal matrices Xpert and Ypert are the ith-order perturbed eigenvectors. Theoretically, if the perturbed modal matrices are exact, then the orthogonal matrix corresponds to the identity matrix. Since the perturbed modal matrices are approximate, the magnitudes of the nonzero off-diagonal terms of the orthogonal matrix can be used to infer the accuracy of the perturbed eigenvectors. As such, an error parameter for the modal matrices is defined, which is given by the average of the magnitudes of the off-diagonal terms, as follows:(39)δ=1pp−1∑i=1p∑j=1,j≠ipI′i,j,where p represents the size of the system matrices and I′i,j denotes the absolute value of the i,jth element of the orthogonal matrix I′. The smaller the value of δ is, the more accurate the perturbed modal matrices are.


### 2.3. Triple Product Method

After the completion of the straightforward perturbation method or the incremental perturbation method, the following triple product is calculated:(40)YpertTBXpert,assuming these perturbed modal matrices have been properly normalized so that the diagonal elements of Equation (38) are all ones. If the perturbed modal matrices were replaced with the exact modal matrices, then the diagonal elements of Equation (40) correspond to the exact eigenvalues and all the off-diagonal terms will be identically zero. Since the perturbation analysis is only an approximation, the diagonal terms from the triple product can also be used to approximate the eigenvalues. In evaluating Equation (40), the modal matrices can be obtained via either a first-order or second-order analysis. Interestingly, numerical experiments show that using an ith-order perturbation analysis, the diagonal elements of Equation (40) are significantly more accurate than those obtained using the ith-order straightforward perturbation method. This drastic improvement in accuracy is achieved by performing one additional step consisting of simple matrix multiplication. As such, in all of the examples, the diagonal elements of the triple product of Equation (40) will also be presented, both to validate the accuracy of the perturbed eigenvectors and to offer yet another approximation for the eigenvalues of the modified system.


## 2.1. Straightforward Perturbation Method

Consider a system whose modes of vibration are governed by the solutions of the following generalized eigenvalue problem:(1)B0x0=λ0A0x0,where B0 and A0 denote the nominal system matrices, both of size p×p, λ0 represents a nominal eigenvalue, and x0 is the corresponding nominal right eigenvector. Assume all the unperturbed eigenvalues of the system are distinct, and B0 and A0 are asymmetric. The initial system’s left eigenvectors, denoted by y0, satisfy(2)B0Ty0=λ0A0Ty0.Assume thejth left eigenvector and the ith right eigenvector have been properly normalized such that they satisfy the following orthogonality conditions:(3)y0jTA0x0i=δij,y0jTB0x0i=λ0iδij,where δij denotes the Kronecker delta and i,j=1,…,p.After an initial analysis has been performed, design changes are introduced so that the system matrices become(4)A=A0+δA,B=B0+δB,where A0 and B0 denote the unperturbed matrices, δA and δB the perturbation matrices, and A and B the perturbed matrices. Moreover, elements of δA and δB are assumed to be an order of magnitude smaller than those of A0 and B0. Because the system parameters have been perturbed, the eigensolutions of the modified system will undergo changes as follows:(5)λj=λ0j+δλj+δ2λj,xj=x0j+δxj+δ2xj,yj=y0j+δyj+δ2yj,where δλj, δxj, and δyj denote the first-order eigensolution perturbations, and δ2λj, δ2xj, and δ2yj represent the second-order eigensolution perturbations. Resolving the eigenvalue problem associated with the perturbed system can be time-consuming and computationally taxing. Instead, the perturbation theory is used to find the approximate eigencharacteristics of the modified system. For completeness, the eigenvalue perturbation theory will be presented in detail.Consider thejth right and left eigenvectors of the perturbed system that satisfy(6)Bxj=λjAxj,yjTB=λjyjTA.Substituting Equations (4) and (5) in Equation (6), expanding, and keeping only the first-order and second-order terms, the following equations are obtained:(7)B0x0j+B0δxj+B0δ2xj+δBx0j+δBδxj=λ0jA0x0j+λ0jA0δxj+λ0jA0δ2xj+λ0jδAx0j+λ0jδAδxj+δλjA0x0j+δλjA0δxj+δλjδAx0j+δ2λjA0x0j,(8)y0jTB0+δyjTB0+δ2yjTB0+y0jTδB+δyjTδB=λ0jy0jTA0+λ0jδyjTA0+λ0jδ2yjTA0+λ0jy0jTδA+λ0jδyjTδA+δλjy0jTA0+δλjδyjTA0+δλjy0jTδA+δ2λjy0jTA0.Equating the first-order terms in Equations (7) and (8) gives(9)B0δxj+δBx0j=λ0jA0δxj+λ0jδAx0j+δλjA0x0j,(10)y0jTδB+δyjTB0=λ0jy0jTδA+λ0jδyjTA0+δλjy0jTA0.Similarly, equating the second-order terms in Equations (7) and (8) yields(11)B0δ2xj+δBδxj=λ0jA0δ2xj+λ0jδAδxj+δλjA0δxj+δλjδAx0j+δ2λjA0x0j,(12)δyjTδB+δ2yjTB0=λ0jδyjTδA+λ0jδ2yjTA0+δλjy0jTδA+δλjδyjTA0+δ2λjy0jTA0.Because the unperturbed eigenvalues are distinct, the unperturbed eigenvectors are linearly independent and form a basis. Thus,δxj and δ2xj can be expressed as a linear combination of x0r as follows:(13)δxj=∑r=1pεrjx0r,(14)δ2xj=∑r=1pηrjx0r,where εrj and ηrj are small coefficients that are to be determined.Consider Equation (9), which consists of the first-order terms. Substituting Equation (13) in Equation (9) gives(15)B0∑r=1pεrjx0r+δBx0j=λ0jA0∑r=1pεrjx0r+λ0jδAx0j+δλjA0x0j.Premultiplying Equation (15) by y0jT and noting the orthogonality conditions of Equation (3), one finds an expression for δλj, the first-order eigenvalue perturbation, as follows:(16)δλj=y0jTδB−λ0jδAx0j.Premultiplying Equation (9) by y0iT, where i≠j, one can readily solve for εij:(17)εij=y0iTδB−λ0jδAx0jλ0j−λ0i.Consider now Equation (11), which consists of the second-order terms. Substituting Equation (14) in Equation (11) gives(18)B0∑r=1pηrjx0r+δBδxj=λ0jA0∑r=1pηrjx0r+λ0jδAδxj+δλjA0δxj+δλjδAx0j+δ2λjA0x0j.Premultiplying Equation (18) by y0jT and exploiting the orthogonality conditions, one obtains the following expression for δ2λj, the second-order eigenvalue perturbation:(19)δ2λj=y0jTδB−λ0jδAδxj−δλjy0jTA0δxj+δAx0j.Premultiplying Equation (18) by y0iT, for i≠j, one obtains(20)ηij=y0iTδB−λ0jδAδxj−δλjy0iTA0δxj+δAx0jλ0j−λ0i.To determineεjj and ηjj, consider(21)yjTAxj=1.Substituting Equations (4) and (5) in Equation (21), expanding the resulting expression, and keeping only terms up to the second-order, one has(22)y0jTA0δxj+y0jTA0δ2xj+y0jTδAx0j+y0jTδAδxj+δyjTA0x0j+δyjTA0δxj+δyjTδAx0j+δ2yjTA0x0j=0.Equating the first-order and second-order terms in Equation (22) gives(23)y0jTA0δxj+y0jTδAx0j+δyjTA0x0j=0,(24)y0jTA0δ2xj+y0jTδAδxj+δyjTA0δxj+δyjTδAx0j+δ2yjTA0x0j=0.Expressingδyj and δ2yj as a linear combination of y0r, one obtains(25)δyj=∑r=1pε^rjy0r,(26)δ2yj=∑r=1pη^rjy0r.Substituting Equations (13) and (25) in Equation (23) gives(27)y0jTA0∑r=1pεrjx0r+y0jTδAx0j+∑r=1pε^rjy0rTA0x0j=0,and because of the orthogonality conditions, Equation (27) reduces to(28)εjj+y0jTδAx0j+ε^jj=0.Assumingεjj=ε^jj, then(29)εjj=ε^jj=−12y0jTδAx0j.Substituting Equations (14) and (26) in Equation (24) leads to(30)y0jTA0∑r=1pηrjx0r+y0jTδAδxj+δyjTA0δxj+δyjTδAx0j+∑r=1pη^rjy0rTA0x0j=0,which simplifies to(31)ηjj+y0jTδAδxj+δyjTA0δxj+δyjTδAx0j+η^jj=0,due to the orthogonality conditions. Assumingηjj=η^jj, then(32)ηjj=η^jj=−12y0jTδAδxj+δyjTA0δxj+δyjTδAx0j.In summary, the perturbed eigenvalues are given by(33)λj=λ0j+δλj+δ2λj,δλj=y0jTδB−λ0jδAx0j,δ2λj=y0jTδB−λ0jδAδxj−δλjy0jTA0δxj+δAx0j.The perturbed right eigenvectors are(34)xj=x0j+δxj+δ2xj,δxj=∑r=1pεrjx0r,δ2xj=∑r=1pηrjx0r,εij=y0iTδB−λ0jδAx0jλ0j−λ0ifori≠j,εjj=−12y0jTδAx0j,ηij=y0iTδB−λ0jδAδxj−δλjy0iTA0δxj+δAx0jλ0j−λ0ifori≠j,ηjj=−12y0jTδAδxj+δyjTA0δxj+δyjTδAx0j,and the perturbed left eigenvectors are(35)yj=y0j+δyj+δ2yj,δyj=∑r=1pε^rjy0r,δ2yj=∑r=1pη^rjy0r,ε^ij=y0jTδB−λ0jδAx0iλ0j−λ0ifori≠j,ε^jj=εjj,η^ij=y0jTδB−λ0jδAδxi−δλjy0jTA0δxi+δAx0iλ0j−λ0ifori≠j,η^jj=ηjj.Equations (33)–(35) can be used to obtain the approximate first-order or the second-order eigencharacteristics of an asymmetric system whose system parameters have undergone design changes or modifications. For a first-order analysis, the perturbed eigensolutions are given by(36)λj=λ0j+δλj,xj=x0j+δxj,yj=y0j+δyj,and for a second-order analysis, the perturbed eigensolutions are given by(37)λj=λ0j+δλj+δ2λj,xj=x0j+δxj+δ2xj,yj=y0j+δyj+δ2yj.While these equations appear to be complicated at first glance, they are easy to code and efficient to run because they only involve simple matrix multiplications. Finally, for a symmetric system whereA0=A0T and B0=B0T, the left and right eigenvectors are identical, i.e., yj=xj. Thus, in applying the perturbation equations, only the perturbed right eigenvectors need to be solved.


## 2.2. Incremental Perturbation Method

While the perturbation method is applicable only when elements ofδA and δB are assumed to be an order of magnitude smaller than those of A0 and B0, the method can be easily extended to analyze systems that are undergoing large changes. Specifically, assume the modifications introduced in the perturbation matrices δA and δB are large. A simple iterative procedure is proposed whereby the perturbations are accounted for incrementally. The proposed iteration scheme is applicable to either a first-order or a second-order analysis. The proposed iterative procedure is as follows:(1)
Solve for the eigensolutions,λ0j, x0j, and y0j, of the nominal system whose matrices are given by A0 and B0.(2)
Divide the large perturbation matrices intoN smaller increments δA/N and δB/N, where N is chosen such that δAi,j/N and δBi,j/N are less than 10% of A0i,j and B0i,j, respectively.(3)
Obtain the perturbed eigensolutions of the modified system whose matrices are given byA′=A0+δA/N and B′=B0+δB/N using Equation (36) or Equation (37), depending on if a first-order or a second-order analysis is desired. Denote the perturbed eigenvalues as λj′ and the perturbed right and left eigenvectors as xj′ and yj′, respectively.(4)
Assemble the perturbed right and left modal matricesXpert and Ypert, whose columns correspond to xj′ and yj′, for j=1,…,p.(5)
Normalize the perturbed modal matricesXpert and Ypert such that the diagonal elements of YpertTA′Xpert equal 1.(6)
For the next iteration, let the nominal system matrices beA0=A′ and B0=B′, and let the nominal eigensolutions be λ0j=λj′, x0j=xj′, and y0j=yj′.(7)
Repeat Steps 3 to 6 untilA′=A and B′=B.Using the proposed scheme, the error between the exact and perturbed eigensolutions at each iteration becomes negligible ifN is sufficiently large. The perturbed eigensolutions are obtained iteratively until the smaller perturbations have accumulated to equal to the total modifications. It is also important to note that after each iteration, the perturbed eigenvectors must be properly normalized such that the diagonal elements of YpertTA′Xpert are equal to 1. This step is necessary because the perturbed eigenvectors must satisfy Equation (3).Of great interest is also the accuracy of the perturbed eigenvectors, which can be determined by checking their self-compatibility using the orthogonality conditions. For the generalized eigenvalue problems, the following orthogonality check may be used:(38)I′=YpertTAXpert,where I′ denotes an orthogonal matrix, A is the system matrix, and Xpert and Ypert are normalized such that the diagonal elements of I′ are identical ones. For an ith-order analysis, the column vectors of the perturbed modal matrices Xpert and Ypert are the ith-order perturbed eigenvectors. Theoretically, if the perturbed modal matrices are exact, then the orthogonal matrix corresponds to the identity matrix. Since the perturbed modal matrices are approximate, the magnitudes of the nonzero off-diagonal terms of the orthogonal matrix can be used to infer the accuracy of the perturbed eigenvectors. As such, an error parameter for the modal matrices is defined, which is given by the average of the magnitudes of the off-diagonal terms, as follows:(39)δ=1pp−1∑i=1p∑j=1,j≠ipI′i,j,where p represents the size of the system matrices and I′i,j denotes the absolute value of the i,jth element of the orthogonal matrix I′. The smaller the value of δ is, the more accurate the perturbed modal matrices are.


## 2.3. Triple Product Method

After the completion of the straightforward perturbation method or the incremental perturbation method, the following triple product is calculated:(40)YpertTBXpert,assuming these perturbed modal matrices have been properly normalized so that the diagonal elements of Equation (38) are all ones. If the perturbed modal matrices were replaced with the exact modal matrices, then the diagonal elements of Equation (40) correspond to the exact eigenvalues and all the off-diagonal terms will be identically zero. Since the perturbation analysis is only an approximation, the diagonal terms from the triple product can also be used to approximate the eigenvalues. In evaluating Equation (40), the modal matrices can be obtained via either a first-order or second-order analysis. Interestingly, numerical experiments show that using an ith-order perturbation analysis, the diagonal elements of Equation (40) are significantly more accurate than those obtained using the ith-order straightforward perturbation method. This drastic improvement in accuracy is achieved by performing one additional step consisting of simple matrix multiplication. As such, in all of the examples, the diagonal elements of the triple product of Equation (40) will also be presented, both to validate the accuracy of the perturbed eigenvectors and to offer yet another approximation for the eigenvalues of the modified system.


## 3. Results

Three different methods are proposed to obtain the eigencharacteristics of modified symmetric or asymmetric systems. They consist of the straightforward perturbation method, the incremental perturbation method, and the triple product method. Each of these methods can be applied using either a first-order or second-order perturbation analysis. Using the straightforward perturbation method, the approximate eigensolutions are readily obtained by the direct application of Equations (33)–(35). Using the incremental perturbation method, Equations (33)–(35) are still applicable, except the perturbations are divided into N increments, and the approximate eigencharacteristics of the modified system are obtained iteratively. Lastly, using the triple product method, the accuracy of the approximate eigenvalues given by the diagonal elements of Equation (40) becomes dramatically improved at the cost of only one additional step consisting of simple matrix multiplication. In the following, the eigensolutions for various symmetric and asymmetric systems will be obtained using the methods presented in this paper. For brevity, only the governing equations for the unperturbed system will be presented in each example. The governing equations for the perturbed system are identical to those of the unperturbed system, except that the perturbed system matrices will include the structural modifications that are subsequently introduced.


### 3.1. Torque-Free Rigid Body

Consider a torque-free rigid body as shown in Figure1 [5], which is symmetric about axis ζ. Two masses, each m/2, are located at distances ζ=±a from the fixed center O. Each mass is supported by 4 springs, each with stiffness k/4. The rigid body spins about axis ζ with a constant angular velocity of Ωζ. The mass moments of inertia of the rigid body about axes ξ, η, and ζ are α, β=α, and c, respectively. The mass moment of inertia of the entire body about a transverse axis is given by α′=α+ma2. The system’s unperturbed matrices are given by [5](41)where α0′=α0+m0a02. For definiteness, the original and unperturbed system’s parameters are Ωζ = Ω0 = 5 × 101 rad/s, ω0 = 6 × 101 rad/s, c0/α0 = 1.5, a0 = 2 m, m0 = 2.5 × 10−2 kg, k0 = 9 × 101 N/m, and α0 = 1 kg·m2.Figure 1
A torque-free rigid body.After an initial analysis has been performed, design changes are introduced. Letx denote an arbitrary system parameter. Assume it is perturbed systematically from its initial value of x0 as follows:(42)x=x01+σδx,where δx is arbitrarily chosen from a set of random numbers with zero mean and a standard deviation of one, and σ represents a disorder strength that can be varied. Together, σδx denotes the percentage change in x0. Thus, a large σ value implies a large perturbation from the nominal value. All the modified system parameters are perturbed in this manner.For the torque-free rigid body, the natural frequenciesω are the imaginary part of the eigenvalues λ, i.e., ω=imagλ. Suppose the modified parameters m and a deviate from their nominal values such that m=m01+σδm and a=a01+σδa, where δm = −4.4070 × 10−1 and δa = 8.4650 × 10−1. Figure 2 shows the natural frequencies as a function of σ. These natural frequencies are obtained via five different means: by solving the problem exactly, by using the straightforward perturbation method via a first-order or a second-order analysis, and using the triple product method of Equation (40), where the modal matrices consist of the first-order or the second-order perturbed eigenvectors. When σ=0, the system is unperturbed. Note that for small values of σ, all the approximate natural frequencies track the exact perturbed natural frequencies accurately. As the magnitude of σ increases, while the first-order and second-order straightforward perturbation approximations deviate from the exact values, the approximations obtained from the diagonal elements of the triple product, using either the first-order or second-order perturbed modal matrices, are indistinguishable from the exact natural frequencies. The fact that the diagonal elements track the exact natural frequencies so closely also validates the accuracy of the perturbed matrices. Table 1 tabulates these natural frequencies for σ=0.65. Thus, m and a deviate from their nominal values by approximately σδm≈−28.6% and σδa≈55.0%, respectively. The numerical value in the parentheses denotes the percentage error in the natural frequencies, defined as(43)ε=exact−approximateexact×100%.Figure 2
The exact and approximate natural frequencies (rad/s) of the torque-free rigid body as a function ofσ.
(a)
(b)
(c)Table 1
The exact and approximate natural frequencies for the perturbed torque-free rigid body of Figure1 for σ=0.65.
Exact
First-order straightforward perturbation
Second-order straightforward perturbation
First-order triple product
Second-order triple product
ω
1
1.2592
×
10
−
1
4.9493
×
10
−
1 2.93×102
1.9843
×
10
−
1 5.75×101
1.2801
×
10
−
1 1.66×100
1.2503
×
10
−
1 7.03×10−1
ω
2
4.3014
×
10
1
4.0276
×
10
1 6.36×100
4.2147
×
10
1 2.01×100
4.2925
×
10
1 2.08×10−1
4.3004
×
10
1 2.43×10−2
ω
3
1.1236
×
10
2
1.1203
×
10
2 2.99×10−1
1.1226
×
10
2 9.32×10−2
1.1236
×
10
2 1.98×10−3
1.1236
×
10
2 5.69×10−4Note the drastic improvement in the approximate natural frequencies obtained using a first-order triple product approach compared with using a first-order straightforward expansion, especially for the fundamental natural frequency. This significant improvement in accuracy is the result of carrying out another step consisting of a simple matrix multiplication given by Equation (40). Interestingly, the approximate natural frequencies obtained using a first-order triple product approach are even more accurate that those obtained using a second-order straightforward perturbation method. Note that among the different methods, the error parameters using the triple product with the second-order modal matrices are the smallest.If the error parameterε is deemed too large for a disorder strength of σ=0.65 or if a larger value of σ is applied, which would represent larger perturbations in the system, dividing the perturbations into smaller pieces and computing the approximations through an incremental method can be utilized. Although this incremental procedure involves more computation, it is deemed worthwhile when dealing with large perturbations in order to obtain accurate approximations. Consider now the case where σ=1.8. Figure 3 shows the natural frequencies obtained using the five methods mentioned previously as a function of N. These approximate natural frequencies monotonically converge to their respective exact values though some methods converge much more quickly than others. Note that the approximate fundamental natural frequency converges to the exact ω1 from above, while the other approximate natural frequencies converge to their respective exact values from below. By N=5, the triple product method with the first-order and second-order modal matrices leads to approximations that are nearly exact while the first-order incremental perturbation method is still very inaccurate by N=10. Across all values of N, the results of the triple product approximations via a second-order analysis are the most accurate. Table 2 illustrates the numerical results of Figure 3 for N=8. With eight iterations, the perturbations at each iteration are less than 10% of the original system matrices. At N=8, the errors for all three natural frequencies are reduced by at least an order of magnitude when approximating the new natural frequencies using the diagonal elements of the second-order triple product as opposed to simply the incremental perturbation method via a first-order analysis. With the second-order triple product, the final error percentages of all three natural frequencies are below 0.83%, substantially smaller than the percentage deviations in m0 and a0 that are introduced into the system, which are approximately σδm≈−79.3% and σδa≈152.3%, respectively. To assess the accuracy of the perturbed modal matrices for this torque-free rigid body system, consider the error parameter given by Equation (39), where p=6. Figure 4 shows δ as a function of N for σ=1.8, obtained using a first-order or a second-order incremental perturbation analysis. Note that δ decreases with N and that the error parameter for the second-order modal matrices decreases more quickly than the error parameter for the first-order modal matrices.Figure 3
The approximate natural frequencies of the torque-free rigid body as a function ofN for σ=1.8.
(a)
(b)
(c)Table 2
The exact and approximate natural frequencies for the perturbed torque-free rigid body atN=8 for σ=1.8.
Exact
First-order incremental perturbation
Second-order incremental perturbation
First-order triple product
Second-order triple product
ω
1
1.3511
×
10
0
2.2065
×
10
0 6.33×101
1.6401
×
10
0 2.13×101
1.4190
×
10
0 5.02×100
1.3622
×
10
0 8.22×10−1
ω
2
3.7885
×
10
1
3.4847
×
10
1 8.01×100
3.6785
×
10
1 2.90×100
3.7427
×
10
1 1.20×100
3.7813
×
10
1 1.89×10−1
ω
3
1.1178
×
10
2
1.1145
×
10
2 2.91×10−1
1.1167
×
10
2 9.99×10−2
1.1177
×
10
2 4.38×10−3
1.1178
×
10
2 1.84×10−5Figure 4
The error parameter,δ, for the perturbed modal matrices of the torque-free rigid body as a function of N for σ=1.8.


### 3.2. Two-Degrees-of-Freedom Gyroscopic System

Figure5 shows a two-degrees-of-freedom gyroscopic system analyzed by Meirovitch and Ryland [10]. The governing equations for the unperturbed system are given by(44)M0q¨0+C0+G0q˙0+K0q0=0,where an overdot denotes a time derivative, [M0], [C0], and [K0] are the mass, damping, and stiffness matrices of the nominal system, respectively, [G0] denotes the unperturbed gyroscopic damping matrix, and q0 is the vector of generalized coordinates for the nominal system. The system matrices are all of size 2×2, and q0 is of length 2. The nominal system matrices are(45)M0=m000m0,G0=0−2m0Ω02m0Ω00,C0=c0000,K0=k01−m0Ω0200k02−m0Ω02.Figure 5
A two-degree-of-freedom gyroscopic system.To solve for the eigensolutions of the system shown in Equation (44), the state vector approach [28] is used. Introducing a state vector of length 4 of the form(46)x0=q˙0q0,allows Equation (44) to be rewritten as(47)A0x˙0−B0x0=0,where A0 and B0 are the nominal state matrices of size 4×4 given by(48)A0=M000K0,B0=−C0−G0−K0K00.Figure5 is an example of a damped, asymmetric system due to the presence of the gyroscopic damping matrix G0. For definiteness, the original or unperturbed system parameters are m0 = 1 kg, Ω0 = 1 rad/s, c0 = 1 × 10−1 N·s/m, k01 = 3 N/m, and k02 = 4 N/m. To illustrate that the perturbation methods presented in this paper can also be used when an asymmetric damping matrix is perturbed, the parameters c0 and m0 are subjected to modifications; thus c=c01+σδc and m=m01+σδm. The random numbers used for the perturbations are δc=−7.2540×10−1 and δm=2.8480×10−1. Figures 6 and 7 illustrate the real and imaginary parts of the system’s eigenvalues (obtained exactly, with the straightforward perturbation method via a first-order or a second-order analysis, and with the diagonal elements of Equation (40) using the first-order or second-order perturbed modal matrices) as a function of σ. Note the diagonal elements of both triple products almost lie on top of the exact eigenvalues, which validates the accuracy of the perturbed modal matrices, and they are more accurate than the second-order straightforward perturbation approximations over the entire range of σ. Tables 3 and 4 compare the real and imaginary parts of the exact and approximate natural frequencies for σ=0.8. Note again that the first-order triple product method yields approximate eigenvalues that are more accurate than those obtained using a second-order straightforward perturbation analysis. Moreover, note that using the second-order triple product method, the error parameters for both the real and imaginary parts of the eigenvalues are reduced by at least two orders of magnitude compared to those when using the first-order straightforward perturbation method. Additionally, for σ=0.8, the real parts of the eigenvalues obtained using the second-order triple product method have ε that are all less than 0.16%, and the imaginary parts have ε all less than 0.0095%. These errors are significantly smaller than the percentage changes in c0 and m0, which are approximately σδm≈22.7% and σδc≈−58.0%, respectively.Figure 6
The real parts of the exact and approximate eigenvalues of the two-degree-of-freedom gyroscopic system as a function ofσ.
(a)
(b)Figure 7
The imaginary parts of the exact and approximate eigenvalues of the two-degree-of-freedom gyroscopic system as a function ofσ.
(a)
(b)Table 3
The real parts of the exact and approximate natural frequencies for the perturbed two-degree-of-freedom gyroscopic system of Figure5 for σ=0.8.
Exact
First-order straightforward perturbation
Second-order straightforward perturbation
First-order triple product
Second-order triple product
Re
λ
1
−
4.5435
×
10
−
3
−
9.7877
×
10
−
4 7.84×101
−
5.4541
×
10
−
3 2.00×101
−
4.4447
×
10
−
3 2.17×100
−
4.5366
×
10
−
3 1.50×10−1
Re
λ
2
−
1.2546
×
10
−
2
−
8.6132
×
10
−
3 3.13×101
−
1.3344
×
10
−
2 6.35×100
−
1.2565
×
10
−
2 1.53×10−1
−
1.2540
×
10
−
2 4.60×10−2Table 4
The imaginary parts of the exact and approximate natural frequencies for the perturbed two-degree-of-freedom gyroscopic system of Figure5 for σ=0.8.
Exact
First-order straightforward perturbation
Second-order straightforward perturbation
First-order triple product
Second-order triple product
Im
λ
1
6.7034
×
10
−
1
6.4020
×
10
−
1 4.49×100
6.7603
×
10
−
1 8.49×10−1
6.7000
×
10
−
1 5.01×10−2
6.7027
×
10
−
1 9.43×10−3
Im
λ
2
2.6928
×
10
0
2.6623
×
10
0 1.13×100
2.6985
×
10
0 2.11×10−1
2.6927
×
10
0 2.90×10−3
2.6928
×
10
0 3.95×10−4Assume larger perturbations are introduced to the system by usingσ=1.37. Figures 8 and 9 show the eigenvalues found using the five methods mentioned previously as a function of N. Note that using a first-order incremental perturbation method, the real parts of the eigenvalue converge to their exact values from above, while the imaginary parts converge to their exact values from below. The trend is reversed using the second-order incremental perturbation analysis, i.e., the real and imaginary parts converge to their exact values from below and above, respectively. For the chosen set of structural parameters, the modified system is stable. However, Figure 8 shows that the first-order incremental perturbation method yields approximate eigenvalues that erroneously indicate the system is unstable. In the case, because the real components of the system’s exact eigenvalues are close to 0, it is important to use a sufficiently large value of N to accurately capture the system’s stability. Also, note that the approximate eigenvalues found using the diagonal elements of the first-order and second-order triple products converge rapidly to the exact eigenvalues for N as small as 3, and for N=10, the approximate and exact eigenvalues become indistinguishable, despite the large perturbations that are introduced in the system. Tables 5 and 6 show the corresponding numerical results. Note that ε of the second-order triple product is at least four orders of magnitude smaller than ε of the first-order incremental perturbation method. Figure 10 shows the error parameter δ for the first-order and second-order perturbed modal matrices, obtained using Equation (39) with p=4, as a function of N. It clearly shows a rapid decrease of δ, which implies a high accuracy in the perturbed eigenvectors, especially those obtained using a second-order analysis.Figure 8
The real parts of the approximate eigenvalues of the two-degree-of-freedom gyroscopic system as a function ofN for σ=1.37.
(a)
(b)Figure 9
The imaginary parts of the approximate eigenvalues of the two-degree-of-freedom gyroscopic system as a function ofN for σ=1.37.
(a)
(b)Table 5
The real parts of the exact and approximate natural frequencies for the perturbed two-degree-of-freedom gyroscopic system atN=10 for σ=1.37.
Exact
First-order incremental perturbation
Second-order incremental perturbation
First-order triple product
Second-order triple product
Re
λ
1
−
5.4248
×
10
−
5
5.7165
×
10
−
4 1.15×103
−
7.0836
×
10
−
5 3.05×101
−
5.0406
×
10
−
5 7.08×100
−
5.4177
×
10
−
5 1.30×10−1
Re
λ
2
−
1.6881
×
10
−
4
6.4366
×
10
−
4 4.81×102
−
1.8005
×
10
−
4 6.65×100
−
1.6942
×
10
−
4 3.62×10−1
−
1.6879
×
10
−
4 1.16×10−2Table 6
The imaginary parts of the exact and approximate natural frequencies for the perturbed two-degree-of-freedom gyroscopic system atN=10 for σ=1.37.
Exact
First-order incremental perturbation
Second-order incremental perturbation
First-order triple product
Second-order triple product
Im
λ
1
5.6913
×
10
−
1
5.6296
×
10
−
1 1.08×100
5.6989
×
10
−
1 3.63×10−2
5.6912
×
10
−
1 2.63×10−3
5.6913
×
10
−
1 4.71×10−5
Im
λ
2
2.5906
×
10
0
2.5843
×
10
0 2.43×10−1
2.5907
×
10
0 6.06×10−3
2.5906
×
10
0 1.08×10−4
2.5906
×
10
0 1.59×10−6Figure 10
The error parameter,δ, for the perturbed modal matrices of the two-degree-of-freedom gyroscopic system as a function of N for σ=1.37.


### 3.3. Whirling Beam with Gyroscopic Damping

Figure11 shows a whirling shaft with a lumped mass M attached at its midspan [12]. The shaft’s mass is uniformly distributed with mass per unit length m. The shaft has stiffnesses EIx and EIy about the x and y axis, respectively. The end supports provide bending moments proportional to the angular deflections, where the proportionality constants are K1 at Z=0 and K2 at Z=L. It is assumed that there is uniformly distributed external dissipation proportional by a constant h to the absolute velocity, and an internal dissipation proportional by a constant c to the velocity relative to the rotating axes. The equations of motion of the unperturbed system are given by [12](49)M0q¨0+C0+G0q˙0+K0+H0q0=0,where the coefficient matrices in Equation (49) are partitioned into s×s submatrices(50)M0=Ms00Ms,G0=0Gs−Gs0,C0=Cs00Cs,K0=Ks100Ks2,H0=0Hs−Hs0.Figure 11
A whirling beam.Thei,jth elements of the submatrices are given by [12](51)Mijs=m0L0δij+2M0siniπ2sinjπ2,Gijs=−2Ω0Mij,Cijs=c0+h0L0δij,Kijs1=2K01+K02cosiπcosjπiπL0jπL0+EI0xiπL02jπL02L0δij−Ω02Mij,Kijs2=2K01+K02cosiπcosjπiπL0jπL0+EI0yiπL02jπL02L0δij−Ω02Mij,Hijs=−h0Ω0L0δij,where i,j=1,…,s. For this example, if s=3, then all of the submatrices will be of size 3×3. For definiteness, the original or unperturbed system parameters are m0 = 1 kg/m, M0 = 10 kg, L0 = 1 m, EI0x = 8.10 × 10−2 N·m2, EI0y = 1.82 × 10−1 N·m2, K01 = K02 = 5 × 10−2 N/m, c0 = h0 = 2.5 × 10−1 N·s/m, Ω0 = 1.46 × 101 rad/s. The eigensolutions of the unperturbed system are obtained using the state vector approach, where the system matrices are given by(52)A0=0M0M0C0+G0,B0=M000−K0−H0.Assume now the parameters affecting the damping matricesC0 and G0 are perturbed so that the modified parameters are given by M=M01+σδM, m=m01+σδm, c=c01+σδc and h=h01+σδh, where δM=−1.1079×100, δm=2.5260×10−1, δc=2.4600×10−2, and δh=−1.1270×100. For the given set of system parameters, λ1 and λ2 are strictly real while the remaining eigenvalues are complex. Additionally, the real components of the first and fourth eigenvalues are positive, indicating that the system is unstable. Figures 12 and 13 illustrate the real and imaginary parts of the eigenvalues as a function of σ, obtained using the methods developed in this paper. The diagonal elements of the second-order triple product again track the exact eigenvalues most closely for all values of σ. Tables 7 and 8 tabulate the numerical results for σ=0.2. For the given set of random variations, the parameters M, m, c, and h deviate from their nominal values by approximately −22.1%, 5.05%, 0.49%, and −22.5%, respectively. Note that the percentage errors for the diagonal elements of Equation (40) with the second-order perturbed modal matrices are all under 0.004%, despite the large deviations in some of the selected system parameters.Figure 12
The real parts of the exact and approximate eigenvalues of the whirling beam as a function ofσ.
(a)
(b)
(c)
(d)
(e)Figure 13
The imaginary parts of the exact and approximate eigenvalues of the whirling beam as a function ofσ.
(a)
(b)
(c)Table 7
The real parts of the exact and approximate natural frequencies for the perturbed whirling beam of Figure11 for σ=0.2.
Exact
First-order straightforward perturbation
Second-order straightforward perturbation
First-order triple product
Second-order triple product
Re
λ
1
2.5713
×
10
0
2.5983
×
10
0 1.04×100
2.5679
×
10
0 1.32×10−1
2.5714
×
10
0 3.18×10−3
2.5713
×
10
0 1.38×10−3
Re
λ
2
−
2.7703
×
10
0
−
2.7934
×
10
0 8.33×10−1
−
2.7670
×
10
0 1.17×10−1
−
2.7703
×
10
0 1.62×10−3
−
2.7703
×
10
0 1.48×10−3
Re
λ
3
−
3.3574
×
10
−
2
−
3.1784
×
10
−
2 5.33×100
−
3.3666
×
10
−
2 2.72×10−1
−
3.3579
×
10
−
2 1.41×10−2
−
3.3574
×
10
−
2 6.39×10−4
Re
λ
4
1.4107
×
10
−
1
1.3856
×
10
−
1 1.78×100
1.4064
×
10
−
1 3.06×10−1
1.4110
×
10
−
1 2.34×10−2
1.4107
×
10
−
1 2.92×10−3
Re
λ
5
−
1.6655
×
10
−
1
−
1.6344
×
10
−
1 1.86×100
−
1.6599
×
10
−
1 3.34×10−1
−
1.6660
×
10
−
1 2.97×10−2
−
1.6656
×
10
−
1 3.76×10−3Table 8
The imaginary parts of the exact and approximate natural frequencies for the perturbed whirling beam of Figure11 for σ=0.2.
Exact
First-order straightforward perturbation
Second-order straightforward perturbation
First-order triple product
Second-order triple product
Im
λ
3
4.6661
×
10
0
4.6783
×
10
0 2.61×10−1
4.6654
×
10
0 1.53×10−2
4.6662
×
10
0 4.82×10−4
4.6661
×
10
0 6.68×10−5
Im
λ
4
1.3917
×
10
1
1.3929
×
10
1 8.71×10−2
1.3919
×
10
1 1.53×10−2
1.3916
×
10
1 1.25×10−3
1.3917
×
10
1 1.59×10−4
Im
λ
5
1.5284
×
10
1
1.5272
×
10
1 7.93×10−2
1.5282
×
10
1 1.39×10−2
1.5284
×
10
1 1.14×10−3
1.5284
×
10
1 1.45×10−4Consider now a larger disorder strength ofσ=0.5. Figures 14 and 15 show the real and imaginary parts of the approximate eigenvalues found using the incremental perturbation method and the triple product method as a function of N. Note that as N increases, all the approximate eigenvalues converge to the exact results. Tables 9 and 10 tabulate the real and imaginary parts of the eigenvalues for N=7 so that the perturbations at each iteration are less than 10% of the original system matrices. The percentage errors for all five eigenvalues obtained using the second-order triple product method are below 0.005%. These results are substantially more accurate than those obtained using the second-order incremental perturbation method, at the cost of one additional computation given by Equation (40). Figure 16 shows the error parameter δ for the first-order and second-order perturbed modal matrices (with p=6) as a function of N. Note that as expected, δ decreases with increasing N, and that a second-order analysis is more accurate.Figure 14
The real parts of the approximate eigenvalues of the whirling beam as a function ofN for σ=0.5.
(a)
(b)
(c)
(d)
(e)Figure 15
The imaginary parts of the approximate eigenvalues of the whirling beam as a function ofN for σ=0.5.
(a)
(b)
(c)Table 9
The real parts of the exact and approximate natural frequencies for the perturbed whirling beam atN=7 for σ=0.5.
Exact
First-order incremental perturbation
Second-order incremental perturbation
First-order triple product
Second-order triple product
Re
λ
1
2.6377
×
10
0
2.6470
×
10
0 3.53×10−1
2.6342
×
10
0 1.30×10−1
2.6371
×
10
0 2.07×10−2
2.6376
×
10
0 2.05×10−3
Re
λ
2
−
2.7433
×
10
0
−
2.7499
×
10
0 2.40×10−1
−
2.7399
×
10
0 1.23×10−1
−
2.7432
×
10
0 4.96×10−3
−
2.7433
×
10
0 5.26×10−4
Re
λ
3
1.2603
×
10
−
2
1.4014
×
10
−
2 1.11×101
1.2605
×
10
−
2 1.42×10−2
1.2600
×
10
−
2 2.04×10−2
1.2603
×
10
−
2 1.45×10−4
Re
λ
4
1.8246
×
10
−
1
1.7763
×
10
−
1 2.64×100
1.8228
×
10
−
1 1.02×10−1
1.8257
×
10
−
1 5.89×10−2
1.8247
×
10
−
1 2.45×10−3
Re
λ
5
−
2.1690
×
10
−
1
−
2.1075
×
10
−
1 2.83×100
−
2.1669
×
10
−
1 9.61×10−2
−
2.1710
×
10
−
1 8.96×10−2
−
2.1691
×
10
−
1 4.37×10−3Table 10
The imaginary parts of the exact and approximate natural frequencies for the perturbed whirling beam atN=7 for σ=0.5.
Exact
First-order incremental perturbation
Second-order incremental perturbation
First-order triple product
Second-order triple product
Im
λ
3
5.0307
×
10
0
5.0397
×
10
0 1.79×10−1
5.0304
×
10
0 4.96×10−3
5.0307
×
10
0 1.54×10−4
5.0307
×
10
0 1.60×10−5
Im
λ
4
1.3720
×
10
1
1.3743
×
10
1 1.73×10−1
1.3721
×
10
1 6.34×10−3
1.3719
×
10
1 4.66×10−3
1.3720
×
10
1 2.16×10−4
Im
λ
5
1.5481
×
10
1
1.5458
×
10
1 1.53×10−1
1.5480
×
10
1 5.61×10−3
1.5482
×
10
1 4.13×10−3
1.5481
×
10
1 1.92×10−4Figure 16
The error parameter,δ, for the perturbed modal matrices of the whirling beam as a function of N for σ=0.5.


### 3.4. Nonuniform Rod

The perturbation method can also be applied to analyze the eigencharacteristics of a nonuniform rod. Specifically, the uniform rod will be considered as the unperturbed system, and any deviations from it will be treated as perturbations. Figure17 shows a nonuniform rod of length L, modulus E, density ρ, and cross-sectional area, Ax, that varies as a function of x from the fixed wall as follows:(53)Ax=A01+axLb,where A0, a, and b are constants. Using the assumed modes method [28], the longitudinal displacement wx,t of the rod can be expressed in the form of a finite series as follows:(54)wx,t=∑i=1suixηit,where s is the number of modes used in the expansion, uix are the eigenfunctions of the uniform rod, and ηit are the unknown generalized coordinates. Formulating the total kinetic and potential energies of the system, one obtains(55)T=12∑i=1s∑j=1sMijη˙itη˙jt,U=12∑i=1s∑j=1sKijηitηjt,where the generalized masses and stiffnesses of the nonuniform rod are given by(56)Mij=∫0LρAxuixujxdx,Kij=∫0LEAxduixdxdujxdxdx.Figure 17
A nonuniform rod.Applying Lagrange’s equations yields the following discretized equations of motion:(57)Mη¨+Kη=0,where the i,jth elements of matrices M and K are Mij and Kij, respectively. Moreover, note that M and K are both symmetric, and because s modes are used in the assumed modes expansion, they are of size s×s.For this particular system, the uniform rod whose cross-sectional area isA0 represents the unperturbed system, and the nonuniform rod is the perturbed system. For a uniform fixed-free rod, its normalized eigenfunctions are(58)uix=2ρA0Lsin2i−1πx2L,and its natural frequencies are(59)ω0i=2i−1π2EρL2.Because the unperturbed system consists of the uniform rod, the unperturbed matrices are(60)M0=I,K0=Λ0,where I denotes the identity matrix, and Λ0 is a diagonal matrix whose elements are ω0i2. For this undamped system, λ=ω2. Finally, Equations (33)–(35) still remain valid. In this case, of course, A=M and B=K. Note that unlike the previous examples, A and B are both symmetric.For definiteness, consider the case where the cross-sectional area varies as a function ofx with b=2 in Equation (53). In this case, the exact natural frequencies for the nonuniform rod are given by the solution of the following transcendental equation (see [29] for detailed derivations)(61)a1+atanβi=βi,where(62)ωi=βiEρL2.Figure18 shows the natural frequencies as a function of a, where a spans −0.35≤a≤0.35, obtained using s=10 in the assumed modes expansion. For a=0.35, the rod’s cross-sectional area at its free end is approximately 82% larger than its cross-sectional area at the base. Again, the natural frequencies obtained using the triple product method closely track the exact natural frequencies. The results are tabulated in Table 11. Although a 82% change in the rod’s cross-sectional area is significant, note that the largest error parameter for the first-order triple product method is less than 0.32%, while for the second-order triple product method it is less than 0.06%.Figure 18
The first five normalized exact and approximate natural frequencies of the nonuniform rod as a function ofa. The natural frequencies are nondimensionalized by dividing by E/ρL2.
(a)
(b)
(c)
(d)
(e)Table 11
The exact and approximate natural frequencies for the nonuniform rod of Figure17 for a=0.35. The natural frequencies are nondimensionalized by dividing by E/ρL2.
Exact
First-order straightforward perturbation
Second-order straightforward perturbation
First-order triple product
Second-order triple product
ω
1
1.3858
×
10
0
1.2825
×
10
0 7.45×100
1.4396
×
10
0 3.88×100
1.3889
×
10
0 2.19×10−1
1.3851
×
10
0 5.22×10−2
ω
2
4.6568
×
10
0
4.6242
×
10
0 6.98×10−1
4.6828
×
10
0 5.58×10−1
4.6422
×
10
0 3.12×10−1
4.6586
×
10
0 3.84×10−2
ω
3
7.8209
×
10
0
7.8014
×
10
0 2.48×10−1
7.8369
×
10
0 2.04×10−1
7.8109
×
10
0 1.27×10−1
7.8226
×
10
0 2.20×10−2
ω
4
1.0972
×
10
1
1.0958
×
10
1 1.26×10−1
1.0983
×
10
1 1.04×10−1
1.0964
×
10
1 6.73×10−2
1.0973
×
10
1 1.11×10−2
ω
5
1.4119
×
10
1
1.4108
×
10
1 7.67×10−2
1.4128
×
10
1 6.36×10−2
1.4113
×
10
1 4.13×10−2
1.4120
×
10
1 6.93×10−3Consider now a larger nonuniformity witha=0.55, implying that the cross-sectional area at the free end of the rod is 140% greater than its cross-sectional area at the fixed end. An incremental perturbation approach is used to obtain the approximate natural frequencies. Figure 19 shows the approximate natural frequencies obtained using the same five methods as a function of N. Although the approximate natural frequencies are inaccurate when N=1 or N=2, they quickly converge to the exact natural frequencies as N continues to increase. Table 12 tabulates the perturbed natural frequencies of Figure 19 for N=10 so that the perturbation at each iteration is less than 10% of the original system matrices. Note that the percentage errors, ε, of the first-order and second-order triple product approximations for all five natural frequencies are less than 0.006%. It should be noted that for this example, because the nonuniform rod has been discretized using the assumed modes method, the perturbed natural frequencies will converge to the natural frequencies of the discretized and not the continuum model of the nonuniform rod. As such, ω3 to ω5 obtained using the first-order triple product method with N=10 are more accurate than the ones found using the second-order triple product approach. This unexpected result is purely an artifact of the discretization process. If the number of modes, s, used in the assumed modes method is sufficiently large, then the second-order triple product approach will yield natural frequencies that are more accurate than those obtained using a first-order triple product method. Figure 20 shows the error parameters δ for the corresponding first-order and second-order modal matrices as a function of N. Note that for N=10, δ<0.001 for the second-order incremental perturbation method, validating the accuracy of the resulting modal matrices.Figure 19
The approximate natural frequencies of the nonuniform rod as a function ofN for a=0.55. The natural frequencies are nondimensionalized by dividing by E/ρL2.
(a)
(b)
(c)
(d)
(e)Table 12
The exact and approximate natural frequencies for the nonuniform rod atN=10 for a=0.55. The natural frequencies are nondimensionalized by dividing by E/ρL2.
Exact
First-order incremental perturbation
Second-order incremental perturbation
First-order triple product
Second-order triple product
ω
1
1.3053
×
10
0
1.2910
×
10
0 1.10×100
1.3060
×
10
0 5.31×10−2
1.3054
×
10
0 5.56×10−3
1.3053
×
10
0 1.10×10−3
ω
2
4.6359
×
10
0
4.6333
×
10
0 5.79×10−2
4.6366
×
10
0 1.42×10−2
4.6359
×
10
0 1.40×10−3
4.6360
×
10
0 1.18×10−3
ω
3
7.8085
×
10
0
7.8070
×
10
0 1.95×10−2
7.8090
×
10
0 5.59×10−3
7.8085
×
10
0 2.47×10−4
7.8086
×
10
0 1.07×10−3
ω
4
1.0963
×
10
1
1.0962
×
10
1 9.09×10−3
1.0963
×
10
1 3.65×10−3
1.0963
×
10
1 8.79×10−4
1.0963
×
10
1 1.30×10−3
ω
5
1.4112
×
10
1
1.4111
×
10
1 4.98×10−3
1.4112
×
10
1 2.75×10−3
1.4112
×
10
1 1.07×10−3
1.4112
×
10
1 1.33×10−3Figure 20
The error parameter,δ, for the perturbed modal matrices of the nonuniform rod as a function of N for a=0.55. The natural frequencies are nondimensionalized by dividing by E/ρL2.Finally, because the computational power of software and hardware these days is high, a few words on the advantages of the proposed methods are warranted. Consider a large, complicated system with a dozen or more parameters. In order to conduct a parametric study on the effects of parameter modifications on the eigenvalues of the system, computationally resolving the eigenvalue problem could be costly and time intensive. The perturbation methods proposed in this paper, on the other hand, can be used efficiently to conduct such an analysis because the proposed methods involve only simple matrix multiplication. In addition to their efficiency, these methods can accurately approximate a perturbed system’s eigensolutions. Moreover, the proposed methods are not limited by the magnitude of the perturbations. Using an incremental approach, these methods maintain their efficiency and accuracy even for large modifications. In addition, the perturbation methods yield physical insights into the problem that solving the problem computationally does not. Using the closed-form expressions given in Equations (34) and (35), one can track how changes to an individual parameter affect specific eigenvalues of the system. One can also calculate the rate at which certain eigenvalues change with respect to specific parameters. The perturbation methods can be harnessed to solve challenging inverse problems as well. Specifically, if the system requires specific eigenvalues and eigenvectors, one could determine which parameters must be changed to obtain said specifications. It would be difficult to solve such an inverse problem numerically on a computer because one could not predict how changing certain parameters may affect the system’s eigensolutions a priori without guessing and checking.


## 3.1. Torque-Free Rigid Body

Consider a torque-free rigid body as shown in Figure1 [5], which is symmetric about axis ζ. Two masses, each m/2, are located at distances ζ=±a from the fixed center O. Each mass is supported by 4 springs, each with stiffness k/4. The rigid body spins about axis ζ with a constant angular velocity of Ωζ. The mass moments of inertia of the rigid body about axes ξ, η, and ζ are α, β=α, and c, respectively. The mass moment of inertia of the entire body about a transverse axis is given by α′=α+ma2. The system’s unperturbed matrices are given by [5](41)where α0′=α0+m0a02. For definiteness, the original and unperturbed system’s parameters are Ωζ = Ω0 = 5 × 101 rad/s, ω0 = 6 × 101 rad/s, c0/α0 = 1.5, a0 = 2 m, m0 = 2.5 × 10−2 kg, k0 = 9 × 101 N/m, and α0 = 1 kg·m2.Figure 1
A torque-free rigid body.After an initial analysis has been performed, design changes are introduced. Letx denote an arbitrary system parameter. Assume it is perturbed systematically from its initial value of x0 as follows:(42)x=x01+σδx,where δx is arbitrarily chosen from a set of random numbers with zero mean and a standard deviation of one, and σ represents a disorder strength that can be varied. Together, σδx denotes the percentage change in x0. Thus, a large σ value implies a large perturbation from the nominal value. All the modified system parameters are perturbed in this manner.For the torque-free rigid body, the natural frequenciesω are the imaginary part of the eigenvalues λ, i.e., ω=imagλ. Suppose the modified parameters m and a deviate from their nominal values such that m=m01+σδm and a=a01+σδa, where δm = −4.4070 × 10−1 and δa = 8.4650 × 10−1. Figure 2 shows the natural frequencies as a function of σ. These natural frequencies are obtained via five different means: by solving the problem exactly, by using the straightforward perturbation method via a first-order or a second-order analysis, and using the triple product method of Equation (40), where the modal matrices consist of the first-order or the second-order perturbed eigenvectors. When σ=0, the system is unperturbed. Note that for small values of σ, all the approximate natural frequencies track the exact perturbed natural frequencies accurately. As the magnitude of σ increases, while the first-order and second-order straightforward perturbation approximations deviate from the exact values, the approximations obtained from the diagonal elements of the triple product, using either the first-order or second-order perturbed modal matrices, are indistinguishable from the exact natural frequencies. The fact that the diagonal elements track the exact natural frequencies so closely also validates the accuracy of the perturbed matrices. Table 1 tabulates these natural frequencies for σ=0.65. Thus, m and a deviate from their nominal values by approximately σδm≈−28.6% and σδa≈55.0%, respectively. The numerical value in the parentheses denotes the percentage error in the natural frequencies, defined as(43)ε=exact−approximateexact×100%.Figure 2
The exact and approximate natural frequencies (rad/s) of the torque-free rigid body as a function ofσ.
(a)
(b)
(c)Table 1
The exact and approximate natural frequencies for the perturbed torque-free rigid body of Figure1 for σ=0.65.
Exact
First-order straightforward perturbation
Second-order straightforward perturbation
First-order triple product
Second-order triple product
ω
1
1.2592
×
10
−
1
4.9493
×
10
−
1 2.93×102
1.9843
×
10
−
1 5.75×101
1.2801
×
10
−
1 1.66×100
1.2503
×
10
−
1 7.03×10−1
ω
2
4.3014
×
10
1
4.0276
×
10
1 6.36×100
4.2147
×
10
1 2.01×100
4.2925
×
10
1 2.08×10−1
4.3004
×
10
1 2.43×10−2
ω
3
1.1236
×
10
2
1.1203
×
10
2 2.99×10−1
1.1226
×
10
2 9.32×10−2
1.1236
×
10
2 1.98×10−3
1.1236
×
10
2 5.69×10−4Note the drastic improvement in the approximate natural frequencies obtained using a first-order triple product approach compared with using a first-order straightforward expansion, especially for the fundamental natural frequency. This significant improvement in accuracy is the result of carrying out another step consisting of a simple matrix multiplication given by Equation (40). Interestingly, the approximate natural frequencies obtained using a first-order triple product approach are even more accurate that those obtained using a second-order straightforward perturbation method. Note that among the different methods, the error parameters using the triple product with the second-order modal matrices are the smallest.If the error parameterε is deemed too large for a disorder strength of σ=0.65 or if a larger value of σ is applied, which would represent larger perturbations in the system, dividing the perturbations into smaller pieces and computing the approximations through an incremental method can be utilized. Although this incremental procedure involves more computation, it is deemed worthwhile when dealing with large perturbations in order to obtain accurate approximations. Consider now the case where σ=1.8. Figure 3 shows the natural frequencies obtained using the five methods mentioned previously as a function of N. These approximate natural frequencies monotonically converge to their respective exact values though some methods converge much more quickly than others. Note that the approximate fundamental natural frequency converges to the exact ω1 from above, while the other approximate natural frequencies converge to their respective exact values from below. By N=5, the triple product method with the first-order and second-order modal matrices leads to approximations that are nearly exact while the first-order incremental perturbation method is still very inaccurate by N=10. Across all values of N, the results of the triple product approximations via a second-order analysis are the most accurate. Table 2 illustrates the numerical results of Figure 3 for N=8. With eight iterations, the perturbations at each iteration are less than 10% of the original system matrices. At N=8, the errors for all three natural frequencies are reduced by at least an order of magnitude when approximating the new natural frequencies using the diagonal elements of the second-order triple product as opposed to simply the incremental perturbation method via a first-order analysis. With the second-order triple product, the final error percentages of all three natural frequencies are below 0.83%, substantially smaller than the percentage deviations in m0 and a0 that are introduced into the system, which are approximately σδm≈−79.3% and σδa≈152.3%, respectively. To assess the accuracy of the perturbed modal matrices for this torque-free rigid body system, consider the error parameter given by Equation (39), where p=6. Figure 4 shows δ as a function of N for σ=1.8, obtained using a first-order or a second-order incremental perturbation analysis. Note that δ decreases with N and that the error parameter for the second-order modal matrices decreases more quickly than the error parameter for the first-order modal matrices.Figure 3
The approximate natural frequencies of the torque-free rigid body as a function ofN for σ=1.8.
(a)
(b)
(c)Table 2
The exact and approximate natural frequencies for the perturbed torque-free rigid body atN=8 for σ=1.8.
Exact
First-order incremental perturbation
Second-order incremental perturbation
First-order triple product
Second-order triple product
ω
1
1.3511
×
10
0
2.2065
×
10
0 6.33×101
1.6401
×
10
0 2.13×101
1.4190
×
10
0 5.02×100
1.3622
×
10
0 8.22×10−1
ω
2
3.7885
×
10
1
3.4847
×
10
1 8.01×100
3.6785
×
10
1 2.90×100
3.7427
×
10
1 1.20×100
3.7813
×
10
1 1.89×10−1
ω
3
1.1178
×
10
2
1.1145
×
10
2 2.91×10−1
1.1167
×
10
2 9.99×10−2
1.1177
×
10
2 4.38×10−3
1.1178
×
10
2 1.84×10−5Figure 4
The error parameter,δ, for the perturbed modal matrices of the torque-free rigid body as a function of N for σ=1.8.


## 3.2. Two-Degrees-of-Freedom Gyroscopic System

Figure5 shows a two-degrees-of-freedom gyroscopic system analyzed by Meirovitch and Ryland [10]. The governing equations for the unperturbed system are given by(44)M0q¨0+C0+G0q˙0+K0q0=0,where an overdot denotes a time derivative, [M0], [C0], and [K0] are the mass, damping, and stiffness matrices of the nominal system, respectively, [G0] denotes the unperturbed gyroscopic damping matrix, and q0 is the vector of generalized coordinates for the nominal system. The system matrices are all of size 2×2, and q0 is of length 2. The nominal system matrices are(45)M0=m000m0,G0=0−2m0Ω02m0Ω00,C0=c0000,K0=k01−m0Ω0200k02−m0Ω02.Figure 5
A two-degree-of-freedom gyroscopic system.To solve for the eigensolutions of the system shown in Equation (44), the state vector approach [28] is used. Introducing a state vector of length 4 of the form(46)x0=q˙0q0,allows Equation (44) to be rewritten as(47)A0x˙0−B0x0=0,where A0 and B0 are the nominal state matrices of size 4×4 given by(48)A0=M000K0,B0=−C0−G0−K0K00.Figure5 is an example of a damped, asymmetric system due to the presence of the gyroscopic damping matrix G0. For definiteness, the original or unperturbed system parameters are m0 = 1 kg, Ω0 = 1 rad/s, c0 = 1 × 10−1 N·s/m, k01 = 3 N/m, and k02 = 4 N/m. To illustrate that the perturbation methods presented in this paper can also be used when an asymmetric damping matrix is perturbed, the parameters c0 and m0 are subjected to modifications; thus c=c01+σδc and m=m01+σδm. The random numbers used for the perturbations are δc=−7.2540×10−1 and δm=2.8480×10−1. Figures 6 and 7 illustrate the real and imaginary parts of the system’s eigenvalues (obtained exactly, with the straightforward perturbation method via a first-order or a second-order analysis, and with the diagonal elements of Equation (40) using the first-order or second-order perturbed modal matrices) as a function of σ. Note the diagonal elements of both triple products almost lie on top of the exact eigenvalues, which validates the accuracy of the perturbed modal matrices, and they are more accurate than the second-order straightforward perturbation approximations over the entire range of σ. Tables 3 and 4 compare the real and imaginary parts of the exact and approximate natural frequencies for σ=0.8. Note again that the first-order triple product method yields approximate eigenvalues that are more accurate than those obtained using a second-order straightforward perturbation analysis. Moreover, note that using the second-order triple product method, the error parameters for both the real and imaginary parts of the eigenvalues are reduced by at least two orders of magnitude compared to those when using the first-order straightforward perturbation method. Additionally, for σ=0.8, the real parts of the eigenvalues obtained using the second-order triple product method have ε that are all less than 0.16%, and the imaginary parts have ε all less than 0.0095%. These errors are significantly smaller than the percentage changes in c0 and m0, which are approximately σδm≈22.7% and σδc≈−58.0%, respectively.Figure 6
The real parts of the exact and approximate eigenvalues of the two-degree-of-freedom gyroscopic system as a function ofσ.
(a)
(b)Figure 7
The imaginary parts of the exact and approximate eigenvalues of the two-degree-of-freedom gyroscopic system as a function ofσ.
(a)
(b)Table 3
The real parts of the exact and approximate natural frequencies for the perturbed two-degree-of-freedom gyroscopic system of Figure5 for σ=0.8.
Exact
First-order straightforward perturbation
Second-order straightforward perturbation
First-order triple product
Second-order triple product
Re
λ
1
−
4.5435
×
10
−
3
−
9.7877
×
10
−
4 7.84×101
−
5.4541
×
10
−
3 2.00×101
−
4.4447
×
10
−
3 2.17×100
−
4.5366
×
10
−
3 1.50×10−1
Re
λ
2
−
1.2546
×
10
−
2
−
8.6132
×
10
−
3 3.13×101
−
1.3344
×
10
−
2 6.35×100
−
1.2565
×
10
−
2 1.53×10−1
−
1.2540
×
10
−
2 4.60×10−2Table 4
The imaginary parts of the exact and approximate natural frequencies for the perturbed two-degree-of-freedom gyroscopic system of Figure5 for σ=0.8.
Exact
First-order straightforward perturbation
Second-order straightforward perturbation
First-order triple product
Second-order triple product
Im
λ
1
6.7034
×
10
−
1
6.4020
×
10
−
1 4.49×100
6.7603
×
10
−
1 8.49×10−1
6.7000
×
10
−
1 5.01×10−2
6.7027
×
10
−
1 9.43×10−3
Im
λ
2
2.6928
×
10
0
2.6623
×
10
0 1.13×100
2.6985
×
10
0 2.11×10−1
2.6927
×
10
0 2.90×10−3
2.6928
×
10
0 3.95×10−4Assume larger perturbations are introduced to the system by usingσ=1.37. Figures 8 and 9 show the eigenvalues found using the five methods mentioned previously as a function of N. Note that using a first-order incremental perturbation method, the real parts of the eigenvalue converge to their exact values from above, while the imaginary parts converge to their exact values from below. The trend is reversed using the second-order incremental perturbation analysis, i.e., the real and imaginary parts converge to their exact values from below and above, respectively. For the chosen set of structural parameters, the modified system is stable. However, Figure 8 shows that the first-order incremental perturbation method yields approximate eigenvalues that erroneously indicate the system is unstable. In the case, because the real components of the system’s exact eigenvalues are close to 0, it is important to use a sufficiently large value of N to accurately capture the system’s stability. Also, note that the approximate eigenvalues found using the diagonal elements of the first-order and second-order triple products converge rapidly to the exact eigenvalues for N as small as 3, and for N=10, the approximate and exact eigenvalues become indistinguishable, despite the large perturbations that are introduced in the system. Tables 5 and 6 show the corresponding numerical results. Note that ε of the second-order triple product is at least four orders of magnitude smaller than ε of the first-order incremental perturbation method. Figure 10 shows the error parameter δ for the first-order and second-order perturbed modal matrices, obtained using Equation (39) with p=4, as a function of N. It clearly shows a rapid decrease of δ, which implies a high accuracy in the perturbed eigenvectors, especially those obtained using a second-order analysis.Figure 8
The real parts of the approximate eigenvalues of the two-degree-of-freedom gyroscopic system as a function ofN for σ=1.37.
(a)
(b)Figure 9
The imaginary parts of the approximate eigenvalues of the two-degree-of-freedom gyroscopic system as a function ofN for σ=1.37.
(a)
(b)Table 5
The real parts of the exact and approximate natural frequencies for the perturbed two-degree-of-freedom gyroscopic system atN=10 for σ=1.37.
Exact
First-order incremental perturbation
Second-order incremental perturbation
First-order triple product
Second-order triple product
Re
λ
1
−
5.4248
×
10
−
5
5.7165
×
10
−
4 1.15×103
−
7.0836
×
10
−
5 3.05×101
−
5.0406
×
10
−
5 7.08×100
−
5.4177
×
10
−
5 1.30×10−1
Re
λ
2
−
1.6881
×
10
−
4
6.4366
×
10
−
4 4.81×102
−
1.8005
×
10
−
4 6.65×100
−
1.6942
×
10
−
4 3.62×10−1
−
1.6879
×
10
−
4 1.16×10−2Table 6
The imaginary parts of the exact and approximate natural frequencies for the perturbed two-degree-of-freedom gyroscopic system atN=10 for σ=1.37.
Exact
First-order incremental perturbation
Second-order incremental perturbation
First-order triple product
Second-order triple product
Im
λ
1
5.6913
×
10
−
1
5.6296
×
10
−
1 1.08×100
5.6989
×
10
−
1 3.63×10−2
5.6912
×
10
−
1 2.63×10−3
5.6913
×
10
−
1 4.71×10−5
Im
λ
2
2.5906
×
10
0
2.5843
×
10
0 2.43×10−1
2.5907
×
10
0 6.06×10−3
2.5906
×
10
0 1.08×10−4
2.5906
×
10
0 1.59×10−6Figure 10
The error parameter,δ, for the perturbed modal matrices of the two-degree-of-freedom gyroscopic system as a function of N for σ=1.37.


## 3.3. Whirling Beam with Gyroscopic Damping

Figure11 shows a whirling shaft with a lumped mass M attached at its midspan [12]. The shaft’s mass is uniformly distributed with mass per unit length m. The shaft has stiffnesses EIx and EIy about the x and y axis, respectively. The end supports provide bending moments proportional to the angular deflections, where the proportionality constants are K1 at Z=0 and K2 at Z=L. It is assumed that there is uniformly distributed external dissipation proportional by a constant h to the absolute velocity, and an internal dissipation proportional by a constant c to the velocity relative to the rotating axes. The equations of motion of the unperturbed system are given by [12](49)M0q¨0+C0+G0q˙0+K0+H0q0=0,where the coefficient matrices in Equation (49) are partitioned into s×s submatrices(50)M0=Ms00Ms,G0=0Gs−Gs0,C0=Cs00Cs,K0=Ks100Ks2,H0=0Hs−Hs0.Figure 11
A whirling beam.Thei,jth elements of the submatrices are given by [12](51)Mijs=m0L0δij+2M0siniπ2sinjπ2,Gijs=−2Ω0Mij,Cijs=c0+h0L0δij,Kijs1=2K01+K02cosiπcosjπiπL0jπL0+EI0xiπL02jπL02L0δij−Ω02Mij,Kijs2=2K01+K02cosiπcosjπiπL0jπL0+EI0yiπL02jπL02L0δij−Ω02Mij,Hijs=−h0Ω0L0δij,where i,j=1,…,s. For this example, if s=3, then all of the submatrices will be of size 3×3. For definiteness, the original or unperturbed system parameters are m0 = 1 kg/m, M0 = 10 kg, L0 = 1 m, EI0x = 8.10 × 10−2 N·m2, EI0y = 1.82 × 10−1 N·m2, K01 = K02 = 5 × 10−2 N/m, c0 = h0 = 2.5 × 10−1 N·s/m, Ω0 = 1.46 × 101 rad/s. The eigensolutions of the unperturbed system are obtained using the state vector approach, where the system matrices are given by(52)A0=0M0M0C0+G0,B0=M000−K0−H0.Assume now the parameters affecting the damping matricesC0 and G0 are perturbed so that the modified parameters are given by M=M01+σδM, m=m01+σδm, c=c01+σδc and h=h01+σδh, where δM=−1.1079×100, δm=2.5260×10−1, δc=2.4600×10−2, and δh=−1.1270×100. For the given set of system parameters, λ1 and λ2 are strictly real while the remaining eigenvalues are complex. Additionally, the real components of the first and fourth eigenvalues are positive, indicating that the system is unstable. Figures 12 and 13 illustrate the real and imaginary parts of the eigenvalues as a function of σ, obtained using the methods developed in this paper. The diagonal elements of the second-order triple product again track the exact eigenvalues most closely for all values of σ. Tables 7 and 8 tabulate the numerical results for σ=0.2. For the given set of random variations, the parameters M, m, c, and h deviate from their nominal values by approximately −22.1%, 5.05%, 0.49%, and −22.5%, respectively. Note that the percentage errors for the diagonal elements of Equation (40) with the second-order perturbed modal matrices are all under 0.004%, despite the large deviations in some of the selected system parameters.Figure 12
The real parts of the exact and approximate eigenvalues of the whirling beam as a function ofσ.
(a)
(b)
(c)
(d)
(e)Figure 13
The imaginary parts of the exact and approximate eigenvalues of the whirling beam as a function ofσ.
(a)
(b)
(c)Table 7
The real parts of the exact and approximate natural frequencies for the perturbed whirling beam of Figure11 for σ=0.2.
Exact
First-order straightforward perturbation
Second-order straightforward perturbation
First-order triple product
Second-order triple product
Re
λ
1
2.5713
×
10
0
2.5983
×
10
0 1.04×100
2.5679
×
10
0 1.32×10−1
2.5714
×
10
0 3.18×10−3
2.5713
×
10
0 1.38×10−3
Re
λ
2
−
2.7703
×
10
0
−
2.7934
×
10
0 8.33×10−1
−
2.7670
×
10
0 1.17×10−1
−
2.7703
×
10
0 1.62×10−3
−
2.7703
×
10
0 1.48×10−3
Re
λ
3
−
3.3574
×
10
−
2
−
3.1784
×
10
−
2 5.33×100
−
3.3666
×
10
−
2 2.72×10−1
−
3.3579
×
10
−
2 1.41×10−2
−
3.3574
×
10
−
2 6.39×10−4
Re
λ
4
1.4107
×
10
−
1
1.3856
×
10
−
1 1.78×100
1.4064
×
10
−
1 3.06×10−1
1.4110
×
10
−
1 2.34×10−2
1.4107
×
10
−
1 2.92×10−3
Re
λ
5
−
1.6655
×
10
−
1
−
1.6344
×
10
−
1 1.86×100
−
1.6599
×
10
−
1 3.34×10−1
−
1.6660
×
10
−
1 2.97×10−2
−
1.6656
×
10
−
1 3.76×10−3Table 8
The imaginary parts of the exact and approximate natural frequencies for the perturbed whirling beam of Figure11 for σ=0.2.
Exact
First-order straightforward perturbation
Second-order straightforward perturbation
First-order triple product
Second-order triple product
Im
λ
3
4.6661
×
10
0
4.6783
×
10
0 2.61×10−1
4.6654
×
10
0 1.53×10−2
4.6662
×
10
0 4.82×10−4
4.6661
×
10
0 6.68×10−5
Im
λ
4
1.3917
×
10
1
1.3929
×
10
1 8.71×10−2
1.3919
×
10
1 1.53×10−2
1.3916
×
10
1 1.25×10−3
1.3917
×
10
1 1.59×10−4
Im
λ
5
1.5284
×
10
1
1.5272
×
10
1 7.93×10−2
1.5282
×
10
1 1.39×10−2
1.5284
×
10
1 1.14×10−3
1.5284
×
10
1 1.45×10−4Consider now a larger disorder strength ofσ=0.5. Figures 14 and 15 show the real and imaginary parts of the approximate eigenvalues found using the incremental perturbation method and the triple product method as a function of N. Note that as N increases, all the approximate eigenvalues converge to the exact results. Tables 9 and 10 tabulate the real and imaginary parts of the eigenvalues for N=7 so that the perturbations at each iteration are less than 10% of the original system matrices. The percentage errors for all five eigenvalues obtained using the second-order triple product method are below 0.005%. These results are substantially more accurate than those obtained using the second-order incremental perturbation method, at the cost of one additional computation given by Equation (40). Figure 16 shows the error parameter δ for the first-order and second-order perturbed modal matrices (with p=6) as a function of N. Note that as expected, δ decreases with increasing N, and that a second-order analysis is more accurate.Figure 14
The real parts of the approximate eigenvalues of the whirling beam as a function ofN for σ=0.5.
(a)
(b)
(c)
(d)
(e)Figure 15
The imaginary parts of the approximate eigenvalues of the whirling beam as a function ofN for σ=0.5.
(a)
(b)
(c)Table 9
The real parts of the exact and approximate natural frequencies for the perturbed whirling beam atN=7 for σ=0.5.
Exact
First-order incremental perturbation
Second-order incremental perturbation
First-order triple product
Second-order triple product
Re
λ
1
2.6377
×
10
0
2.6470
×
10
0 3.53×10−1
2.6342
×
10
0 1.30×10−1
2.6371
×
10
0 2.07×10−2
2.6376
×
10
0 2.05×10−3
Re
λ
2
−
2.7433
×
10
0
−
2.7499
×
10
0 2.40×10−1
−
2.7399
×
10
0 1.23×10−1
−
2.7432
×
10
0 4.96×10−3
−
2.7433
×
10
0 5.26×10−4
Re
λ
3
1.2603
×
10
−
2
1.4014
×
10
−
2 1.11×101
1.2605
×
10
−
2 1.42×10−2
1.2600
×
10
−
2 2.04×10−2
1.2603
×
10
−
2 1.45×10−4
Re
λ
4
1.8246
×
10
−
1
1.7763
×
10
−
1 2.64×100
1.8228
×
10
−
1 1.02×10−1
1.8257
×
10
−
1 5.89×10−2
1.8247
×
10
−
1 2.45×10−3
Re
λ
5
−
2.1690
×
10
−
1
−
2.1075
×
10
−
1 2.83×100
−
2.1669
×
10
−
1 9.61×10−2
−
2.1710
×
10
−
1 8.96×10−2
−
2.1691
×
10
−
1 4.37×10−3Table 10
The imaginary parts of the exact and approximate natural frequencies for the perturbed whirling beam atN=7 for σ=0.5.
Exact
First-order incremental perturbation
Second-order incremental perturbation
First-order triple product
Second-order triple product
Im
λ
3
5.0307
×
10
0
5.0397
×
10
0 1.79×10−1
5.0304
×
10
0 4.96×10−3
5.0307
×
10
0 1.54×10−4
5.0307
×
10
0 1.60×10−5
Im
λ
4
1.3720
×
10
1
1.3743
×
10
1 1.73×10−1
1.3721
×
10
1 6.34×10−3
1.3719
×
10
1 4.66×10−3
1.3720
×
10
1 2.16×10−4
Im
λ
5
1.5481
×
10
1
1.5458
×
10
1 1.53×10−1
1.5480
×
10
1 5.61×10−3
1.5482
×
10
1 4.13×10−3
1.5481
×
10
1 1.92×10−4Figure 16
The error parameter,δ, for the perturbed modal matrices of the whirling beam as a function of N for σ=0.5.


## 3.4. Nonuniform Rod

The perturbation method can also be applied to analyze the eigencharacteristics of a nonuniform rod. Specifically, the uniform rod will be considered as the unperturbed system, and any deviations from it will be treated as perturbations. Figure17 shows a nonuniform rod of length L, modulus E, density ρ, and cross-sectional area, Ax, that varies as a function of x from the fixed wall as follows:(53)Ax=A01+axLb,where A0, a, and b are constants. Using the assumed modes method [28], the longitudinal displacement wx,t of the rod can be expressed in the form of a finite series as follows:(54)wx,t=∑i=1suixηit,where s is the number of modes used in the expansion, uix are the eigenfunctions of the uniform rod, and ηit are the unknown generalized coordinates. Formulating the total kinetic and potential energies of the system, one obtains(55)T=12∑i=1s∑j=1sMijη˙itη˙jt,U=12∑i=1s∑j=1sKijηitηjt,where the generalized masses and stiffnesses of the nonuniform rod are given by(56)Mij=∫0LρAxuixujxdx,Kij=∫0LEAxduixdxdujxdxdx.Figure 17
A nonuniform rod.Applying Lagrange’s equations yields the following discretized equations of motion:(57)Mη¨+Kη=0,where the i,jth elements of matrices M and K are Mij and Kij, respectively. Moreover, note that M and K are both symmetric, and because s modes are used in the assumed modes expansion, they are of size s×s.For this particular system, the uniform rod whose cross-sectional area isA0 represents the unperturbed system, and the nonuniform rod is the perturbed system. For a uniform fixed-free rod, its normalized eigenfunctions are(58)uix=2ρA0Lsin2i−1πx2L,and its natural frequencies are(59)ω0i=2i−1π2EρL2.Because the unperturbed system consists of the uniform rod, the unperturbed matrices are(60)M0=I,K0=Λ0,where I denotes the identity matrix, and Λ0 is a diagonal matrix whose elements are ω0i2. For this undamped system, λ=ω2. Finally, Equations (33)–(35) still remain valid. In this case, of course, A=M and B=K. Note that unlike the previous examples, A and B are both symmetric.For definiteness, consider the case where the cross-sectional area varies as a function ofx with b=2 in Equation (53). In this case, the exact natural frequencies for the nonuniform rod are given by the solution of the following transcendental equation (see [29] for detailed derivations)(61)a1+atanβi=βi,where(62)ωi=βiEρL2.Figure18 shows the natural frequencies as a function of a, where a spans −0.35≤a≤0.35, obtained using s=10 in the assumed modes expansion. For a=0.35, the rod’s cross-sectional area at its free end is approximately 82% larger than its cross-sectional area at the base. Again, the natural frequencies obtained using the triple product method closely track the exact natural frequencies. The results are tabulated in Table 11. Although a 82% change in the rod’s cross-sectional area is significant, note that the largest error parameter for the first-order triple product method is less than 0.32%, while for the second-order triple product method it is less than 0.06%.Figure 18
The first five normalized exact and approximate natural frequencies of the nonuniform rod as a function ofa. The natural frequencies are nondimensionalized by dividing by E/ρL2.
(a)
(b)
(c)
(d)
(e)Table 11
The exact and approximate natural frequencies for the nonuniform rod of Figure17 for a=0.35. The natural frequencies are nondimensionalized by dividing by E/ρL2.
Exact
First-order straightforward perturbation
Second-order straightforward perturbation
First-order triple product
Second-order triple product
ω
1
1.3858
×
10
0
1.2825
×
10
0 7.45×100
1.4396
×
10
0 3.88×100
1.3889
×
10
0 2.19×10−1
1.3851
×
10
0 5.22×10−2
ω
2
4.6568
×
10
0
4.6242
×
10
0 6.98×10−1
4.6828
×
10
0 5.58×10−1
4.6422
×
10
0 3.12×10−1
4.6586
×
10
0 3.84×10−2
ω
3
7.8209
×
10
0
7.8014
×
10
0 2.48×10−1
7.8369
×
10
0 2.04×10−1
7.8109
×
10
0 1.27×10−1
7.8226
×
10
0 2.20×10−2
ω
4
1.0972
×
10
1
1.0958
×
10
1 1.26×10−1
1.0983
×
10
1 1.04×10−1
1.0964
×
10
1 6.73×10−2
1.0973
×
10
1 1.11×10−2
ω
5
1.4119
×
10
1
1.4108
×
10
1 7.67×10−2
1.4128
×
10
1 6.36×10−2
1.4113
×
10
1 4.13×10−2
1.4120
×
10
1 6.93×10−3Consider now a larger nonuniformity witha=0.55, implying that the cross-sectional area at the free end of the rod is 140% greater than its cross-sectional area at the fixed end. An incremental perturbation approach is used to obtain the approximate natural frequencies. Figure 19 shows the approximate natural frequencies obtained using the same five methods as a function of N. Although the approximate natural frequencies are inaccurate when N=1 or N=2, they quickly converge to the exact natural frequencies as N continues to increase. Table 12 tabulates the perturbed natural frequencies of Figure 19 for N=10 so that the perturbation at each iteration is less than 10% of the original system matrices. Note that the percentage errors, ε, of the first-order and second-order triple product approximations for all five natural frequencies are less than 0.006%. It should be noted that for this example, because the nonuniform rod has been discretized using the assumed modes method, the perturbed natural frequencies will converge to the natural frequencies of the discretized and not the continuum model of the nonuniform rod. As such, ω3 to ω5 obtained using the first-order triple product method with N=10 are more accurate than the ones found using the second-order triple product approach. This unexpected result is purely an artifact of the discretization process. If the number of modes, s, used in the assumed modes method is sufficiently large, then the second-order triple product approach will yield natural frequencies that are more accurate than those obtained using a first-order triple product method. Figure 20 shows the error parameters δ for the corresponding first-order and second-order modal matrices as a function of N. Note that for N=10, δ<0.001 for the second-order incremental perturbation method, validating the accuracy of the resulting modal matrices.Figure 19
The approximate natural frequencies of the nonuniform rod as a function ofN for a=0.55. The natural frequencies are nondimensionalized by dividing by E/ρL2.
(a)
(b)
(c)
(d)
(e)Table 12
The exact and approximate natural frequencies for the nonuniform rod atN=10 for a=0.55. The natural frequencies are nondimensionalized by dividing by E/ρL2.
Exact
First-order incremental perturbation
Second-order incremental perturbation
First-order triple product
Second-order triple product
ω
1
1.3053
×
10
0
1.2910
×
10
0 1.10×100
1.3060
×
10
0 5.31×10−2
1.3054
×
10
0 5.56×10−3
1.3053
×
10
0 1.10×10−3
ω
2
4.6359
×
10
0
4.6333
×
10
0 5.79×10−2
4.6366
×
10
0 1.42×10−2
4.6359
×
10
0 1.40×10−3
4.6360
×
10
0 1.18×10−3
ω
3
7.8085
×
10
0
7.8070
×
10
0 1.95×10−2
7.8090
×
10
0 5.59×10−3
7.8085
×
10
0 2.47×10−4
7.8086
×
10
0 1.07×10−3
ω
4
1.0963
×
10
1
1.0962
×
10
1 9.09×10−3
1.0963
×
10
1 3.65×10−3
1.0963
×
10
1 8.79×10−4
1.0963
×
10
1 1.30×10−3
ω
5
1.4112
×
10
1
1.4111
×
10
1 4.98×10−3
1.4112
×
10
1 2.75×10−3
1.4112
×
10
1 1.07×10−3
1.4112
×
10
1 1.33×10−3Figure 20
The error parameter,δ, for the perturbed modal matrices of the nonuniform rod as a function of N for a=0.55. The natural frequencies are nondimensionalized by dividing by E/ρL2.Finally, because the computational power of software and hardware these days is high, a few words on the advantages of the proposed methods are warranted. Consider a large, complicated system with a dozen or more parameters. In order to conduct a parametric study on the effects of parameter modifications on the eigenvalues of the system, computationally resolving the eigenvalue problem could be costly and time intensive. The perturbation methods proposed in this paper, on the other hand, can be used efficiently to conduct such an analysis because the proposed methods involve only simple matrix multiplication. In addition to their efficiency, these methods can accurately approximate a perturbed system’s eigensolutions. Moreover, the proposed methods are not limited by the magnitude of the perturbations. Using an incremental approach, these methods maintain their efficiency and accuracy even for large modifications. In addition, the perturbation methods yield physical insights into the problem that solving the problem computationally does not. Using the closed-form expressions given in Equations (34) and (35), one can track how changes to an individual parameter affect specific eigenvalues of the system. One can also calculate the rate at which certain eigenvalues change with respect to specific parameters. The perturbation methods can be harnessed to solve challenging inverse problems as well. Specifically, if the system requires specific eigenvalues and eigenvectors, one could determine which parameters must be changed to obtain said specifications. It would be difficult to solve such an inverse problem numerically on a computer because one could not predict how changing certain parameters may affect the system’s eigensolutions a priori without guessing and checking.


## 4. Conclusion

In this paper, three methods are proposed based on the eigenvalue perturbation theory to approximate the eigencharacteristics of symmetric and asymmetric systems that have been subjected to modifications. These methods only require the eigensolutions of the nominal system, and they involve only simple matrix multiplications. The straightforward perturbation method is easy to implement, simple to code, efficient to apply, and it yields accurate results when the modifications are small. When large design changes are introduced, an incremental perturbation scheme is proposed whereby the large modifications are perturbed incrementally. A triple product method based on the orthogonality conditions of the perturbed eigenvectors offers yet another means to approximate the new eigenvalues of the modified system. It also serves to validate the accuracy of the perturbed eigenvectors. Numerical experiments are performed on various symmetric and asymmetric systems, and they show that the incremental perturbation method and the triple product method can accurately track the exact eigensolutions even when large design changes are introduced. The eigenvalue perturbation theory yields accurate eigencharacteristics and, more importantly, allows the designer to gain physical insights regarding the effects of modifications on the dynamic characteristics of the system.



---

*Source: 8609138-2018-10-17.xml*