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

数值分析2015/11/10准备工作 算法设计矩阵特征值的求法有幂法、Jacobi法、QR法等，其中幂法可求得矩阵按模最大的特征值（反幂法可求得按模最小特征值），Jacobi法则可以求得对称阵的所有特征值。分析一：由题目中所给条件12n，可得出1、n按模并不一定严格小于或大于其他特征值，且即使按模严格小于或大于其他特征值，也极有可能出现|s|1|n |或|s|n|1 |的情况，导致按幂法和反幂法无法求解1或n二者中的一者；分析二：题目要求求解与数k =1+k(n-1)/40最接近的特征值ik(k=1,2,339)，这个问题其实可以转换为求A-k按模最小的特征值的问题，但因为在第一个问题中无法确定能肯定的求得1和n，所以第二个问题暂先搁浅；分析三：cond(A) 2 = |A| * |A-1| =|max * |min,这可以用幂法和反幂法求得，det(A) =1 *2 * *n，这需要求得矩阵A的所有特征值。由以上分析可知，用幂法和反幂法无法完成所有问题的求解，而用Jacobi法求得矩阵所有特征值后可以求解题目中所给的各个问题。所以该题可以用Jacobi法求解。 模块设计由Jacobi法的经典4步骤，设计以下模块：矩阵初始化模块void init(double *m)迭代模块查找非对角线最大模void find_pq(double *m)计算旋转角度void calCosSin(double *m)更新矩阵void refreshMextrix(double *m)计算最终结果模块void calFinal(double *m) 数据结构设计由于矩阵是对称阵，上下带宽均为2，所以可以考虑用二维数组压缩存储矩阵上半带或下半带。但由于Jacobi法在迭代过程中会破坏矩阵的形态，所以原来为零的元素可能会变为非零，这就导致原来的二维数组无法存储迭代后的矩阵。基于此的考虑，决定采用一维数组存储整个下三角阵，以此保证迭代的正确进行。完整代码如下（编译环境windows10 + visual studio2010）：完整代码/ math.cpp : 定义控制台应用程序的入口点。/#include stdafx.h#include#include#include#define N 501#define V (N+1)*N/2+1#define e 2.718281828459045235360287471352662497757247093699959574966967627724076630353#define a(i) (1.64 - 0.024 * (i) * sin(0.2 * (i) - 0.64 * pow(e , 0.1 / (i)#define b 0.16#define c -0.064#define eps pow(double)10.0,-12)#define PFbits %10.5f #define PFrols 5#define PFe %.11e#define FK 39int p;int q;double cosz;double sinz;double MAX;int kk;/#define PTS pts#ifdef PTSvoid PTS(double *m)printf(-n);printf(迭代第%d次n,kk);for(int i = 1 ; i = PFrols ; i+)for( int j = (i-1)*i/2+1 ; j = (i+1)*i/2 ; j+)printf(PFbits,mj);putchar(10);for(int i = 1 ; i = PFrols+1 ; i+)printf( . );putchar(10);printf( .n);printf( .n);printf( .n);for(int i = 1 ; i = PFrols+2 ; i+)printf( . );putchar(10);#elsevoid PTS(double *m)#endifvoid recounti(int i , int *pp, int *qq)for(int j = 0 ; j = N-1 ; j+)if( (i - (j+1)*j/2) = j+1)*pp = j+1;*qq = i - (j+1)*j/2;break;void refreshMetrix(double *m)int ipr,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(int i = 1; i p)ipr = i;ipc = p;elseipr = p;ipc = i;if(i q)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);/void calCosSin(double *m)double app = m(p+1)*p/2;double aqq = m(q+1)*q/2;double apq = m(p-1)*p/2+q;cosz = cos(atan(2 * apq / (app - aqq)/2);sinz = sin(atan(2 * apq / (app - aqq)/2);/void find_pq(double *m)double max = 0.0;int pp = 0;int qq = 0;for(int i = 1 ; i max)recounti(i,&pp,&qq);if(pp != qq)max = fabs(mi);p = pp;q = qq;MAX = max;void init(double *m)for(int i = 1 ; i = N ;i+)m(i+1)*i/2 = a(i);for(int i = 2 ; i = N ; i+)m(i-1)*i/2+i-1 = b;for(int i = 3 ; i = N ; i+)m(i-1)*i/2+i-2 = c;PTS(m);void calFinal(double *m)printf(-n);printf(结果输出:nn);double conda;double deta = 1.0;double minlumda = pow(double)10.0,12);double maxlumda = pow(double)10.0,-12);double absminlumda = pow(double)10.0,12);for(int i = 1 ; i 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)=%.11e Lumda(%d)=%.11e Lumda(s)=%.11en,minlumda,N,maxlumda,absminlumda);printf( Cond(A)=%.11en,conda);printf( Det(A)=%.11enn,deta);for(int i = 1 ; i = FK ; i+)double muk = minlumda + i * (maxlumda - minlumda) / 40;double lumdak = 0.0;double tempabsmin = pow(double)10.0,12);for(int j = 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(int argc, _TCHAR* argv)double m(N+1)*N/2+1 = 0.0;kk=0;MAX=1.0;time_t t0,t1;t0 = time( &t0);init(m);#ifndef PTSprintf(正在计算.nn);#endifwhile(true)kk+;find_pq(m);if(MAXeps)break;#ifdef PTSprintf( p=%d q=%d |max|=%en,p,q,MAX);printf(-nn);#endifcalCosSin(m);refreshMetrix(m);#ifdef PTSprintf( p=%d q=%d |max|=%en,p,q,MAX);printf(-nn);#endifprintf(矩阵最终形态.n);for(int i = 1 ; i = PFrols ; i+)for( int j = (i-1)*i/2+1 ; j = (i+1)*i/2 ; j+)printf(PFbits,mj);putchar(10);for(int i = 1 ; i = PFrols+1 ; i+)printf( . );putchar(10);printf( .n);printf( .n);printf( .n);for(int i = 1 ; i = PFrols+2 ; i+)printf( . );putchar(10);t1 = time(&t1);#ifdef PTSprintf(计算并输出用时%.2f秒nn,difftime(t1,t0);#elseprintf(迭代次数%d,计算用时%.2f秒nn,kk,difftime(t1,t0);#endifcalFinal(m);return 0;运行结果如下：中间运行状态如下：结果分析 有效性分析1. 由输出结果可见矩阵经过21840次迭代后，非对角元全部为零或接近于零；2. 代码中有定义预编译宏/#define PTS控制程序运行过程是否输出中间迭代结果，如果输出中间迭代结果，可以发现对角元素在迭代的后期变化非常小，达到收敛的效果；3. 算法在多次运行中基本可以在45秒左右完成计算（酷睿i5双核处理器，10G内存，64位windows10操作系统）。由此三条可得出结论：算法是有效的。 时空效率分析1. 算法的空间复杂度为o(N2)，与矩阵的维数相关，当矩阵维数较高时，算法可能在空间的开销上比较大，目前没有想到更好的数据结构来节省空间开销。2. 算法的时间复杂度主要取决于矩阵的收敛速度，