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

1、北航数值分析1Jacobi法计算矩阵特征值准备工作 算法设计矩阵特征值的求法有幂法、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 按模最小的特征值的问题，但因为在第

2、一个问题中无法确定能肯定的求得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 cal

3、CosSin(double *m)更新矩阵void refreshMextrix(double *m)计算最终结果模块void calFinal(double *m) 数据结构设计由于矩阵是对称阵，上下带宽均为2，所以可以考虑用二维数组压缩存储矩阵上半带或下半带。但由于Jacobi法在迭代过程中会破坏矩阵的形态，所以原来为零的元素可能会变为非零，这就导致原来的二维数组无法存储迭代后的矩阵。基于此的考虑，决定采用一维数组存储整个下三角阵，以此保证迭代的正确进行。完整代码如下（编译环境windows10 + visual studio2010）：完整代码/ math.cpp : 定义控制台应用程序

4、的入口点。/#include stdafx.h#include#include#include#define N 501#define V (N+1)*N/2+1#define e 2.6630353#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

5、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( . ); pu

6、tchar(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

7、(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; else ipr = p; ipc = i; if(i q)

8、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; 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; doub

9、le 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

10、 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.

11、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) / absminlum

12、da; else conda = 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.

13、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, _T

14、CHAR* 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 PTS printf(正在计算.nn);#endif while(true) kk+; find_pq(m); if(MAXeps)break;#ifdef PTS printf( p=%d q=%d |max|=%en,p,q,MAX); printf(-nn);#endif calCosSin(m); refreshMetrix(m); #ifdef PTS printf( p=%d q=%

15、d |max|=%en,p,q,MAX); printf(-nn);#endif printf(矩阵最终形态.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 = PFr

16、ols+2 ; i+) printf( . ); putchar(10); t1 = time(&t1);#ifdef PTS printf(计算并输出用时%.2f秒nn,difftime(t1,t0);#else printf(迭代次数%d,计算用时%.2f秒nn,kk,difftime(t1,t0);#endif calFinal(m); return 0;运行结果如下：中间运行状态如下：结果分析 有效性分析1. 由输出结果可见矩阵经过21840次迭代后，非对角元全部为零或接近于零；2. 代码中有定义预编译宏/#define PTS控制程序运行过程是否输出中间迭代结果，如果输出中间迭代结果

17、，可以发现对角元素在迭代的后期变化非常小，达到收敛的效果；3. 算法在多次运行中基本可以在45秒左右完成计算（酷睿i5双核处理器，10G存，64位windows10操作系统）。由此三条可得出结论：算法是有效的。 时空效率分析1. 算法的空间复杂度为o(N2)，与矩阵的维数相关，当矩阵维数较高时，算法可能在空间的开销上比较大，目前没有想到更好的数据结构来节省空间开销。2. 算法的时间复杂度主要取决于矩阵的收敛速度，而每次迭代的时间复杂度为o(N2)，主要耗费在顺序查找非对角线最大元和数组下标分解为矩阵坐标这两个子算法上。目前可以考虑在查找非对角线最大元的算法上做一些改进：在每次查找的过程中记录下上次查找到的第二大元素（第一大元素会被清零），然后在矩阵刷新之后的查找过程中，用上次查找到的第二大元素和刷新过的元素做比较得出本次迭代的最大元和第二大元。 数值角度分析1. 计算过程中，涉及到矩阵变元的量均采用double型变量存储，满足精度要求；2. 由计算结果可见矩阵的条件数较大，可能会由于舍入误差的影响对结果造成较大扭曲。