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

1、數(shù)值計(jì)算解矩陣的按模最大最小特征值及對應(yīng)的特征向量一.哥法.幕法簡介：當(dāng)矩陣A滿足一定條件時，在工程中可用幕法計(jì)算其主特征值(按模最大)及其特征向量。矩陣A需要滿足的條件為：|儲卜|九2戶以九n戶0,%為A的特征值存在n個線性無關(guān)的特征向量，設(shè)為XX2,，4計(jì)算過程：n對任意向量X(0),有x(0)=£%5,%不全為0,則有i1X(k1)=AX的=.=Ak1X(0)nnk1k-1=£A國木uii=1i=1小書1L/裊2k+1上上/'-n'k+I=1忤M+(一)a?U2+()anUnJ九1九11k11lUi2可見，當(dāng)|一|越小時，收斂越快；且當(dāng)k充分大時，有,

2、i(k+1)nk1X()=匕«1U1(k1)X(k)qkX(-'1-1U1X(k).(k1)人1,對應(yīng)的特征向量即是X2算法實(shí)現(xiàn).輸入矩陣A,初始向量X,誤差限"最大迭代次數(shù)N/).k,1,-0;y(k)=芯maX(abs(X(k).計(jì)算X，Ay,；max(x);(4).若|九P|<巴輸由Ky,否則，轉(zhuǎn)(5).若k<N,置kk+1J轉(zhuǎn)3,否則輸由失敗信息，停機(jī).3matlab程序代碼functiont,y=lpowerA,x0,eps,N)%t為所求特征值，y是對應(yīng)特征向量k=1;z=0;%z相當(dāng)于黑y=x0./max(abs(x0);%規(guī)范化初始1可量

3、x=A*y;%迭代格式b=max(x);%b相當(dāng)于Pifabs(z-b)<eps%判斷第一次迭代后是否滿足要求t=max(x);return;endwhileabs(z-b)>eps&&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);%是原值，而非其絕對值。end4舉例驗(yàn)證選取一個矩陣A,代入程序，得到結(jié)果，并與eig(A)的得到結(jié)果比較，再計(jì)算A*y-t*y,驗(yàn)證y是否是對應(yīng)的特征向量。結(jié)果如下：»A=l10.

4、5:110.25;0.3,0.25,2A=1.00001.00000.5000L0000L00000.25000.50000.25002.0000»xO=l;l;lJ?0ps=le-5;N=20;>>t,(A,xOjeps,、)2.53650.74820.6497L0000>>eig(A)ans=-0.0166L48012.5365»A*y-t*yans-1.0e-004+-0.1603-0,16840結(jié)果正確，表明算法和代碼正確，然后利用此程序計(jì)算15階Hilb矩陣，與eig(A)的得到結(jié)果比較，再計(jì)算A*y-t*y,驗(yàn)證y是否是對應(yīng)的特征向量。設(shè)

5、置初始向量為x0=ones(15,1),結(jié)果顯示如下»A=hilb(15);»x0=ones(15,1):eps-le6：»N=30;>>Lt,y-1power(A,x0,eps,t=L8459V-»eig(A)»A*y-t*y1.0000ans=ans=0.62280.47060.3838-0.00DO0.00001.0-007拿00.32640.0000-0.55160.0000-0.64360.28510.0000-0.64650.25370.0000-0.62450.22900.0000-0.59530.20890.oooo

6、-0.565。0.19220.oooo-0.53590.17800,0000-0.50860.16590,0004-0.48340.0056-0.46030.15540.0572-0.43900.14620.4266-0.41950:13811.8439-0.401bjn.GG-可見，結(jié)果正確。得到了15階Hilb矩陣的按模最大特征值和對應(yīng)的特征向二.反哥法.反幕法簡介及其理論在工程計(jì)算中，可以利用反幕法計(jì)算矩陣按模最小特征值及其對應(yīng)特征向量。其基本理論如下，與幕法基本相同：一一i一一一i1.由Ax=Kx=x=A(九x),則Ax=x,可知，A和A1的特征值互為倒數(shù)，求A按模最小特征值即求A-1

7、的按模最大特征值，取倒數(shù)即為A的按模最小特征值所以算法基本相同，區(qū)別就是在計(jì)算x(k+)時,不是令x(k+)=Ay(k),而是x(k+)=A-1y(k)具體計(jì)算時,變換為Ax(k*)=y(k);對A做LU分解,來計(jì)算x(k+).算法實(shí)現(xiàn).輸入矩陣A,初始向量x,誤差限與最大迭代次數(shù)N,.置k1,0,y,max(abs(x).作三角分解A二LU.解方程組LUx=y(Lz=y,Ux=z),.max(x),一1.右|九-腦|<8,輸出1,y,停機(jī)，否則轉(zhuǎn)(7),xj.右k<N,置k-k+diy-mx詢而，轉(zhuǎn)(4);否則輸出失敗信息，停機(jī).3matlab程序代碼functions,y=in

8、vpower(A,x0,eps,n)%s為按模最小特征值，y是對應(yīng)特征向量k=1;r=0;%r相當(dāng)于-oy=x0./max(abs(x0);%規(guī)范化初始向量L,U=lu(A);z=Ly;x=Uz;u=max(x);s=1/u;%按模最小為A-1按模最大的倒數(shù).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

9、=1/x(index);%是原值，而非其絕對值。4舉例驗(yàn)證同事法一樣，選取一個矩陣endA,代入程序，得到結(jié)果，并與eig(A)的得到結(jié)果比較，再計(jì)算A*y-t*y,驗(yàn)證y是否是對應(yīng)的特征向量»A=L1,0.5；1,1,0.210.5,0.25,2;»eig(A)eps=le-6;口=30;_s,y=invpower(A,k。，eps,n)ans二0166L48012.5365s=3*-0.0166ans=y-1.0e-016*L0000-0.1388-0.0694-0.9517-0.1300-0.1908可見結(jié)果正確，然后利用此程序計(jì)算15階Hilb矩陣，eig(A)的得

10、到結(jié)果比較，再計(jì)算A*y-s*y,驗(yàn)證y是否是對應(yīng)的特征向量。設(shè)置初始向量為x0=ones(15,1),結(jié)果顯示如下»A=liilb(15):xO=ones(1571);eps=le-6;n=30;>>s,y=invpower(A?epsn)-5.9673e-0180.0000-0.00000.0000-0.00060.0050-0.02450.0699-0.0968-0.04210.4497-0.9084L0000-0.65270.23790.0375»eig(A)ans=-0,00000.00000.00000.00000,00000.00000.00000

11、.00000.00000.00000.00040.00560.05720.42661.8459>>A*y-s*yans=1.05017*0.3903-0.56380.0000-0.5208-0.47410.35400.4537-0.27460.9724-0.5556-0.62880.66180.06390.1419可見，結(jié)果真確。得到了15階Hilb矩陣的按模最大特征值和對應(yīng)的特征向量。三.計(jì)算條件數(shù)矩陣A的條件數(shù)等于A的范數(shù)與A的逆的范數(shù)的乘積，即cond(A)=IIAll11AA(-1)II,對應(yīng)矩陣的3種范數(shù)，可以定義3種條件數(shù)。函數(shù)cond(A,1)、cond(A)或con

12、d(Ainf)是判斷矩陣病態(tài)與否的一種度量，條件數(shù)越大表明矩陣的病態(tài)程度越大.這里我們選擇矩陣的2范數(shù)，即cond(A)=',%,%分別為矩陣aTA的最大和最小特征值而如果A為對稱矩陣，如Hilb矩陣，ATA的最大最小特征值，分別為A的最大最小特征值的平方。所以cond(A)為A的最大最小特征值得比值。對于本例中的15階Hilb矩陣來說，利用上面計(jì)算結(jié)果得其條件數(shù)(選擇第二種條件數(shù))為：3.0934e+017;這與直接利用cond(A)得到的結(jié)果：2.5083e+017在同一數(shù)量級，再次表明了上述算得得最大最小特征值的正確性，同時又表明陣是病態(tài)矩陣。四.Aitken商加速法1,簡介與原

13、理若瓜收斂與a,且lim電工二a=c#。;即ak蹴性收斂，j'ak-a當(dāng)k充分大時，有虹二世二ak-aak1-a、，、，(xn2-xn1)yn2;xn2xn2-2xn1xn一(ak1-ak)小-a：ak尸akak2-2ak1'ak用ak逼近a,這種方法稱為Aitken加速法.同事法和反幕法計(jì)算最大和最小特征值類似,如果計(jì)算最大特征值,為x(k*)=Ay(k)；計(jì)算最小特征值時，迭代格式為Hilb矩則迭代格式x(s=A,y(kfWk)=Ax(k41)2,算法實(shí)現(xiàn)計(jì)算按模最大特征值算法如下：(1),輸入A=(a。),初始向量x,誤差限以最大迭代次數(shù)N,x.置。=。，=。，<&

14、quot;.。,廣加,計(jì)算x=Ay,置max(x)=a2,2(-?)(4),計(jì)算。-一一匚:2-2：。<3輸出y停機(jī)，否則轉(zhuǎn)(6),(6),若卜<N,置a1na。，a2n%,九二九。,k+i，ma=y轉(zhuǎn),否則，輸出失敗信息，停機(jī).類似幕法和反幕法可以寫出按模最小特征值算法，此處不再贅述3.matlab程序代碼functionr,y=aitken(A,x0,eps,n)%r按模最大特征值，丫為對應(yīng)特征向量k=1;a0=0;%a相當(dāng)于«0a1=1;%a1相當(dāng)于«ir0=1;%相當(dāng)于2中的八0y=x0./max(abs(x0);%規(guī)范化初始向量x=A*y;a2=max

15、(abs(x);%a2相當(dāng)于二2r=a0-(a1-a0)A2/(a2-2*a1+a0);%相當(dāng)于九if(a2-2*a1+a0)=0%若上式中分母為0,則迭代失敗，返回disp"初始向量迭代失敗"return;endifabs(r-r0)<eps%判斷第一次迭代后是否滿足要求，如滿足，則返回結(jié)果returnendwhileabs(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%若分母為0,則迭

16、代失敗，返回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類似可得按模最小特征值和特征向量的代碼如下：與上面類似，所不同的只是迭代格式不同.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=m

17、ax(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&&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

18、=min(abs(eig(A);%以下代碼保證取出來的按模最大特征值aa=eig(A);%是原值，而非其絕對值。ifaa(index)>0|aa(index)=0r=1/r;elser=-1/r;endend4.計(jì)算Hilb矩陣特征值此處不再舉例，而是直接應(yīng)用于15階Hilb矩陣，初始向量選為ones(15,1)結(jié)果如下，并將結(jié)果與幕法和反幕法得到結(jié)果比較L00000.6229»A=hilb(15);xO=ones(13,l):»eps=le-6;»n=30;y_=aitken(AJxOeps,n)1.84590.47060,38380.32640.2851

19、0.25370.22900.20890.19220.17800.16590.15540.14620.1381這與幕法得到的特征值和特征向量一致，表明算法和代碼正確；同理，最小特征值結(jié)果如下：>>r,y=invaitken(A,x0fepSjn)5.9673bo18S0000-0.00000,0000-0.00060.0050-0.02450.0699-0.09687.04210.4497-Q9084L0000-0.65270.23790.0375這與反幕法得到的結(jié)果一致，表明結(jié)果正確五，對稱矩陣的Rayleigh商加速法1,簡介與原理A為對稱矩陣,設(shè)x#0,則稱R(x)x?x為關(guān)于

20、A白Rayleigh商xx(k)yx(k+)R(y(k)I原理如下:設(shè)Xj#0為的特征向量，即AXi=iXi用Xi左乘上式有：(k)xmax(x(k)這稱為Rayleigh商加速法。=Ay(k)_(y(k)Tx(k1)/_(k"T/(kU(y)(y)其中R(y(k)，,y(k)，(k)xmax(x(k)2.算法實(shí)現(xiàn).輸入矩陣A,初始向量x,誤差限叫最大迭代次數(shù)N,x2).置卜*1,r0，0,y，,max(abs(x)T.x,A*y,r若|r-r。卜：；,yxTyy輸出r,y,停機(jī),否則轉(zhuǎn)(5),x.右k<N,置k-k+1,r0-r,y-,轉(zhuǎn)(3);max(abs(x)否則輸出失敗信息，停機(jī).Matlab程序代碼functionr,y=rayleigh(A,x0,eps,n)%r是特征值，y是特征向量k=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類似得計(jì)算按模最小特征值的Rayleigh商