QR方法計算中小型矩陣的全部特征值

1、計算方法課程設計報告學生姓名:學 號：學 院:班 級:題 目:QR方法計算中小型矩陣的全部特征值指導教師： 職稱： 教 授 講 師 實驗師 2015年12月31日目 錄目錄I一、選題背景11.1 QR方法11.2 矩陣的特征值 1二、算法設計12.1 QR方法的理論12.2 基本QR方法 22.3 Householder QR分解 22.4 帶原點位移的QR方法 3三、程序設計及功能說明43.1主要程序以及主要功能43.1.1矩陣約化為上海森伯格矩陣43.1.1 QR方法求實矩陣全部特征值4四、結果分析. 6五、總結及心得體會.19參考文獻.20源程序.21一、選題背景1.1 QR方法矩陣是高

2、等數(shù)學中的常用工具，在很多方面都有重要運用，而矩陣特征值問題在許多領域的研究中有重要的地位，矩陣特征值的一些基本計算方法，研究不同種類矩陣的計算方法和最優(yōu)計算方法.其中求解矩陣的普通方法包括傳統(tǒng)的求法以及初等變換求矩陣的特征值方法；其他的一些優(yōu)化方法包括冪法、反冪法、Jacobi方法、QR方法.在實際的求解矩陣特征值的問題，根據(jù)矩陣的不同特點，選擇最快速的方法求解，從而達到最優(yōu)化解決實際問題。QR分解法是目前求一般矩陣全部特征值的最有效并廣泛應用的方法，一般矩陣先經(jīng)過正交相似變化成為Hessenberg矩陣，然后再應用QR方法求特征值和特征向量。它是將矩陣分解成一個正規(guī)正交矩陣Q與上三角形矩陣

3、R，所以稱為QR分解法，與此正規(guī)正交矩陣的通用符號Q有關。目前QR方法主要用來計算：(1) Hessenberg矩陣的全部特征值問題；(2)計算對稱三對角矩陣的全部特征值問題。1.2 矩陣的特征值關于計算矩陣A的特征值問題，當n=2,3時，我們還可按行列式展開的辦法求，我們還可按行列式展開的辦法求det()=0的根，但當n較大時，如果按展開行列式的辦法，首先求出det()的系數(shù)，再求det()的根，工作量就非常大，用這種辦法求矩陣用這種辦法求矩陣的特征值是不切實際的的特征值是不切實際的，由此需要研究求A的特征值 及特征向量的數(shù)值解法。二、算法設計2.1 QR方法的理論QR方法的理論：對任意一個

4、非奇異矩陣（可逆矩陣）A，可以把它分解成一個正交陣Q和一個上三角陣R的乘積，稱為對矩陣A的QR分解，即A=QR。如果規(guī)定R的對角元取正實數(shù)，這種分解是唯一的。若A是奇異的，則A有零特征值。任取一個不等于A的特征值的實數(shù)，則A-I是非奇異的。只要求出A-I的特征值和特征向量就容易求出矩陣A的特征值和特征向量，所以假設A是非奇異的，不失一般性。2.2 基本QR方法 設A=A1，對A1作QR分解，得A1= Q1R1，交換該乘積的次序，得，由于Q1正交矩陣，A1到A2的變換為正交相似變換，于是A1和A2就有相同的特征值。一般的令A1=A，對k=1，2，3，. A(k)=QkRk （QR分解） A(k-

5、1)=RkQk (迭代定義) 這樣，可得到一個迭代序列Ak，這就是QR方法的基本過程。2.3 Householder QR分解從QR 算法的構造過程可以看到算法的主要計算量出現(xiàn)在QR分解上，如果直接對矩陣A用QR方法求全部特征值，那麼涉及的計算量是很大的，因此應該先對A作預處理。應用中常先對 做正交相似變換將其化為上Hessenberg矩陣H，然后再對H采用QR方法，可以大大減少計算量，這里Hessenberg矩陣也稱為擬三角矩陣，它的非零元素比三角矩陣多了一條次對角線，其形式為： 上Hessenberg矩陣 下Hessenberg矩陣一般矩陣相似約化到Hessenberg矩陣的方法。定理:

6、任取非零向量X=(X1,X2,.,Xn)TRn，可以選擇一個Householder矩陣P，使Px=1式中e1=(1,0,.,0)T是Rn的單位向量， 證： 作u=x+1，用u做一個Householder矩陣P=I-1uuT, 因為 而證畢定理中=(+x1)，為避免+x1出現(xiàn)兩個相似數(shù)相減引起有效數(shù)字損失，實用中常選 這樣可得Householder矩陣的計算公式為u1=x1+=u1 P=I-1uuT2.4 帶原點位移的QR方法設A=A1，對A1-s1I進行QR分解A1-s1I=Q1R1；形成矩陣 A2=R1Q1+s1I=1(A-s1I) Q1+s1I=1A1Q1；求得Ak后，將Ak-skI進行Q

7、R分解Ak-skI=QkRk，k=3，4， 形成矩陣A(k+1)=RkQk+skI=kAkQk.如果令k=Q1Q2Qk, k=RkR2R1,則有A(k+1)= kAk，并且矩陣(A-s1I)(A-s2I) (A-snI) (A)有QR分解式(A)= kk.也可以首先用正交變換(左變換)將Ak-skI化為上三角矩陣，即PnP2P1(Ak-skI)=Rk當A為Householder矩陣或對稱三對角矩陣，Pi可為平面旋轉矩陣，則Ak= P(n-1)P2P1(Ak-skI) 12(n-1)+skI.三、程序及功能說明3.1主要程序以及主要功能3.1.1 矩陣約化為上海森伯格矩陣function k,S

8、k,uk,ck,Pk,Uk,Ak=Householdrer1(A)n=size(A); Ak=A;for k=1:n-2k,Sk=norm(Ak(k+1:n,k)*sign(Ak(k+1,k), uk= Ak(k+1:n,k)+ Sk*eye(n-k,1),ck=(norm(uk,2)2)/2,Pk= eye(n-k,n-k)-uk*uk/ck,Uk=eye(k,k),zeros(k,n-k);zeros(n-k, k),Pk,A1=Uk*Ak;Ak=A1,end 程序功能：通過初等反射矩陣正交相似約化實矩陣為上海森伯格矩陣Ak。來實現(xiàn)矩陣約化為上海森伯格矩陣，以方便使用QR方法求矩陣A的全部

9、特征值。3.1.2 QR方法求實矩陣全部特征值function tzg=qr4(A,t,max1)n,n=size(A); k=0;Ak=A;tzg=zeros(n); state=1;for i=1:n;while(k1)b1=abs(Ak(n,n-1); b2=abs(Ak(n,n); b3=abs(Ak(n-1,n-1);b4=min(b2, b3); jd=10(-t); jd1=jd*b4;if(b1=jd1) sk=Ak(n,n); Bk=Ak-sk*eye(n); Qk,Rk=qr(Bk);At=Rk*Qk+sk*eye(n); k=k+1;tzgk=Ak(n,n);disp(請

10、注意：下面的i表示求第i個特征值，k是迭代次數(shù)，sk是原點位移量，)disp( Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解，At=Rk*Qk+sk*eye(n)迭代矩陣：)i,state=1;k,sk,Bk,Qk,Rk,At,Ak=At;else disp(請注意：i表示求第i個特征值，tzgk是矩陣A的特征值的近似值，k是迭代次數(shù)，) disp( 下面的矩陣B是m階矩陣At的（m-1）階主子矩陣，即At降一階.) i,tzgk=Ak(n,n),tzg(n,1)=tzgk;k=k,sk,Ak;B=Ak(1:n-1,1:n-1),Ak=B;n=n-1;state=1; i=i+1

11、;endendendtzg(1,1)=Ak;tzg=sort(tzg(:,1);tzgk=Akdisp(請注意：n階實對稱矩陣A的全部真特征值lamoda和至少含t個有效數(shù)字的近似特征值tzg如下：)lamoda=sort(eig(A)程序功能：QR方法來計算Hessenberg矩陣或對稱的矩陣的全部特征值至少含t個有效數(shù)字的近似值。四、結果分析實驗結果分析：1.用QR方法求下列矩陣的全部近似特征值，其中精度為；解：（1）先將矩陣轉化為上海森伯格矩陣，MATLAB程序如下： A=1 2 3;1 4 5;0 2 6; k,Sk,uk,ck,Pk,Uk,Ak=Householdrer1(A)結果如

12、下：k = 1Sk = 1uk = 2 0ck = 2Pk = -1 0 0 1Uk = 1 0 0 0 -1 0 0 0 1Ak = 1 2 3 -1 -4 -5 0 2 6然后用最末元位移QR方法求實矩陣k全部特征值.MATLAB程序如下： A=1 2 3;-1 -4 -5;0 2 6; tzg=qr4(A,5,100)結果如下：請注意：下面的i表示求第i個特征值，k是迭代次數(shù)，sk是原點位移量， Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解，At=Rk*Qk+sk*eye(n)迭代矩陣：i = 1k = 4sk = 4.7505Bk = -7.2733 0.1020 -5.

13、3467 -2.5877 -3.9784 -5.5837 0 0.0000 0Qk = -0.9421 0.3352 0.0000 -0.3352 -0.9421 -0.0000 0 0.0000 -1.0000Rk = 7.7199 1.2374 6.9091 0 3.7824 3.4684 0 0 0.0000At = -2.9375 1.4219 -6.9091 -1.2679 1.1870 -3.4684 0 0.0000 4.7505請注意：i表示求第i個特征值，tzgk是矩陣A的特征值的近似值，k是迭代次數(shù)， 下面的矩陣B是m階矩陣At的（m-1）階主子矩陣，即At降一階.i =

14、1tzgk = 4.7505k = 4sk = 4.7505B = -2.9375 1.4219 -1.2679 1.1870請注意：下面的i表示求第i個特征值，k是迭代次數(shù)，sk是原點位移量， Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解，At=Rk*Qk+sk*eye(n)迭代矩陣：i = 2k = 7sk = 0.6930Bk = -3.1366 2.6863 -0.0035 0Qk = -1.0000 0.0011 -0.0011 -1.0000Rk = 3.1366 -2.6863 0 0.0030At = -2.4406 2.6897 -0.0000 0.6900請注

15、意：i表示求第i個特征值，tzgk是矩陣A的特征值的近似值，k是迭代次數(shù)， 下面的矩陣B是m階矩陣At的（m-1）階主子矩陣，即At降一階.i = 2tzgk = 0.6900k = 7sk = 0.6930B = -2.4406tzgk = -2.4406請注意：n階實對稱矩陣A的全部真特征值lamoda和至少含t個有效數(shù)字的近似特征值tzg如下：lamoda = -2.4406 0.6900 4.7505tzg = -2.4406 0.6900 4.75052.用QR方法求下列矩陣的全部近似特征值，其中精度為；解：MATLAB程序如下： A=1 2 2;2 1 2;2 2 1; tzg=q

16、r4(A,5,100)結果如下：請注意：下面的i表示求第i個特征值，k是迭代次數(shù)，sk是原點位移量， Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解，At=Rk*Qk+sk*eye(n)迭代矩陣：i = 1k = 3sk = -0.9883Bk = 5.9418 0.4553 0.2639 0.4553 0.0231 0.0202 0.2639 0.0202 0Qk = 0.9961 0.0764 0.0441 0.0763 -0.9971 0.0034 0.0442 0.0000 -0.9990Rk = 5.9650 0.4562 0.2644 0 0.0117 0.0000 0

17、 0 0.0117At = 5.0000 0.0009 0.0005 0.0009 -1.0000 0.0000 0.0005 0.0000 -1.0000請注意：i表示求第i個特征值，tzgk是矩陣A的特征值的近似值，k是迭代次數(shù)， 下面的矩陣B是m階矩陣At的（m-1）階主子矩陣，即At降一階.i = 1tzgk = -1.0000k = 3sk = -0.9883B = 5.0000 0.0009 0.0009 -1.0000請注意：下面的i表示求第i個特征值，k是迭代次數(shù)，sk是原點位移量， Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解，At=Rk*Qk+sk*eye(n

18、)迭代矩陣：i = 2k = 4sk = -1.0000Bk = 6.0000 0.0009 0.0009 0Qk = 1.0000 0.0001 0.0001 -1.0000Rk = 6.0000 0.0009 0 0.0000At = 5.0000 0.0000 0.0000 -1.0000請注意：i表示求第i個特征值，tzgk是矩陣A的特征值的近似值，k是迭代次數(shù)， 下面的矩陣B是m階矩陣At的（m-1）階主子矩陣，即At降一階.i = 2tzgk = -1.0000k = 4sk = -1.0000B = 5.0000tzgk = 5.0000請注意：n階實對稱矩陣A的全部真特征值la

19、moda和至少含t個有效數(shù)字的近似特征值tzg如下：lamoda = -1.0000 -1.0000 5.0000tzg = -1.0000 -1.00005.00003.用QR方法求下列矩陣的全部近似特征值，其中精度為；解：（1）先將矩陣轉化為上海森伯格矩陣，MATLAB程序如下： A=1 2 3 4;2 3 4 5;3 4 5 6;4 5 6 7; k,Sk,uk,ck,Pk,Uk,Ak=Householdrer1(A)結果如下：k = 2Sk = -0.4549uk = -0.5400 -0.4468ck = 0.2456Pk = -0.1871 -0.9823 -0.9823 0.18

20、71Uk = 1.0000 0 0 0 0 1.0000 0 0 0 0 -0.1871 -0.9823 0 0 -0.9823 0.1871Ak = 1.0000 2.0000 3.0000 4.0000 -5.3852 -7.0564 -8.7277 -10.3989 0.0000 0.4549 0.9097 1.36460.0000 -0.0000 -0.0000 -0.0000（2）然后用最末元位移QR方法求實對稱矩陣全部特征值，因為矩陣非奇異，故MATLAB程序如下： A=1 2 3;-5.3852 -7.0564 -8.7277;0.0000 0.4549 0.9097; tzg=

21、qr4(A,5,100)結果如下：請注意：下面的i表示求第i個特征值，k是迭代次數(shù)，sk是原點位移量， Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解，At=Rk*Qk+sk*eye(n)迭代矩陣：i = 1k = 6sk = -1.4987e-004Bk = -4.6826 11.4681 -3.9762 -0.0002 -0.4637 -0.7271 0 0.0000 0Qk = -1.0000 0.0000 0.0000 -0.0000 -1.0000 -0.0000 0 0.0000 -1.0000Rk = 4.6826 -11.4681 3.9763 0 0.4641 0

22、.7269 0 0 0.0000At = -4.6823 11.4683 -3.9763 -0.0000 -0.4642 -0.7269 0 0.0000 -0.0001請注意：i表示求第i個特征值，tzgk是矩陣A的特征值的近似值，k是迭代次數(shù)， 下面的矩陣B是m階矩陣At的（m-1）階主子矩陣，即At降一階.i = 1tzgk = -1.4992e-004k = 6sk = -1.4987e-004B = -4.6823 11.4683 -0.0000 -0.4642請注意：下面的i表示求第i個特征值，k是迭代次數(shù)，sk是原點位移量， Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR

23、分解，At=Rk*Qk+sk*eye(n)迭代矩陣：i = 2k = 7sk = -0.4642Bk = -4.2181 11.4683 -0.0000 0Qk = -1.0000 0.0000 -0.0000 -1.0000Rk = 4.2181 -11.4683 0 0.0000At = -4.6823 11.4683 -0.0000 -0.4643請注意：i表示求第i個特征值，tzgk是矩陣A的特征值的近似值，k是迭代次數(shù)， 下面的矩陣B是m階矩陣At的（m-1）階主子矩陣，即At降一階.i = 2tzgk = -0.4643k = 7sk = -0.4642B = -4.6823tzg

24、k = -4.6823請注意：n階矩陣A的全部真特征值lamoda和至少含t個有效數(shù)字的近似特征值tzg如下：lamoda = -4.6823 -0.4643 -0.0001tzg = -4.6823 -0.4643 -0.0001五、總結及心得體會在本次課程設計中，雖然我們在完成過程中，遇到這樣那樣問題，通過老師和同學幫助，我們算是基本熟悉QR方法的基本原理；得到QR方法求解矩陣特征值的通用程序完成本次課程設計，使我們了解各種求解矩陣特征值計算方法的之間優(yōu)缺點。參考文獻1 李慶揚，王能超，易大義數(shù)值分析（第5版）北京：清華大學出版社，2008.2 李慶揚，關治，白峰杉 數(shù)值計算原理北京：清華大學出版社，20003 李海，鄧櫻MATLAB程序設計教程北京：高等教育出版社，20024 王萼芳,石生明高等代數(shù)M北京:高等教育出版社2003源程序1 矩陣約化為上海森伯格矩陣function k,Sk,uk,ck,Pk,Uk,Ak=Householdrer1(A)n=size(A); Ak=A;for k=1:n-2k,Sk=norm(Ak(k+1:n,k)*sign(Ak(k+1,k), uk= Ak(k+1:n,k)+ Sk*eye(