反幂法求矩阵特征值

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

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

3、对A作LU分解，即A=LU(3)解线性方程组Ly(k)=u(k?Uv(k)=y(k)(4)计算mk=max(v(k),u(k)=v(k)/mk(5)若|mk-mk|名，则停止计算(1/mk作为绝对值最小特征值n,u(k)作为相应的特征向量)；否则置k=k+1,转(3).三.程序框图1 .主程序2 .子程序(1) .哥法迭代程序框图(2).反哥法迭代程序框图四.结果显示计算结果如下：矩阵A的按模最大特征值为:-1.070011361487e+001矩阵A的按模最小特征值为：-5.557910794230e-003矩阵A最大的特征值为724634101479e+000矩阵A最小的特征值为:-1.0

4、70011361487e+001与各院(k=1,2,.,39)最接近的4k(用Vk表示)的值如下:v1=-1.018293403315e+001u1=-1.018949492196e+001v2=-9.585707425068e+000u2=-9.678876229054e+000v3=-9.172672423928e+000u3=-9.168257536145e+000v4=-8.652284007898e+000v5=-8.093483808675e+000v6=-7.659405407692e+000v7=-7.119684648691e+000v8=-6.611764339397e+0

5、00v9=-6.066103226595e+000v10=-5.585101052628e+000v11=-5.114083529812e+000v12=-4.578872176865e+000v13=-4.096470926260e+000v14=-3.554211215751e+000v15=-3.041090018133e+000v16=-2.533970311130e+000v17=-2.003230769563e+000v18=-1.503557611227e+000v19=-9.935586060075e-001v20=-4.870426738850e-001v21=2.23173

6、6249575e-002v22=5.324174742069e-001v23=1.052898962693e+000v24=1.589445881881e+000v25=2.060330460274e+000v26=2.558075597073e+000v27=3.080240509307e+000v28=3.613620867692e+000v29=4.091378510451e+000v30=4.603035378279e+000v31=5.132924283898e+000v32=5.594906348083e+000v33=6.080933857027e+000v34=6.680354

7、092112e+000v35=7.293877448127e+000v36=7.717111714236e+000v37=8.225220014050e+000v38=8.648666065193e+000v39=9.254200344575e+000u4=-8.657638843237e+000u5=-8.147020150328e+000u6=-7.636401457419e+000u7=-7.125782764510e+000u8=-6.615164071601e+000u9=-6.104545378693e+000u10=-5.593926685784e+000u11=-5.08330

8、7992875e+000u12=-4.572689299966e+000u13=-4.062070607058e+000u14=-3.551451914149e+000u15=-3.040833221240e+000u16=-2.530214528331e+000u17=-2.019595835422e+000u18=-1.508977142514e+000u19=-9.983584496049e-001u20=-4.877397566962e-001u21=2.287893621262e-002u22=5.334976291214e-001u23=1.044116322030e+000u24

9、=1.554735014939e+000u25=2.065353707848e+000u26=2.575972400756e+000u27=3.086591093665e+000u28=3.597209786574e+000u29=4.107828479483e+000u30=4.618447172392e+000u31=5.129065865300e+000u32=5.639684558209e+000u33=6.150303251118e+000u34=6.660921944027e+000u35=7.171540636935e+000u36=7.682159329844e+000u37=

10、8.192778022753e+000u38=8.703396715662e+000u39=9.214015408571e+000MlC;MndQ网式“巧temmcmd.exe1回求RETOu141=-3ulS1=-3ulGJ=-2U(17-2u18J-1u19J7iH丽Jf*090164545378693*0060833C7?728751?*096828892H66e+腼0237Q&S705e&400S551-J51?l-fJ.49e*000940B33221240e*00Q51021429?31e400Qa；l?595S35422t;*006508J771425f

11、4e*0009e3S449&B9e-09iB773975C6V629-B01119e4t4SC91e+3006117643373?7e*BB00B610322&5?5e+9O085101852i28e+900114D83529812e*080578B72i768&5e+00009647092G2GE呼+BO0554211215751ft*0届a41B9Bai8133e+80422970311120e+43O0a83230767563e*000卜：-1-0700119tl487&+001卜M791079429-牺724634101479?.71711171423&e+B0B537A2.22522

12、03100*008uL38J=.64866065193e+00日uL391-fJ2.54aaU44575e+H30请按任竟锹姓媛.五.程序#include#include#defineN501voidmain()doubleQ5501;doublemifa(doubleA5501);doublefanmifa(doubleA5501);doublelm,lmax,lmin,ls,delta,u39,v39;inti,j,k;doubleA5501;A00=A01=A10=A3500=A4499=A4500=0.0;/输入*501矩阵for(i=2;iN;i+)A0i=-0.064;for(i=

13、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+)A3i=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;最大特征值或最小特征值

14、if(lmaxlm)lmin=lmax;lmax=lm;elselmin=lm;for(i=0;iN;i+)/还原Afor(j=0;j5;j+)Aji=Qji;ls=fanmifa(A);/按模最小特征值for(i=0;iN;i+)/还原Afor(j=0;j5;j+)Aji=Qji;for(k=0;k39;k+)/计算u1-u39uk=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+)/还原Afor(j=0;j5;j+)Aji=Qji;printf(

15、矩阵的按模最大特征值为：.12e,lm);printf(n);printf(矩阵的按模最小特征值为：.12e,ls);printf(n);printf(矩阵最大的特征值为：.12e,lmax);printf(n);printf(矩阵最小的特征值为：.12e,lmin);printf(n);for(k=0;k0)return1;elseif(a=0)return0;elsereturn-1;intmax2(inta,intb)returnab?a:b;intmax3(inta,intb,intc)returnmax2(a,b)c?max2(a,b):c;)intmin(inta,intb)(re

16、turnab?a:b;)voidLU(doubleA5501,doubleu501,doubleB501)/LU分解法(doubleX501;inti,j,k,t,l;doublem=0,n=0;for(k=1;k=N;k+)/求L,U(for(j=k;j=min(N,k+2);j+)/U(m=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+)(

17、n+=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;doublemifa(doubleA5501)/嘉法(inti,j,l=0;doubleu501,t501;doubley501;doubleh,b,c;c=0;for(i=0;iN;i

18、+)/嘉法初始向量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*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;)returnb;)doublefanmifa(doubleA5501)/反嘉法doubleu501,y501;doubleP5501,Y501;/LU分解前用于保存A口y的值doublem=0,n=0,b=0,c=0;inti,j;for(i=0;iN;i+)/反募法初始向量u0=1;while(1)b=0;