第4章矩陣特征值和特征向量的計(jì)算

1、第四章矩陣的特征值與特征向量問題物理、力學(xué)和工程技術(shù)中的許多問題在數(shù)學(xué)上都?xì)w結(jié)為求矩陣的特征值和特征向量問題.計(jì)算方陣A的特征值，就是求特征方程nnn _2PlP2Pn 7的根.求出特征值后，再求相應(yīng)的齊次線性方程組A - I x = 0的非零解，即是對應(yīng)于的特征向量.這對于階數(shù)較小的矩陣是可以的 ，但對于階數(shù)較大的 矩陣來說，求解是十分困難，所以用這種方法求矩陣的特征值是不切實(shí)際的我們知道,如果矩陣A與B相似，則A與B有相同的特征值.因此人們就希望在相似變換 下,把A化為最簡單的形式.一般矩陣的最簡單的形式是約當(dāng)標(biāo)準(zhǔn)形 由于在一般情況下，用 相似變換把矩陣A化為約當(dāng)標(biāo)準(zhǔn)形是很困難的，于是人們

2、就設(shè)法對矩陣 A依次進(jìn)行相似變換 使其逐步趨向于一個(gè)約當(dāng)標(biāo)準(zhǔn)形，從而求出A的特征值.本章介紹求部分特征值和特征向量的幕法，反幕法；求實(shí)對稱矩陣全部特征值和特征向量的雅可比方法；求特征值的多項(xiàng)式方法；求任意矩陣全部特征值的QR方法第一節(jié)幕法與反幕法一 冪法幕法是一種求任意矩陣 A的按模最大特征值及其對應(yīng)特征向量的迭代算法該方法最大的優(yōu)點(diǎn)是計(jì)算簡單，容易在計(jì)算機(jī)上實(shí)現(xiàn)，對稀疏矩陣較為合適，但有時(shí)收斂速度很慢為了討論簡單，我們假設(shè)(1) n階方陣A的特征值i，'2,，n按模的大小排列為十2徉屮n|(1)(2) Vi是對應(yīng)于特征值i的特征向量i T，2，n ；(3) Vl，V2，,Vn線性無關(guān)

3、.任取一個(gè)非零的初始向量Xo，由矩陣A構(gòu)造一個(gè)向量序列x1 = Ax0x2 = Ax1AXz稱為迭代向量由于Vl，V2,Vn線性無關(guān)，構(gòu)成n維向量空間的一組基，所以，初始向量X。 可唯一表示成于是x0 二 a1v1a2v2anVn因?yàn)楸戎担? i =2,3, ,n ,所以an2k wk 1k 2o-k丄k亠丄_ ka1 1V1a22V2an n£k 1a-jV-ia2.Xka1v1當(dāng)k充分大時(shí)有kxk1 a1v1從而這說明當(dāng)k充分大時(shí)，兩個(gè)相鄰迭代向量 Xk 1與Xk近似地相差一個(gè)倍數(shù)，這個(gè)倍數(shù)便是矩陣 A的按模最大的特征值 J若用Xk i表示向量Xk的第i個(gè)分量，則.(X"

4、;) 1 :(xk )(8)也就是說兩個(gè)相鄰迭代向量對應(yīng)分量的比值近似地作為矩陣A的按模最大的特征值因?yàn)閤k 1 " ' 1xk,又Xk 1二AXk,所以有AXk 1 Xk,因此向量Xk可近似地作為對應(yīng)于1的特征向量這種由已知的非零向量 X。和矩陣a的乘幕構(gòu)造向量序列'Xk以計(jì)算矩陣a的按模最大 特征值及其相應(yīng)特征向量的方法稱為冪法由(4)式知，幕法的收斂速度取決于比值2'1的大小.比值越小，收斂越快，但當(dāng)比值接近于1時(shí),收斂十分緩慢人>1,則迭代向量Xk的各個(gè)不為零的分量將隨著k無限增用幕法進(jìn)行計(jì)算時(shí)，如果九1 <1,則Xk的各分量將趨于零這樣在

5、有限字長的計(jì)算機(jī)上計(jì)大而趨于無窮反之,如果算時(shí)就可能溢出停機(jī)為了避免這一點(diǎn)，在計(jì)算過程中，常采用把每步迭代的向量Xk進(jìn)行規(guī)范化,即用Xk乘以一個(gè)常數(shù)，使得其分量的模最大為 1.這樣，迭代公式變?yōu)閗 = Ax(mk = max(yk)( k = 1,2；八必=yjmk其中mk是yk模最大的第一個(gè)分量相應(yīng)地取(10)2-10A 二-12-1（精確到小數(shù)點(diǎn)后三位）。_ 0 - 1 2取 X。= 1,1,1 T用幕法求其模為最大的特征值及其相應(yīng)的特征向量kTYkmkXk1101110122-2221-1133-43-4-0.751-0.754-2.53.5-2.53.5-0.7141-0.7145-2

6、.4283.428-2.4283.428-0.7081-0.7086-2.4163.416-2.4163.416-0.7071-0.7077-2.4143.414-2.4143.414-0.7071-0.707解計(jì)算結(jié)果如表4-1所示。表4-1當(dāng)k=7時(shí)，Xk已經(jīng)穩(wěn)定，于是得到m7 二 3.414及其相應(yīng)的特征向量 V1為M %x7 =（-0.707,1,-0.707$應(yīng)用幕法時(shí)，應(yīng)注意以下兩點(diǎn)：（1）應(yīng)用幕法時(shí)，困難在于事先不知道特征值是否滿足（1）式，以及方陣A是否有n個(gè)線性無關(guān)的特征向量.克服上述困難的方法是： 先用幕法進(jìn)行計(jì)算，在計(jì)算過程中檢查是否出現(xiàn) 了預(yù)期的結(jié)果.如果出現(xiàn)了預(yù)期的結(jié)

7、果，就得到特征值及其相應(yīng)特征向量的近似值；否則，只能用其它方法來求特征值及其相應(yīng)的特征向量（2）如果初始向量X。選擇不當(dāng)，將導(dǎo)致公式（3）中V的系數(shù)1等于零.但是，由于舍入誤k差的影響，經(jīng)若干步迭代后，Xk = A X0 .按照基向量v1，v2，,Vn展開時(shí)，v1的系數(shù)可能不等于零。把這一向量 ?？醋鞒跏枷蛄浚媚环ɡ^續(xù)求向量序列 Xk 1，Xk 2，.,仍然會得出預(yù) 期的結(jié)果，不過收斂速度較慢.如果收斂很慢，可改換初始向量.原點(diǎn)平移法的大小當(dāng)比值接近于1時(shí)，收斂可能由前面討論知道，幕法的收斂速度取決于比值 很慢.這時(shí),一個(gè)補(bǔ)救的方法是采用原點(diǎn)平移法.設(shè)矩陣B 二 A - pl(11)其中p為

8、要選擇的常數(shù).我們知道A與B除了對角線元素外，其它元素都相同，而A的特征值'i與B的特征J i 之間有關(guān)系一 p,并且相應(yīng)的特征向量相同.這樣，要計(jì)算A的按模最大的特征值，就是適當(dāng)選擇參數(shù) p,使得1 一 p仍然是B的按模最大的特征值,且使人2 - P爲(wèi)|<九1 一 P|對B應(yīng)用幕法，使得在計(jì)算B的按模最大的特征值1 一 p的過程中得到加速，這種方法稱為原點(diǎn)平移法.例2設(shè)4階方陣A有特征值 i = 15 i, i = 1,2 3 4r比值= 0.9p =12作變換B 二 A - pl則B的特征值為J12,-1,0,-1應(yīng)用幕法計(jì)算B的按模最大的特征值 S時(shí)，確定收斂速度的比值為所

9、以對B應(yīng)用幕法時(shí)，可使幕法得到加速。p值，可以使得幕法得到加速，但由于矩陣的特征值的分布情況事先并不知道，所以在計(jì)算時(shí)，用原點(diǎn)平移法有一定的困難下面考慮當(dāng) A的特征值為實(shí)數(shù)時(shí)，如何選擇參數(shù) P，以使得用幕法計(jì)算1時(shí)得到加速 的方法.設(shè)A的特征值滿足 則對于任意實(shí)數(shù)p, B = A - pl的按模最大的特征值1 一 P或n - P。如果需要計(jì)算 、及V1時(shí)，應(yīng)選擇p使人-P =人- P且確定的收斂速度的比值r = max九+ k2nP =p ,即2 時(shí)，r為最小.這時(shí)用幕法計(jì)算如果需要計(jì)算 n及Vn時(shí)，應(yīng)選擇P使且確定收斂速度的比值> =min丸n - P J+ Z上 1f'-nA

10、P =當(dāng)r-pr-rji-p即2 時(shí)，r為最小這時(shí)用幕法計(jì)算'1及v時(shí)得到加速原點(diǎn)平移的加速方法，是一種矩陣變換方法這種變換容易計(jì)算，又不破壞A的稀疏性， 但參數(shù)p的選擇依賴于對A的特征值的分布有大致了解 三反幕法反幕法用于求矩陣 A的按模最小的特征值和對應(yīng)的特征向量，及其求對應(yīng)于一個(gè)給定的近似特征值的特征向量設(shè)n階方陣A的特征值按模的大小排列為相應(yīng)的特征向量為 v M，£.則A的特征值為對應(yīng)的特征向量仍然為 V1,V2,必.因此，計(jì)算矩陣A的按模最小的特征值，就是計(jì)算 A 的按模最大的特征值.這種把幕法用到 A 上,就是反幕法的基本思想.任取一個(gè)非零的初始向量 X。，由矩陣

11、A，構(gòu)造向量序列Xk=AXk4,k = 1,2,(12)用(12)式計(jì)算向量序列 '兀'時(shí)，首先要計(jì)算逆矩陣 A由于計(jì)算 A時(shí)，一方面計(jì)算 麻煩，另一方面當(dāng) A為稀疏陣時(shí)，A 不一定是稀疏陣，所以利用 A 在實(shí)際計(jì)算時(shí)，常采用解線性方程組的方法求Xk.( 12)式等價(jià)于AXk = Xk，k = 1，2為了防止溢出，計(jì)算公式為7 = Xkm max yN = yjmkk =1,2,(14)相應(yīng)地取Vn1& mJ:yk或 Vn ” Xk進(jìn)行計(jì)算會造成困難(13)(13)式中方程組有相同的系數(shù)矩陣A = LUA,為了節(jié)省工作量，可先對矩陣A進(jìn)行三角分解再解三角形方程組(15)

12、(16)工，“1,2,.(17)當(dāng)A是三對角方陣，或是非零元素較少且分布規(guī)律的方陣時(shí)，無論存儲或計(jì)算都比較便根據(jù)幕法的討論，我們知道，在一定條件下，可求得A 的按模最大的特征值和相應(yīng)的特征向 量，從而得到A的按模最小的特征值和對應(yīng)的特征向量 ，稱這種方法為反幕法反幕法也是 種迭代算法，每一步都要解一個(gè)系數(shù)矩陣相同的線性方程組設(shè)p為任一實(shí)數(shù)，如果矩陣A pI可逆，則A pI的特征值為1 1 11P 2 - P nP 對應(yīng)的特征向量仍為V，V2，,Vn .如果P是矩陣A的特征值'i的一個(gè)近似值，且九i - p v 糾 - p , i h j1則i 一 P是矩陣 A 一 pI 的按模最大的特

13、征值.因此，當(dāng)給出特征值i的一個(gè)近似 值P時(shí)，可對矩陣A 一 pI應(yīng)用反幕法，求出對應(yīng)于i的特征向量反幕法迭代公式中的 yk 通過方程組(A pI 必=Xjk = 1,2，求得 例3 用反幕法求矩陣2-100-12-10A =0 -1 2 -1P 0-12一的對應(yīng)于特征值=0.4的特征向量解取X。一 1111解方程組A - 0.4I 屮二 X。得二 _ 40 廠65廠65，一40 丁，g = max % - -65，1X1* 二 813,1,1,8 13 T .再解方程組A - 0.4I y2 二 X1得y2 F-445.13,- 72013，-720 13，-44513丁，mb 二 max

14、討2 二一 720 13,x2= y2 = （89144,1,1,89144）.Xi與x2的對應(yīng)分量大體上成比例，所以對應(yīng)于 °.4的特征向量為v = 89 144,1,1,89 144】第二節(jié) 雅可比方法雅可比方法是用來計(jì)算實(shí)對稱矩陣A的全部特征值及其相應(yīng)特征向量的一種變換方法.在介紹雅可比方法之前，先介紹方法中需要用到的線性代數(shù)知識與平面上的旋轉(zhuǎn)變換一預(yù)備知識（1） 如果n階方陣A滿足ata= I （即 A = A）則稱A為正交陣.（2）設(shè)A是n階實(shí)對稱矩陣，則A的特征值都是實(shí)數(shù)，并且有互相正交的n個(gè)特征向量（3）相似矩陣具有相同的特征值.（4）設(shè)A是n階實(shí)對稱矩陣，P為n階正交

15、陣，則B - PTAP也是對稱矩陣.（5）n階正交矩陣的乘積是正交矩陣.（6）設(shè)A是n階實(shí)對稱矩陣，則必有正交矩陣P，使PtAP 二(1)其中上的對角線元素的是A的n個(gè)特征值，向量.由（6）可知，對于任意的n階實(shí)對稱矩陣正交陣P的第i列是A的對應(yīng)于特征值i的特征A，只要能求得一個(gè)正交陣P,使PtAP （上的理論基礎(chǔ).為對角陣）,則可得到A的全部特征值及其相應(yīng)的特征向量,這就是雅可比方法旋轉(zhuǎn)變換 設(shè)為二階實(shí)對稱矩陣A|a1131211-321a22,即日12 = 321.因?yàn)閷?shí)對稱矩陣與二次型是一一2 2f x! ,x2 ' a11x12a12x1x2a22x2對應(yīng)的,設(shè)A對應(yīng)的二次型為

16、由解析幾何知識知道， 方程f X1,X2 =C表示在X1 ,X2平面上的一條二次曲線.如果將坐標(biāo)軸 0X1 ,°X2旋轉(zhuǎn)一個(gè)角度,使得旋轉(zhuǎn)后的坐標(biāo)軸 Oy1 ,°y2與該二次曲線的主軸重合 ，如圖4-1 所示，2 y2 = Cn21%這個(gè)變換就是cost -siny1sin cosTy2 一變換(4)把坐標(biāo)軸進(jìn)行旋轉(zhuǎn)COS 日P 二si nxJI,所以稱為旋轉(zhuǎn)變換-sinCOST.其中稱為平面旋轉(zhuǎn)矩陣。顯然有ptapPTP二I ,所以i(5)P是正交矩陣.上面的變換過程即-一2 一.由于a22 sin2 ja12sin2寸Ta11 cos2 二PTAP 二 IL0.5 a2

17、 a11 sin2) a12 cos2所以只要選擇"滿足0.5 a22 一 a11 sin2a11 sin a22cos2 寸a12 cos2a12 sin2日 _1a22 一 a sin2)a12cos2)- 022ai2tan2)=a11 a22(6)316 =(當(dāng)a11=a22時(shí)，可選取 4 ) PtAP就成對角陣，這時(shí) A的特征值為 a11 cosa22sin2ja12 sin2-2 = a sin2 丁a22 cos2 丁 - a12 sin2-相應(yīng)的特征向量為COSV1Sin日- si nV2IL cost雅可比方法雅可比方法的基本思想是通過一系列的由平面旋轉(zhuǎn)矩陣構(gòu)成的正

18、交變換將實(shí)對稱矩陣逐步n.首先引進(jìn)R化為對角陣，從而得到 變換A的全部特征值及其相應(yīng)的特征向量中的平面旋轉(zhuǎn)變換.py,其中X = yi cos： - yi sin=y sin 日 + yi cos日Xj二 ykk = i,jCOSTsincostX1,X2,.,Xny =( y1,y2,.,yn fP y n x xij為R中' j平面內(nèi)的一個(gè)平面旋轉(zhuǎn)變換Px則稱P X Xij稱為 2 j平面內(nèi)的平面旋轉(zhuǎn)矩陣.容易證明ij具有如下簡單性質(zhì)： P為正交矩陣. P的主對角線元素中除第 i個(gè)與第j個(gè)元素為 COS 中除第i行第j列元素為一 sin ,第j行第i列元素為sin外,其它元素均為零

19、.外，其它元素均為 1;非對角線元素(8)pj 二 PtA只改變A的第i行與第j行元素，AP只改變A的第i列與第j列元素，所以 PtAP只改變A的第i行、第j行、第i列、第j列元素.A 二勺 i I n 3a 二 a = 0設(shè) 'ij nn為n階實(shí)對稱矩陣，円 ji為一對非對角線元素.令A(yù) =PTAP = (af)h則A為實(shí)對稱矩陣，且A與A有相同的特征值.通過直接計(jì)算知(i)aH22=aH cos rajj sin 二a0 sin2ajj(1)2 2=aii sin 二 ajj cos 二-aij sin2a“ij(1) aik(1) a ik I jk(1) aki(1)1 /二

20、ajiJ ajj - aii(1).二 aki 二 aik cos(1).二 ajk 二 _ajkSinnii)sin2)cos2ajk sinnajk cost當(dāng)取二滿足關(guān)系式二 akltan 22ajaiajjk",jk",j(9)(10)時(shí)，引1 =a =0，且+蕊ifai： 2a, 2 * aj 22aik2aii2k k = i,ja行 2a2ak1 2a： k,lF,j由于在正交相似變換下，矩陣元素的平方和不變，所以若用平方和，用 SA表示 A的非對角線元素平方和，則由(11)式得；D(A)=D(A) + 2ai2S(A) = S(A)-2ai2D A表示矩陣(

21、11)A的對角線元素(12)這說明用ijPj對A作正交相似變換化為 A后，A1的對角線元素平方和比A的對角線元素平方和增加了 2aij , A1的非對角線元素平方和比A的非對角線元素平方和減少了2aij ,且將事先選)=0A定的非對角線元素消去了 (即ij ).因此，只要我們逐次地用這種變換，就可以使得矩陣 A的非對角線元素平方和趨于零，也即使得矩陣 A逐步化為對角陣.這里需要說明一點(diǎn)：并不是對矩陣A的每一對非對角線非零元素進(jìn)行一次這樣的變換就能得到對角陣.因?yàn)樵谟米儞Q消去 aij的時(shí)候，只有第i行、第j行、第i列、第j列元素在變化， 如果aik或Pj為零，經(jīng)變換后又往往不是零了.雅可比方法就

22、是逐步對矩陣A進(jìn)行正交相似變換，消去非對角線上的非零元素，直到將A的非對角線元素化為接近于零為止，從而求得 A的全部特征值，把逐次的正交相似變換矩陣乘起來，便是所要求的特征向量.雅可比方法的計(jì)算步驟歸納如下：第一步 在矩陣A的非對角線元素中選取一個(gè)非零元素aij. 一般說來，取絕對值最大的非對角線元素；第二步第三步tan2 氧a 一 a”求岀二，從而得平面旋轉(zhuǎn)矩陣 R = Fj ;由公式A = F ARA1的元素由公式（9）計(jì)算.以A代替A,重復(fù)第一、二、三步求岀A2及F，繼續(xù)重復(fù)這一過程，直到A（即小于允許誤差）時(shí)為止.第四步的非對角線元素全化為充分小第五步Am的對角線元素為A的全部特征值的

23、近似值，P=只卩2.巳的第j列為對應(yīng)于特征值j （ ' j為Am的對角線上第j 例1用雅可比方法求矩陣個(gè)元素）的特征向量.-10 12-1-12 一2_ 1_0的特征值與特征向量解首先取i =1，j =2,由于a11=a22 = 2,故取1 2 -1/2R = P；2316 =4 ,所以00-1 .2-1 血A =pTAF=.T 八S-1 2再取 i =1,j =3由2tan 2門1-2:0.45969 ,cos : 0.88808所以0.88808-0.459690459690.888080.63398-0.32505繼續(xù)做下去-0.32505A2 二 F2tAF2=IL,直到非對角

24、線元素趨于零 ，進(jìn)行九次變換后，得-0.62797-0.627972.366030.5875810.000000.0000011A9 =0.000002.000000.000000.000000.000003.41421 jA的對角線元素就是 A的特征值，即鮎賂 0.58758 ,人22.00000 *3止 3.41421相應(yīng)的特征向量為0.500000.707100.50000Vi = 0.70710 ,2 = 0.00000,v3 = -0.70710 .0.50000 一0.70710 一0.50000 一相應(yīng)的特征值的精確值站=2 - V2，人 2 = 2，九 3 = 2 + V2相應(yīng)

25、的特征向量為1 2v3 = | - 1/12 2 一由此可見，雅可比方法變換九次的結(jié)果已經(jīng)相當(dāng)精確了',特別是求得的特征向量正交性很好，所以雅可比方法是求實(shí)對稱矩陣的全部特征值及其對應(yīng)特征向量的一個(gè)較好的方法.但由于上面介紹的雅可比方法，每次迭代都選取絕對值最大的非對角線元素作為消去對象，花費(fèi)很多機(jī)器時(shí)間.另外當(dāng)矩陣是稀疏矩陣時(shí)，進(jìn)行正交相似變換后并不能保證其稀疏的性質(zhì)，所以對階數(shù)較高的矩陣不宜采用這種方法.目前常采用一種過關(guān)雅可比方法.這種方法是選取一個(gè)單調(diào)減小而趨于零的數(shù)列：a即二 0）作為限值，這些限值稱為”關(guān)”,對矩陣的非對角線元素規(guī)定一個(gè)順序（例如先行后列、自左至右的順序）.

26、首先對限值 日1按規(guī)定的順序逐個(gè)檢查矩陣的非對角線元 素,碰到絕對值小于 a1的元素就跳過去，否則就作變換將其化為零 .重復(fù)上述過程，直到所有的非 對角元素的絕對值都小于日1為止.再取日2 ,a3 ,.類似處理，直到所有的非對角線元素的絕對值都小于am時(shí),迭代停止.這時(shí)的am應(yīng)小于給定的誤差限實(shí)際運(yùn)算中常用如下的辦法取限值：對于矩陣A，計(jì)算A = A的非對角線元素平方和=仝 k=1,2,akS A0，任取N - n,取N4. 3多項(xiàng)式方法求特征值問題4.3.1 F-L方法求多項(xiàng)式系數(shù)我們知道，求n階方陣A的特征值就是求代數(shù)方程毋（九）=| a-丸1=0的根。（）稱為A的特征多項(xiàng)式。上式展開為(

27、432)()-' n - Pl ' n J - P2 ' n _2-Pn其中pi,p2,Pn為多項(xiàng)式()的系數(shù)。從理論上講，求 A的特征值可分為兩步：第一步 直接展開行列式| A - I求出多項(xiàng)式：(，)；第二步 求代數(shù)方程'(X)=0的根，即特征值。對于低階矩陣，這種方法是可行的。但對于高階矩陣，計(jì)算量則很大，這種方法是不適用的。這里我們介紹用 F-L ( Faddeev-Leverrier)方法求特征方程(4.3.2)中多項(xiàng)式，()的 系數(shù)。由于代數(shù)方程求根問題在第2章中已經(jīng)介紹，所以本節(jié)中解決特征值問題的關(guān)鍵是確定矩陣A的特征多項(xiàng)式 (-)，所以稱這種方法

28、為多項(xiàng)式方法求特征值問題。記矩陣a= (aj )n n的對角線元素之和為(4.3.3)trA =引1 a 22ann利用遞歸的概念定義以下n個(gè)矩陣Bk(k =1,2,，n):P1=trB1B1 = A,B2 = AQ - pj ),P2P3B3 = A(B2 - P2I),1trB221trB33Pk= -trBkkBk = A(Bk 4- P),Bn =A(Bn-PnJ),Pn1 trBn n(4.3.4)可以證明,(4.3.4)式中Pk,k -1,2,，n,即是所求式求矩陣的特征多項(xiàng)式系數(shù)的方法稱為A F-L方法。的特征多項(xiàng)式")相應(yīng)特征方程為：的各系數(shù)。用(4.3.4)n nn

29、-1n-2(T) (P1P2-Pn ) = 0而且可證矩陣A的逆矩陣可表示為A=丄(Bn4 - Pn)Pn(4.3.6)例1求矩陣324A= 202】423 一的特征值與a1解 用F-L方法求得3 24B!=A=2024 23jPi=trB= 61124B2 = A( Bi Pi I ) = 28242 111P2trB 2 =1528 0 0B3 = A( B2 p21) = 0 8 0 0 0 8 一 1“cp3trB 3=8所以A的特征方程為332(-1) ( -6 '-15'； -8)此方程的根，即特征值為J1 二 8,占2= -1, g = _1214 12 1A(B

30、2-P2l) =P3=0147814121412從例1中的計(jì)算結(jié)果可知 B3二P3l. Faddeev曾經(jīng)證明：對n階矩陣A,按(4.3.4) 式計(jì)算出的Bn總有Bn = Pn I(4.3.7)4.3.2特征向量求法當(dāng)矩陣A的特征向量確定以后，將這些特征值逐個(gè)代入齊次線性程組(A - I )x=0中,由于系數(shù)矩陣A -訂的秩小于矩陣 A - J的階數(shù)n,因此雖然有n個(gè)方程n個(gè)未知數(shù)，但實(shí)際上 是解有n個(gè)未知數(shù)的相互獨(dú)立的r個(gè)方程(r<n).當(dāng)矩陣A的所有特征值互不相同時(shí)，這樣的問題中要解的齊次方程組中有n-1個(gè)獨(dú)立方程，其中含有n個(gè)特征向量分量，因此特征向量分量中至少有一個(gè)需要任意假設(shè)其

31、值，才能求出其他特征分量.在計(jì)算機(jī)中解這樣的齊次線性程組，可用高斯-若當(dāng)消去法，以便把一組n個(gè)方程簡化為等價(jià)的一組n-1個(gè)方程的方程組.然而,用高斯-若當(dāng)消去法簡化一個(gè)齊次線性程組時(shí)，方程之間不都是獨(dú)立的，在消去過程中系數(shù)為零的情況較多.必需交換方程中未知數(shù)的次序，以避免主元素位置上為零的情況.因此，為了提高精度和避免零元素的可能性，我們總是用主元素措施把絕對值最大的系數(shù)放于主元素位置.例如，假設(shè)矩陣A為_42-21A=-532、一241 一其特征方程為4 丸-5展開后為-21 九=0（九一1）（人一2）（九一5） =0故特征值分別為：j = 1, 2 = 2,.3 = 5下面求特征向量，將1

32、代入方程組（A-J）x=°中，得3x1 2x2 - 2x3 = 0一5x1 +2x2 +2x3 = 0_2xr +4x2 +0x3 =0以-5為主元素,交換上式第一與第二個(gè)方程得匚5x2x2 2x3 = 0<3為 +2x2 _2x3 = 0_2xr +4x2 _0x3 =0(438)(439)用高斯-若當(dāng)消去法消去-5所在列中的Xi,并把主元素所在行調(diào)到最后，得1640x1+ X2一 X3 =055“ 0x1164心寫x2x3 = 0522X3 = 0X1 -X?55（4.3.10）再以16/5為主元素，消去它所在列中的X2 ,并把主元素所在的行調(diào)到最后，得0捲 +0x2 +0

33、x3 =01 % +0x2 x3 =0210x1x2x3 = 0L4（4.3.11）這就是用高斯若當(dāng)消去法實(shí)現(xiàn)把一組三個(gè)方程簡化為等價(jià)的一組兩個(gè)獨(dú)立方程的情形 這個(gè)等價(jià)的方程組包含兩個(gè)獨(dú)立的方程，而有三個(gè)未知數(shù)，所以只要假定其中一個(gè)值兩個(gè)值就可以通過兩個(gè)獨(dú)立方程解出比如,令x3 = T,則得到矩陣a的對應(yīng)于'1特征向量為.因?yàn)?貝U其它=1的一個(gè)2_14_ 1對另外兩個(gè)特征值的對應(yīng)特征向量求法與上述對'1 =1的推導(dǎo)過程相同計(jì)算機(jī)中實(shí)現(xiàn)求解這樣的齊次線性方程組的消去步驟是，用第3章討論過的高斯-若當(dāng)消415(4312)去法的公式，方程組(439)的系數(shù)矩陣經(jīng)過第一次消去后的矩陣

34、B為1651652'"5以矩陣為方程組(4.3.10)的系數(shù)矩陣，其中省略了有0和1元素的第一列在進(jìn)行第二次消元之前，要應(yīng)用完全主元素措施對前兩行進(jìn)行最大主元素選擇 ，然后再進(jìn) 行必要的行或列交換每完成一次消元過程，總省略只有0和1元素的第一列，并且計(jì)算機(jī)僅尋 找矩陣的前n-k行中的最大主元素，其中k是消元過程應(yīng)用的次數(shù)對(4.3.12)式再進(jìn)行一次消 元過程，則得到列矩陣-1B1此矩陣是對應(yīng)于方程組02，般來說，最后-4-(4.3.13)(4.3.11)的系數(shù)矩陣，不過省略了含0和1元素的前兩列矩陣列的數(shù)目等于矩陣A-'l的階數(shù)和秩的差值.由于方程組(4.3.8)有

35、三個(gè)未知數(shù)，兩個(gè)獨(dú)立方程，所以計(jì)算機(jī)必須任意給定一個(gè)未知數(shù)的 值，以便可以從其他兩個(gè)獨(dú)立方程中解出另外兩個(gè)未知數(shù).為方便，在計(jì)算機(jī)決定特征向量時(shí)要恰當(dāng)?shù)卦O(shè)定任意選取的未知數(shù)的值.例如，令X3 =由方程組(4.3.11)知道，其他兩個(gè)分量的值正好能從含x3的非零系數(shù)項(xiàng)得出.為此，從計(jì)算機(jī)所存儲的最終矩陣中，令B1最上面的0元素為-1,并把它順次調(diào)到最下面第三行的位置上，就得到所求的特征向量在工程問題中，從特征方程所求出的特征值，少數(shù)情形也有相同的.一般地，當(dāng)一個(gè)特征方程有k重根時(shí)，矩陣A- I的秩可能比其階數(shù)少1,或2,或3,或k，當(dāng)然對應(yīng)于的線性無 關(guān)的特征向量的個(gè)數(shù)也就是1,或2,或3,，或

36、k,下面通過一個(gè)特征值對應(yīng)兩個(gè)線性無關(guān)特征向量的例子進(jìn)一步說明計(jì)算機(jī)求特征向量的方法設(shè)矩陣A為_3 2 4A=2 0 24 2 3j其特征方程為3 -人242-九2=0423 -九展開后得('1)2( ' -8)-0所以特征值為/ .1 1 2為了決定，=-1的特征向量424/J212X2L42，得=0(4314)應(yīng)用一次高斯-若當(dāng)消去法J 1/20門1 10 X?= 01颯J(4.3.15)寫成矩陣形式,(4.3.15)式的系數(shù)矩陣為000101.1/2因?yàn)榉匠探M(4.3.15)的系數(shù)矩陣的秩為 關(guān)的特征向量，必須給兩個(gè)未知數(shù)任意規(guī)定值(4.3.15)式可看出，一般總是選擇x

37、2二，x3另一個(gè)特征向量；這樣有兩個(gè)線性無關(guān)的特征向量1 10(4.3.16)1,它比矩陣階數(shù)少2,因此對應(yīng)于1有兩個(gè)線性無,才能確定這兩個(gè)線性無關(guān)的特征向量，由=0求一個(gè)特征向量；選擇X2 = 0,X3 = -1求1/2：-1=-1,沁=8，將怎-1代入方程組(A - J )x=0,得00，在(4.3.16)式的B中，把第一列中第一個(gè)0計(jì)算機(jī)中求兩個(gè)線性無關(guān)的特征向量的辦法是元素用-1代替,第二列中第二個(gè)0元素也用-1代替，然后把第一、第二行順次調(diào)到最下面一行 的位置上，第三行自然就成了第一行，如此調(diào)換后矩陣的第一列和第二列就是所求的兩個(gè)線性無關(guān)的特征向量。對應(yīng)于=-1的全部特征向量為3/2

38、1I 一1+ k20.0 一L1J其中k1與k2是任意常數(shù)，且不同時(shí)為零。為了說明列交換的必要性，避免主元素為零，再舉一個(gè)例子，設(shè)矩陣| 2 8 121_0其特征方程為( -2)' (，-1)=0特征值為 2""2 0,兒3= 1對應(yīng)于，-2的特征向量可由解下列方程組而求得0-4-8-12 %124 k = 000-01x3 j用一次高斯-若當(dāng)消去法，得0 01 卜0 01 |x2 = 01 2 3眼(4318)若不進(jìn)行列交換，則下一個(gè)消元過程只能在第一行的第二個(gè)元素與第二行的第二個(gè)元素中找 最大主元素，而它們都是零，我們不得不對(4317)式進(jìn)行列交換，即交換未知

39、數(shù)之間的次序，之后再進(jìn)行消去過程對(4.3.17)式進(jìn)行列交換，即把絕對值最大系數(shù)放在主元素位置，顯然是第一列與第三 列的交換，交換后成為一12-8-4 x3 I421 怡=0100 衛(wèi)xj其中未知數(shù)列矩陣中x1與x3也進(jìn)行了交換，這樣才能保證(4.3.17)式與(4.3.19)式等價(jià)，對(4.3.19)式進(jìn)行一次高斯-若當(dāng)消去法，得-001-2/32/32/3-1 / 3 | x1/31/3 一X2.X1 一=0(4.3.20)再進(jìn)行一次消去過程,彳-得001000X20I0 11/2一lx(4.3.21)L1 1在計(jì)算機(jī)中計(jì)算，剩下一個(gè)最終的列矩陣0B =0(4.3.22)】1/2將(4.

40、3.22)式中的列矩陣B中第一個(gè)0元素用-1代替，并隨即調(diào)到最下面一行，便得到1/2LT -( 4.3.23)這就是對應(yīng)于方程組(4.3.19)的解，在計(jì)算機(jī)程序中應(yīng)把原來進(jìn)行列交換的列號次序記住,重新把(4.3.23)式中各分量排列一下，即交換第一行和第三行的元素，就得到對應(yīng)于怎=2的特征向量-111/20 一對應(yīng)于的全部的特征向量為-111/20 J其中k為不等于零的任意常數(shù)4.4 QR算法QR算法也是一種迭代算法，是目前計(jì)算任意實(shí)的非奇異矩陣全部特征值問題的最有效的方法之一.該方法的基礎(chǔ)是構(gòu)造矩陣序列'Ak并對它進(jìn)行QR分解.由線性代數(shù)知識知道，若A為非奇異方陣，則A可以分解為正

41、交矩陣 Q與上三角形矩陣R 的乘積，即A=QR，而且當(dāng)R的對角線元素符號取定時(shí)，分解式是唯一的若A為奇異方陣，則零為A的特征值任取一數(shù)p不是A的特征值，則A-pI為非奇異方陣 只要求出A-pI的特征值，就很容易求出 A的特征值，所以假設(shè)A為非奇異方陣，并不妨礙討論 的一般性形矩陣設(shè)A為非奇異方陣，令A(yù)1 = A，對A1進(jìn)行QR分解，即把A分解為正交矩陣 Q1與上三角 R1的乘積A1 = Q1 R1做矩陣A2 = RiQiA1Q1繼續(xù)對A2進(jìn)行QR分解并定義 一般地，遞推公式為A = A = Q1R1小心=RkQkA2 = Q2 R2A3 = R2Q2 = Q2 A2Q2-Qk AkQk , k

42、 = 2,3,QR算法就是利用矩陣的方陣，則由QR算法就完全確定所有人都相似，它們具有相同的特征值 因?yàn)镼R分解,按上述遞推公式構(gòu)造矩陣序列 %只要A為非奇異 f.這個(gè)矩陣序列'A"具有下列性質(zhì)性質(zhì)1證明Ak 1 二 RkQk 二 QT AkQk= QTQT4Ak4Qk4Qk =-Qk Qk 4 - Q1 AQ1Q2 Qk若令Qk =Q1Q2.Qk，則Qk為正交陣，且有TAk - Qk AkQk因此Ak與A相似，它們具有相同的特征值性質(zhì)2 Ak的QR分解式為kA =QkRk其中 Qk 二Q1Q2Qk,Rk 二 RkRk4Ri證明 用歸納法顯然當(dāng)k=1時(shí)，有A = A, = (

43、1= Q1R1k _1 A二Q2R2假設(shè)Ak 4有分解式Qk Rk - Q1Q2 .Q kd (Qk Rk) Rk.RiT因?yàn)锳k =Qk丄AQk，所以kQk R = Ak Qk 1 Rk 1 = A因?yàn)镼i,Q2,Qk都是正交陣，所以Qk也是正交陣，同樣Rk也是上三角形陣，從而Ak的qr分 解式為k A gRk由前面的討論知 Ak+ = Q AQk .這說明QR算法的收斂性有正交矩陣序列 k >的性質(zhì) 決定定理1如果Qk '收斂于非奇異矩陣 Q :, Rk為上三角形矩陣，則EH :Ak存在并且是上 三角形矩陣證明 因?yàn)镼k "收斂，故下面極限存在kmQk = kmQTQk = q；q 臨=iRM = kimRk kmAk 卅QkT Tlim Qk AQk a = Qy；AQ： k_.由于Rk（k刊2）為上三角形矩陣，所以R:為上三角形矩陣.又因?yàn)锳.- = lim Ak = lim QkRk = R-kJpC所以k'm :Ak存在,并且是上三角形矩陣.定理2 （QR算法的收斂性）設(shè)A為n階實(shí)矩陣，且1)A的特征值滿足入八2 > ' ' '>>0；2)A = Ai =XDX ：其中D =d