(完整word版)第9章矩陣特征值的數(shù)值解法

1、1第 9 章 矩陣特征值的數(shù)值解法9.1 引言矩陣特征值問題有廣泛的應(yīng)用背景.例如動(dòng)力系統(tǒng)和結(jié)構(gòu)系統(tǒng)中的振動(dòng)問題、電力系統(tǒng)的靜態(tài)穩(wěn)定分析上、工程設(shè)計(jì)中的某些臨界值的確定等，都?xì)w結(jié)為矩陣特征值問題.數(shù)學(xué)中諸如方陣的對(duì)角化及解微分方程組等問題，都要用到特征值的理論.本章介紹n階實(shí)矩陣A Rn n的特征值與特征向量的數(shù)值解法 .定義 9.1.1 已知n階實(shí)矩陣A二(aij)Rn n，如果存在常數(shù)和非零向量x,使AxVx或(A I )x=0(9.1.1)那么稱為A的特征值(eigenvalue)，x為A的相應(yīng)于的特征向量(eigenvector).多項(xiàng)式稱為特征多項(xiàng)式(characteristic p

2、olynomial)，det(A - I )=0稱為特征方程( (characteristic equation).注式(9.1.3)是以，為未知量的一元n次代數(shù)方程, 多項(xiàng)式.顯然，A的特征值就是特征方程(9.1.3)的根.特征方程(9.1.3)在復(fù)數(shù)范圍內(nèi)恒有解, 其個(gè)數(shù)為方程的次數(shù)(重根按重?cái)?shù)計(jì)算)，因此n階矩陣A在復(fù)數(shù)范圍內(nèi)有n個(gè)特征值.除特殊情況(如n =2,3或A為上(下)三角矩陣)外，一般不通過直接求解特征方程(9.1.3)來求A的特征值，原因是這樣的算法往往不穩(wěn)定.在計(jì)算上常用的方法是幕法與反幕法和相似變換方法.本章只介紹求矩陣特征值與特征向量的這兩種基本方法.為此將一些特征值

3、和特征向量的性質(zhì)列在此處定理 9.1.2 設(shè)n階方陣A =(aj)n n的特征值為1,2,川，n，那么(1)12Vn二知a?271 - ann；(2)i 2丨In= detA.定理 9.1.3 如果是方陣A的特征值，那么(1)k是Ak的特征值，其中k是正整數(shù)；(2)當(dāng)A是非奇異陣時(shí)，丄是A 的特征值.(3)Pn( )是Pn( A)的特征值，其中Pn(x)是多項(xiàng)式Pn(x)二a盼a?x2丨1|anXn.定義 9.1.4 設(shè)A, B都是n階方陣.若有n階非奇異陣P，使得P二AP = B，則稱矩 陣A與Bai2Pn( )二det(A- I)二a2ia22 -IIIIIIainIa2n(9.1.2)a

4、n1an2III(9.1.3)pn( )二det( A釘)是的n次2相似(similar)，P二AP稱為對(duì)A進(jìn)行相似變換(similarity transformation) ，P稱 為相似變換矩陣( (similarity transformation matrix).定理 9.1.5 若矩陣A與B相似，則A與B的特征值相同.定理 9.1.6 如果A是n階正交矩陣，那么AJ= AT，且det A= 1或-1;(2)若y = Ax，貝 Vy2=x2,即xT= yTy.定理 9.1.7 設(shè)A是任意n階實(shí)對(duì)稱矩陣，則(1)A的特征值都是實(shí)數(shù)；(2)A有n個(gè)線性無關(guān)的特征向量定理 9.1.8 設(shè)A是

5、任意n階實(shí)對(duì)稱矩陣，則必存在n階正交矩陣P，使得P AP =PTAP=上，其中：二diag(n)是以A的n個(gè)特征值仆匕,川，n為對(duì)角元素的對(duì)角矩陣定理 9.1.9 (圓盤定理)矩陣A二(aj)n n的任意一個(gè)特征值至少位于復(fù)平面上的幾個(gè)圓 盤nDi=z|z-aii,i=1,2，山,n,j=, j HJ中的一個(gè)圓盤上。9.2 冪法與反冪法9.2.1幕法及其加速9.2.1.1 幕法幕法是計(jì)算矩陣按模最大特征值(largest eigenvalue in magnitude)及相應(yīng)特征向量的迭代法該方法稍加修改，也可用來確定其他特征值幕法的一個(gè)很有用的特性是：它不僅可以求特征值，而且可以求相應(yīng)的特征

6、向量 實(shí)際上，幕法經(jīng)常用來求通過其他方法確定的特征 值的特征向量下面探討幕法的具體過程設(shè)矩陣A Rnn的n個(gè)特征值滿足人|丸2即乜|蘭川公卩 n|啟0，(921)且有相應(yīng)的n個(gè)線性無關(guān)的特征向量x2l(,xn，則x2l(,xn構(gòu)成n維向量空間Rnfn1的一組基,因此Rn= ?zz = H%兒，ER, i=1,2|,nj.3在Rn中選取某個(gè)滿足- 0的非零向量nZo - 7 JXi.i4用矩陣A同時(shí)左乘上式兩邊，得nnAZo:-i AXi八：iXi.i=1i=1再用矩陣A左乘上式兩邊，得n2、2A Zo二 :i iXi.i丄這樣繼續(xù)下去，一般地有nAkZo - 7、，k =d,2,|）|.i二k

7、記Zk= AZk=A Zo, k =1,2,川，則由式(922)得n kZk1 Xi.式（9.2.4）表明隨著k的增大，序列:豐越來越接近A的對(duì)應(yīng)于特征值1的特征向量Xi的：1倍，由此可確定對(duì)應(yīng)于1的特征向量X1.當(dāng)k充分大時(shí)，可得X1的近似值.上述收斂速度取決于比值.事實(shí)上，由式（9.2.3）知,再由式（9.2.1）得(9.2.2)nn=E昭,=雪1X,吃i人ii=2l入k=1,2,川.(9.2.3)由假設(shè)（921），結(jié)合式（923），得(9.2.4)于是對(duì)充分大的k有(9.2.5)Zkik(9.2.6)3(9.2.7)4結(jié)合式（9.2.6）和式（9.2.7）知，序列 Q/扎：收斂速度取決于

8、比值|婦/人5F 面計(jì)算i.由式（923）知k1k1Zk 1 AZk= AZQ=1當(dāng) k 充分大時(shí),zk+&霜妝X.結(jié)合式（9.2.5），得這表明兩個(gè)相鄰向量大體上只差一個(gè)常數(shù)倍，這個(gè)倍數(shù)就是A的按模最大特征值.記/和(n)；血若Zk, Zk, Zk,則有即兩個(gè)相鄰的迭代向量所有對(duì)應(yīng)分量的比值收斂到1.定義 9.2.1 上述由已知非零向量ZQ及矩陣A的乘幕Ak構(gòu)造向量序列zk來計(jì)算A的 按模最大特征值1及相應(yīng)特征向量的方法稱為幕法（power method），其收斂速度由比值丫 = |為/人|來確定，戈越小，收斂越快.注 由幕法的迭代過程（9.2.3）容易看出，如果 卩 q：1 （或1

9、），那么迭代向量zk的 各個(gè)非零的分量將隨著k-:而趨于無窮（或趨于零），這樣在計(jì)算機(jī)上實(shí)現(xiàn)時(shí)就可能上溢（或下溢）.為了克服這個(gè)缺點(diǎn)，需將每步迭代向量進(jìn)行規(guī)范化：記yk二Azk廠yk1), yk2),|, ykn)若存在yk的某個(gè)分量ykjQ），滿足y：jQ）二徑務(wù)ykj），則記max（ yk） =ykjQ）.將yk規(guī)范化yk；max( yj,這樣就把zk的分量全部控制在-1,1中.例如，設(shè)yk=（-2, 3, Q, -5, 1）T，因?yàn)閥k的所有分量中，絕對(duì)值最大的的是-5，所以max（yQ _ -5，故Zk二y max（yQ =（Q.4, -Q.6, Q, 1, -Q.2）T.綜上所述，得

10、到下列算法：算法 9.2.1 （幕法）設(shè)A是n階實(shí)矩陣，取初始向量ZQ，Rn，通常取ZQ=（1,1,| l|,1）T， 其迭代過程是：對(duì)k =1,2,川，有yk二Amk=max( yQ,Xii =21,2,川，n，(9.2.8)(9.2.9)(j)6定理 9.2.1 對(duì)式（9.2.9）中的序列 和有7其收斂速度由 Y=再/人確定. 證明由迭代過程(9.2.9)知Zk二y/mk二AZk j/mAykjmkmk二A=|( = AkZo | min其中Z0冷Xi.若1 =0，則由(9.2.3)知：AkZ0二1kid(9.2.11)得nna-x1 7aii =2于是4Hkmax x1(9.2.10)n

11、j 、a1x書ai-1k1Ximaxmk二max( yQ二ii =2pmm/nk1 Xi“丿Jaii =2na1x1(9.2.15)Xi(9.2.16)4. mg-nk1ct1%+ 遲aixi,代入式-i-宀丿J(nZn、/k、max口兀+瓦iXi =2嚴(yán)J(9.2.12)Ximax x1(9213)AykjA2Z2kA Zoyk= Azk 1kmax( y) max(AZk/) max(A %)二Ak磯max( Akz).(9.2.14)xi)(9.2.11)Xir 為亠二i=28由式(926)和式(927)知：上述收斂速度由 Y =卩吃/州確定.證畢!例 9.2.1 用幕法求方陣12 3A

12、= 213335_按模最大特征值及相應(yīng)的特征向量，要求mk-mkc10.解選取初始向量zo=(1,1,1)T，按式(929)迭代，結(jié)果見表 9.2.1.表 9.2.1kZkmkg -叫1(0.545455, 0.545455, 1)T112(0.560441, 0.560441, 1)T8.27272.72733(0.559787, 0.559787, 1)T8.36279匯1。24(0.559818, 0.559818, 1)T8.358740乂5(0.559817, 0.559817, 1)T8.35890.2X10，因此，所求按模最大特征值8.3589，相應(yīng)特征向量X - (0.5598

13、17, 0.559817, 1)T.事實(shí)上，A的按模最大特征值.198.358898(，相應(yīng)特征向量為=(0.55981649川，0.55981649 I, 1)T,故所得結(jié)果具有較高的精度.9.2.1.2 幕法的加速從上面的討論可知，由幕法求按模最大特征值，可歸結(jié)為求數(shù)列mk的極限值，其收斂速度由Y=卩2/卅確定.當(dāng)丫=彼/扎1接近1時(shí)，收斂速度相當(dāng)緩慢.為了提高收斂速度，可以利用外推法進(jìn)行加速因?yàn)樾蛄衜k的收斂速度由丫=卩2/州確定，所以若mk收斂，當(dāng)k充分大時(shí)，則有mk-打=0，或改寫為-1、kmu*,其中c是與k無關(guān)的常數(shù)由此可得(9.2.17)這表明幕法是線性收斂的由式(9217)知

14、9g 1九1mk九1mk 2 1mk 1 - 1由上式解出1，并記為mk 2，即機(jī)22mk 2mk - mk 1(mk1k)2,(9.2.18)mk 2-2mk 1mkmk 2- 2叫1mk這就是計(jì)算按模最大特征值的加速公式將上面的分析歸結(jié)為如下算法：算法 9.2.2 (幕法的加速) )設(shè)A是n階實(shí)矩陣，給定非零初始向量ZoRn，通常取Z =(1,1,川,1T.對(duì)k=1,2，用迭代式y(tǒng)k= Azkj, mk=max( yk), Zk= yk/mk求出mi, m2及Z|, Z2.再對(duì)k二3,4, |，迭代過程為Yk=Azki,mk=max( y)(mik-m2)2(9219)削 _mk _2 -

15、,耳二yjk.當(dāng)mk-mk；：，(；是預(yù)先給定的精度)時(shí)，迭代結(jié)束，并計(jì)算Zk；否則繼續(xù)迭代，直至滿足迭代停止條件mk- mk有關(guān)加速收斂技術(shù)，讀者請(qǐng)參考文獻(xiàn)11.8.2.2反幕法及其加速反幕法是計(jì)算矩陣按模最小特征值及相應(yīng)特征向量的迭代法，其基本思想是：設(shè)矩陣A Rnn非奇異，用其逆矩陣A代替A,矩陣A的按模最小特征值就是矩陣A 的按模最大特征值這樣用A4代替A做幕法，即可求出A的按模最大特征值，也就是矩陣A的 按模最小特征值；這種方法稱為反幕法(inverse power method).、d1因?yàn)榫仃嘇非奇異，所以由AXj = Xj可知：A為xi.這說明：如果A的特征值滿足人屮2忙川屮n

16、-1,|人0，(9.2.20)10那么A的特征值滿足且對(duì)A的應(yīng)于特征值 i的特征向量xi也是A的對(duì)應(yīng)于特征值 1 i的特征向量.由上述分析知：對(duì)AJ應(yīng)用幕法求按模最大的特征值1 n及相應(yīng)的特征向量xn，就是求A的按模最小的特征值，n及相應(yīng)的特征向量Xn.算法 9.2.3 (反幕法) )任取初始非零向量Rn，通常取Zo=(1,1,|,1)T.為了避免求A，對(duì)k =1,21，將迭代過程(929)改寫為：Ayk = Zj“ mk=max( yj(9.2.22)k=ymk.仿定理 9.2.1 的證明，可得：定理 9.2.2 對(duì)式(9.2.22)中的序列Uyk 廠v.9.2.3原點(diǎn)平移法為了提高收斂速度

17、，下面介紹加速收斂的原點(diǎn)平移法設(shè)矩陣B = A - pl，其中p是一個(gè)待定的常數(shù)，A與B除主對(duì)角線上的元素外，其他兀素相同.設(shè)A的特征值為1, 2I, n，則B的特征值為r-P，2Pl(,n-P， 且A與B的特征向量相同.9.2.3.1 原點(diǎn)平移法的幕法設(shè)1是A的按模最大的特征值，選擇p，使I2(9221)11-p -p X-p , i =3,4,|,n,及12對(duì)B應(yīng)用反幕法，得：算法 9.2.5 對(duì)k=1,2,Hl，有,一2-P入一P-2 (9.2.24)對(duì)B應(yīng)用幕法，得算法 9.2.4 對(duì)k =1,2川|，有yk=(A- pi)zkJ, mk=max( yjZk= yk；m,(9.2.25

18、)lim m% - p,lim Zk-,k_：:- max x1其收斂速度由(2- P)：(1-P)確定由式(9.2.24)知：在計(jì)算B的按模最大特征值，-p 的過程(9.2.25)中,加速；算法 9.2.4 又稱為原點(diǎn)平移下的幕法 (power method with shift).923.2 原點(diǎn)平移下的反幕法設(shè)n是A的按模最小的特征值，選擇p，使(9.2.26)收斂速度得到入 -P|V人4-P|半廠P，iT，2,|H, n2.及(9.2.27)人- Pi.人n J一P%(9.2.28)若矩陣B = A - pI可逆，貝V B4的特征值為1- P且有1,川，-,幾2P州-P1,21,2川,

19、n2P盒i一p(9.2.29)13(A - pl)yk二zkJ,5k=max(yk),呂=yjmk,其收斂速度由(- p). ( nd - P)確定.的過程(9.2.30)中，收斂速度得n_ P到加速.算法 9.2.5 又稱為 原點(diǎn)平移下的反幕法( (inverse power method with shift).定義 9.2.8 原點(diǎn)平移下的幕法與原點(diǎn)平移下的反幕法統(tǒng)稱為原點(diǎn)平移法.注有的資料上的原點(diǎn)平移法專指原點(diǎn)平移下的反幕法；而有的資料上的反幕法指的就是原點(diǎn)平移下的反幕法.原點(diǎn)平移下的反幕法(算法 9.2.7)的主要應(yīng)用是：已知矩陣的近似特征值后，求矩陣的特征向量.其主要思想是：若已知

20、A的特征值 韋的近似特征值為 叱，則A-*ml的特征值就是- .m(i =1,2j|,m)，其中i(i =1,2| ,m)是A的特征值.而按模最小的特征 值是- $m，相應(yīng)的特征向量與A的特征向量相同.利用公式(9.2.30)可求出曲，皐,并且有只要近似值 叱適當(dāng)好，收斂速度就很快，往往迭代幾次就能得到滿意的結(jié)果.例 9.2.2 已知特征值的近似值：Y 二-0.3589，用原點(diǎn)平移下的反幕法求方陣12 3A= 213335一得對(duì)應(yīng)特征值的特征向量.解取p = -0.3589，對(duì)矩陣(9230)kmmkkimzk二Ximax x-i (9.2.31)由式(9.2.28)知：在計(jì)算B 的按模最大特

21、征值nfmHki休)確定.于是當(dāng)k充分大時(shí)，可取其收斂速度由1141.358923A pl = A+ 0.35891 =21.35893335.3589迭代公式(9.2.30)中的yk是通過解方程組(A - pI ) yk二Zkj100330.6667100-0.64110.4530-11 -00P( A -1.26791) y Pz，即LUyPz.令Lvk二Pzk4, Uyk二Vkj選取Z，使Uy LPz=(1,1,1)T，得y1-(290929.45, 290927.56, -325732.90)T，mb = max( yj = -325732.90，乙=比伽=(-0.8932, -0.8

22、931, 1)T.1由Uy?= LPR得y2=(-845418.49, -845418.49, 946558.42)T，m2= max( y2) = 946558.42，Z2二y2.mi2=(-0.8932, -0.8932, 1)T.由于z,與Z2的對(duì)應(yīng)分量幾乎相等，故A的特征值為1 10.35890.35889894354117,m2946558.42相應(yīng)的特征向量為x =z2=(-0.8932, -0.8932, 1)T.而矩陣A的一個(gè)特征值為=4 - 19二0.358898943540674 |，相應(yīng)的特征向量為(-0.89315,求得的.為了節(jié)省工作量，可先將A _ pI進(jìn)行 LU

23、分解.在 LU 分解中盡量避免較小的Urr當(dāng)除數(shù)，通?？梢韵葘?duì)矩陣A-pI的行進(jìn)行調(diào)換后再分解.為此，我們可用0 0 1P = 010乘A- pl后再進(jìn)行 LU 分解，即10 00P (A pl )= 0L111.358921.358933*_33=25.35891.358935.35891.35893=LU5.3589-0.5726-3.07 10上15-0.89315, 1)T，由此可見得到的結(jié)果具有較高的精度.9.3 QR 算法上一節(jié)我們介紹了求矩陣特征值的幕法和反幕法.幕法主要用來求矩陣的按模最大特征值，而反幕法主要用于求特征向量本節(jié)將介紹幕法的推廣和變形一一 QR 算法，它是求一般中

24、小型矩陣全部特征值的最有效的方法之一，其基本思想就是利用矩陣的QR 分解矩陣A：= Rn n的QR分解就是：用 Householder 變換將矩陣A分解成正交矩陣Q與上三角矩 陣R的乘積，即A二QR.下面首先介紹 Householder 變換.9.3.1豪斯霍爾德變換定義 9.3.1 設(shè)B = (bj)nn是n階方陣，若當(dāng)i1時(shí)，bj=0，則稱矩陣B為上Hessenberg 矩陣(Hessenberg matrix)，又稱為準(zhǔn)上三角矩陣，它的一般形式為川bm!川b2nIIIb3n.(9.3.1)：:時(shí)，Af 收斂到一個(gè)以 A 的特征值i(i=1,2l(, n)為主對(duì)角線元素的上 三 角矩陣.首

25、先介紹用 QR 分解( (QR decomposition 或 QR factorization);即 用Householder 變換將矩陣A分解成正交矩陣Q與上三角矩陣R的乘積，即A = QR.9.3.2.1 QR 分解H2= I- JU2U2011由此得正交矩陣10 00 1I2101 00H2 -00-0.99010.1415J000.14150.9899-1-300-32.3333-0.47140=Q2A2Q20-0.47141.5733-1.5000-00-1.50000.50001000 -0.66670.23570.7070使得QTAQ.02.3333-0.4714-0.4714

26、1.1667-1.5000-1.50000.50000-0.66670.2357-0.70710-0.3333-0.94290.000122算法 9.3.2(QR 分解) )23第 1 步 記A的第 1 列為x(a1(l),a211)|,ani1)T,A = A=佝(叫n利用式(934):fn2=sg詢(1) Q (a )Z 丿3 = %其中e =(1,0川,0)TERn,二嚴(yán)1(匚1事),H1= Ir1Ju1u1r,構(gòu)造出的H1是n階對(duì)稱正交矩陣，使得H1 &“丁待，從而有于是有二Hn4代廠Hn4Hn代二HnnJ!( (H1A二H.斗HH 4令R =An,Q二H1H2山Hn4，其中Q

27、是對(duì)稱正交矩陣，則R = QA.因?yàn)镼是對(duì)稱正交1汀10A2= H1A-)= *-0af?III盤a22)FIIIa22nIFan2IIIr(2)ann第 2 步記X2=(a2?, a32),|,an?)T，同理可構(gòu)造出n -1階對(duì)稱正交矩陣H2，使得H2x2=2e2，其中二2- sgn(a22) |二. 二n皿2，ai2i=2丿e2=(1,0川|,0)TR巴1)若記H2 =山,它仍是對(duì)稱正交矩陣，于是有IH2丿-CT1a；? anIH06a23)IHa23)=H2A2=0+0(3)a33IHa3n+0r0(3)an3IHf(3)ann如此繼續(xù)下去，直到完成第n -1步后，得到上三角矩陣24矩

28、陣，所以得注 1 若A非奇異，則上三角矩陣R也非奇異，從而R的主對(duì)角線元素不為零例 9.3.3 求矩陣1 0 -1A= 214-2 30的 QR 分解A =QR，并使矩陣R的主對(duì)角線上的元素都是正數(shù)。解對(duì)A運(yùn)用算法 9.3.2.第 1 步記A二A,洛=(1,2, -2)T，則 sig n( 1)一1222(-2)2= 3,P1=r1(耳十玄們)=3(3十1) =12,(3, 0, 0)T=(4, 2,-2)T,e(1,0, 0)T,-3A2=H1A1=0.0第 2 步 記x(5 3, 73)T，二2二sign(5 3K (5 3)2(7 3)2二2.86744,鶴二6(6 a2?) =2.86

29、744(2.86744 - 53) =13.0013,業(yè)二X2二2& =(5.3, 7 3)T(2.86744, 0)(4.53411, 2.33333)T,-1!|4.53411(4.53411,2.33333)13.00132.33333/-0.58124 -0.81373-0.813730.58124-注 2 若要求R的主對(duì)角線元素取正數(shù)，則A的QR分解是唯一的5 二 二恂=(1, 2,-2)T1TH I u1u1100101112412-2 J(4, 2, -2)=1-1-2_2-221211,2-7 310.3.2 31 010 1一2二I -次Uj253 1.33333A3

30、= H2A2= 0-2.867440 0-10 0為了使R的主對(duì)角線上的元素都是正數(shù)，取H3=0-10，顯然H3是正交陣,0 0 -1且3-1.33333 2.33333A4= H3A3=02.867442.47995.衛(wèi)02.32494一令3-1.333332.33333R =傀=02.867442.47995，002.32494 _0.333330.15499Q=H1H2H3=0.666670.65874-0.666670.736239.322 QR 算法算法 9.3.3（QR 算法）第 1 步 令A(yù)二A，利用算法 9.3.2 將A1進(jìn)行QR分解，得A二QR，其中Q“為正交矩陣，&

31、為上三角矩陣；然后將Q1與R1逆序相乘，得A?二R1Q1.1 1因?yàn)轸薅1 A，所以有A：二RQ二Q1AA，即A2與A相似.第 2 步 以A?代替A1，再作QR分解，得A?二Q2R2，其中Q2為正交矩陣，R2為.-100 *10!00-0.58124-0.81373H2一0-0.813730.58124H2-2.33333-2.47995-2.32494-0.929980.34874-0.11625了解了QR分解后，下面介紹QR算法(QR algorithm).26上三角矩陣；再將Q2與R2逆序相乘，并記27A3= R2Q2.1 1因?yàn)镽2=Q2_A2，所以A3= R2Q2= Q2A2Q2，

32、即A3與A2相似.依此類推，可得QR算法公式：對(duì)k =1,21,(9.3.10)因?yàn)镽k 1- Qk丄Ak 1,所以Ak= RkjQk J- Qk 1AkjQk,即Ak與Ak J相似.故序列：Ak !這里，我們不加證明地給出 QR 算法收斂的充分條件:定理 9.3.9 (QR 算法的收斂性) )設(shè)A二(aj- Rnn,是由 QR 算法產(chǎn)生的矩陣序 列，其中代=(a(k).若(1) A 二A的特征值打(i =1,2,|(, n)滿足入卜|碼|卄九 n| 0 ；A =P二DP,其中D二diag1,，2ln)，且P有三角分解P = LU(L是單位F 三角矩陣，U 是上三角矩陣)，則即代攵斂到一個(gè)以A

33、的特征值打(i h,2l(, n)為主對(duì)角線元素的上三角矩陣推論 若矩陣A Rn n是對(duì)稱矩陣，且滿足定理9.3.9 中的條件，則由 QR 算法產(chǎn)生的矩陣序列 沖收斂到對(duì)角矩陣D = diag( 1,M,.9.3.3帶原點(diǎn)位移的QR算法QR 算法收斂速度是線性的，需要提高收斂速度.經(jīng)分析知：定理 9.3.9 中l(wèi)im an?kSC1，當(dāng)n越小，收斂速度越快這啟發(fā)我們把 9.2 節(jié)介紹的原點(diǎn)平移法加速幕法和反幕法的思想運(yùn)用到QR 算法的收斂加速，這就是下面將要介紹的帶原點(diǎn)位移的 QR 算法(QR algorithm with shift).因?yàn)榍笊?Hessenberg 矩陣的特征值比一般的方陣

34、簡(jiǎn)單，所以在下面的算法中首先將所考察的矩陣化成上 用位移加速方法進(jìn)行加速lim a(k)k_ -0, i j,i, i = j,(9311)的速度依賴于比值Yn=丸nHessenberg 矩陣，然后再2829其中；是預(yù)先給定的精度.算法 9.3.4(帶原點(diǎn)位移的 QR 算法) )第 1 步 將矩陣A化成上 Hessenberg 矩陣A1.第 2 步 對(duì)k =1,2川 I，給定原點(diǎn)位移數(shù)列d，進(jìn)行迭代加速A f I gRk(QR分解),.Ak 1 =RkQ -Qk 1 & 1Sk1，(9312)并稱由Ak到Ak“的變換為帶原點(diǎn)位移的 QR 變換( (QR transformation w

35、ith shift)，且記=(a(k)Rn.其中g(shù) 的選取方法主要有(k)(a)取Sk= ann；(b)取Sk為Ak的右下角主子矩陣a(k)n d, nVa(k)n,nd(k)的 2 個(gè)特征值中最靠近ann的那一個(gè)注 1 由算法 9.3.4 產(chǎn)生的Ak 1與Ak相似，即Ak 1二QkAkQk二 QjARQk;且它們都 是上HeSSenberg 矩陣.注 2 計(jì)算實(shí)踐表明，取方法(b)中SR的算法比取方法(a)中SR的算法收斂速度更快.特別 是對(duì)稱矩陣，帶方法(b)中SR位移的 QR 算法是無條件收斂的，且收斂階為3.注 3 在迭代過程中,當(dāng)ankn：-o時(shí)，由定理 9.3.4 知：anR)的近

36、似特征值.而lim a? = n的速度依賴于比值nk)二k_ i:收斂速度越快；因此平移值取sk=ank)顯然是一注 4 在迭代過程中，當(dāng)ani/肚0時(shí)，取a(k)an,na(k)：-n，因此可取(n-Sk)；ann)作為Ank)越小，的 2 個(gè)特征值作為A的特征值-n4, n的近似值.注 5 判別an仁和諜心約等于 0 的標(biāo)準(zhǔn)可取為,antz (k)min、an 4,nr,antz30注 6 當(dāng)求得矩陣A的特征值n 或n4, n 后，可以將AR作降階處理：對(duì) 入左上角的31n _1階或n - 2階主子矩陣?yán)^續(xù)運(yùn)用帶原點(diǎn)位移的QR 算法，以求A的其他特征值，這樣可以大大節(jié)省運(yùn)算量9.4 Jac

37、obi 方法上一節(jié)我們介紹了 Householder 變換將矩陣化為上 Hessenberg 矩陣的方法.如果矩陣是 實(shí)對(duì)稱矩陣，用平面旋轉(zhuǎn)變換將矩陣化為上Hessenberg 矩陣比用 Householder 變換更好.本節(jié)介紹的 Jacobi 方法就是這種方法.它主要用于求實(shí)對(duì)稱矩陣的全部特征值和特征向量.Jacobi 方法(Jacobi method)的基本思想是：對(duì)實(shí)對(duì)稱矩陣A = (a人n，定存在正交矩陣Q，使1TQ AQ二Q AQ = D,(941)其中D二diag( id川，n)，j (j =1,2|,n)就是矩陣A的特征值，而正交矩陣Q的 第j列就是對(duì)應(yīng)于j的特征向量.由此可見

38、，Jacobi 方法的實(shí)質(zhì)和關(guān)鍵就是找一個(gè)正交矩陣Q，將A化為對(duì)角矩陣.9.4.1 Jacobi方法定義 9.4.1 設(shè)A= (aj)n n是n階實(shí)對(duì)稱矩陣,(i )G(i, j門)二為旋轉(zhuǎn)矩陣( (rotation matrix)，或 Givens 矩陣(Givens matrix)，簡(jiǎn)記為Gij.對(duì)A進(jìn)行的變換(943)稱為 Givens 旋轉(zhuǎn)變換(Givens rotation).注 1 Give ns 矩陣是在n階單位矩陣I的第i行第i列、第i行第j列、第j行第i列、n階矩陣cos：III1HIIII(i)(9.4.2)IIIHI1IIIcos：(j)GijAGT32第j行第j列的交叉

39、的位置上分別換上rii二COST、rij二sin、5 二-sin、5 二cos而33成的注 2 Give ns 矩陣是正交矩陣，變換(943)是正交相似變換注 3 Jacobi 方法就是通過一系列Give ns 旋轉(zhuǎn)變換，把A化為對(duì)角矩陣，從而求得特征值及相應(yīng)的特征向量的方法.因此，Jacobi 方法也稱為平面旋轉(zhuǎn)法(plane rotationmethod).下面介紹具體介紹將n階實(shí)對(duì)稱矩陣化為對(duì)角矩陣的Jacobi 方法.設(shè)A= (aj)nn是n階實(shí)對(duì)稱矩陣，記(1)TA1=(aij)n n= GijAGij.(944)因?yàn)門TTTAipGjAGj = GjAGj =A1,所以Ai仍是對(duì)稱

40、矩陣.通過直接計(jì)算可得基文斯(J. W. Givens 1910 年 12 月 14 日-993 年 3 月 5 日)是美國數(shù)學(xué)家，計(jì)算機(jī)領(lǐng)域的先驅(qū)之一a(1)=ancos2日 +ajjsin2日+ 2丙cos日sin, a(1)=aiisin2v a” cos2二-2aijcos：sin二a(1)=a(1)=a“ cos日+ajisin日，i, j,a(1)二a二-aHsin二-aicos,l G, j,alm =aml =aml,m,_ i,j,a =a(i1-2(ajj-aii)sin 28+aq(cos28-sin28).不難看出，A經(jīng)過Gj作用后，與A相比，只有A1的第i行、第i列、

41、第j行、第j列的 元素發(fā)生了變化，而其他元素與A的相同.特別地，若aj =0，由(9.4.5)式中最后的式子知：若取二滿足關(guān)系式可使=a=0，也就是說，用Gij對(duì)A進(jìn)行變換，可將A的 2 個(gè)非主對(duì)角線元素aij和aji化為零.Jacobi 方法的一般過程是：記A0二A，選擇A的一對(duì)最大的非零非主對(duì)角線元素aj和aji，使用 Give ns 矩陣Gj對(duì)A作正交相似變換得A1，可使A1的這對(duì)非零非主對(duì)角線元 素a(1)a(1)二0.再選擇A1的一對(duì)最大的非零非主對(duì)角線元素作上述正交相似變換得A2，可使A2的這對(duì)非零非主對(duì)角線元素化為零如此不斷地做下去，可產(chǎn)生一個(gè)矩陣序列A = A0, A1JII,

42、 Ak,UI.注 雖然A至多只有n(n_ 1)2寸非零非主對(duì)角線元素，但是不能期望通過n(n -1) 2次變換使A對(duì)角化.因?yàn)槊看巫儞Q能使一對(duì)非零非主對(duì)角線元素化為零，例如，(9.4.5)cotR引21 -tan22ta nTt _ Tt,44(9.4.6)34aj和aji化為零但在下一次變換時(shí)，它們又可能由零變?yōu)榉橇悴贿^可以證明，如此產(chǎn)生的矩陣序列A = A0, A1J|, Ak,|將趨向于對(duì)角矩陣，即Jacobi 方法是收斂的而這個(gè)對(duì)角矩陣的主對(duì)角線元素就是矩陣A的特征值用 Jacobi 方法求矩陣A的特征值的 步驟為：(1)記A二A，在矩陣A中找出按模最大的非主對(duì)角線元素aij，取相應(yīng)的

43、 Give ns 矩陣Gj，記為Gi= Gj；2 2(2)由條件佝)-aii)sin 2門-2aj(cos二- sin R = 0,定出sin,cos.為避免使 用三角函數(shù)，令daii _ajjd= 二，2aj,匕川)，n) =limQ：AQk二QTAQ.k_sc因?yàn)镼是正交矩陣，所以QT= Q，若記Q = (Q(YQ(川)Q，其中Q(j)(j -1.2JH, n)是Q的第j列，則有A Q Q, D即A(Q(Q,)川Q(n 2=,Q Q) f|(QnD,，得，)(AQAQ，川,AQ(n)=GQ崩Q(2),川，九nQ(n)，亦即AQ一jQ(j), j =1,2,川，n.這說明正交矩陣Q的第j列就

44、是對(duì)應(yīng)于 j 的特征向量,例 9.4.1 利用 Jacobi 方法求矩陣-_1-201A = -2-11013一的全部特征值和特征向量，要求A.的所有非主對(duì)角線元素的絕對(duì)值：:：；=0.1.解 記A =A= (aj)3 3，因?yàn)閍# =-2是A中所有非主對(duì)角線兀素中絕對(duì)值最大的元素，取相應(yīng)的 Give ns 矩陣G12.因?yàn)閍 =1, a20)= T，所以d =(a1a20)/2a1；)= -0.5，t = t an = sg d/(| )d ) = - 0. 6 1 8 0 34cos ) - (1 t2)2=0.850651, sin -1 co - -0.525731,于是因?yàn)?： ：

45、：2n(n1)k 1E(Ao)十n)：1，所以當(dāng)k.E(Ak.i) 0.故k 1E( A).(9.4.15)j =1,2,)|l,n.證畢!為E(A爲(wèi)(A山37則A2二QTAQ2.中絕對(duì)值最大的元素，取相應(yīng)的Give ns 矩陣G13.因?yàn)閍：?=2.236068, a32?=3.134730，則AAQ1.cos日sin 0010.850651-0.52573101-s in日cos 60=0.5257310.8506510001一002.2360680-2.2360680.850651A1 = (aj)3 3.L-0.5257310.850651Q1= Q0G：= IG0.850651-0.

46、5257310.5257310.850651010 ,1用A1代替A因?yàn)閍2；）=0.850651是A1中所有非主對(duì)角線元素中絕對(duì)值最大的元素，取相應(yīng)的Give ns 矩陣G23.因?yàn)閍2； =-2.236068, a)=3,所以d =(a22)（1）-a33 2a：3)= -3.077683；t =tan - - sgncos：-(1 t2)2=0.987688,sin v - tcos：- -0.156434,_100 1_1000cos日sin日=00.987688-0.156434衛(wèi)-sin Bcos日00.1564340.987688一2.2360680.082241-2.37079

47、8.0.519258記=A2=(aj)3x3,3.134730Q2=Q1GT二0.850651-0.5257310.5192580.840178-0.1564340.0822420.133071 ,0.987688用A2代替A1，重復(fù)上述過程:因?yàn)閍；?-0.519258是A2中所有非主對(duì)角線元素0-0.525731記記=G1,12二.1 d2=-0.158384,0.082241-0.519258G23二記G238所以390.0747990記-2.3707880.034190 = A (aj)3.-0.0341903.3720780.8078510.519258Q3=Q2GJ =-0.422

48、828 0.8401780.410602-0.156434則A3二Q3AQ3.因?yàn)锳3的所有非主對(duì)角線元素的絕對(duì)值：:：；=0.1，所以迭代停止.此時(shí)A的特征值 的近似值分別為A3的對(duì)角元素： 1.998721,2：-2.370788,3:3.37208.相應(yīng)的特征向量的近似值分別為x1(0.807851, -0.422828, 0.410602) k1(1.967479, -1.029776, 1)T,x2(0.519258, 0.840178, -0.156434)T= k2(-3.319342, -5.370814, 1)T,X3：(-0.278834, 0.339584, 0.8982

49、95)T=&(-0.310404, -0.378032, 1)T,其中匕=0.410602, k2=-0.156434, k3=0.898295.A的特征值的精確值為2,2=丄-旦二-2.37228132川，3=133= 3.37228132川.22 2 2相應(yīng)的特征向量為花=(2, -1, 1)T,x2=(-3.186140|1|, -5.372281川,1)T,x3=(-0.313859丨1|,0.372281川,1)T.從例 9.4.1 可以看出，即使迭代的次數(shù)不是很多，迭代矩陣Ak的所有非對(duì)角元素的絕對(duì)值并不是很小時(shí)，用Jacobi 方法求得的結(jié)果精度都比較高，因此它是求實(shí)對(duì)稱

50、矩陣的全部特征值和特征向量的一個(gè)較好的方法.d =(a1?-a3l)r/2a1(2)=0.865333；t =tan v - sgn(d)d +Jl+d2)=0.457089,COST - (1 t2) J2=0.909493, sin v -tcosv - 0.415720,于是COST0G13=01-sin日0sin0.9094930 = 0cos日-0.4157200.4157200.9094931.998721G3A2GI =0.074799 0-0.2788340.3395840.898295409.4.3改進(jìn)的Jacobi方法由于每次旋轉(zhuǎn)變換前選非零非主對(duì)角線元素的最大值很費(fèi)時(shí)間，

51、為此介紹如下改進(jìn)方 法.第一種方法：把非主對(duì)角線元素按照行的次序印23,川，ain,a23, a24, I u ,a2n,川，ann依次化為零，稱為一次掃描.一次掃描后，前面已化為零的元素可能成為非零元素，需要再次掃描.這一方法稱為循環(huán) Jacobi 方法，這種方法的缺點(diǎn)是：一些已經(jīng)足夠小的元素也 作化零處理，浪費(fèi)了時(shí)間.第二種方法：首先對(duì)實(shí)對(duì)稱矩陣A計(jì)算n二nI2V=|2ZZ ai2，(9.4.16)I yj4)設(shè)置閥值Vi=v.n，按 盹仔|(,ain,a23,a24,l(),a2n,)l(,an,n的順序進(jìn)行掃描.若a0|v,，則選取旋轉(zhuǎn)矩陣Gj作旋轉(zhuǎn)相似變換將aj和aji化為零；否則讓

52、a0過關(guān)， 即不進(jìn)行旋轉(zhuǎn)相似變換將其化為零.因?yàn)槟承┙^對(duì)值小于V的元素的絕對(duì)值可能在后面的旋轉(zhuǎn)變換中增長(zhǎng)，所以應(yīng)進(jìn)行多次掃描，直到Ai的所有非零非主對(duì)角線元素的絕對(duì)值都小于Vi.再設(shè)置閥值v2vi；n，重復(fù)上述過程，直到達(dá)到精度要求，即V為止(其中；是指定的精度).這種方法稱為 限值 Jacobi 方法.9.5 算法程序9.5.1幕法%幕法% A是方陣， e是精度， 輸出向量 V是A的最大特征值 MaxEig對(duì)應(yīng)的特征向量 fun ctio n PowerMethod(A,e)V = on es(size(A,1), 1);for i=1:2V = A * V;M(i) = norm(Y, I

53、 nf);V = Y / M(i);endwhile abs(M(i)-M(i-1) ei = i + 1;Y = A * V;M(i) = norm(Y , Inf);V = Y / M(i);endMaxEig = M(i);disp(sprintf(最大特征值 MaxEig %.6fn 相應(yīng)的特征向量,MaxEig)V=V ,end%反冪法41例 9.5.1 用幕法求方陣_123A= 213.335_按模最大特征值及相應(yīng)特征向量，要求|mk：10J.解 在 MATLAB 命令窗口中輸入：A=1 2 3;2 1 3; 3 3 5; e=10A-3; PowerMethod(A,e) 回車，

54、輸出結(jié)果：最大特征值 MaxEig =8.358906相應(yīng)的特征向量V =0.55980.55981.00009.5.2幕法的加速%幕法的加速%A 是方陣，e 是精度，輸出向量V 是 A 的最大特征值 MaxEig 對(duì)應(yīng)的特征向量fun ctio n Acc_ PowerMethod (A, e)V = on es(size(A,1), 1);for i=1:2V = A * V;M(i) = norm(Y, I nf);V = Y / M(i);endfor i=3:4Y = A * V;M(i) = norm(Y , Inf);M2(i) = M(i-2) - (M(i-1) - M(i-

55、2)F2 / ( M(i)-2*M(i-1)+M(i-2);V=Y/M2(i);endwhile abs(M2(i)-M2(i-1) ei = i + 1;Y = A * V;M(i) = norm(Y , Inf);M2(i) = M(i-2) - (M(i-1) - M(i-2)F2 / ( M(i)-2*M(i-1)+M(i-2);Y = Y / M2(i);endMaxEig = M2(i);disp(sprintf(最大特征值 MaxEig %.6fn 相應(yīng)的特征向量,MaxEig)V,end9.5.3反幕法% A 是方陣，e 是精度，輸出向量 V 是 A 的最大特征值 MaxEig

56、 對(duì)應(yīng)的特征向量 fun ctio n InversePower(A, e)V = on es(size(A,1), 1);L U = lu(A);for i=1:242Y = U (LV);M(i) = norm(Y, I nf);Vector = Y / M(i);endwhile abs(M(i)-M(i-1) ei = i + 1;Y = U (LV);M(i) = norm(Y , Inf);V = Y / M(i);endMi nEig = 1/M(i);disp(sprintf(該矩陣按模最小的特征值是：%.6fn 相應(yīng)的特征向量是：,MinEig)V,end9.5.4原點(diǎn)平移下

57、的反幕法%原點(diǎn)平移下的反幕法% A 是方陣，e 是精度，%輸出向量 V 是 A 的特征值 Eig 對(duì)應(yīng)的特征向量，Appr 是特征值 Eig 的近似值fun ctio n Tranln versePower(A, Appr, e)V = on es(size(A,1), 1);B = A - Appr * eye( size(A,1);L U = lu(B);for i=1:2V = U (LV);M(i) = norm(Y, I nf);V = Y / M(i);endwhile abs(M(i)-M(i-1) ei = i + 1;Y = U (LV);M(i) = norm(Y , In

58、f);V = Y / M(i);endEig = Appr + 1/M(i);disp(sprintf(該矩陣的特征值是：%.6fn 相應(yīng)的特征向量是：,Eig)V=V,end例 9.5.2 已知特征值的近似值：Y 二0.3589，用原點(diǎn)平移下的反幕法求方陣123A= 21333 5的對(duì)應(yīng)特征值的特征向量.解 在 MATLAB 命令窗口中輸入：A=1 2 3;2 1 3; 3 3 5; Appr=-0.3589; e=10A-3; TranlnversePower(A,Appr,e) 輸出的結(jié)果是：該矩陣的特征值是：-0.358899相應(yīng)的特征向量是：例 9.5.3 設(shè)矩陣43V =0.893

59、10.8931-1.00009.5.5 Householder變換%Householder 變換%參數(shù) X 是列向量function Result = Householder ( X ) d = sig n( X(1) * norm(X, 2);e1 = eye( size(X,1), 1 );u = X + d*e1;p = d * ( d+X(1);H = eye(size(X,1) - pA-1 * (u*u);Result = H;end9.5.6用Householder變換將一個(gè)n階方陣化為上Hessenberg矩陣%用 Householder 變換將一個(gè) n 階方陣化為上 Hess

60、enberg 矩陣%輸出的 Q 是所求的正交矩陣，% 輸出的UpH是用Householder變換將n階方陣 A化成的上 Hessenberg矩陣 function UpHesse nberg( A )N = size(A, 1); |Ai = A;Q = eye(N);for i=2:N-1A = Ai(i-1:N, i-1:N);X = A(2:size(A,1), 1); HH = Householder ( X );Qi = eye(i-1) zeros(i-1, N-i+1); zeros(N-i+1, i-1) H;Q = Q * Qi;Ai = Qi * Ai * Qi;endQ, _UpH=Ai,end4412 122 211A