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

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

2、量即是X.2算法實(shí)現(xiàn).輸入矩陣A,初始向量.k1,0;y(k)x,誤差限,最大迭代次數(shù)Nx(k).計(jì)算X.假設(shè)|Ay,1|,輸由max(abs(x(k)max(x);,y,否那么,轉(zhuǎn)(5)假設(shè)kN,置kk1,轉(zhuǎn)3,否那么輸由失敗信息,停機(jī).3matlab程序代碼functiont,y=lpowerA,x0,eps,N)%t為所求特征值,y是對(duì)應(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)于判斷第一次迭代后是否滿足要求這兩步保證取出來(lái)的按模最大特征值是原值,而非其絕對(duì)值.4舉例驗(yàn)證選取一個(gè)矩陣A,代入程序,得到結(jié)果,并與eig(A)的得到結(jié)果比擬,再計(jì)算A*y-t*y,驗(yàn)證y是否是對(duì)應(yīng)的特征向量.結(jié)果如下：結(jié)果正確,說(shuō)明算法和代碼正確,然后利用此程序計(jì)算15階Hilb矩陣,與eig(A)的得到結(jié)果比擬,再計(jì)算A*y-t*y,驗(yàn)證y是否是對(duì)應(yīng)的特征向量.設(shè)置初始向量為x0=ones(15,1),結(jié)果

4、顯示如下可見(jiàn),結(jié)果正確.得到了15階Hilb矩陣的按模最大特征值和對(duì)應(yīng)的特征向量.二.反哥法1 .反幕法簡(jiǎn)介及其理論在工程計(jì)算中,可以利用反幕法計(jì)算矩陣按模最小特征值及其對(duì)應(yīng)特征向量.其根本理論如下,與幕法根本相同：ii1由AxxxA(x),那么Axx,可知,A和A1的特征值互為倒數(shù),求A按模最小特征值即求A1的按模最大特征值,取倒數(shù)即為A的按模最小特征值所以算法根本相同,區(qū)別就是在計(jì)算x(k1)時(shí),不是令x(k1)Ay(k),而是x(k1)A-1y(k)具體計(jì)算時(shí),變換為Ax(k1)y(k);對(duì)A做LU分解,來(lái)計(jì)算x(k1)2 .算法實(shí)現(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,停機(jī),否那么轉(zhuǎn),(7) .假設(shè)kN,置kk1,0,y/、,轉(zhuǎn)(4);max(abs(x)否那么輸出失敗信息,停機(jī).3matlab程序代碼functions,y=invpower(A,x0,eps,n)%s為按模最小特征值,y是對(duì)應(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);%終止條件.這兩步保證取出來(lái)的按模最大特征值是原值,而非其絕對(duì)值.end4舉例驗(yàn)證同事法一樣,選取一個(gè)矩陣A,代入程序,得到結(jié)果,并與eig(A)的得到結(jié)果比較,再計(jì)算A*y-t*y,驗(yàn)證y是否是對(duì)應(yīng)的特征向量.可見(jiàn)結(jié)果正確,然后利用此程序計(jì)算15階Hilb矩陣,eig(A)的

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

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

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

10、atlab程序代碼functionr,y=aitken(A,x0,eps,n)%r按模最大特征值,y為對(duì)應(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);%以下代碼保證取出來(lái)的按模最大特征值aa=eig(A);%是原值,而非其絕對(duì)值.ifaa(index)>0|aa(index)=0r=r;elser=-r;endendend類(lèi)似可得按

12、模最小特征值和特征向量的代碼如下：與上面類(lèi)似,所不同的只是迭代格式不同.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);%以下代碼保證取出來(lái)的按模最大特征值aa=eig(A);%是原值,而非其絕對(duì)值.ifaa(index)>0|aa(index)=0r=1/r;elser=-1/r;endend4 .計(jì)算Hilb矩陣特征值此處不再舉例,而是直接應(yīng)用于15階Hilb矩陣,初始向量選為ones(15,1)

14、,結(jié)果如下,并將結(jié)果與幕法和反幕法得到結(jié)果比擬這與幕法得到的特征值和特征向量一致,說(shuō)明算法和代碼正確；同理,最小特征值結(jié)果如下：這與反幕法得到的結(jié)果一致,說(shuō)明結(jié)果正確.五,對(duì)稱矩陣的Rayleigh商加速法1,簡(jiǎn)介與原理A為對(duì)稱矩陣,設(shè)x0,那么稱R(x)#x為關(guān)于A白R(shí)ayleigh商xx原理如下：(k)y(k1)XR(y(k)其中R(y(k)1,y(k)(k)Xmax(x(k)2.算法實(shí)現(xiàn).輸入矩陣A,初始向量x,誤差限,最大迭代次數(shù)N,.置k1,r0.xA*y,r假設(shè)|rr.|0,yTyxTyyxmax(abs(x),輸出r,y,停機(jī),否那么轉(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否那么輸出失敗信息,停機(jī).3 .Matlab程序代碼是特征值,y是特征向量functionr,y=rayleigh(A,x0,eps,n)%rk=1;r0=0;y=x0./max(abs(x0);x=A*y;%迭代格式計(jì)算新的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類(lèi)似得計(jì)算按模最小特征值的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)<