幂法,反幂法求解矩阵最大最小特征值及其对应的特征向量

1、数值计算解矩阵的按模最大最小特征值及对应的特征向量一.哥法1.1.幕法简介：当矩阵A满足一定条件时，在工程中可用幕法计算其主特征值(按模最大)及其特征向量。矩阵A需要满足的条件为：|卜|九2户旬K户0,%为A的特征值存在n个线性无关的特征向量，设为Xi，X2，.,Xn1.11.1计算过程：对任意向量x(0),x(k1)=Ax(k)=.二nnk1=A码=Zi1i1n有x(0)=ZctiUi,%不全为0,则有i1Ak1x(0)、k1%Au产科+$口+e严可见，当|上1越小时，收敛越快；且当k充分大时，有,1x(k1)”U1x(k)”UI(k1)1,对应的特征向量即是x。2 2算法实现.输入矩阵A,

2、初始向量X,.k，1,0;y(k)误差限巴最大迭代次数Nx(k).计算x Ay,；max(abs(x(k)max(x);(4)若|BH，输由P,y,否叫转(5)若kN,置k-k+1/转3,否则输由失败信息，停机.3matlab3matlab程序代码functiont,y=lpowerA,x0,eps,N)%t为所求特征值，y是对应特征向量k=1;z=0;%z相当于黑y=x0./max(abs(x0);%规范化初始1可量x=A*y;%迭代格式b=max(x);%b相当于Pifabs(z-b)eps&kt,(A,xOjeps,、)eig(A)ans=-0.0166L48012.5365A*

3、y-t*yans-1.0e-004+-0.1603-0,16840结果正确，表明算法和代码正确，然后利用此程序计算15阶Hilb矩阵，与eig(A)的得到结果比较，再计算 A*y-t*y,验证y是否是对应的特征向量。设置初始向量为 x0=ones(15,1),结果显示如下A=hilb(15);x0=ones(15,1):eps-le6：N=30;Lt,y-1power(A,x0,eps,t=L84592.53650.74820.6497L0000V-eig(A)A*y-t*y1.0000ans=ans=0.62280.47060.3838-0.00DO0.00001.0-007拿00.3264

4、0.0000-0.55160.0000-0.64360.28510.0000-0.64650.25370.0000-0.62450.22900.0000-0.59530.20890.oooo-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-可见，结果正确。得到了15阶Hilb矩阵的按模最大特征值和对应的特征向二.反哥法1 1.反幕法简介及其理论在工程计算中，可以利用反

5、幕法计算矩阵按模最小特征值及其对应特征向量。其基本理论如下，与幕法基本相同：ii1由Ax=，x=x=A(九x),则Ax=x,可知，A和A的特征值互为倒数，儿求A按模最小特征值即求A-1的按模最大特征值，取倒数即为A的按模最小特征值所以算法基本相同，区别就是在计算x(k*时，不是令x(kW=Ay(k),而是x(E=A-1y(k)具体计算时，变换为Ax(k+)=y(k);对A彳故LU分解,来计算x(k41)(2)法实现(7).若kcN,置k-k+1,%K,y-max(abs(x)，转(4)；否则输出失败信息，停机.3 3matlabmatlab程序代码functions,y=invpower(A,

6、x0,eps,n)%s为按模最小特征值，y是对应特征向量k=1;r=0;%r相当于-0y=x0./max(abs(x0);%规范化初始向量L,U=lu(A);z=Ly;x=Uz;u=max(x);s=1/u;%按模最小为A-1按模最大的倒数.ifabs(u-r)eps&knk=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);end%终止条件.%这两步保证取出来的按模最大特征值%是原值，而非其绝对值。4 4举例验证同事法一样，选取一个矩阵A,代入程序，得到结果，并与eig(A)的得

7、到结果比较，再计算 A*y-t*y,验证 y 是否是对应的特征向量.输入矩阵A,初始向量x,误差限以最大迭代次数N,(3).作三角分解A=LU(4).解方程组LUx=y(Lz=y,Ux=z),(5) .max(x),:一1(6).右|九-%|s,y=invpower(A?epsn)-5.9673e-018里。三.计算条件数矩阵A的条件数等于A的范数与A的逆的范数的乘积，即cond(A)=IIAll11AA(-1)II,对应矩阵的3种范数，可以定义3种条件数。函数cond(A,1)、cond(A)或cond(Ainf)是判断矩阵病态与否的一种度量，条件数越大表明矩阵的病态程度越大.这里我们选择矩

8、阵的2范数，即cond(A)=不?,%,%分别为矩阵ATA的最大和最小特征值而如果A为对称矩阵，如Hilb矩阵，ATA的最大最小特征值，分别为A的最大最小特征值的平方。所以cond(A)为A的最大最小特征值得比值。对于本例中的15阶Hilb矩阵来说， 利用上面计算结果得其条件数(选择第二种条件数)为： 3.0934e+017;这与直接利用cond(A)得到的结果：2.5083e+017 在同一0.0000-0.00000.0000-0.00060.0050-0.02450.0699-0.0968-0.04210.4497-0.9084L0000-0.65270.2379-0.0375eig(A

9、)ans=-0,00000.00000.00000.00000,00000.00000.00000.00000.00000.00000.00040.00560.05720.4266A*y-s*y二L0e-017*0.3903-0.56380.0000-0.5208-0.47410.35400.4537-0.27460.9724-0.6556-0.62880.66180.0639可见，结果真确。得到了15阶Hilb矩阵的按模最大特征值和对应的特征向数量级，再次表明了上述算得得最大最小特征值的正确性，同时又表明阵是病态矩阵。Hilb矩四.Aitken.Aitken商加速法1,简介与原理若ak收敛与

10、a,且lim叫二a=c#0;即瓜践性U敛，jak-a当k充分大时，有如二ak-a(X,yn2=Xn2Xn2.烝2-aak1-a)2一2xn1xn二 a.ak-( (4 4ak2-2ak1ak用ak逼近a,这种方法称为Aitken加速法.同事法和反幕法计算最大和最小特征值类似,如果计算最大特征值,为x(k*)=Ay(k)；计算最小特征值时，迭代格式为则迭代格式x(E=A*y(k),即y(k)=Ax”)。2,算法实现计算按模最大特征值算法如下：,输入A=(aj),初始向量x,误差限以最大迭代次数N,x.置k=1,ao=00。.。,丫二加,计算x=Ay,置max(x)=a2,2(-?)(4),计算。

11、-一一屋，:2一21.%巴输出七y停机，否则转(6),(6),若卜N,置。1n0,32na1,九二九0,1工.y转,否则， 输出失败信息，停机.类似幕法和反幕法可以写出按模最小特征值算法，此处不再赘述3.matlab程序代码functionr,y=aitken(A,x0,eps,n)%r按模最大特征值，丫为对应特征向量k=1;a0=0;%a相当于0a1=1;%a1相当于Ir0=1;%相当于2中的八0y=x0./max(abs(x0);%规范化初始向量x=A*y;a2=max(abs(x);%a2相当于二2r=a0-(a1-a0)A2/(a2-2*a1+a0);%相当于九if(a2-2*a1+a

12、0)=0%若上式中分母为0,则迭代失败，返回disp初始向量迭代失败return;endifabs(r-r0)eps&k0|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=max(abs(x);r=a0-(a1-a0)A2/(a2-2*a1+a0);if(a2-2*a1+a0)=0disp

13、初始向量迭代失败return;endifabs(r-r0)eps&k0|aa(index)=0r=1/r;elser=-1/r;endend4.4.计算HilbHilb矩阵特征值此处不再举例， 而是直接应用于15阶Hilb矩阵， 初始向量选为ones(15,1)ones(15,1)结果如下，并将结果与幕法和反幕法得到结果比较L00000.62290.47060,38380.32640.28510.25370.22900.20890.19220.17800.16590.15540.14620.1381这与幕法得到的特征值和特征向量一致，表明算法和代码正确；同理，最小特征值结果如下：r,y

14、=invaitken(A,x0fepSjn)5.9673bo18A=hilb(15);xO=ones(13,l):eps=le-6;n=30;y_=aitken(AJxOeps,n)1.8459这与反幕法得到的结果一致，表明结果正确。五，对称矩阵的 Rayleigh 商加速法1,简介与原理xTAx.、.一、A为对称矩阵，设x#0,则称WxXAx为关于A白Rayleigh冏xx原理如下:设 Xj#0 为的特征向量，即 AXj=%Xj用 Xj左乘上式有：0.00000.00000,00000.00060.0050-0.02450.0699-0.096S-0,04210.44977.9084L000

15、0-0.65270.23790.0375(k)y(k1)xR(y(k)(k)xmax(x(k)=Ay(k)_(y(k)Tx(kP(k)T/(k)(y)(y)这称为 Rayleigh 商加速法。其中 R(y(k)，Z,y(k)，(k)Xmax(x(k)2.算法实现.输入矩阵A,初始向量x,误差限名,最大迭代次数N,x(2).置k，1,r0，0,y，,max(abs(x)TyxTyy输出r,y,停机,否则转(5),x(5).右kN,置k.k+1,r0 r,y.,转(3)max(abs(x)否则输出失败信息，停机.3 3.Matlab程序代码functionr,y=rayleigh(A,x0,eps

16、,n)%r是特征值，y是特征向量k=1;r0=0;y=x0./max(abs(x0);x=A*y;%迭代格式计算新的xr=dot(y,x)/dot(y,y);%Reyleigh商ifabs(r-r0)eps&knk=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)eps&k_r,y=rayleigh(Ax0Pepsn)1,8459L00000.62290.47060.383