北航数值分析1-Jacobi法计算矩阵特征值

数值分析2015/11/101准备工作算法设计矩阵特征值的求法有幂法、Jacobi法、QR法等，其中幂法可求得矩阵按模最大的特征值（反幂法可求得按模最小特征值），Jacobi法则可以求得对称阵的所有特征值。分析一：由题目中所给条件λ1≤λ2≤…≤λn，可得出λ1、λn按模并不一定严格小于或大于其他特征值，且即使按模严格小于或大于其他特征值，也极有可能出现|λs|λ1||λn|或|λs|λn||λ1|的情况，导致按幂法和反幂法无法求解λ1或λn二者中的一者；分析二：题目要求求解与数μk=λ1+k(λn-λ1)/40最接近的特征值λik(k=1,2,3…39)，这个问题其实可以转换为求A-μk按模最小的特征值的问题，但因为在第一个问题中无法确定能肯定的求得λ1和λn，所以第二个问题暂先搁浅；分析三：cond(A)2=||A||*||A-1||=|λ|max*|λ|min,这可以用幂法和反幂法求得，det(A)=λ1*λ2*…*λn，这需要求得矩阵A的所有特征值。由以上分析可知，用幂法和反幂法无法完成所有问题的求解，而用Jacobi法求得矩阵所有特征值后可以求解题目中所给的各个问题。所以该题可以用Jacobi法求解。模块设计由Jacobi法的经典4步骤，设计以下模块：矩阵初始化模块voidinit(double*m)迭代模块查找非对角线最大模voidfind_pq(double*m)计算旋转角度voidcalCosSin(double*m)更新矩阵voidrefreshMextrix(double*m)计算最终结果模块voidcalFinal(double*m)数据结构设计由于矩阵是对称阵，上下带宽均为2，所以可以考虑用二维数组压缩存储矩阵上半带或下半带。但由于Jacobi法在迭代过程中会破坏矩阵的形态，所以原来为零的元素可能会变为非零，这就导致原来的二维数组无法存储迭代后的矩阵。基于此的考虑，决定采用一维数组存储整个下三角阵，以此保证迭代的正确进行。完整代码如下（编译环境windows10+visualstudio2010）：数值分析2015/11/102完整代码//math.cpp:定义控制台应用程序的入口点。//#includestdafx.h#includestdio.h#includemath.h#includetime.h#defineN501#defineV(N+1)*N/2+1#definee2.718281828459045235360287471352662497757247093699959574966967627724076630353#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数值分析2015/11/103#ifdefPTSvoidPTS(double*m){printf(-----------------------------------------------------------------------\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,m[j]);}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*qq){for(intj=0;j=N-1;j++){if((i-(j+1)*j/2)=j+1){*pp=j+1;数值分析2015/11/104*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;i=N;i++){if(i!=p&&i!=q){if(ip){ipr=i;ipc=p;}else{ipr=p;ipc=i;}if(iq){iqr=i;iqc=q;}else{iqr=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;}}数值分析2015/11/105m[(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/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;i=V;i++){if(fabs(m[i])max){recounti(i,&pp,&qq);if(pp!=qq){max=fabs(m[i]);p=pp;q=qq;}}}MAX=max;}voidinit(double*m){for(inti=1;i=N;i++)m[(i+1)*i/2]=a(i);数值分析2015/11/106for(inti=2;i=N;i++)m[(i-1)*i/2+i-1]=b;for(inti=3;i=N;i++)m[(i-1)*i/2+i-2]=c;PTS(m);}voidcalFinal(double*m){printf(---------------------------------------------------------------------------------------------------\n);printf(结果输出:\n\n);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;i=N;i++){if(m[(i+1)*i/2]maxlumda)maxlumda=m[(i+1)*i/2];if(m[(i+1)*i/2]minlumda)minlumda=m[(i+1)*i/2];if(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(%d)=%.11eLumda(s)=%.11e\n,minlumda,N,maxlumda,absminlumda);printf(Cond(A)=%.11e\n,conda);printf(Det(A)=%.11e\n\n,deta);for(inti=1;i=FK;i++)数值分析2015/11/107{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(10);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(正在计算...\n\n);#endifwhile(true){kk++;find_pq(m);if(MAXeps)break;数值分析2015/11/108#ifdefPTSprintf(p=%dq=%d|max|=%e\n,p,q,MAX);printf(-----------------------------------------------------------------------\n\n);#endifcalCosSin(m);refreshMetrix(m);}#ifdefPTSprintf(p=%dq=%d|max|=%e\n,p,q,MAX);printf(-----------------------------------------------------------------------\n\n);#endifprintf(矩阵最终形态...\n);for(inti=1;i=PFrols;i++){for(intj=(i-1)*i/2+1;j=(i+1)*i/2;j++){printf(PFbits,m[j]);}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=time(&t1);#ifdefPTSprintf(计算并输出用时%.2f秒\n\n,difftime(t1,t0));#elseprintf(迭代次数%d,计算用时%.2f秒\n\n,kk,difftime(t1,t0));数值分析2015/11/