矩陣特征問題的求解

1、第五章 矩陣特征問題的求解5.1 引言在科學(xué)技術(shù)的應(yīng)用領(lǐng)域中，許多問題都歸為求解一個特征系統(tǒng)。如動力學(xué)系統(tǒng)和結(jié)構(gòu)系統(tǒng)中的振動問題，求系統(tǒng)的頻率與振型；物理學(xué)中的某些臨界值的確定等等。設(shè)A為n階方陣，若，有數(shù)l使Ax= lx （5.1）則稱l為A的特征值，x為相應(yīng)于l的特征向量。因此，特征問題的求解包括兩方面：1求特征值l，滿足（5.2）2求特征向量，滿足齊方程組（5.3）稱j（l）為A的特征多項式，它是關(guān)于l的n次代數(shù)方程。關(guān)于矩陣的特征值，有下列代數(shù)理論，定義1 設(shè)矩陣A, BÎR n´n，若有可逆陣P，使則稱A與B相似。定理1 若矩陣A, BÎR n´

2、;n且相似，則（1）A與B的特征值完全相同；（2）若x是B的特征向量，則Px便為A的特征向量。定理2 設(shè)AÎR n´n具有完全的特征向量系，即存在n個線性無關(guān)的特征向量構(gòu)成Rn的一組基底，則經(jīng)相似變換可化A為對角陣，即有可逆陣P，使其中l(wèi)i為A的特征值，P的各列為相應(yīng)于li的特征向量。定理3 AÎR n´n，l1, , ln為A的特征值，則（1）A的跡數(shù)等于特征值之積，即（2）A的行列式值等于全體特征值之積，即定理4 設(shè)AÎR n´n為對稱矩陣，其特征值l1l2ln，則（1）對任AÎR n，x0，（2）（3）定理5 （Ger

3、schgorin圓盤定理） 設(shè)AÎR n´n，則（1）A的每一個特征值必屬于下述某個圓盤之中，（5.4）（5.4）式表示以aii為中心，以半徑為的復(fù)平面上的n個圓盤。（2）如果矩陣A的m個圓盤組成的并集S（連通的）與其余n m個圓盤不連接，則S內(nèi)恰包含m個A的特征值。定理4及定理5給出了矩陣特征值的估計方法及界。例1 設(shè)有估計A的特征值的范圍。解 由圓盤定理，A的3個圓盤為圖5.1D1: D2: D3: 見圖5.1。D2為弧立圓盤且包含A的一個實特征值l1（因為虛根成對出現(xiàn)的原理），則3l15。而l2，l3ÎD1D2，則，即5.2 乘冪法與反冪法在實際工程應(yīng)用中，

4、如大型結(jié)構(gòu)的振動系統(tǒng)中，往往要計算振動系統(tǒng)的最低頻率（或前幾個最低頻率）及相應(yīng)的振型，相應(yīng)的數(shù)學(xué)問題便為求解矩陣的按模最大或前幾個按模最大特征值及相應(yīng)的特征向量問題，或稱為求主特征值問題。乘冪法是用于求大型稀疏矩陣的主特征值的迭代方法，其特點是公式簡單，易于上機實現(xiàn)。乘冪法的計算公式為：設(shè)AÎR n´n，取初始向量x(0)ÎR n，令x(1) = Ax(0)，x(2) = Ax(1)，一般有（5.5）形成迭代向量序列x(k)。由遞推公式（5.5），有（5.6）這表明x(k)是用A的k次冪左乘x(0)得到的，因此稱此方法為乘冪法，（5.5）或（5.6）式稱為乘冪公式

5、，x(k)稱為迭代序列。下面分析乘冪過程，即討論當(dāng)k時，x(k)與矩陣A的主特征值及相應(yīng)特征向量的關(guān)系。設(shè)A = (aij)n´n有完全的特征向量系，且l1, l2, ln為A的n個特征值，滿足v1, v2, vn為相應(yīng)的特征向量且線性無關(guān)，從而構(gòu)成Rn上的一組基底。對任取初始向量x(0) Î Rn，可由這組基底展開表示為 （5.7）其中a1, a2, an為展開系數(shù)。將x(0)的展開式（5.7）代入乘冪公式（5.6）中，得 （5.8）利用（5.8）式為 （5.9）（1）如果A有唯一的主特征值，即，設(shè)l1 ¹ 0，且由（5.9）式，有其中，由于，故當(dāng)k充分大時，e

6、 k » 0，此時（5.10）對i = 1, 2, , n，若(a1v1)i ¹ 0，考慮相鄰迭代向量的對應(yīng)分量比值，（5.11）即對i = 1, , n（5.12）這表明主特征值l1可由（5.11）或（5.12）式得到。由于迭代序列x(k)，當(dāng)k充分大時，（5.10）式成立，x(k)與v1只相差一個常數(shù)因子，故可取x(k)作為相應(yīng)于主特征值l1的特征向量的近似值。迭代序列x(k)的收斂速度取決于的大小。（2）如果A的主特征值不唯一，且可分三種情況討論：a）l1 = l2；b）l1 = - l2；c）情況a）當(dāng)l1 = l2時，A的主特征值為二重根，根據(jù)（5.9）式當(dāng)k充分

7、小時，由于，j = 3, n，e k » 0，則對i = 1, 2, , n，如果，則（主特征值）且x(k)收斂到相應(yīng)于l1 (=l2)的特征向量的近似值。這種重主特征值的情況，可推廣到A的r重主特征值的情況，即當(dāng) 且時，上述討論的結(jié)論仍然成立。情況b）當(dāng)l1 = - l2時，A的主特征值為相反數(shù)，（5.9）式為當(dāng)k充分大時，j = 3, 4, , n, e k » 0，則（5.13）由于（5.13）式中出現(xiàn)因子(-1) k，則當(dāng)k變化時，x(k)出現(xiàn)振蕩、擺動現(xiàn)象，不收斂，利用(-1) k的特點，連續(xù)迭代兩步，得從而，對i = 1, 2, , n，若，則（5.14）開方之

8、后，便得到A的以上主特征值l1, l2 = - l1。為計算相應(yīng)于l1, l2的特征向量，采取組合方式，（5.15）（5.16）可見分別為相應(yīng)于l1與l2的特征向量。情況c）當(dāng)時，A的主特征值為共軛復(fù)根。因A為實矩陣，于是由有即（v1與v2為互為共軛向量）。設(shè)，對任取x(0) ÎR n，展開式（5.7）可為（5.17）將（5.17）式代入（5.9）式，同理，當(dāng)k充分大時（5.18）對j = 1, 2, n，設(shè)復(fù)數(shù)表示則（5.18）式的復(fù)數(shù)表示可為連續(xù)迭代，得（5.19）利用三角函數(shù)運算性質(zhì)及l(fā)1、l2的復(fù)數(shù)表示，不難驗證。令（5.20）解方程（j = 1, 2, , n）（5.21）

9、求出p, q后，再解出主特征值l1、l2，得 （5.22）同樣，采取組合方式求相應(yīng)于l1、l2的特征向量。由于（5.23）（5.24）則可分別?。?.23）、(5.24）左端的組合表達式作為相應(yīng)于l1、l2的特征向量的近似值。通過上述分析，有定理6 設(shè)A ÎRn´n有完全特征向量系，若l1, l2, ln為A的n個特征值且滿足對任取初始向量x(0) ÎRn，對乘冪公式確定的迭代序列xk，有下述結(jié)論：（1）當(dāng)時，對i = 1, 2, , n收斂速度取決于的程度，r << 1收斂快，r » 1收斂慢，且x(k)（當(dāng)k充分大時）為相應(yīng)于l1的特征向

10、量的近似值。（2）當(dāng)時a）若l1 = l2，則主特征值l1及相應(yīng)特征向量的求法同（1）；b）若l1 = -l2，對i = 1, 2, , n收斂速度取決于的程度。向量、分別為主特征值l1、l2相應(yīng)的特征向量的近似值。c）若，則連續(xù)迭代兩次，計算出x(k+1)，x(k+2)，然后對j = 1, 2, , n解方程求出、后，由公式解出主特征值l1、l2。此時收斂速度取決于的程度。向量、分別為相應(yīng)于l1，l2的特征向量的近似值。從分析乘冪過程可見，乘冪法可用于求矩陣按模最大的一個（或幾個）特征值及相應(yīng)的特征向量，當(dāng)比值時，收斂速度快，r » 1時，收斂速度慢，且計算公式簡便，便于上機實現(xiàn)。

11、分析中的假設(shè)、，在計算時可不用考慮，如果此條件不滿足，則可通過迭代誤差自行調(diào)整。在用乘冪法求矩陣的主特征值l1及對應(yīng)的特征向量時，迭代向量的分量可能會出現(xiàn)絕對值非常大的現(xiàn)象，從而造成計算中溢出的可能。為此，需對迭代向量x(k)進行規(guī)范化。令max(x)表示向量x分量中絕對值最大者。即如果有某i0，使則max (x) = xi對任取初始向量x(0)，記則一般地，若已知x(k)，稱公式（5.25）為規(guī)范化的乘冪法公式或改進乘冪法公式，這里，乘冪迭代序列y(k)的分量絕對值最大者1。類似前面的分析乘冪過程，有定理7 設(shè)AÎRn´n具有完全特征向量系，l1, l2, , ln為A的

12、n個特征值，且滿足則對任初始向量x(0)，由規(guī)范化的乘冪法公式（5.25）確定的向量序列y(k)，x(k)滿足（1）（5.26）（2）y(k)為相應(yīng)于主特征值l1的特征向量近似值（5.27）例2 用規(guī)范化乘冪法計算矩陣A的主特征值及相應(yīng)特征向量解 A的特征值l1 = 6，l2 = 3，l3 = 2取初始值x(0) = (1, 1, 1)T，用規(guī)范化乘冪法公式（5.25）計算其它結(jié)果見表5.1（表中的向量均為轉(zhuǎn)置向量）。表5.1kmax(y(k)x(k) = y(k)/max(x(k)x(k+1) = Ay(k)01(1, 1, 1)(10, 8, 1)110(1, 0.8, 0.1)(7.2,

13、 5.4, -0.8)27.2(1, 0.75, -0.111111)(6.5, 4.75, -1.222222)36.57(1, 0.730769, -0.203704)(6.230766, 4.499997, -1.407408)46.230766(1, 0.722222, -0.225880)(6.111108, 4.388886, -1.1451767)56.111108(1, 0.718182, -0.237561)(6.054548, 4.336336, -1.475122)66.054548(1, 0.716216, -0.243639)(6.027024, 4.310808,

14、-1.487278)76.027024(1, 0.715247, -0.246768)(6.013458, 4.298211, -1.483536)86.013458(1, 0.714765, -0.248366)(6.00671, 4.291945, -1.496732)96.00671(1, 0.714525, -0.249177)(6.00335, 4.28825, -1.496354)106.00335(1, 0.714405, -0.249586)(6.00167, 4.287265, -1.499172)116.00167(1, 0.714345, -0.239792)(6.000

15、83, 4.286485, -1.499584)126.00083(1, 0.714315, -0.249896)取max(x(12) = 6.00083作為主特征值l1的近似值，與真值l1 = 6相比，有較好的近似程度，相應(yīng)于l1的特征向量的近似值取為y(2) = (1, 0.714315, -0.249896)T。當(dāng)li ( i = 1, 2, , n)為矩陣A Î Rn´n的n個特征值，且 時，乘冪法的收斂速度由決定，r << 1收斂得快。因此為提高收斂速度或改善r » 1的狀況，可以采取原點移位的方法，改變原矩陣A的狀態(tài)。取l0（常數(shù)），用矩陣

16、B = A - l0I來代替A進行乘冪迭代。設(shè)mi (i = 1, 2, , n)為矩陣B的特征值，則B與A特征值之間應(yīng)有關(guān)系式： (i = 1, 2, , n)且若vi是A相應(yīng)于li的特征向量，則vi亦是mi的特征值，即對i = 1, 2, , n因此，對任取x(0)ÎRn，關(guān)于矩陣B的乘冪公式（5.6）可為為加快收斂速度，適當(dāng)選擇參數(shù)l0，使（5.28）達到最小值。如當(dāng)li (i = 1, 2, , n)為實數(shù)，且l1>l2 ln時，取則為w(l0)的極小值點。這時原點移位法是一個矩陣變換過程，變換簡單且不破壞原矩陣的稀疏性。但由于預(yù)先不知道特征值的分布，所以應(yīng)用起來有一定

17、困難，通常對特征值的分布有一個大略估計，設(shè)定一個參數(shù)l0進行試算，當(dāng)所取l0對迭代有明顯加速效應(yīng)以后再進行確定計算。例3 計算A的主特征值解 先用規(guī)范化乘冪法計算，得表5.2表5.2ky(k) = x(k)/max(x(k)l1 » max (x(k)0(1, 1, 1)1(0.9091, 0.8182, 1)T2.7519(0.7482, 0.6497, 1)T2.536537420(0.7482, 0.6497, 1)T2.5365323主特征值l1及特征向量v1為（8位有效數(shù)字）l1 = 2.5362258v1 = (0.74822116, 0.64966116, 1)T而用規(guī)

18、范化乘冪法計算的相應(yīng)近似值為：l1 » max (x(20) = 2.5365323v1 » y(20) = (0.7482, 0.6497, 1)T如果采用原點位移的加速法求解，取l0 = 0.75，矩陣b = A - l0I對矩陣B應(yīng)用規(guī)范化乘冪法公式見表5.3。表5.3ky(k) = x(k)/max(x(k)l1 » max (x(k)0(1, 1, 1)T9(0.7483, 0.6497, 1)T1.786658710(0.7483, 0.6497, 1)T1.7865914可見此結(jié)果與未加速的規(guī)范化乘冪法公式計算結(jié)果相比，收斂速度要快得多。在已經(jīng)求出主

19、特征值l1及特征向量v1以后，可將原矩陣進行修改，使修改后的矩陣按模最大特征值是原矩陣的按模次大特征值，再用乘冪法去求按模次大特征值及特征向量，此方法稱為降階過程。為使問題簡單，設(shè)A ÎRn´n為對稱矩陣。假定已求出主特征值l1及特征向量v1，記A(1) =A，構(gòu)造矩陣（5.29）因A對稱，則具完全特征向量系，且特征向量vi可兩兩相互正交，即滿足（i = 2,n），于是 （i = 2, 3, , n）這表明矩陣A(2)除一個特征值l1 = 0與矩陣A(1)不一樣之外，其余與A(1)具相同特征值，且A(2)的按模最大特征值是l2，用A(2)代替A(1)進行乘冪迭代即可求得l2

20、及相應(yīng)的特征向量v2，而l2又為A(1)的按模次大特征值。這種做法稱為降階法。注意，降階法實際上只可用少數(shù)幾次，求矩陣的前幾個按模最大特征值及特征向量。因為，每降階一次，計算精度就會損失或降低一些。設(shè)AÎRn´n可逆，則無零特征值，由有即若l為矩陣A的特征值，則必為矩陣A-1的特征值，且特征向量相同。如果A的n個特征值li (i = 1, 2, , n)為則A-1的n個特征值更為因此，若乘冪法可求A的主特征值l1，則用A-1做乘冪矩陣，由乘冪迭代格式（5.30）便可求出A-1的按模最大特征值，取倒數(shù)，即為矩陣A的按模最小特征值。因此，對任取初始向量x(0)ÎRn，

21、稱公式（5.30）為求矩陣A按模最小特征值的反冪法。在應(yīng)用公式（5.30）計算時，由于要計算A的逆矩陣A-1，一方面計算復(fù)雜、麻煩，另一方面，有時會破壞A的稀疏性，故改寫（5.30）式為： （k = 0, 1, ）（5.31）類似于公式(5.25)的規(guī)范化乘冪法公式為（5.32）如果考慮到利用原點移位加速的反冪法，則記B = A - l0I，對任取初始向量x(0)ÎRn，（5.33）由于反冪法的主要工作量是每迭代一步都要解一個線性方程組(5.31)，且系數(shù)矩陣A (或B)是不變的，故可利用矩陣的三角分解A = LU (或B = LU)，則每次迭代只需解二個三角形方程組（5.34）且當(dāng)

22、時（5.35）同時x(k+1)便為所求的特征向量，收斂速度為。反冪法的主要應(yīng)用是已知矩陣的近似特征值后，求矩陣的特征向量，且收斂快，精度高，是目前求特征向量最有效的方法之一。5.3 子空間迭代法子空間迭代法也稱為平行迭代法，是乘冪法的推廣。乘冪法每次只能求出矩陣的一個主特值及特征向量，而子空間迭代法一次可求出矩陣的前幾個按模最大特值及特征向量。首先介紹斯密特（Schmidt）正交化過程。以三維為例，設(shè)a1，a2，a3為R3上的三個線性無關(guān)的向量，令，則b1為單位長度的向量，再令可以驗證（b1, b2）= 0，即b1與b2正交。若令則即與b1, b2正交，將其單位化為于是向量組b1, b2, b

23、3構(gòu)成R3上一組標(biāo)準(zhǔn)正交基，且其中Q = b1, b2, b3為正交矩陣，R是上三角陣。對n維向量空間，設(shè)a1, , an為Rn上n個線性無關(guān)的向量，類似有， ，即Q為正交陣，R為上三角陣，將n個線性無關(guān)向量變換為n個兩兩正交向量的方法稱為斯密特正交化方法。即斯密特正交化過程將可逆陣A分解為正交陣與上三角陣的乘積。當(dāng)向量組a1, , ap的個數(shù)p<n時，上述過程仍然成立。子空間迭代法的做法是，若要求A的前p個特征值及相應(yīng)的特征向量（p<n），先取p個線性無關(guān)的向量x1, x2, , xp構(gòu)成一個n ´ p的初始矩陣用矩陣A左乘上式，Z1 = AY0Z1的各列線性無關(guān)，但不

24、一定正交，采用斯密特正交化過程，相當(dāng)于正交分解Z1為Z1 = Y1R1 或 Y1 = Z1 R1-1Y1的各列正交，且R1為可逆上三角陣，或AY0 = Y1R1若已知第k步迭代近似陣為Yk，則子空間迭代法公式，（5.36）子空間迭代法適用于大型對稱稀疏矩陣特征問題的求解。當(dāng)AÎRn´n對稱時，可按下述正交化步驟進行計算，1計算2計算p階對稱陣 ；3求矩陣使其中對角陣的元素為Gk+1的全部特征值，Uk+1的各列u1, up為相應(yīng)特征向量。4計算可以驗證Yk+1為各列規(guī)范正交，即5若（j = 1, 2, p）則便為A的前p個按模最小特征值，Yk+1便為相應(yīng)的特征向量陣，否則Yk

25、+1®Yk，轉(zhuǎn)（1）。注：當(dāng)p<<n時，可用其它方法（見下節(jié)雅可比方法）求低階對稱矩陣的特征值，特征向量。5.4 對稱矩陣的雅克比 (Jacobi) 旋轉(zhuǎn)法雅克比 (Jacobi) 方法是求實對稱矩陣全部特征值及對應(yīng)的特征向量的方法.它也是一種迭代法，其基本思想是把對稱矩陣A經(jīng)一系列正交相似變換約化為一個近似對角陣，從而該對角陣的對角元就是A的近似特征值，由各個正交變換陣的乘積可得對應(yīng)的特征向量。1預(yù)備知識雅克比方法涉及較多的代數(shù)知識，要承認如下一些主要結(jié)論：1）若B是上（或下）三角陣或?qū)顷嚕瑒tB的主對角元素即是B的特征值。2）若矩陣P滿足PTP = I，則稱P為正交矩

26、陣。顯然PT = P-1，且P1, P2, 是正交陣時，其乘積P = P1P2Pk仍為正交矩陣。3）稱矩陣 (5.37）為旋轉(zhuǎn)矩陣，它是在單位陣I的i行，且j行和i列、j列的四個交叉位置上分別置上cosq ，sinq，-sinq 和cosq 而成的。容易驗證旋轉(zhuǎn)矩陣是正交矩陣，即RT (i, j) = R-1 (i, j)，所以用它作相似變換陣時十分方便。雅克比方法就是用這種旋轉(zhuǎn)矩陣對實對稱陣A作一系列的旋轉(zhuǎn)相似變換，從而將A約化為對角陣的。用Pij(i, j)作旋轉(zhuǎn)變換的幾何意義是：在維空間中，以i, j軸形成的平面上，把i, j軸旋轉(zhuǎn)一個角度q 。2雅克比方法先以二階矩陣為例：旋轉(zhuǎn)矩陣為其中為使A的相似矩陣B成為對角陣，只須適當(dāng)選取q，使即 ，其中，a11 = a22時，取。由此q 可以確定，從而旋轉(zhuǎn)矩陣P確定。A的特征值為：l1 = b11, l2 = b22關(guān)于特征向量的計算設(shè)，其中x1, x2為P的行向量，因為 PART = B APT = PTB或 （Ax1, Ax2）= （l1x1, l2x2），即Axi = liZi (i = 1, 2, 3)，所以對應(yīng)于l1，l2的特征向量是考慮n階矩陣的情況：設(shè)矩陣AÎRn´n是對稱矩陣，記A0 = A，對A作一系列旋轉(zhuǎn)相似變換，即（5.38）其中Ak (k = 1, 2,