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

准备工作➢算法设计矩阵特征值的求法有幂法、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法求解。

➢模块设计由➢数据结构设计由于矩阵是对称阵，上下带宽均为2，所以可以考虑用二维数组压缩存储矩阵上半带或下半带。

但由于Jacobi法在迭代过程中会破坏矩阵的形态，所以原来为零的元素可能会变为非零，这就导致原来的二维数组无法存储迭代后的矩阵。

基于此的考虑，决定采用一维数组存储整个下三角阵，以此保证迭代的正确进行。

完整代码如下（编译环境windows10 + visual studio2010）：完整代码// math.cpp : 定义控制台应用程序的入口点。

//#include "stdafx.h"#include<stdio.h>#include<math.h>#include<time.h>#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,m[j]);}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 <= N ;i++){if(i != p && i != q){if(i > p){ipr = i;ipc = p;}else{ipr = p;ipc = i;}if(i > q){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];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 <= V ; i++){if(fabs(m[i]) > max){recounti(i,&pp,&qq);if(pp != qq){max = fabs(m[i]);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("结果输出:\n\n");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 <=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)=%.11e Lumda(%d)=%.11e Lumda(s)=%.11e\n",minlumda,N,maxlumda,absminlumda);printf(" Cond(A)=%.11e\n",conda);printf(" Det(A)=%.11e\n\n",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("正在计算...\n\n");#endifwhile(true){kk++;find_pq(m);if(MAX<eps)break;#ifdef PTSprintf(" p=%d q=%d |max|=%e\n",p,q,MAX);printf("-----------------------------------------------------------------------\n\n"); #endifcalCosSin(m);refreshMetrix(m);}#ifdef PTSprintf(" p=%d q=%d |max|=%e\n",p,q,MAX);printf("-----------------------------------------------------------------------\n\n");#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,m[j]);}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秒\n\n",difftime(t1,t0));#elseprintf("迭代次数%d,计算用时%.2f秒\n\n",kk,difftime(t1,t0)); #endifcalFinal(m);return 0;}运行结果如下：中间运行状态如下：结果分析➢有效性分析1.由输出结果可见矩阵经过21840次迭代后，非对角元全部为零或接近于零；2.代码中有定义预编译宏//#define PTS控制程序运行过程是否输出中间迭代结果，如果输出中间迭代结果，可以发现对角元素在迭代的后期变化非常小，达到收敛的效果；3.算法在多次运行中基本可以在45秒左右完成计算（酷睿i5双核处理器，10G内存，64位windows10操作系统）。