冪法,反冪法求解矩陣最大最小特征值及其對應地特征向量

1、實用標準文案數(shù)值計算解矩陣的按模最大最小特征值及對應的特征向量一. 冪法1. 冪法簡介：當矩陣 A 滿足一定條件時，在工程中可用冪法計算其主特征值( 按模最大 )及其特征向量。矩陣A 需要滿足的條件為：(1) | 1 | | 2 | .| n | 0 ,i 為 A的特征值(2) 存在 n 個線性無關的特征向量，設為1.1 計算過程 ：對任意向量 x(0 )，有 x( 0)ni ui ,i 1x(k 1)Ax(k).Ak 1 x( 0 )nk 1nAk 1u uii ii ii 1i1x1, x2 ,., xni 不全為 0，則有k1u ( 2 )k 1a u( n )k 1a un1112 2

2、n11k1u111可見，當 |2 | 越小時，收斂越快；且當 k 充分大時，有1x(k1 )k1x(k 1 )11u1，對應的特征向量即是 x(k 1 )。x(k )k1u1x(k )112 算法實現(xiàn)輸入矩陣，初始向量 ，誤差限,最大迭代次數(shù)N(1).Ax(2).k1,0; y( k)x( k )max(abs( x( k ) )計算xAy,max(x);(3).若| ,輸出, y,否則轉(zhuǎn)(5)(4).,若kN ,置kk1,轉(zhuǎn)3 ,否則輸出失敗信息,停機.(5).3 matlab程序代碼精彩文檔實用標準文案functiont,y=lpowerA,x0,eps,N) % t為所求特征值， y是對

3、應特征向量k=1;z=0;% z相當于y=x0./max(abs(x0);%規(guī)范化初始向量x=A*y;%迭代格式b=max(x);% b相當于if abs(z-b)<eps % 判斷第一次迭代后是否滿足要求 t=max(x);return ;endwhile abs(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 舉例驗證選取一個矩陣 A，代入程序，得到結(jié)果，并與e

4、ig(A) 的得到結(jié)果比較，再計算 A*y-t*y ，驗證 y 是否是對應的特征向量。結(jié)果如下：精彩文檔實用標準文案結(jié)果正確，表明算法和代碼正確， 然后利用此程序計算15 階 Hilb 矩陣，與eig(A) 的得到結(jié)果比較，再計算A*y-t*y ，驗證 y 是否是對應的特征向量。設置初始向量為 x0=ones(15,1) ，結(jié)果顯示如下精彩文檔實用標準文案可見，結(jié)果正確。得到了15階 Hilb 矩陣的按模最大特征值和對應的特征向量。二. 反冪法1. 反冪法簡介及其理論在工程計算中，可以利用反冪法計算矩陣按模最小特征值及其對應特征向量。其基本理論如下，與冪法基本相同：由Axx x A1，則A11

5、-1的特征值互為倒數(shù)，( x)xx ，可知，A 和 A-1求 A 按模最小特征值即求 A 的按模最大特征值 , 取倒數(shù)即為 A 的按模最小特征值所以算法基本相同，區(qū)別就是在計算x( k 1)時，不是令 x( k 1)Ay(k ) , 而是 x( k 1)A-1 y(k )具體計算時，變換為Ax( k 1)y( k ) ; 對 A做LU 分解，來計算 x( k 1)2. 算法實現(xiàn)精彩文檔實用標準文案(1).輸入矩陣 A,初始向量 x,誤差限, 最大迭代次數(shù) N ,(2).置 k1, 00, yx,max( abs( x)(3).作三角分解 ALU(4).解方程組 LUxy( Lzy,Uxz),(

6、5).max( x),(6).若 |0 |,輸出 1, y,停機 ,否則轉(zhuǎn) (7),(7).若 kN , 置 kk1,0, yx, 轉(zhuǎn) (4);max(abs( x)否則輸出失敗信息,停機 .3 matlab程序代碼functions,y=invpower(A,x0,eps,n) % s為按模最小特征值， y是對應特征向量k=1; r=0;% r相當于 0y=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 % 判斷第一次迭代后是否滿足終止條件 re

7、turnendwhile abs(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(jié)果正確， 然后利用此程序計算15 階 Hilb 矩陣，eig(A) 的得到結(jié)果比較，再計算A*

8、y-s*y ，驗證 y 是否是對應的特征向量。設置初始向量為x0=ones(15,1) ，結(jié)果顯示如下精彩文檔實用標準文案可見，結(jié)果真確。得到了15階 Hilb 矩陣的按模最大特征值和對應的特征向量。三.計算條件數(shù)矩陣 A 的條件數(shù)等于 A 的范數(shù)與 A 的逆的范數(shù)的乘積，即 cond(A)= A· A(-1) ，對應矩陣的 3種范數(shù)，可以定義 3種條件數(shù)。 函數(shù) cond(A,1) 、cond(A)或 cond(A inf) 是判斷矩陣病態(tài)與否的一種度量，條件數(shù)越大表明矩陣的病態(tài)程度越大 .這里我們選擇矩陣的2范數(shù)，即 cond ( A)1, 1, 2分別為2,矩陣 AT A的最大

9、和最小特征值而如果 A 為對稱矩陣，如Hilb 矩陣， AT A 的最大最小特征值，分別為A 的最大最小特征值的平方。所以 cond(A) 為 A 的最大最小特征值得比值。對于本例中的 15階 Hilb 矩陣來說，利用上面計算結(jié)果得其條件數(shù) ( 選擇第二種條件數(shù)）為：3.0934e+017；這與直接利用 cond(A) 得到的結(jié)果： 2.5083e+017 在同一精彩文檔實用標準文案數(shù)量級，再次表明了上述算得得最大最小特征值的正確性，同時又表明 Hilb 矩陣是病態(tài)矩陣。四 .Aitken 商加速法1. 簡介與原理若 a收斂與 a,且 lim ak 1ac 0;即 a線性收斂 ,kkakak當

10、 k充分大時 ,有ak1aak2aakaaka1yn 2( xn2xn 1 )2xn 22xn 1xnxn 2( a1ak)2aakk:a?k2ak 1ak 2ak用 a?k逼近 a, 這種方法稱為 Aitken 加速法 .同冪法和反冪法計算最大和最小特征值類似，如果計算最大特征值， 則迭代格 式為 x(k 1)Ay( k)；計算最小特征值時，迭代格式為x( k 1)A 1 y( k ) ,即 y( k )Ax (k 1) 。2. 算法實現(xiàn)計算按模最大特征值算法如下：(1).輸入 A(aij ), 初始向量 x,誤差限,最大迭代次數(shù) N ,( 2).置k1,0 0,0,0 1.0,x,1yma

11、x( x)(3).計算 xAy, 置 max( x)2 ,( 4).計算( 10 )202,210(5).若0, 輸出, y停機 , 否則轉(zhuǎn) (6),( 6).若kN,置 10 ,21,0 ,k 1k,xy轉(zhuǎn)(3),max( x)否則 , 輸出失敗信息 , 停機 .類似冪法和反冪法可以寫出按模最小特征值算法，此處不再贅述。精彩文檔實用標準文案3.matlab程序代碼function r,y=aitken(A,x0,eps,n) % r按模最大特征值 ,y 為對應特征向量k=1;a0=0;% a相當于0a1=1; % a1相當于1r0=1; %相當于 2中的0y=x0./max(abs(x0);

12、 %規(guī)范化初始向量x=A*y;a2=max(abs(x); % a2相當于2r=a0-(a1-a0)2/(a2-2*a1+a0); %相當于if (a2-2*a1+a0)=0%若上式中分母為 0，則迭代失敗，返回disp" 初始向量迭代失敗 "return;endif abs(r-r0)<eps% 判斷第一次迭代后是否滿足要求，如滿足，則返回結(jié)果returnendwhile abs(r-r0)>eps && k<n %終止條件k=k+1;a0=a1;a1=a2;r0=r;y=x./max(abs(x);x=A*y;%迭代格式a2=max(a

13、bs(x);if (a2-2*a1+a0)=0 %若分母為 0，則迭代失敗，返回return;endr=a0-(a1-a0)2/(a2-2*a1+a0);m,index=max(abs(eig(A); %以下代碼保證取出來的按模最大特征值aa=eig(A);%是原值，而非其絕對值。if aa(index)>0 |aa(index)=0r=r;elser=-r;endend精彩文檔實用標準文案end類似可得按模最小特征值和特征向量的代碼如下：與上面類似，所不同的只是迭代格式不同 .functionr,y=invaitken(A,x0,eps,n)k=1;a0=0;a1=1;r0=1;y=x

14、0./max(abs(x0);L,U=lu(A);%迭代格式的不同z=Ly;x=Uz;a2=max(abs(x);r=a0-(a1-a0)2/(a2-2*a1+a0);i f (a2-2*a1+a0)=0disp" 初始向量迭代失敗 "return;endifabs(r-r0)<eps% 判斷第一次迭代后是否滿足要求，如滿足，則返回結(jié)果returnendwhile abs(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*a

15、1+a0)=0return;endr=a0-(a1-a0)2/(a2-2*a1+a0);endm,index=min(abs(eig(A);%以下代碼保證取出來的按模最大特征值aa=eig(A);%是原值，而非其絕對值。if aa(index)>0 |aa(index)=0r=1/r;elser=-1/r;endend精彩文檔實用標準文案4. 計算 Hilb 矩陣特征值此處不再舉例，而是直接應用于15階Hilb 矩陣，初始向量選為 ones(15,1) ，結(jié)果如下，并將結(jié)果與冪法和反冪法得到結(jié)果比較這與冪法得到的特征值和特征向量一致， 表明算法和代碼正確； 同理，最小特征值結(jié)果如下：精彩

16、文檔實用標準文案這與反冪法得到的結(jié)果一致，表明結(jié)果正確。五，對稱矩陣的 Rayleigh 商加速法1. 簡介與原理為對稱矩陣，設x0,則稱 ()xT Ax 為關于 A的Rayleigh商AR xxT x原理如下：設Xi0為 i的特征向量，即 AXii Xi用Xi 左乘上式有：y( k)x(k )max( x( k) )x( k 1)Ay(k )這稱為商加速法。RayleighR( y(k ) )( y(k ) )T x( k 1)( y(k ) )T ( y(k ) )其中 R(y(k ) )1, y(k )x(k )max(x( k) )精彩文檔實用標準文案2. 算法實現(xiàn)輸入矩陣A,初始向量

17、x,誤差限,最大迭代次數(shù)N ,(1).置k1, r00, yx,( 2).max(abs(x)(3).xA * y, ryT xyT y若| rr0 |,輸出r , y,停機否則轉(zhuǎn)(5),( 4).,若kN,置kk1, r0r, yx,轉(zhuǎn)(3);(5).max(abs( x)否則輸出失敗信息 , 停機 .3. Matlab 程序代碼function r,y=rayleigh(A,x0,eps,n) % r 是特征值， y是特征向量 k=1; r0=0;y=x0./max(abs(x0);x=A*y;%迭代格式計算新的 xr=dot(y,x)/dot(y,y);% Reyleigh商if abs(r-r0)<epsreturnendwhile abs(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 商加速法，如下：function r,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);if abs(r-r0)<epsreturnendwhile