數(shù)值計算課程設(shè)計矩陣特征值與特征向量計算

1、矩陣的特征值與特征向量的計算摘 要矩陣是高等代數(shù)學(xué)中的常見工具，也常見于統(tǒng)計分析等應(yīng)用數(shù)學(xué)學(xué)科中.在物理學(xué)中，矩陣于電路學(xué)、力學(xué)、光學(xué)和量子物理中都有應(yīng)用；計算機科學(xué)中，三維動畫制作也需要用到矩陣. 矩陣的運算是數(shù)值分析領(lǐng)域的重要問題.將矩陣分解為簡單矩陣的組合可以在理論和實際應(yīng)用上簡化矩陣的運算.在本論文中，我們主要討論矩陣的特征值和特征向量的計算，我們知道，有很多現(xiàn)實中的問題都可以用到矩陣特征值與特征向量計算的知識，比如，在物理、力學(xué)和工程技術(shù)方面有很多的應(yīng)用，并且發(fā)揮著極其重要的作用.因為這些問題都可歸結(jié)為求矩陣特征值的問題，具體到一些具體問題，如振動問題，物理中某些臨界值的確定問題以及

2、一些理論物理中的問題.在本論文中，我們主要介紹求矩陣的特征值與特征向量的一些原理和方法，原理涉及高得代數(shù)中矩陣的相關(guān)定理，方法主要介紹冥法及反冥法，Jacobi方法和QR算法，并利用MATLAB，VC等軟件編寫相關(guān)算法的程序來求解相關(guān)問題，加以驗證.關(guān)鍵詞： 矩陣；特征值；特征向量；冥法；反冥法；Jacobi方法；QR算法；VC軟件；MATLAB軟件THE CALCULATIONS OF EIGENVALUE AND EIGENVECTOR OF MATRIXABSTRACT The matrix is an usual tool in Advanced Algebra, which also

3、 used by applied mathematics such as Statistics Analysis. In Physics, we can see the important usage of matrix including Electric Circuits, Mechanics, Optics and Quantum Physics. Making three dimension needs matrix in Computer. The arithmetic of matrix is a very important part in Numerical Analysis.

4、 It can simplify the calculation of matrix that we decompose the matrix into several simple parts.In this thesis, we mainly talk about the calculation of eigenvalue and eigenvector of matrix. As we all know, there are lots of realistic problems which need the knowledge of the thesis to solve. We can

5、 see the important usage of matrix including Electric Circuits, Mechanics, Optics and Quantum Physics. It play an important role in these problems inferred above. Because these problems can regarded as the calculation of eigenvalue and eigenvector of matrix, like vibrating problems and critical valu

6、e problems and so on.We primarily introduce the principle and approach of the calculation of eigenvalue and eigenvector of matrix that infer the relevant principle in Advanced Algebra. We mainly talk about iteration methods, Jacobi method and QR method by using MATLAB.Key words: Matrix；Eigenvalue；Ei

7、genvector；Iteration methods; Jacobi method; QR method;MATLAB目 錄1 引言.12 相關(guān)定理。.13 符號說明.24 冥法及反冥法.2 4.1冥法.3 4.2反冥法.85 QR算法.14參考文獻(xiàn).18 附錄.19 1 引言在本論文中，我們主要討論矩陣的特征值和特征向量的計算，我們知道，有很多現(xiàn)實中的問題都可以用到矩陣特征值與特征向量計算的知識，比如，在物理、力學(xué)和工程技術(shù)方面有很多的應(yīng)用，并且發(fā)揮著極其重要的作用.因為這些問題都可歸結(jié)為求矩陣特征值的問題，具體到一些具體問題，如振動問題，物理中某些臨界值的確定問題以及一些理論物理中的問題

8、.在本論文中，我們主要介紹求矩陣的特征值與特征向量的一些原理和方法，原理涉及高得代數(shù)中矩陣的相關(guān)定理，方法主要介紹冥法及反冥法，Jacobi方法和QR算法，并利用MATLAB，VC等軟件編寫相關(guān)算法的程序來求解相關(guān)問題，加以驗證.2 相關(guān)定理定理2.1 如果 是矩陣A的特征值，則有定理2.2 設(shè)A與B為相似矩陣，則 A與B有相同的特征值；若是的一個特征向量，則是A的特征向量定理2.3 設(shè)，則A的每一個特征值必屬于下述某個圓盤之中： 定義2.1 設(shè)A是n階是對稱矩陣，對于任意非零向量x，稱為對應(yīng)于向量x的Rayleigh商.定理2.4 設(shè)為對稱矩陣（其特征值次序記作，對應(yīng)的特征向量組成規(guī)范化正交

9、組，即），則 （對于任何非零向量x）；3 符號說明A:n階矩陣B:n階矩陣I：n階單位陣:矩陣特征值x:實數(shù)域上的n維向量：實數(shù)域上的n維向量：實屬上的規(guī)范化向量 4 冥法及反冥法4.1 冥法冪法是一種計算矩陣的主特征值的一種迭代法，它最大優(yōu)點是方法簡單，適合于計算大型稀疏矩陣的主特征值.設(shè)，其特征值為，對應(yīng)特征向量為即 且線性無關(guān).設(shè)特征值滿足：（即為強占優(yōu)） （4.1.1）冪法的基本思想，是任取一個非零初始向量，由矩陣的乘冪構(gòu)造一向量序列 (4.1.2)稱為迭代向量.下面來分折.由設(shè)為中一個基本，于是，有展開式 （且設(shè)）且有(4.1.3 ) 由假設(shè)（4.1.1）式,則即且收斂速度由比值確定

10、.且有（41.4） 這說明，當(dāng)充分大時，有，或越來越接近特征向量.下面考慮主特征值的計算.用表示的第個分量，考慮相鄰迭代向量的分量的比值.從而是 (4.1.5)說明相鄰迭代向量分量的比值收斂到主特征,且收斂速度由比值來度量,越小收斂越快,但越小收斂越快,但,而接近于1時,收斂可能很慢.定理4.1 （1）設(shè)n個線性無關(guān)的特征向量：（2）設(shè)特征值滿足（3）冪法： ）則 （1）；（2） 如果主特征值為實的重根，即有 又設(shè)A有個線性無關(guān)的特征向量，其中對于任意初始向量則由冪法有 且有 （設(shè)不全為零） 由此，當(dāng)充分大時，接近于與對應(yīng)的特征向量的某個線性組合.應(yīng)用冪法計算的主特征值及對應(yīng)的特征向量時，如果

11、），迭代向量的各個不等于零的分量將隨而趨于無究（或趨于零），這樣電算時就可能溢出.為此，就南非要將迭代向量加以規(guī)范化.設(shè)有非零向量其中表示向量絕對值最大的元素，即如果有草藥則其中為所有絕對值最大的分量中最小指標(biāo). 顯然有下面性性質(zhì)： 設(shè)，則 在定理4.1條件下冪法可改進(jìn)為： 任取初始向量. 迭代： 規(guī)范化： ， （4.16） 于是，由上式產(chǎn)生迭代向量序列及規(guī)范化向量且改進(jìn)冪法計算公式為： 設(shè) 對于 （4.1.7） 下面考查與計算的關(guān)系. 由 且有 （4.1.8） 其中 (1) 考查規(guī)范化向量序列:由(4.1.7)及(4.1.8)式,則有 (2) 考查迭代向量序列:于是， 定理 (改進(jìn)冪法)(1

12、) 設(shè)有個線性無關(guān)特征向量;(2) 設(shè)特征值滿足 且 (3)由改進(jìn)冪法得到(4.1.7)式),則有 (a) (b)且收斂速度由比值確定.實現(xiàn)冪法，每迭代一次主要是計算一次矩陣乘向量，可編一個子程序.例1.用MATLAB編寫冥法程序求矩陣主特征值及近似主特征量用冪法計算下列矩陣的主特征值和對應(yīng)的特征向量的近似向量，精度.并把輸出的結(jié)果真實結(jié)果進(jìn)行比較.解 輸入MATLAB程序 B=1 2 3;2 1 3;3 3 6; V0=1,1,1; k,lambda,Vk,Wc=mifa(B,V0,0.00001,100), V,D = eig (B), Dzd=max(diag(D), wuD= abs(

13、Dzd- lambda), wuV=V(:,3)./Vk,運行后屏幕顯示結(jié)果請注意：迭代次數(shù)k,主特征值的近似值lambda,主特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：k = lambda = Wc = Dzd = wuD = 3 9 0 9 0Vk = wuV = 0.00 0.873 0.00 0.873 1.00 0.873V = 0.55 0.63 0.86 -0.55 0.63 0.86 0 -0.63 0.8734.2 反冥法及位移反冥法(1) 反冪法可用來計算矩陣按模最小的特征值及對應(yīng)的特征向量.設(shè)為非廳異矩陣，特征值滿足對應(yīng)特征向量為線性無關(guān)，則特征求值為特征向量為

14、因此計算的按模最小的特征值的部題就是計算按模最大的特征值部題.對于應(yīng)用冪法迭代（稱為反冪法），可求矩陣的主特征值.反冪法迭代公式：任取初始向量， 1，2， （4.2.1）其中迭代向量可通過解方程組求得:如果個線性無關(guān)特征向量且特征值滿足:則由反冪法(2.11)構(gòu)造的向量序列滿足 且收斂速度由比值確定.（2）應(yīng)用反冪法求一個的似特征值對應(yīng)的特征向量.設(shè)已知的特征值的一個近似值（通常是用其它方法得到），現(xiàn)要求對應(yīng)的特征向量（近似），在反冪法中也可用原點平移法來加速收斂.如果存在，顯然，特征值為對應(yīng)的特征向量.現(xiàn)?。ǖ荒苋。?，且設(shè)與其它特征值是分離的，即即 說明是的主特征值.現(xiàn)對應(yīng)用冪法得到反冪法

15、計算公式：取初始向量 （4.2.2）與定理8證明類似,可得下述結(jié)果.定理10 （1）設(shè)有個線性無關(guān)特征向量即.（2）?。樘卣髦狄粋€近似值），設(shè)存在且則由反冪法迭代公式（2，12）構(gòu)造向量序列滿足：或 且收斂速度由比值 確定.由定理可知，反冪法計算公式（4.2.2）可用計算特征向量.選擇是的一個近似且的特征值分離情況較好，一般很小，所以迭代過程收斂較快，同時改進(jìn)特征值.反冪法迭代公式中是以通過解方程組求得.為了節(jié)省計算量，可先將進(jìn)行三角分解.其中為置換陣，于是每次迭代求相當(dāng)于求解兩個三角形方程組可按下述方法取，即選使回代求解即求得.反冪法計算公式：1分解計算，且保存及信息2反冪法迭代（1） （

16、2） 1）求 求 2） 3）對于計算對稱三對角陣，或計算Hessenberg陣對應(yīng)于一個給定的近似特征值的特征向量，反冪 法是一個有效方法.例2 .用反冪法計算對應(yīng)于近似特征值（精確特征值為）的特征向量 解 取，用部分選主元分解法實現(xiàn)，其中（1）求解 （2）求解 （3）求解特征向量（真解）是 由此，相當(dāng)好的近似.例3.用原點位移反冪法的迭代公式，根據(jù)給定的下列矩陣的特征值的初始值，計算與對應(yīng)的特征向量的近似向量,精確到0.000 1.,解 輸入MATLAB程序 A=1 -1 0;-2 4 -2;0 -1 2;V0=1,1,1;k,lambda,Vk,Wc=ydwyfmf(A,V0,0.2,0.

17、0001,10000)運行后屏幕顯示結(jié)果請注意：因為A-aE的各階主子式都不等于零，所以A-aE能進(jìn)行LU分解.A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：k = lambda = Wc = hl = 3 0.2384 1.0213e-007 0.8000 1.0400 0.2720Vk = V = D = 1.0000 -0.2424 -1.0000 -0.5707 5.1249 0 0 0.7616 1.0000 -0.7616 0.3633 0 0.2384 0 0.4323 -0.

18、3200 -0.4323 1.0000 0 0 1.6367例4. 用原點位移反冪法的迭代公式（5.28），計算的分別對應(yīng)于特征值的特征向量的近似向量，相鄰迭代誤差為0.001.將計算結(jié)果與精確特征向量比較.解 計算特征值對應(yīng)的特征向量的近似向量.輸入MATLAB程序 A=0 11 -5;-2 17 -7;-4 26 -10;V0=1,1,1;k,lambda,Vk,Wc= ydwyfmf(A,V0,1.001, 0.001,100),V,D=eig(A);Dzd=min(diag(D), wuD= abs(Dzd- lambda),VD=V(:,1),wuV=V(:,1)./Vk,運行后屏幕

19、顯示結(jié)果請注意：因為A-aE的各階主子式都不等于零，所以A-aE能進(jìn)行LU分解.A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：hl = -1.00 5.00 -0.00k = lambda = RA1 = 5 1.00 -0.00Vk = VD = wuV = -0.00 -0.86 0.873 -0.00 -0.86 0.873 -1.00 -0.873 0.873Wc = Dzd = wuD = 1.562e-009 1.00 0.00 從輸出的結(jié)果可見，迭代5次，特征向量的近似向量的相

20、鄰兩次迭代的誤差Wc1.379 e-009，由wuV可以看出，= Vk與VD 的對應(yīng)分量的比值相等.特征值的近似值lambda 1.002與初始值1.001的絕對誤差為0.001，而與的絕對誤差為0.002，其中,. 5 QR方法用最末元位移QR方法求下列實對稱矩陣的全部近似特征值，并將計算結(jié)果與全部真特征值比較.其中精度為解 首先保存用最末元位移QR方法求實對稱矩陣全部特征值的MATLAB主程序為M文件，取名為qr4.m.在MATLAB工作窗口輸入程序 A=5 2 2 1;2 -4 1 1;2 1 3 1;1 1 1 2; tzg=qr4(A,5,100)運行后屏幕顯示結(jié)果請注意：下面的i表

21、示求第i個特征值，k是迭代次數(shù)，sk是原點位移量， Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解，At=Rk*Qk+sk*eye(n)迭代矩陣：i = 3tzgk = -4.549k = 5sk = -4.549B = 7.142tzgk = 7.142請注意：n階實對稱矩陣A的全部真特征值lamoda和至少含t個有效數(shù)字的近似特征值tzg如下：lamoda = -4.549 1.84 2.00 7.181tzg = -4.549 1.30 2.00 7.142 用求根位移QR方法求實對稱矩陣全部特征值，精度為.并將計算結(jié)果與全部真特征值比較.其中解 首先把用求根位移QR方法求實

22、對稱矩陣全部特征值的MATLAB主程序保存為M文件，命名為qr8.m.然后在工作窗口輸入MATLAB程序 A=5 2 2 1;2 -3 1 1;2 1 3 1;1 1 1 2; tzg=qr8(A,0.0001,100)運行后屏幕顯示結(jié)果如下：請注意：下面的i表示求第i個特征值，如果迭代矩陣Ak的階數(shù)2,且m 階矩陣Ak的m行第m-1列的元近似等于零.則原n階矩陣A的第j個特征值j=skj，j=1,2,.,n-2;下面的矩陣Ak降一階.i = 3tzgk = 2.0004Ak = 4.8235 -2.0282 -2.0282 -5.2085如果迭代矩陣Ak的階數(shù)=2,則原n階矩陣A的最后兩個特

23、征值j=k+xj，k=n-2,j=1,2.x1 = 5.2180x2 = -5.6030tzg1 = 7.2183tzg2 = -3.6027請注意：n階實對稱矩陣A的全部真特征值lamoda和精度為jd的近似特征值tzg如下：lamoda = -3.6027 1.3843 2.0000 7.2183tzg = -3.6027 1.3843 2.0004 7.2183參 考 文 獻(xiàn)1 姜啟源，謝金星，葉俊編數(shù)學(xué)模型（第三版）M北京：高等教育出版社，2005：1-202.2 王建衛(wèi)，曲中水 凌濱編著. MATLAB 7.X 程序設(shè)計M. 北京：中國水利水電出版社，2007：55-80.3 李慶揚

24、，王能超，易大義編著.數(shù)值分析（第四版）M. 武漢：華中科技大學(xué)出版社，2006：219-245.4 王萼芳編著.高等代數(shù)M. 北京：高等教育出版社,2009.161-210.5 張軍編著.數(shù)值計算M.北京：清華大學(xué)出版社，2008.7.1-200.6 莫勒編著.MATLAB數(shù)值計算M.北京：機械工業(yè)出版社，2006.6.1-150.附 錄程序1： 用冪法計算矩陣的主特征值和對應(yīng)的特征向量的MATLAB主程序function k,lambda,Vk,Wc=mifa(A,V0,jd,max1)lambda=0;k=1;Wc =1; ,jd=jd*0.1;state=1; V=V0;while(k

25、jd)state=1;endk=k+1;Wc=Wc;endif(Wc=jd)disp(請注意：迭代次數(shù)k,主特征值的近似值lambda,主特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：) elsedisp(請注意：迭代次數(shù)k已經(jīng)達(dá)到最大迭代次數(shù)max1,主特征值的迭代值lambda,主特征向量的迭代向量Vk,相鄰兩次迭代的誤差Wc如下：) end Vk=V;k=k-1;Wc;程序2.用原點位移反冪法計算矩陣的特征值和對應(yīng)的特征向量的MATLAB主程序1function k,lambdan,Vk,Wc=ydwyfmf(A,V0,jlamb,jd,max1)n,n=size(A); A1=A

26、-jlamb*eye(n); jd= jd*0.1;RA1=det(A1); if RA1=0disp(請注意：因為A-aE的n階行列式hl等于零，所以A-aE不能進(jìn)行LU分解.)returnendlambda=0;if RA1=0 for p=1:nh(p)=det(A1(1:p, 1:p);endhl=h(1:n);for i=1:nif h(1,i)=0disp(請注意：因為A-aE的r階主子式等于零，所以A-aE不能進(jìn)行LU分解.) returnendend if h(1,i)=0 disp(請注意：因為A-aE的各階主子式都不等于零，所以A-aE能進(jìn)行LU分解.)k=1;Wc =1;

27、state=1; Vk=V0;while(kjd)state=1;endk=k+1;%Vk=Vk2,mk=mk1,endif(Wc=jd)disp(A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：) elsedisp(A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k已經(jīng)達(dá)到最大迭代次數(shù)max1,按模最小特征值的迭代值lambda,特征向量的迭代向量Vk,相鄰兩次迭代的誤差Wc如下：) endhl,RA1endendV,D=eig(A,nobalance),Vk;k=k-1;Wc;l

28、ambdan=jlamb+1/mk1;程序3.用原點位移反冪法計算矩陣的特征值和對應(yīng)的特征向量的MATLAB主程序function k,lambdan,Vk,Wc=wfmifa1(A,V0,jlamb,jd,max1)n,n=size(A); jd= jd*0.1;A1=A-jlamb*eye(n);nA1=inv(A1); lambda1=0;k=1;Wc =1;state=1; U=V0;while(kjd)state=1;endk=k+1;endif(Wc=jd) disp(請注意迭代次數(shù)k,特征值的近似值lambda,對應(yīng)的特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：) els

29、edisp(請注意迭代次數(shù)k已經(jīng)達(dá)到最大迭代次數(shù)max1, 特征值的近似值lambda,對應(yīng)的特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：) endV,D =eig(A,nobalance), Vk=U;k=k-1;Wc;lambdan=jlamb+1/mk;程序4.用最末元位移QR方法求實對稱矩陣全部特征值的MATLAB主程序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(A

30、k(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(請注意：下面的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