aaaaaaeded
aaaaaaeded
采纳率100%
2020-05-20 21:04

C语言实现古典雅可比法求十阶矩阵的特征值和特征向量,为什么结果全是0?

已采纳
#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]);
        }
    }
}

参考资料图片说明图片说明图片说明

  • 点赞
  • 写回答
  • 关注问题
  • 收藏
  • 复制链接分享
  • 邀请回答

1条回答