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

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

2、一些理論物理中的問(wèn)題.在本論文中，我們主要介紹求矩陣的特征值與特征向量的一些原理和方法，原理涉及高得代數(shù)中矩陣的相關(guān)定理，方法主要介紹冥法及反冥法，Jacobi方法和QR算法，并利用MATLAB，VC等軟件編寫相關(guān)算法的程序來(lái)求解相關(guān)問(wèn)題，加以驗(yà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 符號(hào)說(shuō)明.24 冥法及反冥法.2 4.1冥法.3 4.2反冥法.85 QR算法.14參考文獻(xiàn).18 附錄.19 1 引言在本論文中，我們主要討論矩陣的特征值和特征向量的計(jì)算，我們知道，有很多現(xiàn)實(shí)中的問(wèn)題都可以用到矩陣特征值與特征向量計(jì)算的知識(shí)，比如，在物理、力學(xué)和工程技術(shù)方面有很多的應(yīng)用，并且發(fā)揮著極其重要的作用.因?yàn)檫@些問(wèn)題都可歸結(jié)為求矩陣特征值的問(wèn)題，具體到一些具體問(wèn)題，如振動(dòng)問(wèn)題，物理中某些臨界值的確定問(wèn)題以及一些理論物理中的問(wèn)題

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

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

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

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

12、) 設(shè)有個(gè)線性無(wú)關(guān)特征向量;(2) 設(shè)特征值滿足 且 (3)由改進(jìn)冪法得到(4.1.7)式),則有 (a) (b)且收斂速度由比值確定.實(shí)現(xiàn)冪法，每迭代一次主要是計(jì)算一次矩陣乘向量，可編一個(gè)子程序.例1.用MATLAB編寫冥法程序求矩陣主特征值及近似主特征量用冪法計(jì)算下列矩陣的主特征值和對(duì)應(yīng)的特征向量的近似向量，精度.并把輸出的結(jié)果真實(shí)結(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,運(yùn)行后屏幕顯示結(jié)果請(qǐng)注意：迭代次數(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) 反冪法可用來(lái)計(jì)算矩陣按模最小的特征值及對(duì)應(yīng)的特征向量.設(shè)為非廳異矩陣，特征值滿足對(duì)應(yīng)特征向量為線性無(wú)關(guān)，則特征求值為特征向量為

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

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

16、2） 1）求 求 2） 3）對(duì)于計(jì)算對(duì)稱三對(duì)角陣，或計(jì)算Hessenberg陣對(duì)應(yīng)于一個(gè)給定的近似特征值的特征向量，反冪 法是一個(gè)有效方法.例2 .用反冪法計(jì)算對(duì)應(yīng)于近似特征值（精確特征值為）的特征向量 解 取，用部分選主元分解法實(shí)現(xiàn)，其中（1）求解 （2）求解 （3）求解特征向量（真解）是 由此，相當(dāng)好的近似.例3.用原點(diǎn)位移反冪法的迭代公式，根據(jù)給定的下列矩陣的特征值的初始值，計(jì)算與對(duì)應(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)運(yùn)行后屏幕顯示結(jié)果請(qǐng)注意：因?yàn)锳-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. 用原點(diǎn)位移反冪法的迭代公式（5.28），計(jì)算的分別對(duì)應(yīng)于特征值的特征向量的近似向量，相鄰迭代誤差為0.001.將計(jì)算結(jié)果與精確特征向量比較.解 計(jì)算特征值對(duì)應(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,運(yùn)行后屏幕

19、顯示結(jié)果請(qǐng)注意：因?yàn)锳-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 的對(duì)應(yīng)分量的比值相等.特征值的近似值lambda 1.002與初始值1.001的絕對(duì)誤差為0.001，而與的絕對(duì)誤差為0.002，其中,. 5 QR方法用最末元位移QR方法求下列實(shí)對(duì)稱矩陣的全部近似特征值，并將計(jì)算結(jié)果與全部真特征值比較.其中精度為解 首先保存用最末元位移QR方法求實(shí)對(duì)稱矩陣全部特征值的MATLAB主程序?yàn)镸文件，取名為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)運(yùn)行后屏幕顯示結(jié)果請(qǐng)注意：下面的i表

21、示求第i個(gè)特征值，k是迭代次數(shù)，sk是原點(diǎn)位移量， 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請(qǐng)注意：n階實(shí)對(duì)稱矩陣A的全部真特征值lamoda和至少含t個(gè)有效數(shù)字的近似特征值tzg如下：lamoda = -4.549 1.84 2.00 7.181tzg = -4.549 1.30 2.00 7.142 用求根位移QR方法求實(shí)對(duì)稱矩陣全部特征值，精度為.并將計(jì)算結(jié)果與全部真特征值比較.其中解 首先把用求根位移QR方法求實(shí)

22、對(duì)稱矩陣全部特征值的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)運(yùn)行后屏幕顯示結(jié)果如下：請(qǐng)注意：下面的i表示求第i個(gè)特征值，如果迭代矩陣Ak的階數(shù)2,且m 階矩陣Ak的m行第m-1列的元近似等于零.則原n階矩陣A的第j個(gè)特征值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的最后兩個(gè)特

23、征值j=k+xj，k=n-2,j=1,2.x1 = 5.2180x2 = -5.6030tzg1 = 7.2183tzg2 = -3.6027請(qǐng)注意：n階實(shí)對(duì)稱矩陣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è)計(jì)M. 北京：中國(guó)水利水電出版社，2007：55-80.3 李慶揚(yáng)

24、，王能超，易大義編著.數(shù)值分析（第四版）M. 武漢：華中科技大學(xué)出版社，2006：219-245.4 王萼芳編著.高等代數(shù)M. 北京：高等教育出版社,2009.161-210.5 張軍編著.數(shù)值計(jì)算M.北京：清華大學(xué)出版社，2008.7.1-200.6 莫勒編著.MATLAB數(shù)值計(jì)算M.北京：機(jī)械工業(yè)出版社，2006.6.1-150.附 錄程序1： 用冪法計(jì)算矩陣的主特征值和對(duì)應(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(請(qǐng)注意：迭代次數(shù)k,主特征值的近似值lambda,主特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：) elsedisp(請(qǐng)注意：迭代次數(shù)k已經(jīng)達(dá)到最大迭代次數(shù)max1,主特征值的迭代值lambda,主特征向量的迭代向量Vk,相鄰兩次迭代的誤差Wc如下：) end Vk=V;k=k-1;Wc;程序2.用原點(diǎn)位移反冪法計(jì)算矩陣的特征值和對(duì)應(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(請(qǐng)注意：因?yàn)锳-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(請(qǐng)注意：因?yàn)锳-aE的r階主子式等于零，所以A-aE不能進(jìn)行LU分解.) returnendend if h(1,i)=0 disp(請(qǐng)注意：因?yàn)锳-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.用原點(diǎn)位移反冪法計(jì)算矩陣的特征值和對(duì)應(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(請(qǐng)注意迭代次數(shù)k,特征值的近似值lambda,對(duì)應(yīng)的特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：) els

29、edisp(請(qǐng)注意迭代次數(shù)k已經(jīng)達(dá)到最大迭代次數(shù)max1, 特征值的近似值lambda,對(duì)應(yīng)的特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：) endV,D =eig(A,nobalance), Vk=U;k=k-1;Wc;lambdan=jlamb+1/mk;程序4.用最末元位移QR方法求實(shí)對(duì)稱矩陣全部特征值的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(請(qǐng)注意：下面的i表示求第i個(gè)特征值，k是迭代次數(shù)，sk是原點(diǎn)位移量，)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(請(qǐng)注意：i表示求第i個(gè)特征值，tzgk是矩陣A的特征值的近似值，k是迭代次數(shù)，) disp( 下面的矩陣B是m