王平心,吳頡爾,楊習貝

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

線性阻尼系統的自由振動方程可以表示為：

(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的鄰域內解析的矩陣值函數．

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

1 特征值的導數

假定二次特征值問題(3)在p=p0處有r重半單特征值λ0,當參數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.

考慮計算特征值的導數：

特征向量的導數：

利用求解二次特征值問題的數值方法可計算二次特征值問題(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)的特征向量時,幾乎不可能正好找到對參數p解析的特征向量U,而只能計算Φ,再由Φ尋找矩陣Γ,從而得到U=ΦΓ.記

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

(9)

(10)

2 特征向量的導數

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)

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

(18)

討論如何計算系數矩陣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)

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

算法: 等導情況下特征對導數的算法

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

步驟：

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

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

3 數值例子

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

圖1 彈簧-質點-阻尼系統Fig.1 Mass-damper-spring system

則系統的質量、剛度和阻尼矩陣分別為：

設參數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.因而,

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

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

利用算法1可得：

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

由算法1得：

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

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

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

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

4 總結

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

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