有理特征值問題兩類改進的數(shù)值求解算法

陳小平, 魏 偉, 潘小明, 史雪瑩, 羅 安

(1. 泰州學(xué)院數(shù)理學(xué)院, 江蘇 泰州 225300; 2. 東南大學(xué)數(shù)學(xué)學(xué)院, 南京 210009;3. 南京林業(yè)大學(xué)理學(xué)院, 南京 210037)

考慮有理特征值問題: 尋找特征值λ∈R, 非零向量x∈Rn, 使得

T(λ)x=0,

(1)

其中T(λ)是關(guān)于參數(shù)λ的有理矩陣函數(shù)且T(λ)∈Rn×n, 并稱(λ,x)為式(1)的一個特征對．上述問題廣泛應(yīng)用于高速列車自動發(fā)射器的優(yōu)化[1]、流-固結(jié)構(gòu)振動[2]和量子點結(jié)構(gòu)計算[3-4]等．然而, 實際應(yīng)用過程中往往無須求出有理特征值問題(1)的所有特征對, 研究者真正感興趣的主要是部分特征對(λ,x)∈I×H, 式中I?R為某一區(qū)間,H?Rn是Hilbert空間．目前, 求解有理特征值問題(1)的第一類常見方法是通過線性化等方法[5-6]將原有理特征值問題(1)轉(zhuǎn)化為等價的多項式特征值問題, 然后進行求解．第二類常用方法是直接將有理特征問題(1)視為一般的非線性特征值問題, 這類方法主要包括矩陣分解方法[7-10]、 插值類方法[11]、 投影類方法[12-14]和逐次近似方法[15-16]等．這些方法雖能較好地求解有理特征值問題(1), 但在求解時不能利用原有理特征值問題(1)的某些特殊結(jié)構(gòu)或性質(zhì)．最近, Chen等[17]研究了具有特殊結(jié)構(gòu)的有理特征值問題(1)的譜性質(zhì), 并結(jié)合這些譜性質(zhì)發(fā)展了求解有理特征值問題的二分迭代算法和Rayleigh函數(shù)迭代算法．然而, 該算法在實際計算時依賴初值的選取, 不合適的初值可能會導(dǎo)致迭代不收斂, 故本文擬通過區(qū)間變換來改進有理特征值問題的數(shù)值求解方法．

眾所周知,有理特征值問題(1)中的有理矩陣函數(shù)T(λ)通常可表示成如下標準形式:

(2)

其中P(λ)=λmAm+λm-1Am-1+…+λA1+A0為矩陣多項式;fi(λ),gi(λ)分別是關(guān)于參數(shù)λ的階數(shù)為mi,ni互質(zhì)的標量多項式函數(shù), 且mi

(3)

其中fi(λ)=ai,mi-1λmi-1+…+ai,1λ+ai,0,gi(λ)=λni+bi,ni-1λni-1+…+bi,1λ+bi,0．

類似文獻[17], 本文僅探討有理特征值問題(3)在正半實軸上的特征值分布情況, 并假設(shè)如下[17]:

假設(shè)1矩陣Ai為對稱正定矩陣, 對于任意i=1,2,…,m, 且A0為對稱矩陣．

假設(shè)2矩陣Bi為低秩對稱半正定矩陣, 對于任意i=1,2,…,j．

1 兩種數(shù)值方法求解有理特征值問題

有理特征值問題(3)的譜理論詳見文獻[17]．根據(jù)假設(shè)2和低秩擾動理論[18], 有理特征值問題(3)可近似轉(zhuǎn)化為多項式特征值問題

(4)

2 兩種改進的數(shù)值方法求解有理特征值問題

2.1 兩種改進的數(shù)值方法的構(gòu)造

一般事先無法預(yù)知有理特征值問題(3)的λ, 故若選定的σ接近于區(qū)間Il的邊界, 則會導(dǎo)致算法1和算法2的計算效率較差,甚至算法失效, 即最接近σ的μ不在Il內(nèi)．為改善上述問題, 本文基于區(qū)間變換法提出求解有理特征值問題的兩種改進的數(shù)值方法．

T1(ξ)x=0．

(5)

不難發(fā)現(xiàn),λ為問題(3)特征值的充分必要條件是:ξ為有理特征值問題(5)的特征值, 且上述問題具有相同的特征向量．

另外, 若計算有理特征值問題(3)在Il=(σl-1,σl)內(nèi)的特征值λ．不妨令ζ=(λ-σl-1)/(σl-σl-1), 則λ=ζ(σl-σl-1)+σl-1, 式中ζ∈(0,1)．將此式代入有理特征值問題(3), 則問題(3)變?yōu)樾碌挠欣硖卣髦祮栴}, 記作

T2(ζ)x=0,

(6)

注1類似算法3, 可得求解有理特征值問題(3)的改進Rayleigh泛函加速迭代法(算法4), 即將算法3第8步中“ω=(ω+μ)/2”改為“ω=pl(x)”即可．

2.2 改進算法的收斂性

定理1若假設(shè)1～3均成立, 則算法3和算法4產(chǎn)生的迭代序列分別具有局部線性收斂性和二階收斂性．

證明 由二分法以及Rayleigh商迭代收斂性[17]可知,算法3具有局部線性收斂性, 且算法4具有局部二階收斂性．

3 數(shù)值實驗

通過與算法1和算法2的對比, 說明本文改進算法的有效性, 其中λk表示各算法運行第k步的近似特征值．類似文獻[17], 數(shù)值分析時僅考慮有理特征值問題在I1=(0,η)和I2=(η,100)上特征值的分布, 并取參數(shù)η=1和精度τ=10-12．

例1[5]有理特征值問題

(7)

取n=1 000, 容易驗證λ*=0.457 318 5為有理特征值問題(7)的精確特征值, 因為T(λ*)的最小奇異值小于10-9．利用MATLAB R2019a編程, 置σ=0.5, 運行算法3和算法4, 計算結(jié)果如表1所示．由表1可見, 算法3和算法4分別具有線性和二階收斂性．

表1 算法3和算法4的數(shù)值結(jié)果

取n=2 000, 容易驗證λ*=0.457 318 5和λ*=63.723 821 1分別是有理特征值問題(7)在區(qū)間I1=(0,1)和區(qū)間I2=(50,100)內(nèi)的2個精確特征值, 因為它們的最小奇異值均小于10-9．分別置σ=0.5和σ=75, 不同算法的計算結(jié)果見表2．由表2可見, 算法3和算法4的數(shù)值效果分別優(yōu)于算法1和算法2, 因此區(qū)間變換是有必要的．此外, 算法1和算法2能否計算出有效結(jié)果主要取決于σ取值是否合適．數(shù)值計算結(jié)果表明, 若σ取值接近各區(qū)間的邊界,則計算會出現(xiàn)不收斂的情形．

表2 四種算法在區(qū)間I1和I2內(nèi)的數(shù)值結(jié)果

取n=1 000,σ=1．ω取0.5和13時, 算法4的計算耗時分別為0.03和0.04 s, 而投影方法[19]在ω1=0.3,ω2=0.4,ω3=0.5,ω4=0.6和ω1=5,ω2=9,ω3=13,ω4=17條件下的計算耗時均達0.08 s．計算結(jié)果表明, 本文提出的算法4比文獻[19]中的投影方法更快．

4 結(jié)論

本文基于區(qū)間變換法提出兩種求解有理特征值問題的改進的迭代算法以及相應(yīng)的理論結(jié)果,新方法的收斂因子有待進一步研究.雖然數(shù)值實驗驗證了結(jié)論的正確性以及新算法的有效性, 但是對于更一般的有理特征值問題譜性質(zhì)及其相應(yīng)的數(shù)值求解算法有待深入研究．