#include<stdio.h>
#include<math.h>
int main()
{
double Rb[1000][10][10]={0},R[10][10]={0},A[1000][10][10],max,tanth,costh,sinth,sum=0,pfh,c,d;
int p,q,i,j,k=1,m,n;
for(i=0;i<10;i++)
{
for(j=0;j<10;j++)
{
if(i==j) Rb[0][i][j]=1;
else Rb[0][i][j]=0;
}
}
printf("请输入矩阵A:\n");
for(i=0;i<10;i++)
{
scanf("%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",&A[0][i][0],&A[0][i][1],&A[0][i][2],&A[0][i][3],&A[0][i][4],&A[0][i][5],&A[0][i][6],&A[0][i][7],&A[0][i][8],&A[0][i][9]);
}
while(pfh>0.001)
{
for(m=0,max=0;m<9;m++)
{
for(n=m+1;n<10;n++)
{
if(fabs(A[k-1][m][n])>max)
{
max=fabs(A[k-1][m][n]);
p=m,q=n;
}
}
}
if(A[k-1][p][p]-A[k-1][q][q]==0)
{
if(A[k-1][p][q]>0) costh=0.5*sqrt(2),sinth=0.5*sqrt(2);
else if(A[k-1][p][q]<0) costh=0.5*sqrt(2),sinth=(-1)*0.5*sqrt(2);
}
else if(A[k-1][p][p]-A[k-1][q][q]!=0)
{
c=(A[k-1][p][p]-A[k-1][q][q])/(2*A[k-1][p][q]);
d=1/c;
if(fabs(A[k-1][p][q])>0.005*fabs(A[k-1][p][p]-A[k-1][q][q]))
{
if(c>0) tanth=1/(c+sqrt(1+c*c));
else if(c<0) tanth=(-1)/(fabs(c)+sqrt(1+c*c));
}
else tanth=d/(1+sqrt(1+d*d));
costh=1/(sqrt(1+tanth*tanth));
sinth=tanth*costh;
}
for(i=0;i<10;i++)
{
for(j=0;j<10;j++)
{
if(i!=p&&i!=q&&j!=p&&j!=q) A[k][i][j]=A[k-1][i][j];
else if(i==p&&j!=p&&j!=q) A[k][i][j]=A[k-1][p][j]*costh+A[k-1][q][j]*sinth;
else if(i==q&&j!=p&&j!=q) A[k][i][j]=(-1)*A[k-1][p][j]*sinth+A[k-1][q][j]*costh;
else if(j==p&&i!=p&&i!=q) A[k][i][j]=A[k-1][p][i]*costh+A[k-1][q][i]*sinth;
else if(j==q&&i!=p&&i!=q) A[k][i][j]=(-1)*A[k-1][p][i]*sinth+A[k-1][q][i]*costh;
else if(i==p&&j==p) A[k][i][j]=A[k-1][p][p]*costh*costh+2*A[k-1][p][q]*sinth*costh+A[k-1][q][q]*sinth*sinth;
else if(i==q&&j==q) A[k][i][j]=A[k-1][p][p]*sinth*sinth-2*A[k-1][p][q]*sinth*costh+A[k-1][q][q]*costh*costh;
else if((i==p&&j==q)||(i==q&&j==p)) A[k][i][j]=0;
}
}
for(i=0;i<10;i++)
{
for(j=0;j<10;j++)
{
if(i==j)
{
if(i==p||i==q) R[i][j]=costh;
else R[i][j]=1;
}
else
{
if(i==p&&j==q) R[i][j]=sinth;
else if(i==q&&j==p) R[i][j]=(-1)*sinth;
else R[i][j]=0;
}
}
}
for(i=0;i<10;i++)
{
for(j=0;j<10;j++)
{
for(m=0;m<10;m++)
{
sum=sum+Rb[k-1][i][m]*R[m][j];
}
Rb[k][i][j]=sum;
sum=0;
}
}
for(i=0,pfh=0;i<10;i++)
{
for(j=0;j<10;j++)
{
if(i!=j) pfh=pfh+A[k][i][j]*A[k][i][j];
}
}
k++;
}
for(i=0;i<10;i++)
{
printf("特征值%d为:%lf\n",i+1,A[k][i][i]);
printf("对应特征向量的转置为:\n");
for(j=0;j<10;j++)
{
printf("%lf ",Rb[k][j][i]);
}
}
}
参考资料