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

Jacobi计算过程如下：

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

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

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

#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;
}

©️2019 CSDN 皮肤主题: 大白 设计师: CSDN官方博客