# An Improved Method for Computing Eigenpair Derivatives of Damped System

Authors: Pingxin Wang; Jieer Wu; Xibei Yang
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

The calculation of eigenpair derivatives plays an important role in vibroengineering. This paper presents an improved algorithm for the eigenvector derivative of the damped systems by dividing it into a particular solution and general solution of the corresponding homogeneous equation. Compared with the existing methods, the proposed algorithm can significantly reduce the condition number of the equation for particular solution. Therefore, the relative errors of the calculated solutions are notably cut down. The results on two numerical examples show that such strategy is effective in reducing the condition numbers for both distinct and repeated eigenvalues.


---

## Body

## 1. Introduction

It is well-known that eigenvalues and eigenvectors of a structure represent the dynamic characteristics of the structure. The eigenvalues are the natural frequencies of vibration, and the eigenvectors are the shapes of these vibrational modes. According to different vibrational models, eigenvalues problems can be divided into four categories.Definition 1 (standard eigenvalue problem).
LetA∈Cn×n; the standard eigenvalue problem is to find scalars λ and nonzero vectors x∈Cn satisfying(1)Ax=λx.λ is called the eigenvalue of A and x is the eigenvectors corresponding to the eigenvalue λ of A, the set of the eigenvalues of A.Definition 2 (generalized eigenvalue problem).
LetA,B∈Cn×n, and the generalized eigenvalue problem is to find scalars λ and nonzero vectors x∈Cn satisfying(2)Ax=λBx.λ is called the eigenvalue of matrix pencil (A,B) and x is the eigenvectors corresponding to the eigenvalue of (A,B).Definition 3 (quadratic eigenvalue problem).
LetM,C,K∈Cn×n, and the quadratic eigenvalue problem is to find scalars λ and nonzero vectors u∈Cn satisfying(3)λ2M+λC+Ku=0.λ is called the eigenvalue of (M,C,K) and u is the eigenvectors corresponding to the eigenvalue λ of (M,C,K).Definition 4 (nonlinear eigenvalue problem).
LetD(λ)∈Cn×n is a matrix function with parameter λ; the nonlinear eigenvalue problem is to find scalars λ and nonzero vectors x∈Cn satisfying(4)Dλx=0.λ is called the eigenvalue of D(λ) a nd x is the eigenvectors corresponding to the eigenvalue λ of D(λ).In particular, whenD(λ)=A-λI,D(λ)=A-λB and D(λ)=λ2M+λC+K, nonlinear eigenvalue problem (4) degenerates into standard problem (1), generalized eigenvalue problem (2), and quadratic eigenvalue (3) respectively.Eigensensitivity analysis of eigenvalue and eigenvector plays an important role in a variety of problems, such as structural optimal design [1, 2], finite model updating [3], structural damage detection [4], and system identification [5]. One of the main tasks of eigensensitivity analysis is to compute the eigenpair derivatives with respect to structural design parameters. Methods for computing the derivatives of eigenvalues and eigenvectors have been studied by many researchers over the past four decades. Although calculating eigenvalue derivative is straightforward, determining eigenvector sensitivity raises several challenges, due in part to the singularity of the coefficient matrix of the linear equation of the eigenvector derivatives. Different numerical methods have been developed for computing derivatives of eigenvectors. Most of the methods for computing eigenpair derivatives of standard eigenvalues, generalized eigenvalues, and quadratic eigenvalue can be divided into five categories: modal method [6–12], algebraic method [6, 13–24], finite difference method [25], iterative method [26–28], and Nelson’s method [6, 29–39]. As for nonlinear eigenvalue problem, Mehrmann and Voss [40, 41] summarize the numerical solution of the general nonlinear eigenvalue problem, which can be divided into Newton method based on matrix decomposition [42–44], Newton method based on nonlinear equations [45], successive approximation method [46], subspace projection method [47–49], perimeter integration method [50–52], and so on.Damped system arises frequently in many areas such as applied mechanics, electrical oscillation, vibroacoustics, fluid dynamics, and signal processing. The equation of motion for the free vibration of a linear damped discrete system withn degrees of freedom can be expressed by(5)Mpq¨t+Cpq˙t+Kpqt=0, where p is a design parameter and matrices M(p),C(p), and K(p)∈Cn×n are symmetric, respectively, mass, damping, and stiffness matrices which are analytically dependent on parameter p. Suppose that q(t)=u(p)eλ(p)t (u(p) does not depend on time t) is a solution of (5). Substituting it into (5) will lead to the following quadratic eigenvalue problem:(6)λ2pMp+λpCp+Kpup=0. If λ(p) and nonzero vectors u(p) satisfy (6), λ(p) is called the eigenvalue and u(p) is the eigenvectors corresponding to the eigenvalue λ(p).Some numerical methods for computing eigenpair derivatives of quadratic eigenvalue problem are presented in [17–24]. Recently, Wang and Dai [39] extended Nelson’s method to compute the eigenvector derivatives of damped system by expressing the eigenvector derivatives as addition of particular solutions and homogeneous solutions. However, when we solve a set of linear equations, we hope the corresponding coefficient matrix has small condition number [34, 53]. The equations of finding particular solutions in [39] are in different levels in their values and it thus has large condition number. In order to overcome this problem, we present an improved algorithm for computing the derivatives of eigenvectors.The remainder of this paper is arranged as follows. In Section2, we derive the derivatives of eigenvalues and review some related methods for computing the derivatives eigenvectors. In Section 3, we focus on dealing with the derivatives of eigenvectors and proposed an algorithm which can significantly reduce the condition number of the equation for particular solution. In Section 4, two numerical examples are performed by our proposed method. Finally, we make some concluding remarks in Section 5.Throughout this paper, we use the following notation.Cm×n denotes the set of complex m×n matrices, Cn=Cn×1, C=C1. In is the identity matrix of order n. diag(a1,…,an) stands for the diagonal matrix with diagonal elements a1,…,an. AT denotes the transpose of a matrix A.


## 2. Methods for Sensitivity Analysis of Damped System

We assume thatλi(p) is the ith eigenvalue of (6) and ui(p) is the eigenvector corresponding to λi(p); that is,(7)λi2pMp+λipCp+Kpuip=0.In this paper, we consider eigenpair derivatives of damped system with distinct or repeated eigenvalues. Without loss of generality, we suppose thatλ1(p0)=λ2(p0)=⋯=λr(p0) are a semisimple eigenvalue with multiplicity r>1 of (6). Note that the case of r=1 corresponds to a distinct eigenvalue and the following results are also applicable.Sun [54] has pointed out that the derivatives of repeated eigenvalues are only directionally differentiable and the eigenvectors associated with multiple eigenvalues are not necessarily continuous functions of the design parameters and are generally not derivable. In order to ensure the existence of derivatives, we assume that λi(p) and ui(p) are sufficiently differentiable at p=p0 and eigenvalues are distinct in the neighborhood of p0. Under these assumptions, the eigenvectors of the quadratic eigenvalue problem (6) are uniquely determined, to within a normalizing condition, for all p≠p0 in the neighborhood of p0. Often the following normalization is adopted to ensure uniqueness of eigenvectors for simple eigenvalues:(8)uiTp2λipMp+Cpuip=1. As for repeated eigenvalues λ1(p0), we impose another normalization:(9)UTp02λ1p0Mp0+Cp0Up0=Ir, where U(p0)=u1p0,…,urp0.When an eigenvalue is repeated, the associated eigenvectors are only defined up to a subspace with dimension equal to the geometric multiplicity of the eigenvalue. The computed eigenvalues of (6) at p=p0 are λ1(p0)=λ2(p0)=⋯=λr(p0), while the corresponding eigenvectors may not be u1(p0),u2(p0),…,ur(p0) described as above. To simplify notation, here and henceforth, “p0” is omitted for variables evaluated at p=p0. Let the computed linearly independent eigenvectors corresponding to λ1 be the columns of the n×r matrix X. It is natural that the eigenvector matrix X satisfies (7) and the normalizing condition (9) at p=p0 too. Then the adjacent eigenvectors U=(u1,u2,…,ur) can be expressed in terms of X as follows:(10)U=XΓ, where Γ=[γ1,…,γr]∈Cr×r is a transformation matrix to be determined and should satisfy ΓTΓ=Ir.For convenience, the following notation is adopted in this study.(•)p≡∂(•)/∂p,D=λ12M+λ1C+K,Dp=λ12Mp+λ1Cp+Kp,H=2λ1M+C. Differentiating (7) with respect to design parameter p and letting p=p0 yield(11)Duip=-λipHui-Dpui. Set i=1,2,…,r in (11) and assemble them into a linear system of algebraic equations:(12)DUp=-HUΛp-DpU, where Λp=diag(λ1p,…,λrp). Premultiplying each side of (12) by UT, the eigenvalue derivatives can be obtained by(13)Λp=-UTDpU. Substituting (10) into (13) will derive a following eigenvalue problem:(14)ΞΓ=ΓΛp, where Ξ=-XTDpX. Equation (14) shows that diagonal elements of matrix Λp are n eigenvalues of Ξ and matrix Γ are corresponding to right eigenvectors matrix. Hence, solving eigenproblem ΞΓ=ΓΛp, yields eigenvalues derivative Λp and matrix Γ.The derivatives of the eigenvectors cannot be solved directly by (11) since the coefficient matrix is singular. There are many different methods to address the singularity problem such as modal method, algebraic method, Nelson’s method, iterative method, and finite difference method. Here we briefly introduce modal method, algebraic method, and Nelson’s method.


### 2.1. Modal Method

The basic idea of modal method is to express each eigenvector derivative as a linear combination of all the eigenvectors:(15)uip=∑k=12Ncikuk, where uk are the complete eigenvalues space and cik are sets of complex constants which are given by(16)cik=1λk-λiuiTλi2Mp+λiCp+Kpuik≠i-0.5uiT2λiMp+Cpuik=i.However, it is difficult to obtain all eigenvectors, in practice, especially for those large and complex engineering structures. So we can only use a part of the low-order modal as the basis vector and ignore the contribution of those unknown high-order modes, which is the so-called truncated modal method. The corrections to the problem of modal truncation error for damped system have been studied by several authors [14, 15]. However, most of the existing methods obtain the derivatives of complex mode shapes of viscously damped systems by considering the state-space forms which linearize the quadratic eigenvalue problem to standard or generalized eigenvalue problem. These procedures will complicate the derivations and the practical computations since the size of the problem is doubled.


### 2.2. Algebraic Method

In algebraic method, the eigenvector derivatives are calculated by assembling the derivatives of eigenproblems and the additional constraints obtained from the derivative of normalization condition into a linear system of algebraic equations. In this subsection, we briefly introduce the process of algebraic method for distinct eigenvalue [6].Suppose thatλi(p) is a distinct eigenvalue of (6) and ui(p) is the eigenvector corresponding to λi(p). Differentiating (8) with respect to p yields(17)uiTHuip=-λ1uiTMpui-uiTMuiΛp-0.5uiTCpui. Assembling (17) and the derivatives of (6) into a linear system of algebraic equation,(18)DHuiuiTH0uip0=-HuiΛp-Dpui-λ1uiTMpui-uiTMuiΛp-0.5uiTCpui, where the order of coefficient matrix on the left side of (18) is (N+1)×(N+1) and the matrix on the right side is (N+1)×1. The derivatives uip can be found by solving (18). The coefficient matrix can be decomposed into upper and lower triangular forms and then a forward and backward substitution scheme may be used to compute uip. This method has been extended to repeated eigenvalue of symmetric systems with viscous damping by [13]. However, Wu et al. [34] pointed out that there was a mistake in Lee’s extension. Later, Xu et al. [19] extended the algebraic method to compute the eigensolution sensitivities of asymmetric viscously damped systems. Li et al. [20, 21] extended the algebraic method to the more general nonviscous damped systems.


### 2.3. Nelson’s Method

By expressing the eigenvector derivatives as a particular solution and a homogeneous solution, Nelson [29] proposed an efficient method for computing the eigenpair derivatives for undamped systems, where only eigenpair under consideration is required. Friswell and Adhikari [31] extended Nelson’s method to symmetric and asymmetric systems with viscous damping by expressing the derivative of each eigenvector as a particular solution and a homogeneous solution of the singularity problem. Recently, Wang and Dai [39] extended Nelson’s method to compute the eigenvector derivatives of damped system with repeated eigenvalues. In [39], the solution of (12) can be expressed by(19)U′=V+UA, where U is the general solution of the corresponding homogeneous equations. A=(aij)∈Cr×r is a constant matrix. V is a particular solution of (12), that is,(20)DV=-HUΛp-DpU. The particular solution V may be required to be H-orthogonal with respect to U, that is,(21)UTHV=0. The particular solution V can be obtained by combining (20) and (21):(22)DHUUTH0V0=-HUΛp-DpU0.Using nonsingularity of the coefficient matrix in (22), we may solve (22) to achieve the particular solution V to (12). The computation of the constant matrix A in (19) employs the second-order differential information of (7). The method is an exact method and only requires the eigenvector of interest. Similar to algebraic method, Nelson’s method also needs matrix decomposition to obtain the particular solution for each eigensolution and calculate the particular solution by forward and backward substitutions.


## 2.1. Modal Method

The basic idea of modal method is to express each eigenvector derivative as a linear combination of all the eigenvectors:(15)uip=∑k=12Ncikuk, where uk are the complete eigenvalues space and cik are sets of complex constants which are given by(16)cik=1λk-λiuiTλi2Mp+λiCp+Kpuik≠i-0.5uiT2λiMp+Cpuik=i.However, it is difficult to obtain all eigenvectors, in practice, especially for those large and complex engineering structures. So we can only use a part of the low-order modal as the basis vector and ignore the contribution of those unknown high-order modes, which is the so-called truncated modal method. The corrections to the problem of modal truncation error for damped system have been studied by several authors [14, 15]. However, most of the existing methods obtain the derivatives of complex mode shapes of viscously damped systems by considering the state-space forms which linearize the quadratic eigenvalue problem to standard or generalized eigenvalue problem. These procedures will complicate the derivations and the practical computations since the size of the problem is doubled.


## 2.2. Algebraic Method

In algebraic method, the eigenvector derivatives are calculated by assembling the derivatives of eigenproblems and the additional constraints obtained from the derivative of normalization condition into a linear system of algebraic equations. In this subsection, we briefly introduce the process of algebraic method for distinct eigenvalue [6].Suppose thatλi(p) is a distinct eigenvalue of (6) and ui(p) is the eigenvector corresponding to λi(p). Differentiating (8) with respect to p yields(17)uiTHuip=-λ1uiTMpui-uiTMuiΛp-0.5uiTCpui. Assembling (17) and the derivatives of (6) into a linear system of algebraic equation,(18)DHuiuiTH0uip0=-HuiΛp-Dpui-λ1uiTMpui-uiTMuiΛp-0.5uiTCpui, where the order of coefficient matrix on the left side of (18) is (N+1)×(N+1) and the matrix on the right side is (N+1)×1. The derivatives uip can be found by solving (18). The coefficient matrix can be decomposed into upper and lower triangular forms and then a forward and backward substitution scheme may be used to compute uip. This method has been extended to repeated eigenvalue of symmetric systems with viscous damping by [13]. However, Wu et al. [34] pointed out that there was a mistake in Lee’s extension. Later, Xu et al. [19] extended the algebraic method to compute the eigensolution sensitivities of asymmetric viscously damped systems. Li et al. [20, 21] extended the algebraic method to the more general nonviscous damped systems.


## 2.3. Nelson’s Method

By expressing the eigenvector derivatives as a particular solution and a homogeneous solution, Nelson [29] proposed an efficient method for computing the eigenpair derivatives for undamped systems, where only eigenpair under consideration is required. Friswell and Adhikari [31] extended Nelson’s method to symmetric and asymmetric systems with viscous damping by expressing the derivative of each eigenvector as a particular solution and a homogeneous solution of the singularity problem. Recently, Wang and Dai [39] extended Nelson’s method to compute the eigenvector derivatives of damped system with repeated eigenvalues. In [39], the solution of (12) can be expressed by(19)U′=V+UA, where U is the general solution of the corresponding homogeneous equations. A=(aij)∈Cr×r is a constant matrix. V is a particular solution of (12), that is,(20)DV=-HUΛp-DpU. The particular solution V may be required to be H-orthogonal with respect to U, that is,(21)UTHV=0. The particular solution V can be obtained by combining (20) and (21):(22)DHUUTH0V0=-HUΛp-DpU0.Using nonsingularity of the coefficient matrix in (22), we may solve (22) to achieve the particular solution V to (12). The computation of the constant matrix A in (19) employs the second-order differential information of (7). The method is an exact method and only requires the eigenvector of interest. Similar to algebraic method, Nelson’s method also needs matrix decomposition to obtain the particular solution for each eigensolution and calculate the particular solution by forward and backward substitutions.


## 3. The Proposed Algorithm for Eigenvector Derivatives

The algorithm presented in [39] gives a simple method for computing the eigenvector derivatives of damped system with repeated eigenvalues. As we know, we hope that the calculated solution is a close representation of the true solution when we solve a set of linear equations. This requires that the corresponding coefficient matrix has small condition number [34, 53]. However, elements in the coefficient matrix of (22) are in different levels in their values and it thus has large condition number. In order to overcome this problem, we will propose an improved method to reduce condition number of the coefficient matrix in this section.Note that the coefficient matrixD in (11) is of rank N-r; it cannot be inverted. We assume that the eigenvector derivative uip have the following form:(23)uip=vi+Uci, where U is the general solution of the corresponding homogeneous equations and vi is a particular solution of (11); that is, vi satisfies(24)Dvi=-λipHui-Dpui. Particular solutions vi may be required to be H-orthogonal with respect to uj(j=1,…,r); that is,(25)ujTHvi=0,j=1,…,r. Combining (24) and (25), we obtain the following linear system:(26)DHu1Hu2⋯Huru1TH00⋯0u2TH00⋯0⋮⋮⋮⋱⋮urTH00⋯0vi00⋮0=-λipHui-Dpui00⋮0.The nonsingularity of the coefficient matrix in (26) has been proved in [39]. From above discussions, we know that the part vi of the solution to (26) is the particular solution to (11). Using nonsingularity of the coefficient matrix in (26), we can obtain the particular solution vi. When we solve a set of linear equations, we hope that the calculated solution is a close representation of the true solution. This requires that the corresponding coefficient matrix has small condition number [34, 53]. However, elements in the coefficient matrix of (26) are in different levels in their values and it thus has large condition number. In order to overcome this problem, we revise (26) into the following form:(27)Dl1Hu1l2Hu2⋯lrHurl1u1TH00⋯0l2u2TH00⋯0⋮⋮⋮⋱⋮lrurTH00⋯0vi00⋮0=-λipHui-Dpui00⋮0, where li=max1≤j≤N(hjj)/Hui∞,hjj is the jth diagonal element of H.It is easy to show that (26) and (27) have the same solution and vi given by that solution is the particular solution of (11). However, coefficient matrices in (26) and (27) have different condition numbers. We will show that condition number of (27) is much smaller than that of (26) in numerical examples in Section 4.Once the particular solutionvi is obtained, coefficient ci can be determined by the following [39]:(28)cii=-0.5uiTHpui-λipuiTMui,cij=gijλjp-λip,i≠j,j=1,…,r, where gij=-0.5uiTDppuj-λjpuiTHpuj-λjp2uiTMuj-uiTDpvj,Dpp=λ12Mpp+λ1Cpp+Kpp,Hp=2λ1Mp+Cp,(•)pp≡∂2(•)/∂p2.


## 4. Numerical Examples

In order to test the effectiveness of the proposed method, two numerical examples are employed in this section. The first one is a three-degree-of-freedom damped system with distinct eigenvalues and the second one is four-degree-of-freedom mass-spring-damper system with repeated eigenvalues. We use maximum row sum matrix norm·∞ to compute the condition numbers of (26) and (27). All codes are run in Matlab R2016a with machine precision 10-16 on a personal computer.Example 1.
Consider a three-degree-of-freedom damped system shown in Figure1. From Figure 1, we have the mass, stiffness, and damping matrices as (29)M=m1000m2000m3,C=c+c2+c3-c-c2-c2c+c1-c-c2-c2c+c2,K=k1+k2-k20-k2k2+k3+k5-k30-k3k3+k4.Figure 1
A three-degree-of-freedom damped system.Setm1=m2=m3=1kg,k1=k4=k5=1000N/m,k2=k3=0N/m,c1=10Ns/m,c2=c3=5Ns/m and c is chosen as the parameter. The system has six distinct eigenvalues at c=0Ns/m, which are -2.5000±31.5238i,-1.4645±31.5888i and -8.5355±30.4491i and the corresponding normalized eigenvectors are (30)u1=00.0891+0.0891i0,u2=0.0340+0.0340i00.0822+0.0822i,u3=0.0837+0.0837i0-0.0347-0.0347i. Substituting the above eigenvectors to (26) and (27), respectively, we compute the condition numbers of (26) and (27) and list the results in Table 1.Table 1
Condition numbers of coefficient matrices in (26) and (27).
Eigenvalue Condition number of (26) Condition number of (27) - 2.500 0  - 3 1.523 8 i 59.7389 9.6172 - 2.5000  + 31.523 8 i 59.7389 9.6172 - 1.464 5  - 3 1.588 8 i 89.2379 11.2299 - 1.4645  + 31.588 8 i 89.2379 11.2299 - 8.535 5  - 3 0.449 1 i 90.8835 2.5304 - 8.5355  + 30.449 1 i 90.8835 2.5304From Table1, we can see that the condition numbers of (27) are considerably reduced. Therefore, the relative errors of the calculated solutions in the proposed method are notably cut down. Example 1 shows that the proposed algorithm is effective when the eigenvalues of damped are distinct. Next, we will give another example with repeated eigenvalues.Example 2.
Figure2 shows a four-degree-of-freedom mass-spring-damper system, which is an example of [39]. We recalculate the example by using our proposed method and compare condition number of (27) with the result of (26). Let m1=m2=m3=m4=1kg,k1=k3=0,k2=k4=103N/m,c1=c3=0Ns/m,c2=1Ns/m. The damping c4 is chosen as the parameter. It is easy to verify that the system has two 2-repeated eigenvalues: -1.0000+44.7102i,-1.0000-44.7102i and one 4-repeated eigenvalue: 0 at c4=1Ns/m. We do not consider eigenvalue 0 because it is defective in this problem. For the eigenvalues -1.0000+44.7102i and -1.0000-44.7102i, the corresponding normalized eigenvectors are (31)U=0.0529-0.0529i0-0.0529+0.0529i000.0529-0.0529i0-0.0529+0.0529i. Substituting the above eigenvectors to (26) and (27), we can get the condition numbers of coefficient matrices as in Table 2.Table 2
Condition numbers of coefficient matrices in (26) and (27).
Eigenvalue Condition number of (26) Condition number of (27) - 1.0000  + 44.7102 i 223.3214 22.3665 - 1.0000  - 44.7102 i 223.3214 22.3665Figure 2
A four-degree-of-freedom mass-spring-damper system.From Table2, we can see that the condition numbers of coefficient matrices in (27) decrease significantly. Example 2 shows that our proposed method is effective in reducing the condition numbers for repeated eigenvalues.


## 5. Concluding Remarks

This paper has outlined a numerical method for computing the derivatives of eigenvalues and corresponding eigenvectors of symmetric damped system by dividing eigenvector derivatives into a particular solution and general solution of the corresponding homogeneous equation. The proposed algorithm requires only the information of those eigenvectors corresponding to the repeated eigenvalues to extend the governing equations of particular solutions. Compared with the existing methods, the proposed algorithm can significantly reduce the condition number of the governing equations of particular solutions. Therefore, the relative errors of the calculated solutions are notably cut down. Numerical examples have demonstrated the validity of the proposed method for distinct eigenvalues and repeated eigenvalues. The present study is the first step for the research of eigensensitivity. The high-order derivatives of eigenvalues and eigenvectors are specially interesting topic that remains to be further investigated.



---

*Source: 8050132-2018-02-21.xml*