王平心,吳頡爾,楊習貝

(1.江蘇科技大學 理學院, 鎮(zhèn)江 212003) (2.江蘇科技大學 計算機學院, 鎮(zhèn)江 212003) (3.河北師范大學 數(shù)學與信息科學理學院, 石家莊 050024)

線性阻尼系統(tǒng)的自由振動方程可以表示為：

(1)

式中：M,C,K分別為質量、阻尼和剛度矩陣；u(t)為位移向量．令u(t)=ueλt(其中u是不依賴時間的常向量),代入上式則可得到如下二次特征值問題：

(λ2M+λC+K)u=0

(2)

當M,C,K均為實對稱矩陣,則稱問題(2)為對稱二次特征值問題．如果質量、阻尼和剛度矩陣依賴于設計參量p1,p2,…,pN,則得如下二次特征值問題：

[λ2(p)M(p)+λ(p)C(p)+K(p)]u(p)=0

(3)

式中：p=(p1,…,pN)T；M(p),C(p),K(p)是在p=p0的鄰域內解析的矩陣值函數(shù)．

阻尼系統(tǒng)的結構動力特性主要表現(xiàn)在其特征對上,即特征值和特征向量．因此,特征對導數(shù)的計算是理解并確定參數(shù)變化對系統(tǒng)影響必不可少的方法．特征靈敏度分析的主要任務之一就是計算特征對的導數(shù),其結果在結構故障診斷[1]、結構優(yōu)化設計[2]、結構分析與識別[3-4]、結構模型修正[5]和代數(shù)特征值反問題[6-7]等領域都有重要的應用．在特征對導數(shù)的計算中,特征值導數(shù)的計算比較簡單,已有一些有效的方法．特征向量導數(shù)的計算需要求解系數(shù)矩陣奇異的線性方程組,是特征對導數(shù)計算中的主要難點．為了解決這一問題,許多學者在這一領域進行了大量研究,提出了多種特征對導數(shù)的計算方法[8-18]．但是,目前存在的算法大多針對特征值是單根或者特征值是重根，但特征值的導數(shù)是互異的情形．文中利用矩陣廣義逆理論,推導一種計算對稱二次特征問題重特征值在等導情況下特征對導數(shù)的新算法．

1 特征值的導數(shù)

假定二次特征值問題(3)在p=p0處有r重半單特征值λ0,當參數(shù)p在p0的某鄰域內變化時,假設r重特征值λ0變成r個單特征值λ1(p),…,λr(p),u1(p),…,ur(p)分別是二次特征值問題(3)對應于特征值λ1(p),…,λr(p)的特征向量,并且λ1(p0)=…=λr(p0)=λ0≠λi(p0) (i>r).為便于討論,記

Λ(p)=diag(λ1(p),…,λr(p))

U(p)=[u1(p),…,ur(p)]

由式(3),有

M(p)U(p)Λ2(p)+C(p)U(P)Λ(p)+K(p)U(p)=0

(4)

為了保證特征向量的唯一性,使用如下規(guī)范化條件：

(5)

當p=p0時,因為Λ(p0)=λ0Ir,此時規(guī)范化條件轉化為：

UT(p0)(2λ0M(p0)+C(p0))U(p0)=Ir

(6)

為方便起見,如無特別聲明，矩陣或向量在p=p0處取值時將p0省略,例如M(p)在p=p0處的取值M(p0)就簡記為M.

考慮計算特征值的導數(shù)：

特征向量的導數(shù)：

利用求解二次特征值問題的數(shù)值方法可計算二次特征值問題(3)在p=p0處的特征值λ0和相應滿足規(guī)范化條件的特征向量φ1,…,φr.記Φ=[φ1,…,φr],則Φ滿足：

ΦT(2λ0M+C)Φ=Ir

(7)

注意計算所得的特征向量Φ不一定恰好是特征向量U(p)在p=p0的值U.但span(Φ)=span(U),則存在非奇異矩陣Γ,使得：

U=ΦΓ

(8)

并且Γ滿足ΓTΓ=Ir.

通常在計算二次特征值問題(3)的特征向量時,幾乎不可能正好找到對參數(shù)p解析的特征向量U,而只能計算Φ,再由Φ尋找矩陣Γ,從而得到U=ΦΓ.記

對式(3)兩邊關于pk求導并取p→p0有：

(9)

(10)

2 特征向量的導數(shù)

U[U1,U2,…,Uh]=ΦΓ=Φ(Γ1,Γ2,…,Γh)

(11)

Γs=ρsβs,s=1,…,h

(12)

式中，ρs滿足方程：

(13)

對任意ms階正交矩陣βs,Γs=ρss也是特征值問題 (13) 的解, 因此Γs和Us=ΦΓs不是唯一確定的．

(14)

式中：G∈Cn×n為矩陣D的一個廣義逆,即G滿足DGD=D；α1為任意的r×r矩陣．事實上,GQ1為方程組(9)的一個特解．記H=2λ0M+C,取α1=-ΦTHGQ1,代入(13),可得方程組(9)的另一個特解,

(15)

式中，P=I-HΦΦT.由于：

則

(16)

(UTH-ΓTΦTHΦΦTH)GQ1=0

(17)

因此,特征向量的導數(shù)可以表示為：

(18)

討論如何計算系數(shù)矩陣c1.記

對式(3)兩邊關于pk求二階導，并取p→p0,可得：

(19)

上式兩邊左乘UT,可得：

(20)

將式(16)代入式(20),并利用規(guī)范化條件(6)得

(21)

把方程(21)寫成ms×mt的分塊形式,可得：

(22)

式中：

在方程(21)中令s=t，并在兩邊左乘βs,可得：

(23)

在方程(21)中取s≠t,可得矩陣c1的非對角子塊：

(24)

類似于方程(3)的討論,將式(19)的解表示為：

(25)

(26)

對式(3)兩邊關于pk求三階導，p→p0,并利用關系(9，21，25)有：

(27)

在式(27)兩邊取出相應的對角子塊可得：

(28)

它的純量形式為：

(29)

(30)

矩陣c1的對角元可由規(guī)范化條件求得．對方程(5)兩邊求導并令p→p0有：

(31)

將式(6，7)代入上式可得：

(32)

因此可得矩陣c1的對角元為：

(33)

通過特征向量導數(shù)的計算可以得到在特征值的導數(shù)相等的情況下特征對導數(shù)的算法．

算法: 等導情況下特征對導數(shù)的算法

輸入:矩陣M，C，K以及他們關于參數(shù)的一階、二階、三階偏導數(shù)，特征值λ0和相應的特征向量Φ，滿足規(guī)范化條件ΦT(2λ0M+C)Φ=Ir.

步驟：

矩陣D的廣義逆G以及Ge=(I-ΦΦTH)G·(I-HΦΦT);

(3) 若Ξ的特征值有重根,計算

3 數(shù)值例子

考慮三自由度的彈簧質點阻尼系統(tǒng)(圖1).

圖1 彈簧-質點-阻尼系統(tǒng)Fig.1 Mass-damper-spring system

則系統(tǒng)的質量、剛度和阻尼矩陣分別為：

設參數(shù)m1=m2=m3=1 kg,

k1=k4=k5=1 000 N/m，

k2=k3=0 N/m,c1=c2=10 Ns/m,

c3=20 Ns/m.因而,

選取阻尼矩陣中的參數(shù)c為待定參數(shù)．考慮系統(tǒng)在c=c0=0處特征對的導數(shù)．顯然,

當c=c0=0時，系統(tǒng)有二重特征值λ=-5-31.225i其對應的特征向量為:

利用算法1可得：

即特征值的導數(shù)有重根,而特征值二階導數(shù)為:

由算法1得：

為了驗證上述結果的正確性,設參數(shù)c的一個擾動Δc=0.1，然后計算當c=c0+Δc時系統(tǒng)的特征值與特征向量,并與如下近似計算結果：

比較，結果見表1、2．可見,c=c0+Δc時系統(tǒng)的特征值與特征向量，與近似計算的特征值與特征向量相差很小,說明文中的算法是有效的．

表1 特征值的比較Table 1 Comparison of eigenvalues

表2 特征向量的比較Table 2 Comparison of eigenvectors

4 總結

特征對導數(shù)的計算是特征靈敏度研究的主要問題．文中考慮二次特征值問題在特征值導數(shù)相等情況下特征對導數(shù)的計算問題．通過將特征向量的導數(shù)表示為一個用廣義逆矩陣表示的特解和對應齊次方程組的通解之和,給出了一種計算二次特征值問題等導特征對導數(shù)的算法．該算法在n維空間中直接計算對稱二次特征值問題重特征對的導數(shù),利用廣義逆矩陣導出了控制方程的一個特解,從而給出了計算重特征對導數(shù)的算法．和目前存在的算法相比,該算法可以適用于特征值導數(shù)有相等的情形,而且計算過程只需要計算系統(tǒng)當前特征值對應的特征對,計算量和存儲量較?。?/p>

在文中研究工作的基礎之上,下一步研究可從以下兩方面展開:① 文中考慮了二次特征值問題特征值等導情況下特征向量一階導數(shù)的計算,還可以考慮特征向量高階導數(shù)的計算問題．② 非線性常微分方程特征值問題、阻尼結構系統(tǒng)的動力分析、流體力學中的穩(wěn)定性分析、線性時滯系統(tǒng)的穩(wěn)定性分析等領域會出現(xiàn)一般的非線性特征值問題,非線性特征值問題特征對的導數(shù)研究值得關注．