数值分析幂法和反幂法

题目幂法和反幂法求矩阵特征值具体内容随机产生一对称矩阵，对不同的原点位移和初值(至少取3个)分别使用幂法求计算矩阵的主特征值及主特征向量，用反幂法求计算矩阵的按模最小特征值及特征向量，并比较不同的原点位移和初值说明收敛。要求1.认真读题，了解问题的数学原形；2.选择合适问题求解的数值计算方法；3.设计程序并进行计算；4.对结果进行解释说明；采用方法及结果说明对于幂法和反幂法求解矩阵特征值和特征向量的问题将从问题分析，算法设计和流程图，理论依据，程序及结果进行阐述该问题。一．问题的分析：求n阶方阵A的特征值和特征向量，是实际计算中常常碰到的问题，如：机械、结构或电磁振动中的固有值问题等。对于n阶矩阵A，若存在数和n维向量x满足Ax=x（1）则称为矩阵A的特征值，x为相应的特征向量。由高等代数知识可知，特征值是代数方程|I-A|=n+a11n+…+a1n+an=0（2）的根。从表面上看，矩阵特征值与特征向量的求解问题似乎很简单，只需求解方程（2）的根，就能得到特征值，再解齐次方程组（I-A）x=0（3）的解，就可得到相应的特征向量。上述方法对于n很小时是可以的。但当n稍大时，计算工作量将以惊人的速度增大，并且由于计算带有误差，方程（2）未必是精确的特征方程，自然就不必说求解方程（2）与（3）的困难了。幂法是一种计算矩阵主特征值（矩阵按模最大的特征值）及对应特征向量的迭代方法，特别是用于大型稀疏矩阵。反幂法是计算海森伯格阵或三角阵的对应一个给定近似特征值的特征向量的有效方法之一。二．算法设计及流程图1、幂法算法（1）取初始向量u)0(（例如取u)0(=(1,1,…1)T）,置精度要求，置k=1.（2）计算v)(k=Au)1(k,mk=max(v)(k),u)(k=v)(k/mk（3）若|mk=m1k|，则停止计算（mk作为绝对值最大特征值1，u)(k作为相应的特征向量）否则置k=k+1，转（2）2、反幂法算法（1）取初始向量u)0(（例如取u)0(=(1,1,…1)T）,置精度要求，置k=1.（2）对A作LU分解，即A=LU（3）解线性方程组Ly)(k=u)1(k,Uv)(k=y)(k（4）计算mk=max(v)(k),u)(k=v)(k/mk（5）若|mk=m1k|，则停止计算（1/mk作为绝对值最小特征值n，u)(k作为相应的特征向量）;否则置k=k+1，转（3）.幂法流程图：开始输入A；[m,u,index]=pow(A,1e-6)k=0;m1=v=A*u[vmax,i]=max(abs(v))m=v(i);u=v/mabs(m-m1)1e-6index=1;break;输出：m，u，index结束m1=m;k=k+1反幂法流程图开始输入A；[m,u,index]=pow_inv(A,1e-6)k=0;m1=0v=invA*u[vmax,i]=max(abs(v))m=v(i);u=v/mabs(m-m1)1e-6index=1;break;输出：m，u，index结束m1=m;k=k+1输入A；[m,u,index]=pow(A,1e-6)三、算法的理论依据及其推导（一）幂法算法的理论依据及推导幂法是用来确定矩阵的主特征值的一种迭代方法，也即，绝对值最大的特征值。稍微修改该方法，也可以用来确定其他特征值。幂法的一个很有用的特性是它不仅可以生成特征值，而且可以生成相应的特征向量。实际上，幂法经常用来求通过其他方法确定的特征值的特征向量。1、幂法的迭代格式与收敛性质设n阶矩阵A的特征值1，2,…,n是按绝对值大小编号的，xi(i=1,2,…,n)为对应i的特征向量，且1为单根，即|1||2|≥…≥|n|则计算最大特征值与特征向量的迭代格式为v)(k=Au)1(k,mk=max(v)(k),u)(k=v)(k/mk（1）其中max(v)(k)表示向量v)(k绝对值的最大分量。2、对于幂法的定理按式（1）计算出mk和u)(k满足klimmk=1,klimu)(k=)max(11xx（二）反幂法算法的理论依据及推导反幂法是用来计算绝对值最小的特征值忽然相应的特征向量的方法。是对幂法的修改，可以给出更快的收敛性。1、反幂法的迭代格式与收敛性质设A是非奇异矩阵，则零不是特征值，并设特征值为|1|≥|2|≥…≥|1n||n|则按A1的特征值绝对值的大小排序，有|n1||11n|≥…≥|11|对A1实行幂法，就可得A1的绝对值最大的特征值1/n和相应的特征向量，即A的绝对值最小的特征值和相应的特征向量。由于用A1代替A作幂法计算，因此该方法称为反幂法，反幂法的迭代格式为v)(k=A1u)1(k,mk=max(v)(k),u)(k=v)(k/mk（2）2、对于反幂法的定理按式（2）计算出的mk和u)(k满足：klimmk=n1,klimu)(k=)max(nnxx在式（2）中，需要用到A1，这给计算带来很大的不方便，因此，把（2）式的第一式改为求解线性方程组Av)(k=u)1(k（3）但由于在反幂法中，每一步迭代都需求解线性方程组（3）式，迭代做了大量的重复计算，为了节省工作量，可事先把矩阵A作LU分解，即A=LU所以线性方程组（3）改为Ly)(k=u)1(k,Uv)(k=y)(k四、算法程序设计代码幂法程序，在matlab中建立一个M文件并保存。%pow.mfunction[m,u,index,k]=pow(A,u,ep,it_max)ifnargin4it_max=1000;endifnargin3ep=1e-5;endn=length(A);index=0;k=0;m1=0;m0=0;I=eye(n);T=A-m0*I;whilek=it_maxv=T*u;[vmax,i]=max(abs(v));m=v(i);u=v/m;ifabs(m-m1)ep;index=1;break;endm=m+m0;m1=m;k=k+1;end在matlab输入面板，输入A=rand（4）；%产生一个4维随机矩阵B=A+A’；u=[1111]’;%设立初始向量[m,u,index,k]=pow(B，u,ep,it_max)%最多可省略2个参数程序结束。在M文件中可以通过改变m0的值改变原点位移，从而达到原点位移加速。反幂法程序设计代码：在matlab中建立一个M文件并保存。%pow_inv.mfunction[m,u,index,k]=pow_inv(A,u,ep,it_max)ifnargin4it_max=1000;endifnargin3ep=1e-5;endn=length(A);index=0;k=0;m1=0;m0=0;I=eye(n);T=A-m0*I;invT=inv(T);whilek=it_maxv=invT*u;[vmax,i]=max(abs(v));m=v(i);u=v/m;ifabs(m-m1)epindex=1;break;endm1=m;k=k+1;endm=1/m;m=m+m0;在matlab输入面板，输入A=rand（4）；%产生一个4维随机矩阵B=A+A’；u=[1111]’;%设立初始向量[m,u,index,k]=pow_inv(B，u,ep,it_max)%最多可省略2个参数程序结束。在M文件中可以通过改变m0的值改变原点位移，从而达到原点位移加速。【结果显示】%在M0=1e-4B=rand(4);A=B+B’A=0.26750.57760.63441.31300.57761.15030.76410.13670.63440.76410.02570.41931.31300.13670.41931.2248u=[1111]';[m,u,index,k]=pow(A,u)m=2.6813u=0.85760.69340.56231.0000index=1k=49修改M0=1e-3m=2.6814u=0.85760.69340.56231.0000index=0k=1001修改M0=0%此时为幂法m=2.6815u=0.85760.69350.56231.0000index=1k=10修改U=[1234]修改M0=1e-4m=2.6813u=0.85760.69340.56231.0000index=1k=9修改M0=1e-3m=2.6805u=0.85760.69340.56221.0000index=1k=7修改M0=0m=2.6814u=0.85760.69340.56231.0000index=1k=9修改U=[3567]修改M0=1e-4m=2.6819u=0.85770.69370.56241.0000index=1k=7修改M0=1e-3m=2.6814u=0.85760.69340.56231.0000index=0k=1001修改M0=0m=2.6820u=0.85770.69370.56241.0000index=1k=7总结以上,幂法如下：Um0muindexk[1111]0.00012.6813[0.85760.69340.56231.0000]1490.0012.6814[0.58760.69340.56231.0000]0100102.6815[0.85760.69350.56231.0000]110[1234]0.00012.6813[0.85760.69340.56231.0000]190.0012.6805[0.85760.69340.56221.0000]1702.6814[0.85760.69340.56231.0000]19[3567]0.00012.6819[0.85770.69370.56241.0000]170.0012.6914[0.85760.69340.56231.0000]0100102.692[0.85770.69370.56241.0000]17反幂法结果显示：在m0为0时M0=0.001U=[1111]M0=0.1u=[1111]M0=0u=[1357]M0=0.1u=[1357]M0=0.5u=[1357]M0=0u=[2345]M0=0.1u=[2345]M0=0.7u=[2345]综上，反幂法结果如下：um0muindexk[1111]0.10.3847[-0.89961.00000.2726-0.2364]1150.0010.3847[-0.89961.00000.2726-0.2364]11600.3847[-0.89961.00000.2726-0.2364]116[1357]0.50.3847[-0.89951.00000.2726-0.2364]1270.10.3847[-0.89961.00000.2726-0.2364]11700.3847[-0.89961.00000.2726-0.2364]120[2345]0.70.7091[-0.6962-0.44970.21961.0000]150.10.3847[-0.89951.00000.2726-0.2364]11700.3847[-0.89961.00000.2726-0.2364]119五、结果分析采用幂法和反幂法，求矩阵的最大和最小特征值，从原理上看，这两种方法都是迭代法，因此迭代初始向量的选择对计算结果会产生一定影响，主要表现在收敛速度上。同时，原点位移m的选取也影响收敛的速度。但原点位移m0的适当选取依赖于对矩阵A的大致了解。成员1007024104辛志贤1007024107张容1007024108罗言月