- 注册时间
- 2008-6-8
- 最后登录
- 1970-1-1
|
雅可比法求实对称矩阵所有特征值和对应特征向量
程序如下
#define EPSL 1.0e-6 //迭代的精度水平
#define PI 3.1415926
#include <stdio.h>
#include <math.h>
void main() //主程序开始;
{
int size,//协方差矩阵的阶数
row,col,k,
p,q;//记录maxa的行与列
double a[99][99],//存放协方差矩阵
R[99][99]={0},//平面旋转变换矩阵
RT[99][99],//矩阵R的转置
S[99][99],//存放特征向量,赋初值为单位矩阵
tan,//旋转角度的正切
alfa,//旋转角度
maxa;//非主对角线上绝对值最大元素
FILE *fp;
fp=fopen("sj.txt","r");
fscanf(fp,"%d",&size);//读入协方差矩阵的阶数
for(row=0;row<size;row++)
for(col=0;col<size;col++)
fscanf(fp,"%lf",&a[row][col]);
fclose(fp);
for(row=0;row<size;row++)//给矩阵S赋初值为单位矩阵(经验证无误)
for(col=0;col<size;col++)
if(row==col)
S[row][col]=1;
else
S[row][col]=0;
while(fabs(maxa)>EPSL)
{
maxa=a[0][1];
for(row=0;row<size;row++)
for(col=row+1;col<size;col++)
if(fabs(a[row][col])>maxa)
{
maxa=a[row][col];//找出非主对角线上绝对值最大的元素
p=row;q=col;//记录maxa的行与列数
}
if(a[p][p]=a[q][q])
tan=PI/2.0;
else
tan=2*a[p][q]/(a[p][p]-a[q][q]);
alfa=atan(tan);//求出旋转角度
R[p][p]=cos(alfa);
R[q][q]=cos(alfa);
R[p][q]=-sin(alfa);
R[q][p]=sin(alfa);//形成矩阵R
for(row=0;row<size;row++)
for(col=0;col<size;col++)
RT[col][row]=R[row][col];//形成R矩阵的转置矩阵RT
for(row=0;row<size;row++)
for(col=0;col<size;col++)
for(k=0;k<size;k++)
a[row][col]+=RT[row][k]*a[k][col]; //计算RT与a的乘积,并赋给a
for(row=0;row<size;row++)
for(col=0;col<size;col++)
for(k=0;k<size;k++)
a[row][col]+=a[row][k]*R[k][col]; //计算a与R的乘积,并赋给a
for(row=0;row<size;row++)
for(col=0;col<size;col++)
for(k=0;k<size;k++)
S[row][col]+=S[row][k]*R[k][col]; //计算S与R的乘积,并赋给S
}
fp=fopen("sc1.txt","w");
for(row=0;row<size;row++)
{
fprintf(fp,"\n");
for(col=0;col<size;col++)
fprintf(fp,"%10.4f",a[row][col]);
}
fclose(fp);
fp=fopen("sc2.txt","w");
for(row=0;row<size;row++)
{
fprintf(fp,"\n");
for(col=0;col<size;col++)
fprintf(fp,"%10.4f",S[row][col]);
}
fclose(fp);
}
/*for(row=0;row<size;row++)//验证用语句
{
printf("\n");
for(col=0;col<size;col++)
printf("%10.2f",S[row][col]);
}
printf("\n");*/
[ 本帖最后由 bobnie 于 2008-6-21 21:55 编辑 ] |
|