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

1、数值计算解矩阵的按模最大最小特征值及对应的特征向量一.哥法1 .幕法简介：当矩阵A满足一定条件时,在工程中可用幕法计算其主特征值(按模最大)及其特征向量.矩阵A需要满足的条件为：(1) |l|2|n|0一为A的特征值(2)存在n个线性无关的特征向量,设为Xi,X2,.,Xn1.1 计算过程：n对任意向量X(0),有X(0)iUi,i不全为0,那么有i1x(k1)Ax的Ak1x(0)nnAk1码q*1Uii1i1、k12k1n、k111U1()a2U2(-)anUn11k1u11U1可见,当2|越小时,收敛越快；且当k充分大时,有1(k1)XX(k)1U1x(k1)1U1(k1)1,对应的特征向

2、量即是X.2算法实现.输入矩阵A,初始向量.k1,0;y(k)x,误差限,最大迭代次数Nx(k).计算X.假设|Ay,1|,输由max(abs(x(k)max(x);,y,否那么,转(5)假设kN,置kk1,转3,否那么输由失败信息,停机.3matlab程序代码functiont,y=lpowerA,x0,eps,N)%t为所求特征值,y是对应特征向量k=1;z=0;%z相当于y=x0./max(abs(x0);%标准化初始向量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相当于判断第一次迭代后是否满足要求这两步保证取出来的按模最大特征值是原值,而非其绝对值.4举例验证选取一个矩阵A,代入程序,得到结果,并与eig(A)的得到结果比拟,再计算A*y-t*y,验证y是否是对应的特征向量.结果如下：结果正确,说明算法和代码正确,然后利用此程序计算15阶Hilb矩阵,与eig(A)的得到结果比拟,再计算A*y-t*y,验证y是否是对应的特征向量.设置初始向量为x0=ones(15,1),结果

4、显示如下可见,结果正确.得到了15阶Hilb矩阵的按模最大特征值和对应的特征向量.二.反哥法1 .反幕法简介及其理论在工程计算中,可以利用反幕法计算矩阵按模最小特征值及其对应特征向量.其根本理论如下,与幕法根本相同：ii1由AxxxA(x),那么Axx,可知,A和A1的特征值互为倒数,求A按模最小特征值即求A1的按模最大特征值,取倒数即为A的按模最小特征值所以算法根本相同,区别就是在计算x(k1)时,不是令x(k1)Ay(k),而是x(k1)A-1y(k)具体计算时,变换为Ax(k1)y(k);对A做LU分解,来计算x(k1)2 .算法实现.输入矩阵A,初始向量x,误差限,最大迭代次数N,(2

5、) .置k1,00,yx,max(abs(x)(3) .作三角分解ALU(4) .解方程组LUxy(Lzy,Uxz),(5) .max(x),一一,1(6) .右|0|,输出一,y,停机,否那么转,(7) .假设kN,置kk1,0,y/、,转(4);max(abs(x)否那么输出失败信息,停机.3matlab程序代码functions,y=invpower(A,x0,eps,n)%s为按模最小特征值,y是对应特征Ik=1;r=0;%r相当于0y=x0./max(abs(x0);%标准化初始向量L,U=lu(A);z=Ly;x=Uz;u=max(x);s=1/u;%按模最小为A1按模最大的倒数.

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,代入程序,得到结果,并与eig(A)的得到结果比较,再计算A*y-t*y,验证y是否是对应的特征向量.可见结果正确,然后利用此程序计算15阶Hilb矩阵,eig(A)的

7、得到结果比较,再计算A*y-s*y,验证y是否是对应的特征向量.设置初始向量为x0=ones(15,1),结果显示如下可见,结果真确.得到了15阶Hilb矩阵的按模最大特征值和对应的特征向量.三.计算条件数矩阵A的条件数等于A的范数与A的逆的范数的乘积,即cond(A)=IIAll11AA(-1)II,对应矩阵的3种范数,可以定义3种条件数.函数cond(A,1)、cond(A)或cond(Ainf)是判断矩阵病态与否的一种度量,条件数越大说明矩阵的病态程度越大.这里我们选择矩阵的2范数,即cond(A),1,2分别为2矩阵ATA的最大和最小特征值而如果A为对称矩阵,如Hilb矩阵,ATA的最

8、大最小特征值,分别为A的最大最小特征值的平方.所以cond(A)为A的最大最小特征值得比值.对于本例中的15阶Hilb矩阵来说,利用上面计算结果得其条件数(选择第二种条件数)为：3.0934e+017；这与直接利用cond(A)得到的结果：2.5083e+017在同一数量级,再次说明了上述算得得最大最小特征值的正确性,同时又说明Hilb矩阵是病态矩阵.四.Aitken商加速法1 .简介与原理假设ak收敛与a,且imak1aakac0;即ak线性收敛,当k充分大时,有曳工招曳工招akaak1ayn2xn2(xn2xn1)xn22xn1xnaak(ak1aQ2ak22ak1用ak逼近a,这种方法称

9、为一:仇akAitken加速法.同事法和反幕法计算最大和最小特征值类似,如果计算最大特征值,那么迭代格式为x(k1)Ay(k)；计算最小特征值时,迭代格式为x(k1)A1y(k),即y(k)Ax(k1).2 .算法实现计算按模最大特征值算法如下：输入Aa),初始向量x,误差限,最大迭代次数N,3 2).置k1,o0,i0,o1.0,y(3) .计算x(4) .计算max(x)Ay,置max(x)(1o)2022210.假设(6).假设k0,输出,y停机,否那么转(6),N,置10,21,0,k1k,七通否那么,输出失败信息,停机.类似幕法和反幕法可以写出按模最小特征值算法,此处不再赘述3 .m

10、atlab程序代码functionr,y=aitken(A,x0,eps,n)%r按模最大特征值,y为对应特征向量k=1;a0=0;%a相当于0a1=1;%a1相当于1r0=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+a0)=0%假设上式中分母为0,那么迭代失败,返回disp"初始向量迭代失败"return;endifabs(r-r0)<eps%判断第一次迭代后是否满足要求,如满足,那么返回结果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%假设分母为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%判断第一次迭代后是否满足要求,如满足,那么返回结果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矩阵特征值此处不再举例,而是直接应用于15阶Hilb矩阵,初始向量选为ones(15,1)

14、,结果如下,并将结果与幕法和反幕法得到结果比拟这与幕法得到的特征值和特征向量一致,说明算法和代码正确；同理,最小特征值结果如下：这与反幕法得到的结果一致,说明结果正确.五,对称矩阵的Rayleigh商加速法1,简介与原理A为对称矩阵,设x0,那么称R(x)#x为关于A白Rayleigh商xx原理如下：(k)y(k1)XR(y(k)其中R(y(k)1,y(k)(k)Xmax(x(k)2.算法实现.输入矩阵A,初始向量x,误差限,最大迭代次数N,.置k1,r0.xA*y,r假设|rr.|0,yTyxTyyxmax(abs(x),输出r,y,停机,否那么转5,设Xi0为i的特征向量,即AXiiXi用

15、Xj左乘上式有：(k)Xmax(x(k)这称为Rayleigh商加速法.Ay(k)(y(k)Tx(k1)(y(k)T(y(k)假设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)<