eigenvector2.txt

求矩阵的特征值的源码,用于求矩阵的特征值.

#include<stdio.h>
#include<math.h>
#define n 3
int sgn(double x)
{
  int h;
  if(x>0)h=1;
  if(x<0)h=-1;
  return h;
} 
  
void main()/*古典雅可比方法（精确稍慢）求全部特征值*/
{
  double v,t;
  static double a[n][n]={1.0,0.5,0.5,0.5,2,0.5,0.5,0.5,3},u[n][n],a1[n][n],a2[n][n],b1,b2,b3,sin,cos;
  int h,i,j,k,p,q;
  for(p=0;p<n;p++)
      for(q=0;q<n;q++)
          {a2[p][q]=a[p][q];
           if(p!=q&&fabs(a[p][q])>v)v=fabs(a[p][q]); 
          };
  for(h=0;h<50;h++)
{ 
  for(i=0;i<n;i++)
     for(j=0;j<n;j++)  
         if(i!=j&&fabs(a2[i][j])>=v)
            {b1=-a2[i][j];v=0;  
             b2=(a2[i][i]-a2[j][j])/2;
                 if(b2!=0)b3=sgn(b2)*b1/sqrt(b1*b1+b2*b2);
                 if(b2==0)b3=1;
             sin=b3/sqrt(2.0+2.0*sqrt(1.0-b3*b3));
             cos=sqrt(1-sin*sin);  
             for(p=0;p<n;p++)
                for(q=0;q<n;q++)
                    {u[p][q]=0;
                     if(p==i&&q==i)u[p][q]=cos;
                     if(p==j&&q==j)u[p][q]=cos;
                     if(p==i&&q==j)u[p][q]=-sin;
                     if(p==j&&q==i)u[p][q]=sin;
                     if(p==q&&q!=i&&q!=j)u[p][q]=1;
                    };
             for(p=0;p<n;p++)
                 for(q=0;q<n;q++)                       
                     for(k=0;k<n;k++)
                        a1[p][q]=a1[p][q]+u[p][k]*a2[k][q];
             for(p=0;p<n;p++)
                 for(q=0;q<n;q++)
                    a2[p][q]=0; 
             for(p=0;p<n;p++)   /*求转置矩阵*/
                 for(q=0;q<p;q++)
                    {t=u[p][q];u[p][q]=u[q][p];u[q][p]=t;};
             for(p=0;p<n;p++)
                 for(q=0;q<n;q++)                       
                    for(k=0;k<n;k++)
                         a2[p][q]=a2[p][q]+a1[p][k]*u[k][q];
             for(i=0;i<n;i++)/*求最大值*/
                 for(j=0;j<n;j++)
                     if(i!=j&&fabs(a2[i][j])>v)v=fabs(a2[i][j]);  
             for(p=0;p<n;p++)
                 for(q=0;q<n;q++)
                     a1[p][q]=0; 
           for(p=0;p<n;p++)
              for(q=0;q<n;q++)
                  {printf("%f\t",a2[p][q]);
                   if(q==n-1)printf("\n");
                  };  
           printf("h=%d,v=%f\n",h,v);        
           };
};     
}