1.txt

使用带双步位移QR法求解矩阵的特征值和特征向量

#include "stdio.h"
#include "math.h"
int m=10;
#define N 10
double Re[10],Im[10],z[10];
#include<iostream>
#include<iomanip>
double Hessenberg_fun(double A[][N]);
double Householder(double C[][N],double B[][N],int n);
using namespace std;



/*求特征根*/ 
void roots(double b,double c)
       {double disc;
       	disc=b*b-4*c;
       	if(disc==0)
       	  {Re[m-1]=b/2;Re[m-2]=b/2;
       	   Im[m-1]=0;Im[m-2]=0;
       	  }
       	else if(disc>0)
       	   {Re[m-1]=(b+sqrt(disc))/2;
       	   	Re[m-2]=(b-sqrt(disc))/2;
       	   	Im[m-1]=0;Im[m-2]=0;
       	   }
       	else
       	   {Re[m-1]=b/2;Re[m-2]=b/2;
       	   	Im[m-1]=sqrt(-disc)/2;
       	   	Im[m-2]=sqrt(-disc)/2;
       	   }
       }


/*用列主元素高斯消去法解方程组*/        
void solve(double v[10],double a[10][10])
       {int i,j,k,t;
       	double max,e,f,sum;
       	for(k=0;k<=8;k++)
       	   {max=0;t=0;
       	   	for(i=k;i<=9;i++)
       	   	   {if(fabs(a[i][k])>max)
       	   	   	  {max=fabs(a[i][k]);t=i;
       	   	   	  }
       	   	   }
       	   	for(j=k;j<=9;j++)
       	   	   {e=a[t][j];a[t][j]=a[k][j];a[k][j]=e;
       	   	   }
       	   	f=v[t];v[t]=v[k];v[k]=f;
       	   	for(i=k+1;i<=9;i++)
       	   	   {e=a[i][k]/a[k][k];
       	   	   	for(j=k+1;k<=9;k++)
       	   	   	   {a[i][j]=a[i][j]-e*a[k][j];
       	   	   	   }
       	   	   	v[i]=v[i]-e*v[k];
       	   	   }
       	   }
       	 z[9]=v[9]/a[9][9];
       	 for(i=8;i>=0;i--)
       	    {sum=0;
       	     for(j=i+1;j<=9;j++)
       	        {sum=sum+a[i][j]*z[j];
       	        }
       	     z[i]=(v[i]-sum)/a[i][i];
       	    }
       }
       
       
/*主程序*/       
void main()
{
	int i,j,k,r,L;
	double ep,s,t,sum,w,u,g;
    double a[10][10],f[10][10],I[10][10],x[10][10],A[10][10],B[10][10],y[10],o[10];

	  ep=1.0e-12;
	  L=100;
	  w=0.0;
	  k=1;

    printf("矩阵A是：\n\n");
    for(i=0;i<=9;i++)
         {for(j=0;j<=9;j++)
         	 {if(i==j)
         	 	{a[i][j]=1.5*cos(i+1+1.2*(j+1));
         	 	 A[i][j]=1.5*cos(i+1+1.2*(j+1));
         	 	 I[i][j]=1;
         	 	}
         	  else
         	    {a[i][j]=sin(0.5*(i+1)+0.2*(j+1));
         	     A[i][j]=sin(0.5*(i+1)+0.2*(j+1));
         	     I[i][j]=0;
         	    }
         	  printf("%.11e  ",a[i][j]);
         	 }
          printf("\n\n");
         }

      Hessenberg_fun(a);

      printf("A的拟上三角矩阵是：\n\n");

      for(i=0;i<=9;i++)
         {for(j=0;j<=9;j++)
         	 {printf("%.11e  ",a[i][j]);
         	 }
          printf("\n\n");
         }
label_1:if(fabs(a[m-1][m-2])<=ep)
          {Re[m-1]=a[m-1][m-1];Im[m-1]=0;
           m=m-1;
           label_2:if(m==0)
                     {Re[0]=a[0][0];Im[0]=0;
                      goto label_3;
                     }
                   else
                     goto label_1;
          }
        else
          {s=a[m-2][m-2]+a[m-1][m-1];
           t=a[m-2][m-2]*a[m-1][m-1]-a[m-1][m-2]*a[m-2][m-1];
           if(m==1)
             {roots(s,t);goto label_3;
             }
           else
             {if(fabs(a[m-2][m-3])<=ep)
             	{roots(s,t);
             	 m=m-2;
             	 goto label_2;
             	}
              else
                {if(k==L)
                   {printf("未得到A的全部特征值。\n\n");
             	    goto label_3;
                   }
                  else 
                   {for(i=0;i<=9;i++)
                   	   {for(r=0;r<=9;r++)
                   	   	   {sum=0;
                   	   	   	for(j=0;j<=9;j++)
                   	   	   	   {sum=sum+a[i][j]*a[j][r];
                   	   	   	   }
                   	   	   	f[i][r]=sum;
                   	   	   }
                   	   }
                   	 for(i=0;i<=9;i++)
                   	    {for(j=0;j<=9;j++)
                   	    	{x[i][j]=f[i][j]-s*a[i][j]+t*I[i][j];
                   	    	}
                   	    }
                Householder(a,x,m);
                   }
                }
             }
          }
         k=k+1;
         goto label_1;
label_3:;
      printf("A的拟上三角阵经过QR分解后得到的矩阵如下：\n\n");
	cout<<"第1~5列";
	cout<<endl;
	for (i=0;i<N;i++)
	{
		for(j=0;j<5;j++)
			printf("%.11f\t",a[i][j]);
		cout<<endl;
	}
	cout<<endl;
	cout<<"第7~10列";
	cout<<endl;
	for (i=0;i<N;i++)
	{
		for(j=5;j<N;j++)
			printf("%.11f\t",a[i][j]);
		cout<<endl;
	}
      printf("A的全部特征值如下：\n\n");
      for(i=0;i<=9;i++)
      {
		  cout<<Re[i]<<setw(15)<<Im[i];
		  cout<<endl;
	  }
      printf("A的实特征值依次对应的特征向量如下：\n\n");
      for(i=0;i<=9;i++)
         {if(Im[i]==0)
         	{for(j=0;j<=9;j++)
         	   for(k=0;k<=9;k++)
                        B[j][k]=A[j][k]-Re[i]*I[j][k];

         	
          for(k=0;k<=9;k++)
             {z[k]=1;
             }
             w=0;
          do
  label_4:{sum=0;u=w;
             for(k=0;k<=9;k++)
                {sum=sum+z[k]*z[k];
                }
             g=sqrt(sum);
             for(k=0;k<=9;k++)
                {o[k]=y[k]=z[k]/g;
                }
             solve(y,B);
             sum=0;
             for(k=0;k<=9;k++)
                {sum=sum+o[k]*z[k];
                }
             w=sum;
             if(u==0)
               goto label_4;
            }
            while(fabs(1/w-1/u)/fabs(1/w)>0.005);
            for(k=0;k<=9;k++)
               {printf("%.11e  ",o[k]);
			cout<<endl;
               }
            printf("\n\n");
         	}
         }
     }
         	    	
         	    	
double Hessenberg_fun(double A[][N])
{
	int i,j,r;
	double A0[N][N],up[N][N],wu[N][N];
	double flag,c,d,h,t;	
	double u[N],p[N],q[N],w[N];

	for(r=0;r<N-1;r++)//拟上三角化
	{
		for(i=0;i<N;i++)
			for(j=0;j<N;j++)
				A0[i][j]=A[i][j];
		t=0.0;
		flag=0.0;
		for(i=0;i<N;i++)
		{
			u[i]=0.0;
		    p[i]=0.0;
		    q[i]=0.0;
		    w[i]=0.0;
		}

		for(i=r+2;i<N;i++)
			flag=flag+A0[i][r]*A0[i][r];

		if(flag!=0)
		{
			d=sqrt(flag+A0[r+1][r]*A0[r+1][r]);
			if(A[r+1][r]>0)  //若A[r+1][r]<0,c取正值
				c=-d;
			else
				c=d;   
			h=c*c-c*A[r+1][r];
			
			for(i=r+1;i<N;i++)  //给u赋值，与教材的p62的第三条对应
				u[i]=A0[i][r];
			u[r+1]=u[r+1]-c; 

			for(i=0;i<N;i++)
			{
				for(j=0;j<N;j++)
				{
					p[i]+=A0[j][i]*u[j]/h;
					q[i]+=A0[i][j]*u[j]/h;
				}
			}

			for(i=0;i<N;i++)
				t+=p[i]*u[i]/h;				

			for(i=0;i<N;i++)
				w[i]=q[i]-t*u[i];

			for(i=0;i<N;i++)
			{
				for(j=0;j<N;j++)
				{
					wu[i][j]=w[i]*u[j];
					up[i][j]=u[i]*p[j];
				}
			}
			
			for(i=0;i<N;i++)
				for(j=0;j<N;j++)
					A[i][j]=A0[i][j]-wu[i][j]-up[i][j];
		}
		else
			for(i=0;i<N;i++)
				for(j=0;j<N;j++)
					A[i][j]=A0[i][j];
	}
	return 1;
}       	    	
         	    	
       

double Householder(double C[][N],double B[][N],int n)
{
	int i,j,r;
	double c,d,h,t,flag;
	double p[N],q[N],u[N],v[N],w[N];

	for(r=0;r<n-1;r++)
	{
		t=0.0;
		flag=0.0;
		for(i=r+1;i<n;i++)
			flag+=B[i][r]*B[i][r];

		if(flag!=0)
		{
			d=sqrt(flag+B[r][r]*B[r][r]);
			if(B[r][r]<=0)
				c=d;
			else
				c=-d;
			h=c*c-c*B[r][r];
			for(i=0;i<n;i++)
			{
				u[i]=0.0;
				v[i]=0.0;
				w[i]=0.0;
				q[i]=0.0;
				p[i]=0.0;
			}

			u[r]=B[r][r]-c;
			for(i=r+1;i<n;i++)
				u[i]=B[i][r];

			for(i=0;i<n;i++)
				for(j=0;j<n;j++)
					v[i]+=B[j][i]*u[j]/h;

			for(i=0;i<n;i++)
				for(j=0;j<n;j++)
				{
					p[i]+=C[j][i]*u[j]/h;
					q[i]+=C[i][j]*u[j]/h;
				}

			for(i=0;i<n;i++)
				t+=p[i]*u[i]/h;  

			for(i=0;i<n;i++)
				w[i]=q[i]-t*u[i];
  
		    for(i=0;i<n;i++)
				for(j=0;j<n;j++)
				{
					B[i][j]=B[i][j]-u[i]*v[j];
					C[i][j]=C[i][j]-w[i]*u[j]-u[i]*p[j];
				}
		}
	}
	return 1;
}