数值分析报告幂法和反幂法

1、题目幕法和反幕法求矩阵特征值具/、 体 内 容随机产生一对称矩阵，对不同的原点位移和初值 （至少取3个）分别使用幕 法求计算矩阵的主特征值及主特征向量，用反幂法求计算矩阵的按模最小特征 值及特征向量，并比较不同的原点位移和初值说明收敛。要求1. 认真读题，了解问题的数学原形；2. 选择合适冋题求解的数值计算方法；3. 设计程序并进行计算；4. 对结果进行解释说明；采 用 方 法 及 结 果 说 明对于幕法和反幕法求解矩阵特征值和特征向量的问题将从问题分析，算 法设计和流程图，理论依据，程序及结果进行阐述该问题。一.问题的分析：求n阶方阵A的特征值和特征向量，是实际计算中常常碰到的问题，如： 机

2、械、结构或电磁振动中的固有值冋题等。 对于n阶矩阵A，若存在数九和 n维向量x满足Ax=x（1）则称入为矩阵A的特征值，x为相应的特征向量。由高等代数知识可知，特征值是代数方程|丸 l-A|=九 n+a, k n_l+azk+an=O（2）的根。从表面上看，矩阵特征值与特征向量的求解问题似乎很简单，只需 求解方程（2）的根，就能得到特征值人，再解齐次方程组（k I-A ） x=0（3）的解，就可得到相应的特征向量。上述方法对于n很小时是可以的。但当n稍大时，计算工作量将以惊 人的速度增大，并且由于计算带有误差，方程（2）未必是精确的特征方程， 自然就不必说求解方程（2）与（3）的困难了。幕法是

3、一种计算矩阵主特 征值（矩阵按模最大的特征值）及对应特征向量的迭代方法，特别是用于 大型稀疏矩阵。反幕法是计算海森伯格阵或三角阵的对应一个给定近似特 征值的特征向量的有效方法之一。二.算法设计及流程图1、幕法算法(1) 取初始向量u(例如取u(0)=(1,1,1)T),置精度要求S置k=1.(2) 计算v(k)=Au(k),mk =max(v(k), u (k)= v (k)/ m k(3) 若| m k= mlvE，则停止计算(mk作为绝对值最大特征值九1 , u(k)作为 相应的特征向量)否则置k=k+1,转(2)2、反幕法算法(1) 取初始向量u(例如取u(0)=(1,1,1)T),置精

4、度要求s置k=1.(2) 对A作LU分解，即A=LU(3) 解线性方程组 Ly (k)=u(k),Uv(k)=y(k)(4) 计算mk=max(v(k), u (k)= v (k)/ m k(5) 若|mk=qjv e,则停止计算(1/mk作为绝对值最小特征值九n, u(k)作为 相应的特征向量)；否则置k=k+1，转(3).幕法流程图:反幕法流程图三、算法的理论依据及其推导(一) 幕法算法的理论依据及推导幕法是用来确定矩阵的主特征值的一种迭代方法，也即，绝对值最大的特 征值。稍微修改该方法，也可以用来确定其他特征值。幕法的一个很有用的特 性是它不仅可以生成特征值，而且可以生成相应的特征向量。

5、实际上，幕法经 常用来求通过其他方法确定的特征值的特征向量。1、幕法的迭代格式与收敛性质设n阶矩阵A的特征值，七，,n是按绝对值大小编号的，Xj(i=1,2,n)为对应-i的特征向量，且!为单根，即1 l|2| n |则计算最大特征值与特征向量的迭代格式为v(k) =Au(k d) ,m k =max(v(k), u (k) = v (k)/ m k(1)其中max(v(k)表示向量v(k)绝对值的最大分量。2、对于幕法的定理按式(1)计算出mk和u(k)满足lim mk=人，lim u(k) = Xk - ：kmax(x1)(二) 反幕法算法的理论依据及推导反幕法是用来计算绝对值最小的特征值

6、忽然相应的特征向量的方法。是对 幕法的修改，可以给出更快的收敛性。1、反幕法的迭代格式与收敛性质设A是非奇异矩阵，则零不是特征值，并设特征值为I |2| n4| n |则按A 4的特征值绝对值的大小排序，有1 1 1I II I - I In n4 1对A，实行幕法，就可得A，的绝对值最大的特征值1/ n和相应的特征向量， 即A的绝对值最小的特征值和相应的特征向量。由于用A代替A作幕法计算，因此该方法称为反幕法，反幕法的迭代格 式为v(k)= Au(z,mk =max(v(k), u (k)= v (k)/ m k(2)2、对于反幕法的定理按式(2)计算出的mk和u(k)满足：lim mk =

7、 1 ,lim u(k)= Xnk-却ZnkXmaX(Xn)在式(2)中，需要用到A，这给计算带来很大的不方便，因此，把(2)式 的第一式改为求解线性方程组A v(k)= u(k)(3)但由于在反幕法中，每一步迭代都需求解线性方程组(3)式，迭代做了大量的 重复计算，为了节省工作量，可事先把矩阵A作LU分解，即A=LU所以线性方程组(3)改为Ly (k) =u(k JL) ,Uv (k) =y(k)四、算法程序设计代码幕法程序，在matlab中建立一个M文件并保存。%pow.mfunction m,u,index,k=pow(A,u,ep,it_max)if nargin4it_max=100

8、0;endif nargin3ep=1e-5;endn=len gth(A);in dex=0;k=0;m仁0;m0=0;I=eye (n);T=A-m0*l;while k=it_maxv=T*u;vmax,i=max(abs(v);m=v(i);u=v/m;if abs(m-m1)ep;in dex=1;break;endm=m+m0;m1=m;k=k+1;end在matlab输入面板，输入A=rand ( 4)； %产生一个4维随机矩阵B=A+Au=1 1 1 1 ;%设立初始向量m,u,index,k=pow(B， u,ep,it_max)% 最多可省略 2 个参数程序结束。在M文件中

9、可以通过改变m0的值改变原点位移，从而达到原点位移加速。反幕法程序设计代码：在matlab中建立一个M文件并保存。%pow_i nv.mfun ctio n m,u,i ndex,k=pow_i nv(A,u,ep,it_max)if nargin4it_max=1000;endif nargin3ep=1e_5;endn=len gth(A);in dex=0;k=0;m仁0;m0=0;I=eye (n);T=A-m0*l;in vT=i nv(T);while k=it_maxv=in vT*u;vmax,i=max(abs(v);m=v(i);u=v/m;if abs(m-m1)B=ra

10、 nd(4);A=B+BA =0.26750.57760.63441.31300.57761.15030.76410.13670.63440.76410.02570.41931.3130 0.13670.41931.2248 u=1 1 1 1; m,u,i ndex,k=pow(A,u) m =2.6813u =0.85760.69340.56231.0000in dex =1k =49修改 M0=1e-3m =2.6814u =0.85760.69340.56231.0000in dex =0k =1001修改M0=0 %此时为幕法m =2.6815u =0.85760.69350.562

11、31.0000in dex =1k =10修改 U=1 2 3 4 修改 M0=1e-4m =2.6813u =0.85760.69340.56231.0000in dex =1k =9修改 M0=1e-3 m =2.6805u =0.85760.69340.56221.0000in dex =1k =7修改M0=0 m =2.6814u =0.85760.69340.56231.0000in dex =1k =9修改 U=3 5 6 7修改 M0=1e-4m =2.6819u =0.85770.69370.56241.0000in dex =1k =7修改 M0=1e-3m =2.6814u

12、 =0.85760.69340.56231.0000in dex =k =1001修改M0=0m =2.6820u =0.85770.69370.56241.0000in dex =1总结以上，幕法如下:Um0muin dexk0.00012.68130.85760.69340.56231.00001491 1 1 10.0012.68140.58760.69340.56231.00000100102.68150.85760.69350.56231.00001100.00012.68130.85760.69340.56231.0000191 2 3 40.0012.68050.85760.69

13、340.56221.00001702.68140.85760.69340.56231.0000190.00012.68190.85770.69370.56241.0000173 5 6 70.0012.69140.85760.69340.56231.00000100102.6920.85770.69370.56241.000017反幕法结果显示：在mO为0时 A=rand(4A =0, 81800. 34120,83850. 54U60, S6020,53410.56810.4490. 34200. 72710.3704. 69460. 28970-30930.702?0. 6213 B=A+

14、ATB =1,63591.00141.18050. 83631,00141.0S821,29520, 75421. 18051.29520.74081. 39730,83630.?542K39?31.2426 1 1 1 1?u =111I in, 11, index, kJ =pow_iiw (Bj le-4)Jk =0.3037U 二-0. 54751. 0000CL 1266-0. 5040index =1k =16M0=0.001 U=1 1 1 1B =1- 63591.00141. 18050*83331. 00141.06821.29520,57421.13051-20520.7

15、408h39730.83630.57421.3973h2426 u=l 1 1 1Fu =1111 叫 u, index, k =p aw_inv (B* le-4jm -山 3347u =-0. 89961. 00000. 2726T 2364index =1k =16M0=0.1 u=1 1 1 1B =L63E9tT00141- 18050. B3631.01141.06321.29520.57421. 13051.29520.7403L 39730. S3630.57421.3973L 2426 炉1111Fu =1111 itkj u, indeKj klpoW-invCBj le-

16、4)血=0. 384?u =-0. 89951.00000. 2726-0. 2364index =1k =15M0=0 u=1 3 5 7B =L63591.0014L 18050.8363L00141.0S82L 29520.5742L 18051.29520.7408L 39730.33630.5742L 3973L 2426 u=l 3 5 YVu -1357 % vi, indexj k=pov_inv(EJ u)ji -0.3&47u =-0. 8996LOOOQ0,2725-0. 2364index -1k =20M0=0.1 u=1 3 5 7B -L63591.00141,

17、18050.83631,00141.0682L25520.542L 18051.29520,740S1.39730,33630.5742L 35731.2426 il 3 5 71u =1357q indeKj k =pow_inv (Bf u)m -0.3847u =-0.89951.00000. 2726-0.2364inde -1k =17M0=0.5 u=1 3 5 7B 二1.63591.0014L1S050.8363E0014L0682L 29520.5742L18O51.2A520,74081.39730.83630, 5742L39731.2426 u=l 3 5 7Yu =1

18、357 m，u; index, kJpoinvfBj u)m -0.3S47u -0.89951. 00000.272S-0.2364irides =1k =27M0=0 u=2 3 4 5B -1.63591.OOM+ 1GO50.1.00141.06821 + 29520. bf421. 1EO51.29520,74061.39730.83630.57421,3S731.2426 3 2 3 4 5*u =2345 m, u, index, k =pow_inv (B, u)m =0.3847L =-0.3996LOOOO0.27260. 2364index =1k =ISM0=0.1 u

19、=2 3 4 5B =L63591.0014L 18050.E363R00141.0682L2S520.5742L18051.29520.74081.39730.83630.5742L 30731.2426 u=2 3 4 5*u =2345 jr, u, indexj k -pow_inv (Bj u)ji =0.38*17u 二-0,8995LOOOQ0,2725-0.2364index =1k =17M0=0.7 u=2 3 4 5B -1+635S1.0014L 180583631+ 0014LO6S21.29520+5742+ 18051.29520.74081.3973(L 肝備0.&742L 39731.2426 u=2 35u =234E m, u, index, k =pow_