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;