幂法和反幂法求矩阵特征值_课程设计.doc

题目幂法和反幂法求矩阵特征值 课程设计具体内容随机产生一对称矩阵，对不同的原点位移和初值(至少取3个)分别使用幂法求计算矩阵的主特征值及主特征向量，用反幂法求计算矩阵的按模最小特征值及特征向量，并比较不同的原点位移和初值说明收敛。要求1.认真读题，了解问题的数学原形；2.选择合适问题求解的数值计算方法；3.设计程序并进行计算；4.对结果进行解释说明；采用方法及结果说明对于幂法和反幂法求解矩阵特征值和特征向量的问题将从问题分析，算法设计和流程图，理论依据，程序及结果进行阐述该问题。一问题的分析：求n阶方阵a的特征值和特征向量，是实际计算中常常碰到的问题，如：机械、结构或电磁振动中的固有值问题等。对于n阶矩阵a，若存在数和n维向量x满足 ax=x （1）则称为矩阵a的特征值，x为相应的特征向量。由高等代数知识可知，特征值是代数方程 |i-a|=+a+a+a=0 （2）的根。从表面上看，矩阵特征值与特征向量的求解问题似乎很简单，只需求解方程（2）的根，就能得到特征值，再解齐次方程组 （i-a）x=0 （3）的解，就可得到相应的特征向量。上述方法对于n很小时是可以的。但当n稍大时，计算工作量将以惊人的速度增大，并且由于计算带有误差，方程（2）未必是精确的特征方程，自然就不必说求解方程（2）与（3）的困难了。幂法是一种计算矩阵主特征值（矩阵按模最大的特征值）及对应特征向量的迭代方法，特别是用于大型稀疏矩阵。反幂法是计算海森伯格阵或三角阵的对应一个给定近似特征值的特征向量的有效方法之一。二算法设计及流程图1、幂法算法（1）取初始向量u（例如取u=(1,1,1)）,置精度要求，置k=1. （2）计算v=au,m=max(v), u= v/ m（3）若| m= m|，则停止计算（m作为绝对值最大特征值，u作为相应的特征向量）否则置k=k+1，转（2）2、反幂法算法（1）取初始向量u（例如取u=(1,1,1)）,置精度要求，置k=1. （2）对a作lu分解，即a=lu（3）解线性方程组 ly=u,uv=y（4）计算 m=max(v), u= v/ m（5）若|m=m|，则停止计算（1/m作为绝对值最小特征值，u作为相应的特征向量）;否则置k=k+1，转（3）.幂法流程图：开始输入a；m,u,index=pow(a,1e-6)k=0;m1=0v=a*uvmax,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*uvmax,i=max(abs(v)m=v(i);u=v/mabs(m-m1)|则计算最大特征值与特征向量的迭代格式为v=au,m=max(v), u= v/ m （1）其中max(v)表示向量v绝对值的最大分量。2、对于幂法的定理按式（1）计算出m和u满足 m=, u=（二）反幂法算法的理论依据及推导反幂法是用来计算绝对值最小的特征值忽然相应的特征向量的方法。是对幂法的修改，可以给出更快的收敛性。1、反幂法的迭代格式与收敛性质设a是非奇异矩阵，则零不是特征值，并设特征值为|则按a的特征值绝对值的大小排序，有 |对a实行幂法，就可得a的绝对值最大的特征值1/和相应的特征向量，即a的绝对值最小的特征值和相应的特征向量。由于用a代替a作幂法计算，因此该方法称为反幂法，反幂法的迭代格式为 v= au,m=max(v), u= v/ m （2）2、对于反幂法的定理按式（2）计算出的m和u满足： m=, u=在式（2）中，需要用到a，这给计算带来很大的不方便，因此，把（2）式的第一式改为求解线性方程组 a v= u （3）但由于在反幂法中，每一步迭代都需求解线性方程组（3）式，迭代做了大量的重复计算，为了节省工作量，可事先把矩阵a作lu分解，即 a=lu所以线性方程组（3）改为 ly=u,uv=y四、算法程序设计代码幂法程序，在matlab中建立一个m文件并保存。%pow.mfunction m,u,index,k=pow(a,u,ep,it_max)if nargin4 it_max=1000;endif nargin3 ep=1e-5;endn=length(a);index=0;k=0;m1=0;m0=0; i=eye(n);t=a-m0*i;while k=it_max v=t*u;vmax,i=max(abs(v); m=v(i); u=v/m; if abs(m-m1)ep; index=1; break; end m=m+m0; m1=m; k=k+1;end在matlab输入面板，输入a=rand（4）；%产生一个4维随机矩阵b=a+a；u=1 1 1 1;%设立初始向量m,u,index,k=pow(b，u,ep,it_max)%最多可省略2个参数程序结束。在m文件中可以通过改变m0的值改变原点位移，从而达到原点位移加速。反幂法程序设计代码：在matlab中建立一个m文件并保存。%pow_inv.mfunctionm,u,index,k=pow_inv(a,u,ep,it_max)if nargin4 it_max=1000;endif nargin3 ep=1e-5;endn=length(a);index=0;k=0;m1=0;m0=0; i=eye(n);t=a-m0*i;invt=inv(t);while k=it_max v=invt*u; vmax,i=max(abs(v); m=v(i); u=v/m; if abs(m-m1)b=rand(4);a=b+ba =0.2675 0.5776 0.6344 1.3130 0.5776 1.1503 0.7641 0.1367 0.6344 0.7641 0.0257 0.4193 1.3130 0.1367 0.4193 1.2248 u=1 1 1 1; m,u,index,k=pow(a,u)m = 2.6813u = 0.8576 0.6934 0.5623 1.0000index = 1k =49修改m0=1e-3m = 2.6814u = 0.8576 0.6934 0.5623 1.0000index = 0k = 1001修改m0=0 %此时为幂法m = 2.6815u = 0.8576 0.6935 0.5623 1.0000index = 1k = 10修改u=1 2 3 4修改m0=1e-4m = 2.6813u = 0.8576 0.6934 0.5623 1.0000index = 1k = 9修改m0=1e-3m = 2.6805u = 0.8576 0.6934 0.5622 1.0000index = 1k = 7修改m0=0m = 2.6814u = 0.8576 0.6934 0.5623 1.0000index = 1k = 9修改u=3 5 6 7修改m0=1e-4m = 2.6819u = 0.8577 0.6937 0.5624 1.0000index = 1k = 7修改m0=1e-3m = 2.6814u = 0.8576 0.6934 0.5623 1.0000index = 0k = 1001修改m0=0m = 2.6820u = 0.8577 0.6937 0.5624 1.0000index = 1k = 7总结以上,幂法如下：um0muindexk1 1 1 10.00012.68130.8576 0.6934 0.5623 1.00001490.0012.68140.5876 0.6934 0.5623 1.00000100102.68150.8576 0.6935 0.5623 1.00001101 2 3 40.00012.68130.8576 0.6934 0.5623 1.0000190.0012.68050.8576 0.6934 0.5622 1.00001702.68140.8576 0.6934 0.5623 1.0000193 5 6 70.00012.68190.8577 0.6937 0.5624 1.0000170.0012.69140.8576 0.6934 0.5623 1.00000100102.6920.8577 0.6937 0.5624 1.000017反幂法结果显示：在m0为0时m0=0.001 u=1 1 1 1m0=0.1 u=1 1 1 1m0=0 u=1 3 5 7m0=0.1 u=1 3 5 7m0=0.5 u=1 3 5 7m0=0 u=2 3 4 5m0=0.1 u=2 3 4 5m0=0.7 u=2 3 4 5综上，反幂法结果如下：um0muindexk1 1 1 10.10.3847-0.8996 1.0000 0.2726 -0.23641150.0010.3847-0.8996 1.0000 0.2726 -0.236411600.3847-0.8996 1.0000 0.2726 -0.23641161 3 5 70.50.3847-0.8995 1.0000 0.2726 -0.23641270.10.3847-0.8996 1.0000 0.2726 -0.236411700.3847-0.8996 1.0000 0.2726 -0.23641202 3 4 50.70.7091-0.6962 -0.4497 0.2196 1.0000150.10.3847-0.8995 1.0000 0.2726 -0.2