《北航数值分析1-Jacobi法计算矩阵特征值》由会员分享，可在线阅读，更多相关《北航数值分析1-Jacobi法计算矩阵特征值（11页珍藏版）》请在金锄头文库上搜索。
1、如果您需要使用本文档，请点击下载按钮下载!准备工作 算法设计矩阵特征值的求法有幂法、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 calC
3、osSin(double *m)更新矩阵void refreshMextrix(double *m)计算最终结果模块void calFinal(double *m) 数据结构设计由于矩阵是对称阵，上下带宽均为2，所以可以考虑用二维数组压缩存储矩阵上半带或下半带。但由于Jacobi法在迭代过程中会破坏矩阵的形态，所以原来为零的元素可能会变为非零，这就导致原来的二维数组无法存储迭代后的矩阵。基于此的考虑，决定采用一维数组存储整个下三角阵，以此保证迭代的正确进行。完整代码如下（编译环境windows10 + visual studio2010）：如果您需要使用本文档，请点击下载按钮下载!完整代码/
4、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
5、 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
6、; 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+
7、)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
8、* 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
9、(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
10、,&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);doub
11、le 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
12、);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 - minlum
