幂法求矩阵特征值

1、问题描述用幂法与反幂法求解矩阵特征值求n阶方阵A的特征值和特征向量， 是实际计算中常常碰到的问题，如：机械、结构或电磁振动中的固有值问题等。对于n阶矩阵A,若存在数和n维向量x满足Ax=x(1)则称为矩阵A的特征值，x为相应的特征向量。由线性代数知识可知，特征值是代数方程(2)只需求解方程(2)(3)|l-A|= +玄1 + +an+an =0的根。从表面上看，矩阵特征值与特征向量的求解问题似乎很简单, 的根，就能得到特征值 ，再解齐次方程组( I-A ) x=0的解，就可得到相应的特征向量。上述方法对于n很小时是可以的。但当n稍大时，计算工作量将以惊人的速度增大， 并且由于计算带有误差，方程

2、(2)未必是精确的特征方程，自然就不必说求解方程 (2)与(3)的困难了。幕法与反幕法是一种计算矩阵主特征值及对应特征向量的迭代方法 特别是用于大型稀疏矩阵。这里用幕法与反幕法求解带状稀疏矩阵A501501的特征值。二.算法设计1 .幕法(1) 取初始向量u(例如取u(0)=(1,1,1)T ),置精度要求；，置k=1.(2) 计算v(k) =Au(k ), m k =max(v(k), u (k) = v (k) / m k(3) 若| m k-mk则停止计算(mk作为绝对值最大特征值, u (k)作为相应的特征 向量)否则置k=k+1，转(2)2.反幕法(1) 取初始向量u(0)(例如取u

3、 (0) =(1,1，1) T ),置精度要求；，置k=1.(2) 对A作LU分解，即A=LU(3) 解线性方程组Ly (k)=u(kT,Uv (k) =y(k)(4) 计算mk=max(v(k), u (k) = v (k)/ m k(5) 若|mk-mkJ| 111241=1.5547350l193?e+0B0 id 羽 J-2.S&3E3?0?818o*0B0 u(26 J2.575572-100756c t0B0 uI2?=3,08&S?1093665e+0B0 u!2S1-3.E9720?8&74+OB0 uI2? 1=4.10782 87? *483e *060 u(30-4.E1

4、94471723?2e+060 u31I-E.139fl6ES6&30afl*060 ul32J-&.63?8455S?0?e*OB0 uI33J-6.1503032SL118e+000 uE341,t9219402?0*0S0 uC35J-7.17LS4063&935e+00 u36J-7_68215?32?e44e+0G0 ul2?J8.192778022?&3o*0G0 uf3818.70i3?6715662e+O00 u391=9,2141540&?le 也陌五.程序#include #include #define N 501void main()double Q5501;doubl

5、e mifa( double A5501);double fanmifa( double A5501);double lm,lmax,lmin,ls,delta,u39,v39;int i,j,k;double A5501;A00=A01=A10=A3500=A4499=A4500=0.0; / 输入 *501 矩阵 for (i=2;iN;i+)A0i=-0.064;for (i=1;iN;i+) A1i=0.16;for (i=0;iN;i+) A2i=(1.64-0.024*(i+1)*sin(0.2*(i+1)-0.64*exp(0.1/(i+1);for (i=0;i500;i+)A

6、3i=0.16;for (i=0;i499;i+)A4i=-0.064;for (i=0;i5;i+) / 保存 Afor (j=0;j501;j+) Qij=Aij;lm=mifa(A); /按模最大特征值,函数mifa()不会改变矩阵A的值，不需还原for (i=0;iN;i+)/ 平移 AA2i=A2i-lm;lmax=mifa(A); /平移后A的按模最大特征值lmax=lmax+lm; / 最大特征值或最小特征值if (lmaxlm)lmin=lmax;lmax=lm;elselmin=lm;for (i=O;iN;i+)/ 还原 Afor (j=0;j5;j+)Aji=Qji;ls

7、=fanmifa(A); / 按模最小特征值for (i=0;iN;i+)/ 还原 Afor (j=O;j5;j+)Aji=Qji;for (k=0;k39;k+) / 计算 u1-u39 uk=lmin+(k+1)*(lmax-lmin)/40);for (k=0;k39;k+)for (j=0;jN;j+)A2j=A2j-uk; vk=fanmifa(A)+uk; for (i=0;iN;i+) / 还原 A for (j=0;j5;j+) Aji=Qji;printf( 矩阵的按模最大特征值为:%.12e ,lm);printf(n );printf( 矩阵的按模最小特征值为:%.12e

8、 ,ls);printf(n );printf( 矩阵最大的特征值为 :%.12e ,lmax);printf(n );printf( 矩阵最小的特征值为 :%.12e ,lmin);printf(n );for (k=0;k0)return 1;else if (a=0)return 0;else return -1;int max2( int a, int b)return ab?a:b;int max3( int a, int b, int c)return max2(a,b)c?max2(a,b):c;int min( int a, int b)return ab?a:b;/LU 分解

9、法void LU( double A5501, double u501, double B501) double X501;int i,j,k,t,l;double m=0,n=0;for (k=1;k=N;k+) / 求L, Ufor (j=k;j=min(N,k+2);j+)/Um=0;for (t=max3(1,k-2,j-2);t=k-1;t+) m+=Ak-t+2t-1*At-j+2j-1;Ak-j+2j-1=Ak-j+2j-1-m;for (i=k+1;i=min(N,k+2);i+)/Lif (kN)n=0;for (l=max3(1,i-2,k-2);l=k-1;l+) n+=

10、Ai-l+2l-1*Al-k+2k-1; Ai-k+2k-1=(Ai-k+2k-1-n)/A2k-1;for (i=2;i=N;i+)/ 回代过程m=0;for (t=max2(1,i-2);t=1;i-)n=0;for (t=i+1;t=min(N,i+2);t+)n+=Ai-t+2t-1*Xt-1;Xi-1=(Bi-1-n)/A2i-1;for (i=1;i=N;i+) / 输出方程结果ui-1=Xi-1;double mifa( double A5501) / 幂法int i,j,l=0;double u501,t501;double y501;double h,b,c;c=0;for

11、(i=0;iN;i+) / 幂法初始向量ui=1;while (1)for (i=0;iN;i+)ti=0;h=u0;for (i=0;iN;i+) / 无穷范数if (fabs(h)fabs(ui)h=ui;l=i;for (i=0;iN;i+) yi=ui/fabs(h);for (i=2;i499;i+)for (j=i-2;j=i+2;j+)ti=ti+Ai-j+2j*yj;ui=ti; u0=A20*y0+A11*y1+A02*y2; u1=A30*y0+A21*y1+A12*y2+A03*y3; u499=A4497*y497+A3498*y498+A2499*y499+A1N-1

12、*yN-1; uN-1=A4498*y498+A3499*y499+A2N-1*yN-1;b=sgn(h)*ul;if (fabs(b-c)/fabs(b)=1e-12)/printf( 幂法成功！ );/printf(n);break ;c=b;return b;double fanmifa( double A5501) / 反幂法double u501,y501;double P5501,Y501; /LU分解前用于保存A和y的值double m=0,n=0,b=0,c=0;int i,j;for (i=0;iN;i+)/ 反幂法初始向量u0=1;while (1)b=0;n=0;for (i=0;iN;i+)n=n+ui*ui;n=sqrt(n);for (i=0;iN;i+)yi=ui/n;for (i=0;iN;i+)/ 保存 A和yYi=yi;for (j=0;j5;j+)Pji=Aji;LU(A