冪法-反冪法求解矩陣最大最小特征值及其對應(yīng)的特征向量

1、數(shù)值計算解矩陣的按模最大最小特征值及對應(yīng)的特征向量一.哥法1 .幕法簡介：當(dāng)矩陣A滿足一定條件時,在工程中可用幕法計算其主特征值(按模最大)及其特征向量.矩陣A需要滿足的條件為：(1) |l|2|n|0一為A的特征值(2)存在n個線性無關(guān)的特征向量,設(shè)為Xi,X2,.,Xn1.1 計算過程：n對任意向量X(0),有X(0)iUi,i不全為0,那么有i1x(k1)Ax的Ak1x(0)nnAk1碼q*1Uii1i1、k12k1n、k111U1()a2U2(-)anUn11k1u11U1可見,當(dāng)2|越小時,收斂越快；且當(dāng)k充分大時,有1(k1)XX(k)1U1x(k1)1U1(k1)1,對應(yīng)的特征向

2、量即是X.2算法實現(xiàn).輸入矩陣A,初始向量.k1,0;y(k)x,誤差限,最大迭代次數(shù)Nx(k).計算X.假設(shè)|Ay,1|,輸由max(abs(x(k)max(x);,y,否那么,轉(zhuǎn)(5)假設(shè)kN,置kk1,轉(zhuǎn)3,否那么輸由失敗信息,停機.3matlab程序代碼functiont,y=lpowerA,x0,eps,N)%t為所求特征值,y是對應(yīng)特征向量k=1;z=0;%z相當(dāng)于y=x0./max(abs(x0);%標(biāo)準(zhǔn)化初始向量x=A*y;%迭代格式b=max(x);%bifabs(z-b)<eps%t=max(x);return;endwhileabs(z-b)>eps&

3、&k<Nk=k+1;z=b;y=x./max(abs(x);x=A*y;b=max(x);endm,index=max(abs(x);%t=x(index);%end相當(dāng)于判斷第一次迭代后是否滿足要求這兩步保證取出來的按模最大特征值是原值,而非其絕對值.4舉例驗證選取一個矩陣A,代入程序,得到結(jié)果,并與eig(A)的得到結(jié)果比擬,再計算A*y-t*y,驗證y是否是對應(yīng)的特征向量.結(jié)果如下：結(jié)果正確,說明算法和代碼正確,然后利用此程序計算15階Hilb矩陣,與eig(A)的得到結(jié)果比擬,再計算A*y-t*y,驗證y是否是對應(yīng)的特征向量.設(shè)置初始向量為x0=ones(15,1),結(jié)果

4、顯示如下可見,結(jié)果正確.得到了15階Hilb矩陣的按模最大特征值和對應(yīng)的特征向量.二.反哥法1 .反幕法簡介及其理論在工程計算中,可以利用反幕法計算矩陣按模最小特征值及其對應(yīng)特征向量.其根本理論如下,與幕法根本相同：ii1由AxxxA(x),那么Axx,可知,A和A1的特征值互為倒數(shù),求A按模最小特征值即求A1的按模最大特征值,取倒數(shù)即為A的按模最小特征值所以算法根本相同,區(qū)別就是在計算x(k1)時,不是令x(k1)Ay(k),而是x(k1)A-1y(k)具體計算時,變換為Ax(k1)y(k);對A做LU分解,來計算x(k1)2 .算法實現(xiàn).輸入矩陣A,初始向量x,誤差限,最大迭代次數(shù)N,(2

5、) .置k1,00,yx,max(abs(x)(3) .作三角分解ALU(4) .解方程組LUxy(Lzy,Uxz),(5) .max(x),一一,1(6) .右|0|,輸出一,y,停機,否那么轉(zhuǎn),(7) .假設(shè)kN,置kk1,0,y/、,轉(zhuǎn)(4);max(abs(x)否那么輸出失敗信息,停機.3matlab程序代碼functions,y=invpower(A,x0,eps,n)%s為按模最小特征值,y是對應(yīng)特征Ik=1;r=0;%r相當(dāng)于0y=x0./max(abs(x0);%標(biāo)準(zhǔn)化初始向量L,U=lu(A);z=Ly;x=Uz;u=max(x);s=1/u;%按模最小為A1按模最大的倒數(shù).

6、ifabs(u-r)<eps%判斷第一次迭代后是否滿足終止條件returnendwhileabs(u-r)>eps&&k<n%k=k+1;r=u;y=x./max(abs(x);z=Ly;x=Uz;u=max(x);endm,index=max(abs(x);%s=1/x(index);%終止條件.這兩步保證取出來的按模最大特征值是原值,而非其絕對值.end4舉例驗證同事法一樣,選取一個矩陣A,代入程序,得到結(jié)果,并與eig(A)的得到結(jié)果比較,再計算A*y-t*y,驗證y是否是對應(yīng)的特征向量.可見結(jié)果正確,然后利用此程序計算15階Hilb矩陣,eig(A)的

7、得到結(jié)果比較,再計算A*y-s*y,驗證y是否是對應(yīng)的特征向量.設(shè)置初始向量為x0=ones(15,1),結(jié)果顯示如下可見,結(jié)果真確.得到了15階Hilb矩陣的按模最大特征值和對應(yīng)的特征向量.三.計算條件數(shù)矩陣A的條件數(shù)等于A的范數(shù)與A的逆的范數(shù)的乘積,即cond(A)=IIAll11AA(-1)II,對應(yīng)矩陣的3種范數(shù),可以定義3種條件數(shù).函數(shù)cond(A,1)、cond(A)或cond(Ainf)是判斷矩陣病態(tài)與否的一種度量,條件數(shù)越大說明矩陣的病態(tài)程度越大.這里我們選擇矩陣的2范數(shù),即cond(A),1,2分別為2矩陣ATA的最大和最小特征值而如果A為對稱矩陣,如Hilb矩陣,ATA的最

8、大最小特征值,分別為A的最大最小特征值的平方.所以cond(A)為A的最大最小特征值得比值.對于本例中的15階Hilb矩陣來說,利用上面計算結(jié)果得其條件數(shù)(選擇第二種條件數(shù))為：3.0934e+017；這與直接利用cond(A)得到的結(jié)果：2.5083e+017在同一數(shù)量級,再次說明了上述算得得最大最小特征值的正確性,同時又說明Hilb矩陣是病態(tài)矩陣.四.Aitken商加速法1 .簡介與原理假設(shè)ak收斂與a,且imak1aakac0;即ak線性收斂,當(dāng)k充分大時,有曳工招曳工招akaak1ayn2xn2(xn2xn1)xn22xn1xnaak(ak1aQ2ak22ak1用ak逼近a,這種方法稱

9、為一:仇akAitken加速法.同事法和反幕法計算最大和最小特征值類似,如果計算最大特征值,那么迭代格式為x(k1)Ay(k)；計算最小特征值時,迭代格式為x(k1)A1y(k),即y(k)Ax(k1).2 .算法實現(xiàn)計算按模最大特征值算法如下：輸入Aa),初始向量x,誤差限,最大迭代次數(shù)N,3 2).置k1,o0,i0,o1.0,y(3) .計算x(4) .計算max(x)Ay,置max(x)(1o)2022210.假設(shè)(6).假設(shè)k0,輸出,y停機,否那么轉(zhuǎn)(6),N,置10,21,0,k1k,七通否那么,輸出失敗信息,停機.類似幕法和反幕法可以寫出按模最小特征值算法,此處不再贅述3 .m

10、atlab程序代碼functionr,y=aitken(A,x0,eps,n)%r按模最大特征值,y為對應(yīng)特征向量k=1;a0=0;%a相當(dāng)于0a1=1;%a1相當(dāng)于1r0=1;%相當(dāng)于2中的0y=x0./max(abs(x0);%標(biāo)準(zhǔn)化初始向量x=A*y;a2=max(abs(x);%a2相當(dāng)于2r=a0-(a1-a0)A2/(a2-2*a1+a0);%相當(dāng)于if(a2-2*a1+a0)=0%假設(shè)上式中分母為0,那么迭代失敗,返回disp"初始向量迭代失敗"return;endifabs(r-r0)<eps%判斷第一次迭代后是否滿足要求,如滿足,那么返回結(jié)果retu

11、rnendwhileabs(r-r0)>eps&&k<n%終止條件k=k+1;a0=a1;a1=a2;r0=r;y=x./max(abs(x);x=A*y;%迭代格式a2=max(abs(x);if(a2-2*a1+a0)=0%假設(shè)分母為0,那么迭代失敗,返回return;endr=a0-(a1-a0)A2/(a2-2*a1+a0);m,index=max(abs(eig(A);%以下代碼保證取出來的按模最大特征值aa=eig(A);%是原值,而非其絕對值.ifaa(index)>0|aa(index)=0r=r;elser=-r;endendend類似可得按

12、模最小特征值和特征向量的代碼如下：與上面類似,所不同的只是迭代格式不同.functionr,y=invaitken(A,x0,eps,n)k=1;a0=0;a1=1;r0=1;y=x0./max(abs(x0);L,U=lu(A);%迭代格式的不同z=Ly;x=Uz;a2=max(abs(x);r=a0-(a1-a0)A2/(a2-2*a1+a0);if(a2-2*a1+a0)=0disp"初始向量迭代失敗"return;endifabs(r-r0)<eps%判斷第一次迭代后是否滿足要求,如滿足,那么返回結(jié)果returnendwhileabs(r-r0)>eps

13、&&k<nk=k+1;a=b;b=c;r0=r;y=x./max(abs(x);z=Ly;x=Uz;a2=max(abs(x);if(a2-2*a1+a0)=0return;endr=a0-(a1-a0)A2/(a2-2*a1+a0);endm,index=min(abs(eig(A);%以下代碼保證取出來的按模最大特征值aa=eig(A);%是原值,而非其絕對值.ifaa(index)>0|aa(index)=0r=1/r;elser=-1/r;endend4 .計算Hilb矩陣特征值此處不再舉例,而是直接應(yīng)用于15階Hilb矩陣,初始向量選為ones(15,1)

14、,結(jié)果如下,并將結(jié)果與幕法和反幕法得到結(jié)果比擬這與幕法得到的特征值和特征向量一致,說明算法和代碼正確；同理,最小特征值結(jié)果如下：這與反幕法得到的結(jié)果一致,說明結(jié)果正確.五,對稱矩陣的Rayleigh商加速法1,簡介與原理A為對稱矩陣,設(shè)x0,那么稱R(x)#x為關(guān)于A白Rayleigh商xx原理如下：(k)y(k1)XR(y(k)其中R(y(k)1,y(k)(k)Xmax(x(k)2.算法實現(xiàn).輸入矩陣A,初始向量x,誤差限,最大迭代次數(shù)N,.置k1,r0.xA*y,r假設(shè)|rr.|0,yTyxTyyxmax(abs(x),輸出r,y,停機,否那么轉(zhuǎn)5,設(shè)Xi0為i的特征向量,即AXiiXi用

15、Xj左乘上式有：(k)Xmax(x(k)這稱為Rayleigh商加速法.Ay(k)(y(k)Tx(k1)(y(k)T(y(k)假設(shè)kN,置kk1,r0r,y否那么輸出失敗信息,停機.3 .Matlab程序代碼是特征值,y是特征向量functionr,y=rayleigh(A,x0,eps,n)%rk=1;r0=0;y=x0./max(abs(x0);x=A*y;%迭代格式計算新的xr=dot(y,x)/dot(y,y);%Reyleigh商ifabs(r-r0)<epsreturnendwhileabs(r-r0)>eps&&k<nk=k+1;r0=r;y=x./max(abs(x);x=A*y;r=dot(y,x)/dot(y,y);endend類似得計算按模最小特征值的Rayleigh商加速法,如下：functionr,y=invrayleigh(A,x0,eps,n)k=1;r0=0;y=x0./max(abs(x0);L,U=lu(A);%迭代格式不同z=Ly;x=Uz;r=dot(y,x)/dot(y,y);ifabs(r-r0)<