准备工作
? 算法设计
矩阵特征值的求法有幂法、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法求解。
? 模块设计
由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）：
编辑版word
完整代码
// math.cpp : 定义控制台应用程序的入口点。 //
#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 \ #define PFrols 5 #define PFe %.11e
#define FK 39
int p; int q;
double cosz; double sinz;
double MAX; int kk;
//#define PTS pts
编辑版word
#ifdef PTS
void PTS(double *m) {
printf(\ printf(\ 迭代第%d次\\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(\ printf(\ printf(\
. .\\n\ . .\\n\ . .\\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) {
编辑版word
北航数值分析1-Jacobi法计算矩阵特征值
