北航数值分析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

#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 39

int p;

int q;

double cosz;

double sinz;

double MAX;

int kk;

//#define PTS pts

#ifdef PTS

void 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);

}

#else

void PTS(double *m)

{

}

#endif

void 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;

}