JOURNAL OF SOUND AND VIBRATION
Journal of Sound and Vibration 290 (2006) 839–864
Symmetry preserving eigenvalue embedding in finite-element
model updating of vibrating structures
Joao B. Carvalho
a, Biswa N. Datta
b,, Wen-Wei Lin
c, Chern-Shuh Wang
d aDepartment of Pure and Applied Mathematics, Universidade Federal do Rio Grande do Sul, Porto Alegre, Brazil
bDepartment of Mathematical Sciences, Vibration and Acoustic Center of the College of Engineering and Engineering
Technology, Northern Illinois University, DeKalb, IL 60115, USA
cDepartment of Mathematics, National Tsing Hua University, Hsinchu 30043, Taiwan, ROC d
Department of Mathematics, National Cheng Kung University, Tainan 701, Taiwan, ROC Received 15 September 2003; received in revised form 5 April 2005; accepted 26 April 2005
Available online 19 August 2005
Abstract
The eigenvalue embedding problem addressed in this paper is the one of reassigning a few troublesome eigenvalues of a symmetric finite-element model to some suitable chosen ones, in such a way that the updated model remains symmetric and the remaining large number of eigenvalues and eigenvectors of the original model is to remain unchanged. The problem naturally arises in stabilizing a large-scale system or combating dangerous vibrations, which can be responsible for undesired phenomena such as resonance, in large vibrating structures. A new computationally efficient and symmetry preserving method and associated theories are presented in this paper. The model is updated using low-rank symmetric updates and other computational requirements of the method include only simple operations such as matrix multiplications and solutions of low-order algebraic linear systems. These features make the method practical for large-scale applications. The results of numerical experiments on the simulated data obtained from the Boeing company and on some benchmark examples are presented to show the accuracy of the method. Computable error bounds for the updated matrices are also given by means of rigorous mathematical analysis.
r2005 Elsevier Ltd. All rights reserved.
www.elsevier.com/locate/jsvi
0022-460X/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.jsv.2005.04.030
Corresponding author.
E-mail addresses: [email protected] (J.B. Carvalho), [email protected] (B.N. Datta), [email protected] (W.-W. Lin), [email protected] (C.-S. Wang).
1. Introduction
Vibrating structures such as bridges, highways, buildings, automobiles, air and space crafts, etc., are very often modelled by using finite-element methods (FEMs). These methods generate structured systems of matrix second-order differential equations of the form
M €x þ C _x þ Kx ¼ 0, (1)
where the coefficient matrices M, C and K are called, respectively, the mass, damping and stiffness matrices. In most applications, these matrices have very special exploitable properties such as the symmetry, positive definiteness, sparsity and others. The matrix M is often symmetric positive definite and denoted by M40; and K is symmetric positive semi-definite, denoted by KX0. The damping matrix C is hard to determine in practice; however, very often, for the sake of computational convenience and other practical considerations, it is assumed to be symmetric.
It is critical and very important that these properties are preserved while solving a vibration problem or updating a FEM to achieve certain design objectives.
In this paper, we will assume throughout that M40; K40 and C ¼ CT.
The classical approach is to use separation of variables, accounting for a solution xðtÞ ¼ yeltto
(1), where y is a constant vector. This leads to the quadratic matrix eigenvalue problem F ðlkÞyk ¼0; k ¼ 1; 2; . . . ; 2n,
where
F ðlÞ ¼ l2M þ lC þ K (2)
is the so-called associated quadratic matrix pencil. The quantities ðlk; ykÞ, k ¼ 1;. . . ; 2n are the
eigenpairs of the pencil (2).
It is well-known [1] that the dynamical behavior of a vibrating system, which can show
undesired phenomena such as instability and resonance, is determined by their natural frequencies and corresponding mode shapes, that is, the eigenvalues and eigenvectors of the pencil F ðlÞ. It is desirable that such behaviors are altered by making minimal changes in the system and keeping the structural properties invariant, as much as possible. Realistically, while dealing with a large system, it is often found in practice that only a small number of eigenvalues are ‘‘troublesome’’. Thus, it makes sense to reassign to suitable locations, chosen by the designer, only these troublesome eigenvalues, while keeping the remaining large number of eigenvalues unchanged.
Such a problem in control theory is known as the partial pole-placement problem and feedback control is used to solve this problem. For the standard first-order state–space systems of the form
_
xðtÞ ¼ AxðtÞ þ BuðtÞ, though there exist many numerical methods for the complete pole-placement
(see Ref. [2] for details), only two methods have so far been developed for the partial
pole-placement problem: (i) the projection method due to Saad [3], and (ii) the Sylvester equation
method by Datta and Sarkissian [4]. For a matrix second-order system, the choices are either to
transform the latter to a standard first-order form and then use one of the above methods or to
use the Independent Modal Space Control (IMSC) approach [1]. Both have some severe
engineering and computational limitations. The first approach might require an ill-conditioned matrix inversion or solution of a descriptor control problem (no method still exists for the partial pole-placement in descriptor systems). The IMSC approach requires complete knowledge of the
spectrum and the associated eigenvectors of the quadratic pencil (1) for decoupling of the open-loop pencil. Furthermore, the decoupling of the closed-open-loop pencil requires some very stringent
conditions on actuators and sensors[1], which is unpractical for real-life applications.
In several recent papers[5–10]numerically effective methods have been developed for both the
partial pole-placement and eigenstructure assignment problems; they overcome the difficulties associated with the above two approaches. These methods are designed directly in matrix second-order setting without resorting to first-second-order transformations and without requiring complete
knowledge of the spectrum of the pencil F ðlÞ, as needed by the IMSC approach[1]. Although they
satisfy control design requirements and are practical for control applications, unfortunately, they are not capable of preserving the symmetry of the original model.
In this paper, a novel symmetry preserving partial spectrum assignment method for vibrating system (1) is proposed. Specifically, the following problem is solved:
Let flig2ni¼1 and fyig2ni¼1 be, respectively, the spectrum and the eigenvector set of F ðlÞ. Given (i)
symmetric n n matrices M; C; and K of the pencil (2) with M40; KX0, and C ¼ CT, (ii) a part
of the spectrum fl1; . . . ; lrg; rp2n of FðlÞ and the corresponding eigenvectors fy1; . . . ; yrg, and
(iii) a set of r complex conjugate numbers fm;. . . ; mrg. Assuming the both sets fl1; . . . ; lng and
fm1; . . . ; mrgare closed under complex conjugations, find real symmetric matrices Mnew; Cnew; and
Knew such that the spectrum of Fnew¼l2MnewþlCnewþKnew is fm1; . . . ; mr; lrþ1; . . . ; l2ng
andfurthermore, the eigenvectors corresponding to lrþ1; . . . ; l2nremain unchanged. Furthermore,
characterize the eigenvectors of Fnew corresponding to m1; . . . ; mr.
The last property is highly significant from practical applications view points. It says that certain important physical properties of the system are completely preserved by updating. However, the most important benefits obtained by this new method over the existing non-symmetric pole-placement methods for the second-order model are that the updated model remains symmetric and the changes made in the data matrices M; K; and C might be significantly less than those obtained by the pole-placement algorithms usingg feedback control.
To distinguish this problem from the partial pole-placement problem in control theory, we will call this problem ‘‘Eigenvalue Embedding’’ Problem (EEP). Our major contributions to EEP in this paper are as follows:
(i) An algorithm and associated theories are developed, using low-rank symmetric updates. (ii) Computable error bounds are derived by means of rigorous error analysis.
(iii) The accuracy of the algorithm is demonstrated by both an illustrative, and a real-life example with simulated data from the Boeing Company.
(iv) A complete characterization of the eigenvectors of the updated model is also given. It is shown by mathematical proofs that the eigenvectors corresponding to the eigenvalues which are not reassigned also remain invariant.
Finally, it is noted that the EEP addressed in this paper is clearly related to the well-known problem in vibrating engineering, called ‘‘Finite-Element Model Updating Problem’’ (FEMUP). The FEMUP is concerned with updating a symmetric FEM in such a way that the updated model remains symmetric and a set of measured eigenvalues and eigenvectors are incorporated into the updated model, while the other eigenvalues and eigenvectors remain invariant or at least do not spill over the regions of resonance and instability.
The problem has been well-studied: a couple of hundred papers and a book [11] have been published on the problem. For an extensive list of papers on this topic, see the reference list of the
book[11]. The existing so-called ‘‘direct methods’’[11–15]can reproduce the given set of measured
data, but cannot guarantee that the remaining eigenvalues and eigenvectors of the FEM remain unchanged. Furthermore, these methods deal with undamped systems only; thus the underlying eigenvalue problem in this setting is a generalized eigenvalue problem in the liner pencil K lM
[2,16] rather than quadratic eigenvalue problem for the pencil (2). The quadratic eigenvalue
problem is much harder to solve numerically [17].
The solution proposed in this paper for EEP can be considered as a partial but meaningful solution to the FEMUP. In contrast with the existing direct methods for FEMUP, the proposed method deals with the damped second-order model and can guarantee mathematically that the eigenvalues and eigenvectors that do not participate in the updating process remain unchanged.
2. Embedding of a real eigenvalue
In this section, we construct the updated matrices Mnew, Knewand Cnew, such that a distinct real
eigenpair ðl1; y1Þ of the pencil F ðlÞ ¼ l2M þ lC þ K is replaced by ðm1; y1Þ, where m1 is
preassigned; m1al1, and the other eigenvalues and eigenvectors remain invariant. To achieve this
goal, we consider a low-rank transformation, called the non-equivalence transformation for the quadratic matrix pencil F ðlÞ. A non-equivalence transformation for the rational l-matrix
functions has been previously considered in Refs. [18–23]. However, the non-equivalence
transformation reported in this paper cannot be derived by using a straightforward generalization of the results in the above papers.
Since ðl1; y1Þis a real eigenpair of F ðlÞ, we have
F ðl1Þy1 ðl21M þ l1C þ KÞy1¼0. (3)
Since K is positive definite, the eigenvector y1 can be normalized such that y>
1Ky1¼1. Suppose
that l12R is a distinct unwanted eigenvalue that needs to be replaced by a prescribed real
number m1. The following theorem provides a non-equivalence transformation of F ðlÞ such that
the updated matrix pencil, FnewðlÞ, keeps the eigenstructure of F ðlÞ except that m1 replaces l1 to
become an eigenvalue of FnewðlÞ.
Theorem 1 (Real eigenvalue embedding). Let ðl1; y1Þ be a distinct real eigenpair of F ðlÞ with
y> 1Ky1¼1, suppose m1al1, and 1 l1m1y1a0 and 1 l21y1a0 and define y1¼y>1My1, (4) e1¼ l1m1 1 l1m1y1 . (5)
Then the updated matrix pencil
FnewðlÞ ¼ l2MnewþlCnewþKnew, (6)
where Mnew¼M e1l1My1y>1M, Cnew¼C þ e1ðMy1y > 1K þ Ky1y > 1MÞ, Knew¼K e1 l1 Ky1y>1K ð7Þ
is symmetric, and has the following spectral properties:
(a) The number m1 is in the spectrum of FnewðlÞ and the remaining eigenvalues of FnewðlÞ are the
same as those of F ðlÞ.
(b) (i) y1 is also an eigenvector of FnewðlÞ corresponding to the eigenvalue m1. (ii) The remaining
eigenvectors of FnewðlÞ are the same as those of F ðlÞ; that is, if l2al1and ðl2; y2Þis an eigenpair
of F ðlÞ, then it is also an eigenpair of FnewðlÞ.
Proof. (a) Substituting the result of Eq. (3) into F ðlÞ, we obtain
F ðlÞy1¼l2My1þlCy1þKy1
¼l2My1þlCy1l21My1l1Cy1
¼ ðl l1Þððl þ l1ÞM þ CÞy1. ð8Þ
By using the identity
detðInþRSÞ ¼ detðImþSRÞ, (9)
where R 2 Cnm and S 2 Cmn, together with Eq. (8), we have
detðFnewðlÞÞ ¼ detðl2MnewþlCnewþKnewÞ
¼ det l2M þ lC þ K l2e1l1My1y>1M þle1ðMy1y>1K þ Ky1y>1MÞ e1 l1 Ky1y> 1K ¼ detðF ðlÞ þ e1ððl þ l1ÞM þ CÞy1y>1ðK ll1MÞÞ ¼ det F ðlÞ þ e1 l l1 F ðlÞy1y>1ðK ll1MÞ ¼ detðF ðlÞÞ 1 þ e1 l l1 ð1 ll1y1Þ ¼ detðF ðlÞÞ l l1 ðl l1þe1ð1 ll1y1ÞÞ.
Since 1 l21y1a0, we now use Eq. (5) to get
l l1þe1ð1 ll1y1Þ ¼ ðl m1Þ
ð1 l21y1Þ
1 l1m1y1
.
Therefore, we conclude that detðFnewðlÞÞ has the same roots as detðF ðlÞÞ, except that l1is replaced
by m1.
(b) We first prove (b)(i). From Eq. (7), we have
Fnewðm1Þy1¼m21ðM 1l1My1y>1MÞy1þm1ðC þ 1
ðMy1y> 1K þ Ky1y>1MÞÞy1þ K 1 l1 Ky1y> 1K y1 ð10Þ ¼ ðm21m211l1y1þm11ÞMy1þm1Cy1þ m11y1þ1 1 l1 Ky1. ð11Þ
Again using Eq. (5), we have
m11y1þ1 1 l1 ¼1 l1m1y11 l1 þ1 ¼m1 l1 . (12)
Since F ðl1Þy1¼0, we have
Ky1¼ l21My1l1Cy1. (13)
Substituting Eqs. (12) and (13) into Eq. (10), we then obtain
Fnewðm1Þy1 ¼ ðm21m211l1y1þm11l1m1ÞMy1.
Once more, from Eq. (5), we conclude that
m21m211l1y1þm11l1m1¼m1ðm1l1Þ þm11ð1 m1l1y1Þ
¼m1ðm1l1Þ þm1ðl1m1Þ
¼0.
This implies that Fnewðm1Þy1¼0, and so (b)(i) is proven.
To prove (b)(ii), we observe that
F ðl2Þy2¼ ðl22M þ l2C þ KÞy2¼0,
that is, Ky2¼ l22My2l2Cy2. This implies
Using the same arguments as in the proof of (a) and Eq. (14), we obtain Fnewðl2Þy2¼ ðl22Mnewþl2CnewþKnewÞy2
¼F ðl2Þy2þ e1 l2l1 ðF ðl2Þy1y>1ðK l2l1MÞy2Þ ¼ e1 l2l1 ðF ðl2Þy1y>1ðl2ððl1þl2ÞM þ CÞÞy2Þ ¼ l2e1 ðl2l1Þ2 ðF ðl2Þy1y>1F ðl1Þy2Þ ¼0.
Hence, ðl2; y2Þ is also an eigenpair of FnewðlÞ. &
Remarks. (i) Note that if l1¼m1, then 1¼0, and there will be no updating at all. Of course, in
practice, it does not make any sense to reassign an eigenvalue which is not desirable to have in the spectrum.
(ii) An alternative and shorter proof of Theorem 1, using orthogonality relations between the eigenvectors of a symmetric positive semi-definite pencil, appear in the Ph.D. Dissertation of
Carvalho [24](available from the website: www.math.niu.edu/~dattab).
3. Embedding of a complex conjugate pair of eigenvalues
We now develop the results in this section, analogous to those of Theorem 1, to show how to
compute the updated symmetric matrices Mnew, Knew and Cnew, such that a distinct complex
conjugate pair of eigenvalues, m1 and ¯m1 is assigned to the spectrum of FnewðlÞ, while the other
eigenvalues of FnewðlÞ and the corresponding eigenvectors remain the same as those of F ðlÞ. We
also give a characterization of the eigenvectors associated with the complex conjugate pair that is reassigned. For simplicity, a matrix pair ðL; Y Þ satisfying
MY L2þCY L þ KY ¼ 0
will be called an eigenpair of F ðlÞ. The notation spec ðTÞ stands for spectrum of the matrix T.
Let ðl1; y1Þbe a complex eigenpair of F ðlÞ, associated with a distinct eigenvalue l1¼a1þib1,
a1, b12R, b1a0, and y1¼y1rþiy1i, y1r, y1i2Rn. Suppose that y1r and y1i are linearly
independent, then y1 and ¯y1 are linearly independent, and ð¯l1; ¯y1Þ is also an eigenpair of F ðlÞ.
Since ðl1; y1Þ is an eigenpair of F ðlÞ, we have
MZ1L21þCZ1L1þKZ1¼0, (15) where L1¼ a1 b1 b1 a1 " # and Z1¼ ½y1r y1i.
Thus, ðL1; Z1Þ is an eigenpair of F ðlÞ. Since K is positive definite, S1¼Z>1KZ1 is also positive
definite. Thus there exists an orthogonal matrix S12R22, and a positive diagonal matrix
D1¼ d1 0 0 d2 " # , such that S1¼S1D1D1S>1.
Therefore, the definitions
Y1¼Z1S1D11 , (16) L1¼D1S>1L1S1D11 , (17) clearly imply Y> 1KY1¼I2, L1¼ d1 0 0 d2 " # a1 b1 b1 a1 " # 1 d1 0 0 1 d2 2 6 6 4 3 7 7 5 ¼ a1 b1=d db1 a1 " # , where d ¼ d1=d2.
To present our main result, we need the following Lemma.
Lemma 2. Given a complex number, m1¼j1þic1, c1a0, there is a real diagonal matrix, EM, such
that m1 is an eigenvalue of the matrix pair
ðL1L>1 EM; L1>EMY1L>1Þ,
where Y1¼Y>1MY1 and Y1, L1 are given by Eqs. (16) and (17), respectively.
Proof. Let Y1¼Y>1MY1¼ y11 y12 y12 y22 " # and EM ¼ x 0 0 Z " # 2R22,
where x, Z are two unknowns. By expanding the following two conjugated equations: det½m1ðL>1 EMY1L>1Þ ðL1L>1 EMÞ ¼0,
det½ ¯m1ðL>1 EMY1L1>Þ ðL1L>1 EMÞ ¼0, ð18Þ
we conclude that x, Z satisfy a system of two real two degree polynomials p1þp2x þ p3Z þ p4xZ ¼ 0,
where p1¼2ðj1r1a1s1Þ, p2¼s1 a1y11þ a1 r1þdb1y12 2j1 r1 ða 2 1þd 2b2 1Þ, p3¼s1 a1 r1þa1y22 b1y12 d 2j1 r1 a 2 1þ b21 d2 , p4¼s1 r1 db1y12 b1y12 d a1y11a1y22 þ2j1 r1 , q1¼s1r1, q2¼ 1 r1ða 2 1þd 2b2 1Þ s1y11, q3¼ 1 r1 a 2 1þ b21 d2 s1y22, q4¼s1ðy11y22y212Þ 1 r1.
Here, yj;k is the ðj; kÞth entry of Y1, j, k ¼ 1; 2; r1¼a21þb21 and s1¼j21þc21. Hence, from Eq.
(19), we can find EM by setting x ¼ q1þq3Z q2þq4Z and Z ¼ ‘2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ‘224‘1‘3 q 2‘1 , (20)
where ‘1¼p3q4p4q3, ‘2 ¼p1q4p2q3þp3q2p4q1 and ‘3¼p1q2p2q1. &
Remark 3.1. (i) It is easily seen from above that x and Z are real provided that ‘224‘1‘3X0. This
will always happen whenever the assumptions of Lemma 2 hold. (ii) Formula (20) usually will give
two possibly solution pairs ðx; ZÞ. The pair ðx; ZÞ that gives the smaller matrix norm kEMkshould
be chosen in a numerical implementation.
The next theorem provides a low-rank transformation of the matrix pencil F ðlÞ, such that the
eigenvalues of the updated symmetric pencil FnewðlÞ are the same as those of F ðlÞ, except for the
complex pair of eigenvalues ðl1; ¯l1Þ of F ðlÞ that is replaced by a prescribed complex pair of
numbers ðm1; ¯m1Þ.
Theorem 3 (Embedding of a pair of complex conjugate eigenvalues). Let Y1 and L1be the same as
those defined in Eqs. (16) and (17). Let EM be the same as in Lemma 2. Define
Mnew¼M MY1EMY>1M,
Cnew¼C þ MY1ECY>1K þ KY1E>CY > 1M,
where
EK ¼L11 EML>1 and EC ¼EML>1 . (22)
Then the real symmetric pencil FnewðlÞ ¼ l2MnewþlCnewþKnew has the following properties:
(i) The eigenvalues of the matrix pencil FnewðlÞ are the same as those of F ðlÞ except that the pair of
complex conjugate eigenvalues l1; ¯l1 of F ðlÞ are replaced by the complex conjugate numbers
m1; ¯m1.
(ii) The eigenvectors associated with the other eigenvalues remain the same as those of the original pencil.
(iii) The eigenvector associated with m1 is given by ¯y1¼Y1X1e1, where X1 is a non-singular matrix
that diagonalizes the matrix f1
¯c1
c1
f1
, and e1is the first unit vector. (Note that m1¼f1þic1.)
Proof. (i) From Eq. (15) and the definitions of Y1and L1, we see that ðL1; Y1Þis an eigenpair of
F ðlÞ and therefore
MY1L21þCY1L1þKY1¼0.
Now, letting L ¼ lI2, we have
F ðlÞY1¼ ðl2M þ lC þ KÞY1
¼ ðMY1ðL þ L1Þ þCY1ÞðL L1Þ. ð23Þ
From Eqs. (21)–(23), we obtain
FnewðlÞ ¼ l2MnewþlCnewþKnew
¼F ðlÞ þ lMY1ECY>1K KY1EKY>1K þ lKY1E>CY > 1M l2MY1EMY>1M ¼F ðlÞ þ ðCY1þMY1ðL þ L1ÞÞL1EKðY>1K lL > 1Y > 1MÞ ¼F ðlÞ þ F ðlÞY1ðL L1Þ1L1EKðY>1K lL>1Y>1MÞ. This implies
det½FnewðlÞ ¼ det½F ðlÞ þ F ðlÞY1ðL L1Þ1L1EKðY>1K lL > 1Y
> 1MÞ
¼ det½F ðlÞ det½InþY1ðL L1Þ1L1EKðY>1K lL>1Y>1MÞ
¼ det½F ðlÞ det½I2þ ðL L1Þ1L1EKðI2lL>1Y1Þ
¼ det½F ðlÞ ðl l1Þðl ¯l1Þ det½ðlI2L1Þ þL1EKðI2lL>1Y1Þ ¼ det½F ðlÞ ðl l1Þðl ¯l1Þ det½lðI2L1EKL>1Y1Þ L1ðI2EKÞ.
By Lemma 3.1, the pair fm1; ¯m1gis a complex conjugate pair of the eigenvalues of the matrix pencil ðL1L>1 EM; L1>EMY1L>1Þ. This implies that
det½lðI2L1EKL>1Y1Þ L1ðI2EKÞ ¼
ðl m1Þðl ¯m1Þ
l1¯l1
.
Therefore, FnewðlÞ has the same eigenvalues as F ðlÞ, except that l1, ¯l1are now replaced by m1, ¯m1.
(ii) Let l2¼a2þib2and y2¼y2rþiy2i. Define Y2and L2in the same way as Y1 and L1 have
been defined in Eqs. (15)–(17). Then ðL1; Y1Þand ðL2; Y2Þare eigenpairs of F ðlÞ, with Y>1KY1¼
I2 and Y>2KY2¼I2. Thus Y> 2KY1þY > 2CY1L1þY > 2MY1L 2 1¼0, (24) Y>2KY1þL>2Y > 2CY1þ ðL>1Þ 2Y> 2MY1¼0. (25)
Eliminating the terms involving ‘‘Y>
2CY1’’ from Eqs. (24) and (25), we have
L>
2ðY>2KY1Þ ðY>2KY1ÞL1þL>2ðY2>MY1ÞL21 ðL>2Þ2ðY>2MY1ÞL1 ¼0.
Let KY ¼Y>2KY1; MY ¼Y>2MY1. Let and vecðÞ denote the Kronecker product and
vectorizing operator, respectively. Then vectorizing the last equation, we have ðI L>
2 L>1 I ÞvecðKYÞ ¼vecðL>2ðL>2MY MYL1ÞL1Þ
¼ ðL>1 L>2ÞvecðL2>MY MYL1Þ
¼ ðL>1 L>2ÞðI L2>L>1 IÞvecðMYÞ
¼ ðI L>2 L>1 IÞðL>1 L>2ÞvecðMYÞ.
Suppose l1al2, then specðL1Þ \specðL2Þ ¼ ;. This implies that the matrix ðI L>2 L
> 1 I Þ is
non-singular and hence, ðL>1 L>2ÞðvecðMYÞÞ ¼vecðKYÞ. Thus,
Y>2KY1¼L>2ðY >
2MY1ÞL1. (26)
Since ðL2; Y2Þ is an eigenpair of F ðlÞ, we have
MY2L22þCY2L2þKY2¼0.
From Eq. (21) it then follows that MnewY2L22þCnewY2L2þKnewY2
¼ ðM MY1EMY>1MÞY2L22þ ðC þ MY1ECY>1K þ KY1E>CY>1MÞY2L2 þ ðK KY1EKY>1KÞY2 ¼MY2L22MY1EMY>1MY2L22þCY2L2þMY1ECYT1KY2L2 þKY1E>CY > 1MY2L2þKY2KY1EKYT1KY2 ¼ MY1L1EKL>1Y > 1MY2L22þMY1L1EKY>1KY2L2þKY1EkL>1Y > 1 MY2L2KY1EKY>1KY2 ¼MY1L1EKðY1>KY2L2L>1Y > 1MY2L22Þ þKY1EKðL>1Y > 1MY2L2Y>1KY2Þ.
(iii) Let O1¼ j1
c1
c1 j1
h i
, where m1¼j1þic1 is a complex eigenvalue of FnewðlÞ, with c1a0.
From Eq. (18), there exists a non-singular matrix V12R22 such that
ðI2L1EKL>1Y1ÞV1O1L1ðI2EKÞV1¼0.
By setting O1¼V1O1V11 , we obtain
ðI2L1EKL>1Y1ÞO1L1ðI2EKÞ ¼0. (27)
Now,
MnewY1O21þCnewY1O1þKnewY1¼MY1ðO21EMY1O12þECO1Þ þCY1O1
þKY1ðE>CY1O1þI2EKÞ, ð28Þ
where Y1¼Y>1MY1. Since ðL1; Y1Þis an eigenpair of F ðlÞ, and O1satisfies Eq. (27), we conclude that
CY1O1þKY1ðE>CY1O1þI2EKÞ
¼CY1O1þ ðMY1L21CY1L1ÞðEKL>1Y1O1þI2EKÞ
¼ MY1L21ðEKL>1Y1O1þI2EKÞ þCY1½ðI2L1EKL>1Y1ÞO1L1ðI2EKÞ
¼ MY1L21ðEKL>1Y1O1þI2EKÞ.
Therefore, Eq. (28) becomes
MY1½O21EMY1O21þECO1L21ðEKL>1Y1O1Þ L21ðI2EKÞ
¼MY1½ððI2EMY1ÞO1L1ðI2EKÞÞO1þL1ððI2EMY1Þ
O1L1ðI2EKÞÞ
¼0.
Thus, ðY1; O1Þis an eigenpair of FnewðlÞ. Letting T1O1T11 ¼ ðm01 0
¯m1Þand setting X ¼ T1V1, the (iii) is
proved. &
Based on the above theorem, we present the following algorithm for assigning a pair of complex conjugate numbers to be eigenvalues of the updated symmetric matrix pencil.
Algorithm 3.1 (Assignment of a Pair of Complex Conjugate Eigenvalues). Input:
(i) An unwanted distinct complex eigenvalue, l1¼a1þib1; a1; b12R; b1a0 (and its complex
conjugate), and the corresponding eigenvector, y1¼y1rþiy1i, y1r; y1i2Rn, with y1r, y1i being
linearly independent.
(ii) A pair of complex conjugate numbers, m1 and ¯m1, that needs to be embedded.
(iii) Symmetric matrices, M, C and K, with M, and K positive definite.
Output: Symmetric matrices Mnew, Cnew and Knew such that the updated pencil FnewðlÞ ¼
l2MnewþCnewl þ Knew has the eigenpair fm1; ¯m1g in its spectrum, the remaining eigenvalues and
eigenvectors are the same, and the eigenvector associated with m1 is given by Y1X1e1; where Y1
Step 1: Use Eqs. (16) and (17) to find the eigenpair ðL1; Y1Þ of the original matrix pencil
F ðlÞ ¼ l2M þ lC þ K such that specðL1Þ ¼ fl1; ¯l1g and Y>1KY1¼I2.
Step 2: Determine x and Z by using formula (20). If x or Z is complex then stop and return to Step 1.
Step 3: Set EM ¼ x0 0Z
h i
, EK ¼L11 EML1> and EC ¼EML>1 .
Step 4: Computed the updated matrices
Mnew¼M MY1EMY>1M,
Cnew¼C þ MY1ECY>1K þ KY1E>CY > 1M,
Knew¼K KY1EKY>1K.
Remark 3.2. Above, we have discussed how to replace an unwanted complex conjugate pair
fl1; ¯l1g by a prescribed conjugate pair ðm1; ¯m1Þ, assuming that the associated eigenvector y1¼
y1rþiy1i is such that y1r and y1i are linearly independent.
We now consider the degenerate case where the real and the imaginary parts of the eigenvector,
y1r and y1i are linearly dependent. In this case, the eigenvectors corresponding to l1 and ¯l1, are
also linearly dependent. Hence, the eigenvector y1 can be a real vector, i.e., y12Rn. Since both
ðl1; y1Þand ð¯l1; y1Þ are eigenpairs of F ðlÞ, we have
l21My1þl1Cy1þKy1¼0,
¯l21My1þ ¯l1Cy1þKy1¼0.
Then, we obtain ðl1þ ¯l1ÞMy1þCy1¼0. This implies that Cy1==My1, and thus, Ky1==My1. Let
Q 2 Rnn be an orthogonal matrix such that Q>y1¼e1. Let
e
M ¼ Q>MQ; C ¼ Qe >CQ; K ¼ Qe >KQ,
then the first columns of eM, eC, eK are mutually parallel. Furthermore, since eM, eC, eK are
symmetric, the first row vectors of eM, eC, eK are also mutually parallel. Hence, if we apply an
elementary matrix L to eliminate the second through the nth elements of the column of eM (see
Ref.[16]) then the first columns and rows of the matrices L eML>, L eCL>, L eKL>are parallel to e 1.
Hence, the dimension of the quadratic problem in this case can be reduced to n 1 by removing
the first row and column of matrices L eML>, L eCL>, L eKL>simultaneously. Thus, the unwanted
eigenvalues l1 and ¯l1 are deflated simultaneously, reducing the dimension of the problem by 1.
Algorithm 3.1 now can be applied to the reduced problem.
Remark 3.3. In case l1¼a1þib12C is a multiple eigenvalue, the formula (21) can still be used
to update the matrices M; C and K; however, in this case to construct Y1 we must consider not
only the eigenvector y1, but the associated generalized eigenvector as well. Thus, if l1 is an
eigenvalue with multiplicity 2, then Y1 is computed by normalizing ½y1r; y1i; z1r; z1i with
Y>
1KY1¼I4, and y1r, z1rðy1i; z1iÞare, respectively, the real (and imaginary) parts of y1, z1, and EM
is diagonal and is yet to be determined, EC ¼EMLe
> 1 , EK ¼ eL 1 1 EMLe > 1 . In addition, eL12R44 is similar to L1 0 I L1 h i , and L1¼ ba1 1 b1 a1 h i
system of equations:
detðFnewðm1ÞÞ ¼0; detðFnewð¯m1ÞÞ ¼0,
d
dl½detðFnewðm1ÞÞ ¼0; d
dl½detðFnewð¯m1ÞÞ ¼0.
An Illustrative Example. Consider application of Algorithm 3.1 to a free beam with
I ¼ Moment of inertia ¼ 1:136 109m4,
E ¼ Young’s modulus ¼ 72 Gpa, l ¼ Length of the beam ¼ 0:4005 m. The stiffness matrix has the form
K ¼EI l3 12 6l 12 6l 6l 4l2 6l 2l2 12 6l 12 6l 6l 2l2 6l 4l2 0 B B B @ 1 C C C A. With the above values of I ; E, and l, we have
K ¼ 104 1:5330 0:3066 1:5330 0:3066 0:3066 0:0818 0:3066 0:0409 1:5330 0:3066 1:5330 0:3066 0:3066 0:0409 0:3066 0:0818 0 B B B @ 1 C C C A, M ¼ 0:1349 0:0076 0:0467 0:0045 0:0076 0:0006 0:0045 0:0004 0:0467 0:0045 0:1349 0:0076 0:0045 0:0004 0:0076 0:0006 0 B B B @ 1 C C C A, D ¼ 0.
The eigenvalues of F ðlÞ are: 103ð5:4363i; 1:5916i; 0; 0; 0; 0Þ. The pair of the complex
eigenvalues, 103ð1:5916iÞ were charged to 103ð1:3509iÞ, obtained from an experiment at the
vibration laboratory at Northern Illinois University. The updated stiffness matrix is given by
Knew¼104 1:5330 0:3066 1:5330 0:3066 0:3066 0:0787 0:3066 0:0440 1:5330 0:3066 1:5330 0:3066 0:3066 0:0440 0:3066 0:0787 0 B B B @ 1 C C C A.
The entries of the updated mass matrix Mneware almost the same as those of the original matrix
and the entries of the matrix Dneware of Oð1014Þ. The results on both 1 and 10 elements of the
4. Error analysis for the assignment of a complex conjugate pair of eigenvalues
In this section k k denotes the 2-norm,^ (hat) denotes a computed quantity and the term HOT
stands for ‘‘the higher-order terms.’’
First, we estimate the error bounds for the computed Mnew. From Eq. (21), we have
k bMnewMnewk ¼ kM bY1EbMYb > 1M MY1EMY>1Mk pkMk2k bY 1EbMYb > 1 Y1EMY > 1k. ð29Þ 1 2 3 4 0 0.5 1 1.5 2 2.5 3 3.5 4 Degrees of Freedom
Percentage Change in Stiffness
along diagonal
1 Beam Element
Fig. 2. Percent change in the diagonal entries of the stiffness matrix for one beam element.
0 5 10 15 20 25 0 5 10 15 20 25 30 35 40 Degrees of Freedom
Percentage Change in Stiffness
along diagonal
10 Beam Element
By using the triangular inequality, we obtain k bY1EbMYb > 1 Y1EMY>1k pk bY1EbMYb > 1 bY1EbMY>1k þ k bY1EbMY>1 bY1EMY>1k þ k bY1EMY>1 Y1EMY>1k þ k bY1Y1k kEMY>1k pk bY1EbMkk bY1 bY1Y1k þ k bY1kk bY kk bEM EMk þ k bY1Y1kkEMY>1k p½ðk bY1Y1k þ kY1kÞk bEMEMk þ k bY1Y1k kEMk þ kY1EMkk bY1Y1k þ ðk bY1Y1k þ kY1kÞkY1k k bEM EMk þ k bY1Y1k kEMY>1k. ð30Þ 5 10 15 20 0 5 10 15 20 25 -4 -2 0 2 4 6 × 106 Degrees of Freedom 10 Beam Element Degrees of Freedom Change in Stiffness
Fig. 4. Change in the entries of the stiffness matrix for 10 beam element.
1 2 3 4 1 2 3 4 -40 -20 0 20 40 Degrees of Freedom 1 Beam Element Degrees of Freedom Change in Stiffness
From the definition of Y1 in Eq. (16), we then have k bY1Y1k ¼ k bZ1Sb1Db 1 1 Z1S1D 1 1 k ð31Þ p½ðk bZ1Z1k þ kZ1kÞk bS1S1k þ k bZ1Z1k kS1k þ kZ1S1kk bD 1 1 D11 k þ ðk bZ1Z1k þ kZ1kÞkD11 k k bS1S1k þ k bZ1Z1k kS1D11 k. ð32Þ
It is known (see Ref. [17]) that the error bound for Z1 satisfies
k bZ1Z1kpc1e, (33) where c1¼ X2n k¼2 kzkk jlkl1jð1 þ jlkl1jÞjzHky1j ,
zk and yk are, respectively, the left and right eigenvectors corresponding to the eigenvalue lk of
F ðlÞ, Zk¼ ½ykryki. Similarly,
k bS1S1kpc2e, (34)
where c2 is a constant. Since D1¼ ½d01 d0
2 2R
22, d
14d2, and S12R22 is orthogonal, we have
kD11 k ¼ kS1D11 k ¼ 1 d2 (35) and k bD11 D11 k ¼ kD11 ðD1 bD1Þ bD 1 1 k pkD1 1 k2kD1 bD1k þHOT ¼ d1 d22e þ HOT. ð36Þ
From Eqs. (33)–(36) and (31), we then obtain
k bY1Y1kpc3e þ HOT, (37) where c3¼ d1 d22þ c2 d2 ! kZ1k þ c1 d2 . Since EM ¼ ½x0 0Z, we have k bEM EMk ¼maxfjbx xj; jbZ Zjg. (38)
We now obtain the bounds for j^z zj, and j^Z Zj. From Eq. (20) and relations (17)–(19), we know that
x ¼ xðl1; m1Þ ¼xða1; b1; j1; c1Þ,
where l1¼a1þib1 and m1¼j1þic1. In addition, we have bx ¼ bxða1; b1; j1; c1Þ ¼xðba1; bb1;jb1; bc1Þ þHOT, bZ ¼ bZða1; b1; j1; c1Þ ¼Zðba1; bb1;jb1; bc1Þ þHOT. So, bx x ¼ qx qa1 Da þ qx qb1Db þ qx qj1 Dj þ qx qc1 Dc þ HOT, (39) bZ Z ¼qaqZ 1 Da þ qZ qb1Db þ qZ qj1 Dj þ qZ qc1 Dc þ HOT, (40)
where Da ¼ba1a1, Db ¼ bb1b1, Dj ¼jb1j1, and Dc ¼ bc1c1. Since m1¼j þ ic is a
prescribed number, we need not calculate it. The numbers Dj and Dc are usually much smaller than Da or Db. We can hence ignore terms involving Dj or Dc in Eqs. (39) and (40). Hence, we are only concerned with those terms related to Da or Db in the estimation of the error bounds for x and Z. From Eqs. (39) and (40), we have
jbx xjp qx qa1 jDaj þ qx qb1 jDbj þ HOT p qaqx 1 þ qx qb1 jbl1l1j þHOT, jbZ Zjp qZ qa1 jDaj þ qZ qb1 jDbj þ HOT p qZ qa1 þ qZ qb1 jbl1l1j þHOT.
After performing some tedious calculations, it can be shown that jqx=qa1j, jqx=qb1j, jqZ=qa1jand
jqZ=qb1jare bounded by the relative rational functions in a1, b1and jl1j. More precisely, one can
prove that jbx xjpjz1ða1; b1Þj z2ðjl1jÞ jbl1l1j þHOT, (41) jbZ ZjpjB1ða1; b1Þj B2ðjl1jÞ jbl1l1j þHOT, (42)
where z1, B1 are low degree polynomials in a1, b1, and z2, B2 are low-degree polynomials in jl1j.
Since jl1ja0, z2 and B2 are non-zero, and both bounds in Eqs. (41) and (42) are finite. Again,
jbl1l1jp
1
ð1 þ jl1j2ÞjzH1y1j
e. (43)
Substituting Eqs. (41)–(43) into Eq. (38), we have
where c4¼ zða1; b1Þ Bðjl1jÞ 1 ð1 þ jl1j2ÞjzH1y1j ,
zða1; b1Þ ¼maxfjz1ða1; b1Þj; jB1ða1; b1Þjg,
Bðjl1jÞ ¼ minfz1ðjl1jÞ; B1ðjl1jÞg.
By using Eqs. (30), (37) and (44), we then obtain k bY1EbMYb
>
1 Y1EMY>1kp½2kY1k kEMkc3þ kY1k2c4 e þ HOT. (45)
Using Eq. (45) in Eq. (29), we finally obtain the following error bound for Mnew:
k bMnewMnewkpkMk2½2kY1k kEMkc3þ kY1k2c4. (46)
To estimate the error bounds for Cnew and Knew, we first need to find the error bound for L11 .
From Eq. (17), we have
kbL11 L11 k ¼ k bD1Lb 1 1 Db 1 1 D1L11 D11 kpc5e þ HOT, where c5¼ d1 d2 d1 d2jl1j þ 1 jl1j þ 1 ð1 þ jl1j2ÞjzH1y1j ! .
Hence, by a similar process as above, we obtain the error bounds for Cnew and Knew as given
below
k bCnewCnewkp2ekMk kKk kEMk
d1c3 d2jl1j þc5 þ c4 jl1j , (47)
k bKnewKnewkpekKk2 2kY1k kEMk
d21c3 d22jl1j2 þ kY1k2 2d1 d2jl1j kEMkc5þ d21c4 d22jl1j2 ! " # . (48)
5. Simultaneous assignment of several real eigenvalues
So far, we have considered the problem of assigning either one real or a pair of complex conjugate eigenvalues. In this section, we consider the simultaneous assignment of several real eigenvalues.
It is always possible to embed the sequence of real eigenvalues, fm1; . . . ; mmrg in the updated
symmetric matrix pencil, FnewðlÞ, by using the formula (7) recursively, for s ¼ 1;. . . ; mr
Ms ¼Ms1eslsMs1ysy>s Ms1, Cs¼Cs1þesðMs1ysys>Ks1þKs1ysy>s Ms1Þ, Ks ¼Ks1 es ls Ks1ysy>s Ks1, ð49Þ
where M0 ¼M, C0¼C and K0¼K, and ys and es are given by
ys ¼y>sMs1ys and es¼
lsms
1 lsmsys
. (50)
By doing so, the eigenvalues will be embedded one at a time. However, it is possible to assign several of them at a time as long as the mass and stiffness matrices remain positive definite.
The method proposed below delays the updating of the coefficient matrices until all the real
numbers, fysg and fesg, needed for the multi-assignment, have been computed. After all these
quantities have been computed, the coefficient matrices are updated with only one rank-mr
symmetric update. The process will not only be more efficient than that which assigns one eigenvalue at a time, but it will be rich in Basic Linear Algebra Subroutines Level 3 (BLAS-3), such as matrix–matrix multiplications, rank-r updates, etc., which will make it suitable for high-performance computing.
Given r real numbers, fm1; . . . ; mrg, the following method computes a positive integer, mrpr, the
matrices W and U, and the diagonal matrices DM, DC and DK, such that the updated symmetric
matrix pencil, FnewðlÞ ¼ l2MnewþlCnewþKnew, has the spectrum
specðFnewðlÞÞ ¼ fm1; . . . ; mmr; lmrþ1; . . . ; l2ng,
where
Mnew¼M WDMW>, (51)
Cnew¼C þ UDCW>þWDCU>, (52)
Knew¼K UDKU>. (53)
To develop formula (51), we consider the mrth iteration of Eq. (49) and observe that
Mmr ¼M0 Xmr s¼1 eslsMs1ysy>s Ms1 ¼M0 ½M0y1; . . . ; Mmr1ymr e1l1 . . . emrlmr 2 6 6 6 4 3 7 7 7 5½M0y1; . . . ; Mmr1ymr >. ð54Þ
We also observe that, for s ¼ 1;. . . ; mr,
ys¼y>s Ms1ys ¼y>s ½Ms2es1ls1Ms2ys1y>s1Ms2ys ¼y>s Ms2yses1ls1ðy>sMs2ys1Þðy>s1Ms2ysÞ ð55Þ and Ms1ys¼ ½Ms2es1ls1Ms2ys1y>s1Ms2ys ¼Ms2yses1ls1ðy>s1Ms2ysÞMs2ys1. ð56Þ
Therefore, formula (51) can be derived from Eq. (54) by letting W ¼ ½M0y1; . . . ; Mmr1ymr and DM ¼ e1l1 . . . emrlmr 2 6 6 4 3 7 7 5.
In addition, the matrices DM and W can be determined by using recursions (55) and (56).
Similarly, formulae (52) and (53) can be obtained for the appropriate matrices U, DK and DC.
Our discussions above are summarized in the algorithms below. Algorithm 5.1 (Simultaneous Assignment of Real Eigenvalues). Input:
(i) A set of real numbers fmigri¼1,
(ii) A set of unwanted real eigenpairs fðli; yiÞgri¼1,
(iii) Symmetric matrices M, C and K such that M, and K are positive definite.
Output: Integer mr, and the symmetric matrices Mnew, Cnew and Knew such that the updated
quadratic matrix pencil FnewðlÞ contains the mr eigenvalues in the spectrum ðmrprÞ while the
other eigenvalues and the associated eigenvectors remain unchanged.
Step 1: Compute mi ¼Myi, ki¼Kyi, i ¼ 1;. . . ; r. Step 2: Compute aij ¼y>i mj; bij ¼y>i kj; j ¼ i; . . . ; r; i ¼ 1; . . . ; r. Step 3: Set Z1¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi y> 1Ky1 q Update a1j ¼a1j=Z1; b1j¼b1j=Z1; j ¼ 1; . . . ; r. Step 4: Set e1¼1ll1m1m1 1a11. Step 5: For s ¼ 2; . . . ; r. For i ¼ s;. . . ; r. For j ¼ i;. . . ; r.
Update aij ¼aijes1ls1as1;ias1;j; bij ¼bijles1s1bs1;ibs1;j.
End for j. End for i. Compute es ¼1llsms smsass. If bss40,then Compute Zs ¼pffiffiffiffiffiffibss.
Update ai;s¼ai;s=Zs; bi;s¼bi;s=Zs; i ¼ 1; . . . ; s.
Update as;j ¼as;j=Zs; bs;j ¼bs;j=Zs; j ¼ 1; . . . ; s.
Compute mr¼s.
Else Exit Loop. End for s.
Step 7: For s ¼ 2; . . . ; mr,
For i ¼ s;. . . ; mr,
Update mi ¼mies1ls1as1;ims1; ki¼kiles1s1bs1;iks1.
End for i. End for s.
Step 8: Set W ¼ ½m1; m2; . . . ; mmr; U ¼ ½k1; k2; . . . ; kmr; DM ¼diagðe1l1; e2l2; . . . ; emrlmrÞ,
DC ¼diagðl1; l2; . . . ; lmrÞ and DK ¼diag
e1 l1 ;e2 l2 ; . . . ;emr lmr . Step 9: Update Mnew¼M WDMW>, Cnew¼C þ UDCW>þWDCU>, Knew¼K UDKU>. Return
To show the efficiency of the simultaneous assignment process, we compare the flop counts of Algorithm 5.1, with those of the successive assignment strategy by using non-equivalence
transformation (7). InTable 1, we list the flop counts of these two methods.
From Table 1, we see that the simultaneous assignment method is more efficient than the successive assignment procedure.
6. Numerical results
In this section, we illustrate the efficiency and reliability of the proposed method by using two
examples: The first one is taken from Harwell–Boeing Collections[25]. The data of the second is
the simulated data of a real-life aerospace example provided to us by the Boeing company. All numerical implementations were performed on a IBM Pentium III machine using MATLAB. 6.1. Example 1 (Updating of a statistically condensed oil rig model)
Consider the model ðM; D; KÞ where
The matrices M 2 R6666 and K 2 R6666 come from the statically condensed oil rig model ofthe Harwell–Boeing set BCSSTRUC1[25]. The matrix M is symmetric positive definite and the
matrix K is symmetric positive semi-definite.
The damping matrix C is defined by C ¼ rI66, with r ¼ 1:55.Table 1
Approximate flop counts for embedding r ðr5nÞ real eigenvalues
Strategy Parameters Mnew, Cnew, Knew Total
r sequential assignment 6n2r 7n2r
2
19n2r
2
This model has 132 eigenvalues out of which eight are real eigenvalues fl1; . . . ; l8g, given by
fl1 l2 l3 l4g ¼ f3:4628 3:5709 5:3584 9:2761g,
fl5 l6 l7 l8g ¼ f13:1972 13:4480 27:5536 44:5031g
and 62 pairs of complex conjugate eigenvalues that are not shown here. The set fl1; . . . ; l8g is
changed to the set fm1; . . . ; m8g, where
fm1 m2 m3 m4g ¼ f3:32 3:75 5:05 9:07g,
fm5 m6 m7 m8g ¼ f13:59 13:04 27:31 42:11g.
Algorithm 5.1 is then applied, giving matrices DM, DC, and DK as follows:
DM ¼diagð0:6697 0:9138 3:6368 2:4231 2:6340 2:6111 17:3927 197:1462Þ,
DC ¼diagð0:1934 0:2559 0:6787 0:2612 0:1996 0:1942 0:6312 4:4299Þ,
DK ¼diagð0:0558 0:0717 0:1267 0:0282 0:0151 0:0144 0:0229 0:0995Þ.
The matrices W and U are not shown here. The matrices Mnew; Cnewand Kneware then computed,
using the update formulas, as a single rank-8 update of the matrices M; D, and K. Verification: Define L ¼ diagðl1; . . . ; l132Þ, ~ L ¼ diagðm1; . . . ; m8; l9; . . . l132Þ, Y ¼ ½y1::: y132 then kMnewY ~L 2 þDnewY ~L þ KnewX kF ¼1:7709 107,
which shows that the multiple embedding was successful and produced no spill-over.
Fig. 5shows the bar graphs of the magnitude of the components of the matrix K Knew. Similar
graphs exist for the matrices M Mnew and D Dnew.
6.2. Example 2
The Boeing Simulated Example. The test matrices K, C, M in this example come from an aerospace industry problem in constructing aircraft structural models.
Ten complex conjugate pairs of eigenvalues, which seem to be ‘‘troublesome’’, need to be embedded in the given model. This is done by applying Algorithm 3.1 ten times, assigning one pair at a time. The results of implementation are plotted in the figure below. To understand the error behaviors more clearly, both the absolute and the logarithms of the error matrices have been
of the matrix M, denoted by log M, is defined by
log Mði; jÞ ¼ log10jMnewði; jÞ Mði; jÞj if jMnewði; jÞ Mði; jÞj410
4;
0 otherwise:
(
The results clearly show that our updating with low-rank transformations is successful. Furthermore, the results of the type obtained here provide an insight for the practicing engineers into what rows of the mass, stiffness or damping matrices need modification. For this particular example, our plots show that the largest errors occur around 3rd and 37th rows and columns in all these matrices. These rows and columns, therefore, need most modifications for the application under considerations.
0 10 20 30 40 50 0 5 10 15 20 25 30 35 40 45 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 × 10 9 abs (Knew−K)
Fig. 6. The absolute error of jKnewKj. 1.5
1
0.5
0
Update in stiffness matrix 0
10 20 30 40 50 60 70 10 20 30 40 50 60
7. Conclusion
The symmetric eigenvalue embedding problem addressed in this paper is the one of updating a symmetric finite element generated second-order model in such a way that the updated model remains symmetric, and a small subset of unwanted eigenvalues is replaced by a suitably user-chosen set, while the remaining large number of eigenvalues and eigenvectors do not change. The problem is intimately related to the partial eigenvalue assignment problem in control theory, which is usually solved by using feedback control. Unfortunately, with the use of feedback control, the symmetry of the model is completely destroyed. A novel symmetry preserving algorithm and the associated theories are presented in this paper. The proposed method results in a symmetric low-rank transformation of the original model, with the required properties. The method allows simultaneous assignment of several real eigenvalues; however, complex eigenvalues have to be assigned one at a time. Further research on simultaneous assignment of more than one complex eigenvalues is currently underaway. The results of the paper contribute to the progress in the solution of a well-known problem of immense practical importance in vibration industries: namely, the finite-element model updating problem, which is concerned with updating a symmetric finite-element model such that the updated model is symmetric, a small number of measured eigenvalues and eigenvectors from a practical structure is incorporated into the model, and the remaining large number of eigenvalues and eigenvectors that do not participate in the updating process remain invariant. Furthermore, because the proposed algorithms are rich in Basic Linear Algebra Subroutine-3 (BLAS-3) level operations, they can be implemented using high-performance software packages such as LAPACK on today’s high-speed computers.
Acknowledgements
This work was supported in part by the National Center for Theoretical Sciences, Hsinchu, Taiwan. J.B. Carvalho was also partially supported by Brazilian CAPES and USA-NSF Grant ECS-0074411 and B.N. Datta was partially supported by USA-NSF Grant ECS-0074411 and by US Department of Education Grant #P116Z040180 to the Vibration and Acoustic Center of the College of Engineering and Engineering Technology of Northern Illinois University. The authors would like to thank the reviewers for their constructive suggestions which led to this improved version of the paper. They would also like to thank M. Lagadapati, a graduate student of Professor Datta for his help with running experiments on the proposed algorithm using data from NIU’s vibration and Laboratory.
References
[1] D. Inman, Vibration with Control, Measurement and Stability, Prentice-Hall, Englewood Cliffs, NJ, 1989. [2] B. Datta, Numerical Methods for Linear Control Systems Design and Analysis, Elsevier Academic Press, Boston,
2003.
[3] Y. Saad, Projection and deflation methods for partial pole assignment in linear state feedback, IEEE Transactions on Automatic Control 33 (1988) 290–297.
[4] B. Datta, D. Sarkissian, Partial eigenvalue assignment in linear systems, in: Proceedings of 15th International Symposium on Mathematical Theory of Networks and Systems, University of Notre Dame, Notre Dame, USA, August 12–16, 2002.
[5] B. Datta, S. Elhay, Y.M. Ram, Orthogonality and partial pole assignment for the symmetric definite quadratic pencil, Linear Algebra and its Applications 257 (1997) 29–48.
[6] B. Datta, D. Sarkissian, Multi-input partial eigenvalue assignment for the symmetric quadratic pencil, in: Proceedings of the American Control Conference, San Diego, California, June 2–4, 1999, pp. 2244–2247.
[7] B. Datta, S. Elhay, Y.M. Ram, D.R. Sarkissian, Partial eigenstructure assignment for the quadratic pencil, Journal of Sound and Vibration 230 (2000) 101–110.
[8] B. Datta, D. Sarkissian, Theory and computations of some inverse eigenvalue problems for the quadratic pencil, Contemporary Mathematics 280 (2001) 221–240.
[9] Y.M. Ram, S. Elhay, Pole assignment in vibrating systems by multi-input control, Journal of Sound and Vibration 230 (2001) 101–119.
[10] Y.M. Ram, Pole-zero assignment of vibratory systems by state feedback control, Journal of Vibration and Control 4 (1998) 165–185.
[11] M. Friswell, J. Mottershed, Finite Element Model Updating in Structural Dynamics, Kluwer Academic Publishing, London, 1995.
[12] B. Datta, Finite element model updating, eigenstructure assignment and eigenvalue embedding techniques for vibrating systems, Mechanical Systems and Signal Processing 16 (2002) 83–96.
[13] C. Minas, D. Inman, Correcting finite element models with measured modal results using eigenstructure assignment methods, in: Proceedings of the Sixth International Modal Analysis Conference, Orlando, Florida, February 1988, pp. 583–587.
[14] C. Minas, D. Inman, Matching finite element models to model data, Transactions of the Association of Mechanical Engineering 112 (1990) 84–92.
[15] D. Zimmerman, M. Widengren, Correcting finite element models using a symmetric eigenstructure assignment technique, American Institute of Aeronautics and Astronautics Journal 28 (1990) 1670–1676.
[16] B. Datta, Numerical Linear Algebra and Applications, Brooks/Cole Publishing Company, Pacific Grove, CA, 1995. [17] F. Tisseur, K. Meerbergen, The quadratic eigenvalue problem, SIAM (Society for Industrial and Applied
Mathematics) Review 43 (2001) 253–286.
[18] P. Dewilde, P. Vandewalle, On the factorization of a nonsingular rational matrix, IEEE Transactions Circuits and Systems 22 (1975) 637–645.
[19] W.R. Ferng, W.-W. Lin, D.J. Pierce, C.-S. Wang, Nonequivalence transformation of l-matrix eigenproblems and model embedding approach to model tuning, Numerical Linear Algebra with Applications 8 (2001) 53–70. [20] J.-S. Guo, W.-W. Lin, C.-S. Wang, Numerical solution for large sparse quadratic eigenvalue problem, Linear
Algebra and its Applications 225 (1995) 57–89.
[21] J.-S. Guo, W.-W. Lin, C.-S. Wang, Nonequivalent deflation for the solution of matrix latent value problems, Linear Algebra and its Applications 231 (1995) 15–45.
[22] W.-W. Lin, An algorithm for computing the eigenstructure of a regular matrix polynomial, Linear Algebra and its Applications 91 (1990) 195–211.
[23] Van Dooren, P. Dewilde, The Eigenstructure of an arbitrary polynomial matrix: computational aspects, Linear Algebra and its Applications 50 (1983) 545–579.
[24] J. Carvalho, State Estimation and Finite Element Model Updating for Vibrating Systems, PhD Dissertation, Northern Illinois University, DeKalb, Illinois, USA, 2002.
