Jacobi迭代求矩阵特征值和特征向量+C代码

Jacobi计算过程如下：

1. 选择矩阵A非对角元中最大值A[i][j]，运用公式 tan 2O = 2*A[i][j] / (A[i][i] -A[j][j]) 获得选择平面矩阵J，使J * A *J'=A1

2. 计算得到A1后，重复第1歩，计算A2 =J2 *A1 *J2'，连续这个过程，直到非对角线元素全化为充分小，得到特征值矩阵Am。

3. 特征向量的计算为 D = J1' *J2'* ... *Jm' ，D 的列向量就是对应特征值矩阵 Am 对角线的特征向量。

下面是一个课本上的例子：

计算结果：（A的对角线元素为矩阵的特征值，Q为对应的特征向量(列对应)）

实现代码：（数据为书上的另一个例子）

#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "vector"
using namespace std;


#define N 5
#define E 0.000001   //非对角线数充分小
#define PI 3.141592657
#define MAXITER 10000  //最大迭代次数


typedef vector<double> dim1Vector;
typedef vector<dim1Vector> dim2Vector;
typedef vector<dim2Vector> dim3Vector;


bool QueryArray(dim2Vector Array);   //检查是否满足
dim2Vector matTran(dim2Vector Array);   //矩阵转置
dim2Vector matMul(dim2Vector mat1, dim2Vector mat2);   //矩阵相乘


void main()
{
int i, j;
int count;
bool flag = false;
double A[N][N] = {1267.9, -307.23, 0, 0, 0,  
-307.23, 710.24, -403.01, 0, 0,  
0, -403.01, 927.21, -524.2, 0,  
0, 0, -524.2, 770.07, -245.88,
0, 0, 0, -245.88, 245.88};
dim1Vector tempArray(N, 0);
dim2Vector Array;
dim2Vector charatMat(N, tempArray);
dim2Vector dim2Jac;
dim2Vector dim2JacT;
dim3Vector dim3Jac;
double maxArrayNum;
int laber_j, laber_i;
double theta;

for(i=0; i<N; i++)
{
tempArray.clear();
for(j=0; j<N; j++)
tempArray.push_back(A[i][j]);
Array.push_back(tempArray);
}


//开始迭代
count = 0;
tempArray.clear();
tempArray.resize(N, 0);
while(count<MAXITER && !flag)
{
count++;
dim2Jac.clear();
dim2Jac.resize(N, tempArray);
maxArrayNum = 0;
laber_i = laber_j = 0;


//寻找非对角元中绝对值最大的A[i][j]
for(i=0; i<N; i++)
for(j=0; j<N; j++)
{
if(i==j)
continue;


if(maxArrayNum<fabs(Array[i][j]))
{
maxArrayNum = fabs(Array[i][j]);
laber_i = i;
laber_j = j;
}
}


theta = atanf(Array[laber_i][laber_j]*2/(Array[laber_i][laber_i]-Array[laber_j][laber_j]+E));


//构造雅克比矩阵
for(i=0; i<N; i++)
dim2Jac[i][i] = 1;


dim2Jac[laber_i][laber_i] = dim2Jac[laber_j][laber_j] = cosf(theta/2);
dim2Jac[laber_i][laber_j] = sinf(theta/2);
dim2Jac[laber_j][laber_i] = -sinf(theta/2);


dim2JacT = matTran(dim2Jac);  //矩阵转置
dim3Jac.push_back(dim2JacT);  //保存矩阵


Array = matMul(matMul(dim2Jac, Array), dim2JacT);


if(QueryArray(Array))
flag = true;

}


//初始化特征矩阵
for(i=0; i<N; i++)
charatMat[i][i] = 1;


//计算特征矩阵
for(i=0; i<dim3Jac.size(); i++)
charatMat = matMul(charatMat, dim3Jac[i]);

printf("迭代次数：%d\n", count);


printf("\n特征值：\n");
for(i=0; i<N; i++)
{
for(j=0; j<N; j++)
printf("%lf  ", Array[i][j]);
printf("\n");
}


printf("\n特征矩阵：\n");
for(i=0; i<N; i++)
{
for(j=0; j<N; j++)
printf("%lf  ", charatMat[i][j]);
printf("\n");
}



}




//检查是否满足
bool QueryArray(dim2Vector Array)
{
int i, j;


for(i=0; i<N; i++)
for(j=0; j<N; j++)
{
if(i==j)
continue;


if(fabs(Array[i][j])>E)
return false;
}


return true;
}




//矩阵转置
dim2Vector matTran(dim2Vector Array)
{
int i, j;
dim1Vector temp(N, 0);
dim2Vector dst(N, temp);


for(i=0; i<N; i++)
for(j=0; j<N; j++)
dst[j][i] = Array[i][j];


return dst;
}


//矩阵相乘
dim2Vector matMul(dim2Vector mat1, dim2Vector mat2)
{
  int i, j, k;
  dim1Vector temp(N, 0);
  dim2Vector dst(N, temp);


  for(i=0; i<mat1.size(); i++)
  for(j=0; j<mat2[0].size(); j++)
  for(k=0; k<mat2.size(); k++)
  dst[i][j] += mat1[i][k]*mat2[k][j];


return dst;
}

运行结果：