北航數(shù)值分析1Jacobi法計(jì)算矩陣特征值

1、準(zhǔn)備工作算法設(shè)計(jì)矩陣特點(diǎn)值的求法有冪法、Jacobi法、QR法等，此中冪法可求得矩陣按模最大的特點(diǎn)值（反冪法可求得按模最小特點(diǎn)值），Jacobi法例能夠求得對稱陣的所有特點(diǎn)值。剖析一：由題目中所給條件12n，可得出1、n按模其實(shí)不必定嚴(yán)格小于或大于其余特征值，且即便按模嚴(yán)格小于或大于其余特點(diǎn)值，也極有可能出現(xiàn)|s|1|n|或|s|n|1|的狀況，致使按冪法和反冪法沒法求解1或n兩者中的一者；剖析二：題目要求求解與數(shù)k=1+k(n-1)/40最靠近的特點(diǎn)值ik(k=1,2,339)，這個問題其實(shí)能夠變換為求A-k按模最小的特點(diǎn)值的問題，但因?yàn)樵诘谝粋€問題中沒法確立能必定的求得1和n，因此第二個問

2、題暫先擱淺；剖析三：cond(A)2=|A|*|A-1|=|max*|min,這能夠用冪法和反冪法求得，det(A)=1*2*n，這需要求得矩陣A的所有特點(diǎn)值。由以上剖析可知，用冪法和反冪法沒法達(dá)成所有問題的求解，而用Jacobi法求得矩陣所有特點(diǎn)值后可以求解題目中所給的各個問題。因此該題能夠用Jacobi法求解。模塊設(shè)計(jì)由Jacobi法的經(jīng)典4步驟，設(shè)計(jì)以下模塊：矩陣初始化模塊voidinit(double*m)查找非對角線最大模voidfind_pq(double*m)迭代模塊計(jì)算旋轉(zhuǎn)角度voidcalCosSin(double*m)更新矩陣voidrefreshMextrix(doubl

3、e*m)計(jì)算最后結(jié)果模塊voidcalFinal(double*m)?數(shù)據(jù)構(gòu)造設(shè)計(jì)因?yàn)榫仃囀菍ΨQ陣，上下帶寬均為2，因此能夠考慮用二維數(shù)組壓縮儲存矩陣上半帶或下半帶。但因?yàn)镴acobi法在迭代過程中會損壞矩陣的形態(tài)，因此本來為零的元素可能會變成非零，這就致使本來的二維數(shù)組沒法儲存迭代后的矩陣。鑒于此的考慮，決定采納一維數(shù)組儲存整個下三角陣，以此保證迭代的正確進(jìn)行。完好代碼以下（編譯環(huán)境windows10+visualstudio2010）：完好代碼math.cpp:定義控制臺應(yīng)用程序的進(jìn)口點(diǎn)。/#includestdafx.h#include#include#include#defineN50

4、1#defineV(N+1)*N/2+1#definee2.6630353#definea(i)(1.64-0.024*(i)*sin(0.2*(i)-0.64*pow(e,0.1/(i)#defineb0.16#definec-0.064#defineepspow(double)10.0,-12)#definePFbits%10.5f#definePFrols5#definePFe%.11e#defineFK39intp;intq;doublecosz;doublesinz;doubleMAX;intkk;/#definePTSpts#ifdefPTSvoidPTS(double*m)pri

5、ntf(-n);printf(迭代第%d次n,kk);for(inti=1;i=PFrols;i+)for(intj=(i-1)*i/2+1;j=(i+1)*i/2;j+)printf(PFbits,mj);putchar(10);for(inti=1;i=PFrols+1;i+)printf(.);putchar(10);printf(.n);printf(.n);printf(.n);for(inti=1;i=PFrols+2;i+)printf(.);putchar(10);#elsevoidPTS(double*m)#endifvoidrecounti(inti,int*pp,int*

6、qq)for(intj=0;j=N-1;j+)if(i-(j+1)*j/2)=j+1)*pp=j+1;*qq=i-(j+1)*j/2;break;voidrefreshMetrix(double*m)intipr,ipc,iqr,iqc;m(p+1)*p/2=m(p+1)*p/2*pow(cosz,2)+m(q+1)*q/2*pow(sinz,2)+2*m(p-1)*p/2+q*cosz*sinz;m(q+1)*q/2=m(p+1)*p/2*pow(sinz,2)+m(q+1)*q/2*pow(cosz,2)-2*m(p-1)*p/2+q*cosz*sinz;for(inti=1;ip)ipr

7、=i;ipc=p;elseipr=p;ipc=i;if(iq)iqr=i;iqc=q;elseiqr=q;iqc=i;m(ipr-1)*ipr/2+ipc=m(ipr-1)*ipr/2+ipc*cosz+m(iqr-1)*iqr/2+iqc*sinz;m(iqr-1)*iqr/2+iqc=-m(ipr-1)*ipr/2+ipc*sinz+m(iqr-1)*iqr/2+iqc*cosz;m(p-1)*p/2+q=0;PTS(m);/voidcalCosSin(double*m)doubleapp=m(p+1)*p/2;doubleaqq=m(q+1)*q/2;doubleapq=m(p-1)*p

8、/2+q;cosz=cos(atan(2*apq/(app-aqq)/2);sinz=sin(atan(2*apq/(app-aqq)/2);/voidfind_pq(double*m)doublemax=0.0;intpp=0;intqq=0;for(inti=1;imax)recounti(i,&pp,&qq);if(pp!=qq)max=fabs(mi);p=pp;q=qq;MAX=max;voidinit(double*m)for(inti=1;i=N;i+)m(i+1)*i/2=a(i);for(inti=2;i=N;i+)m(i-1)*i/2+i-1=b;for(inti=3;i=

9、N;i+)m(i-1)*i/2+i-2=c;PTS(m);voidcalFinal(double*m)printf(-n);printf(結(jié)果輸出:nn);doubleconda;doubledeta=1.0;doubleminlumda=pow(double)10.0,12);doublemaxlumda=pow(double)10.0,-12);doubleabsminlumda=pow(double)10.0,12);for(inti=1;imaxlumda)maxlumda=m(i+1)*i/2;if(m(i+1)*i/2minlumda)minlumda=m(i+1)*i/2;if(

10、fabs(m(i+1)*i/2)absminlumda)absminlumda=fabs(m(i+1)*i/2);deta*=m(i+1)*i/2;if(fabs(minlumda)fabs(maxlumda)conda=fabs(maxlumda)/absminlumda;elseconda=fabs(minlumda)/absminlumda;printf(Lumda(1)=%.11eLumda(s)=%.11en,minlumda,N,maxlumda,absminlumda);printf(Cond(A)=%.11en,conda);printf(Det(A)=%.11enn,deta

11、);Lumda(%d)=%.11efor(inti=1;i=FK;i+)doublemuk=minlumda+i*(maxlumda-minlumda)/40;doublelumdak=0.0;doubletempabsmin=pow(double)10.0,12);for(intj=1;j=N;j+)if(fabs(muk-m(j+1)*j/2)tempabsmin)lumdak=m(j+1)*j/2;tempabsmin=fabs(muk-m(j+1)*j/2);printf(Lumda(i%d)=%.11e,i,lumdak);if(i%3=0)putchar(10);putchar(1

12、0);printf(-n);putchar(10);putchar(10);int_tmain(intargc,_TCHAR*argv)doublem(N+1)*N/2+1=0.0;kk=0;MAX=1.0;time_tt0,t1;t0=time(&t0);init(m);#ifndefPTSprintf(正在計(jì)算.nn);#endifwhile(true)kk+;find_pq(m);if(MAXeps)break;#ifdefPTSprintf(p=%dq=%d|max|=%en,p,q,MAX);printf(-nn);#endifcalCosSin(m);refreshMetrix(m

13、);#ifdefPTSprintf(p=%dq=%d|max|=%en,p,q,MAX);printf(-nn);#endifprintf(矩陣最后形態(tài).n);for(inti=1;i=PFrols;i+)for(intj=(i-1)*i/2+1;j=(i+1)*i/2;j+)printf(PFbits,mj);putchar(10);for(inti=1;i=PFrols+1;i+)printf(.);putchar(10);printf(.n);printf(.n);printf(.n);for(inti=1;i=PFrols+2;i+)printf(.);putchar(10);t1=t

14、ime(&t1);#ifdefPTSprintf(計(jì)算并輸出用時%.2f秒nn,difftime(t1,t0);#elseprintf(迭代次數(shù)%d,計(jì)算用時%.2f秒nn,kk,difftime(t1,t0);#endifcalFinal(m);return0;運(yùn)轉(zhuǎn)結(jié)果以下：中間運(yùn)轉(zhuǎn)狀態(tài)以下：結(jié)果剖析有效性剖析1.由輸出結(jié)果可見矩陣經(jīng)過21840次迭代后，非對角元所有為零或靠近于零；2.代碼中有定義預(yù)編譯宏/#definePTS控制程序運(yùn)轉(zhuǎn)過程能否輸出中間迭代結(jié)果，假如輸出中間迭代結(jié)果，能夠發(fā)現(xiàn)對角元素在迭代的后期變化特別小，達(dá)到收斂的成效；算法在多次運(yùn)轉(zhuǎn)中基本能夠在45秒左右達(dá)成計(jì)算（酷睿i5雙核辦理器，10G存，64位windows10操作系統(tǒng)）。由此三條可得出結(jié)論：算法是有效的。時空效率剖析算法的空間復(fù)雜度為o(N2)，與矩陣的維數(shù)有關(guān)，當(dāng)矩陣維數(shù)較高時，算法可能在空間的開支上比較大，當(dāng)前沒有想到更好的數(shù)據(jù)構(gòu)造來節(jié)儉空間開支。2.算法的時間復(fù)雜度主要取決于矩陣的收斂速度，而每次迭代的時間復(fù)雜度為o(N2)，主要耗資在順