t5.cpp

采用双步移位的QR分解算法求解方程组的解

#include "stdio.h"
 #include "math.h"
/*******************************************************************主函数***************************************************/
void main ()
 {

/************************变量初始化****************************************************/
int i,j,tt,r,num,rr,ij;
double a[11][11],aa[11][11],a3[11][11],at[11][11];/* a为原始生成矩阵，aa为做双步QR分解时用到的矩阵，a3求解方程组时用到的矩阵*/
int zero=0;
double dr=0.0;
double ur[11],wr[11],vr[11];
double b[11][11],mkt[11][11],mk[11][11];     /*b为a的转置，mkt为mk的转置*/
double pr[11]={0,0,0,0,0,0,0,0,0,0,0};
double qr[11]={0,0,0,0,0,0,0,0,0,0,0};
double bb[11]={0,0,0,0,0,0,0,0,0,0,0};
double tr=0;
double error=10e-13;                         /*允许的迭代误差*/
double error1=10e-13;
double t11,t22,cr,dk,bk,bdk;
int m=10;
double solve[11],solve3[11],solve4[11];      /*存放特征值*/
double solve2[11]={0,0,0,0,0,0,0,0,0,0,0};
double solves[11]={0,0,0,0,0,0,0,0,0,0,0};
double solvess[11]={0,0,0,0,0,0,0,0,0,0,0};
double u0[11]={1,1,1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1};
int k=1;
double hr,s1,s2,s,t,s3,s4;
double x[11],y[11],l[11];
double u1[11][11],l1[11][11];
double mi=0;
int ik,im;
double max,middle;
double beta1=1;
double beta2=1;;
double beta0,eit,p;

for (i=1;i<11;i++)
 { for (j=1;j<11;j++)                       /*初始化矩阵mk*/
     mk[i][j]=0;
l1[i][j]=0;
u1[1][j]=0;}

for (i=1;i<11;i++)
 { for (j=1;j<11;j++)
 
 {if(i==j) a[i][j]=1.5*cos(i+1.2*j);       /*生成初始矩阵a*/ 
  else a[i][j]=sin(0.5*i+0.2*j);
  }
    }

for (i=1;i<11;i++)
{ for (j=1;j<11;j++)
{printf("%13e  ",a[i][j]);                  /*输出a*/
aa[i][j]=a[i][j];
}
printf("\n");
printf("\n");
printf("\n");}

printf("\n");
printf("\n");
printf("\n");

/************************************变量初始化结束*************************************************/

/************************a拟上三角化*************************************/

 for (r=1;r<9;r++)
 {

    for (i=r+2;i<11;i++)          /*判断是否为零*/
 { 
	 if (abs(a[i][r])<=error)      
	 {zero=zero+1;a[i][r]=0;}
      }
    if (zero==9-r)
    { for (i=1;i<11;i++)
	  { for (j=1;j<11;j++)
		
	    {a[i][j]=a[i][j];}
	}
	     
     continue;
      }
    else
    {
    
    for (i=r+1;i<11;i++)
    {dr=dr+a[i][r]*a[i][r];             /*求解dr*/
    }
    dr=sqrt(dr);

    if (a[r+1][r]<=0)
      cr=dr;
    else cr=-dr;                      /*求解cr*/
      hr=cr*cr-cr*a[r+1][r];           /*求解hr*/

     for (i=1;i<11;i++)
      { if (i<r+1)  ur[i]=0;
	else if(i==r+1)
	ur[i]=a[r+1][r]-cr;
	else ur[i]=a[i][r];               /*得到ur*/
	}

      for (i=1;i<11;i++)
      { for (j=1;j<11;j++)
	{b[i][j]=a[j][i];
	}
      }

      for (i=1;i<11;i++)
      { for (j=1;j<11;j++)
       {
       pr[i]=b[i][j]*ur[j]+pr[i];        /*得到pr,qr*/
       qr[i]=a[i][j]*ur[j]+qr[i];
       }
       pr[i]=pr[i]/hr;
       qr[i]=qr[i]/hr;
       }


       for (i=1;i<11;i++)
       { tr=tr+pr[i]*ur[i];              /*得到tr*/
       }
       tr=tr/hr;


       for (i=1;i<11;i++)
       {wr[i]=qr[i]-tr*ur[i];
       }                                 /*得到wr*/


       for (i=1;i<11;i++)
       {
       for (j=1;j<11;j++)

	{a[i][j]=a[i][j]-wr[i]*ur[j]-ur[i]*pr[j];
	  if (abs(a[i][j])<=error)
	  a[i][j]=0;}                              /*得到新a*/
	}


	}             /*大else循环*/

  zero=0;
  dr=0;
  cr=0;
  hr=0;
  tr=0;
  for (i=1;i<11;i++)
  {ur[i]=0;
  pr[i]=0;
  qr[i]=0;
  wr[i]=0;}

 }/*大for 循环*/

  /**********************************************矩阵拟上三角化结束*********************************************/
	for (i=1;i<11;i++)
    {for (j=1;j<11;j++)
	{printf("%13e  ",a[i][j]);}
        printf("\n");
		printf("\n");
        printf("\n");
    }
/*输出拟上三角化之后的a*/

		
	/**************************************************qr分解**************************************************/
  zero=0;
  dr=0;
  cr=0;
  hr=0;
  tr=0;
  for (i=1;i<11;i++)
  {ur[i]=0;
  pr[i]=0;
  qr[i]=0;
  wr[i]=0;
  vr[i]=0;}
  /***变量初始化结束*****/

	for (m=10;m>1;)

	 {
		
  if (abs(a[m][m-1])<=error1)
       { solve[k]=a[m][m];
        
         k=k+1;
         m=m-1;    /*第一次判断转4*/
		
       }                                            
  /* 否则转5*/
  else
  {
   dk=a[m-1][m-1]*a[m][m]-a[m-1][m]*a[m][m-1];
   bk=-(a[m-1][m-1]+a[m][m]);
   bdk=bk*bk-4*dk;
   if (bdk>=0)
   {s1=(-bk+sqrt(bdk))/2;s3=0;
    s2=(-bk-sqrt(bdk))/2;s4=0;
 
   }    /*得到两个特征值*/
  else 
   {s1=-bk/2;s3=sqrt(-bdk)/2;
    s2=-bk/2;s4=-sqrt(-bdk)/2;
  }
   if (m==2)
  {solve[k]=s1;solve2[k]=s3;
   solve[k+1]=s2;solve2[k+1]=s4;
  
   k=k+1;
   break;}                                /*如果m=2,得到两个特征值，break跳出循环*/
 else if (abs(a[m-1][m-2])<=error1)
   {solve[k]=s1;solve2[k]=s3;
   solve[k+1]=s2;solve2[k+1]=s4;
  
    k=k+2;
    m=m-2;
	
    if (m>1) continue;               /*判断是否为a的特征值*/
    else break;
   }        /*如果是m<=1,则跳出循环*/


/*****************************对mk进行qr分解****************************************/
 else
   {
	 s=a[m-1][m-1]+a[m][m];
     t=a[m-1][m-1]*a[m][m]-a[m][m-1]*a[m-1][m];
     for (i=1;i<m+1;i++)
     {
         for (j=1;j<m+1;j++)
         {
             for (ij=1;ij<m+1;ij++)
             { mk[i][j]=mk[i][j]+a[i][ij]*a[ij][j];}   /*得到mk*/
			 mk[i][j]=mk[i][j]-s*a[i][j];
             if (i==j)
             {mk[i][j]=mk[i][j]+t;}
         }}

          for (rr=1;rr<m;rr++)
     {
       for (i=rr+1;i<m+1;i++)
     { if (abs(mk[i][rr])<=error)
	   {zero=zero+1;                   /*判断是否需要进行分解*/
	   mk[i][rr]=0;}
      }
    if (zero==m-rr)
    {  
		for (i=1;i<m+1;i++)
	{ 
		for (j=1;j<m+1;j++)
	    {mk[i][j]=mk[i][j];
		a[i][j]=a[i][j];}
		}
}
    else
       {
           for (i=rr;i<m+1;i++)
	   { dr=dr+mk[i][rr]*mk[i][rr];           /*得到dr*/
           }
           dr=sqrt(dr);
           if (mk[rr][rr]<=0)
               cr=dr;
           else cr=-dr;                       /*得到cr*/
	   hr=cr*cr-cr*mk[rr][rr];

           for (i=1;i<m+1;i++)
           {
             if (i<rr)   ur[i]=0;
             else if (i==rr)
		   ur[i]=mk[rr][rr]-cr;
	       else ur[i]=mk[i][rr];             /*得到ur*/
           }

           for (i=1;i<m+1;i++)
           {   for (j=1;j<m+1;j++)
           { mkt[i][j]=mk[j][i];}}

           for (i=1;i<m+1;i++)
           {   for (j=1;j<m+1;j++)
           { vr[i]=mkt[i][j]*ur[j]+vr[i];}    /*得到vr*/
		      vr[i]=vr[i]/hr;
           }

           for (i=1;i<m+1;i++)
           {
               for (j=1;j<m+1;j++)
               {mk[i][j]=mk[i][j]-ur[i]*vr[j];    /*得到新的mk*/
               }
           }

           for (i=1;i<m+1;i++)
           {
               for (j=1;j<m+1;j++)
               {at[i][j]=a[j][i];
               }
           }

           for (i=1;i<m+1;i++)
           {
               for (j=1;j<m+1;j++)
               {pr[i]=at[i][j]*ur[j]+pr[i];
                qr[i]=a[i][j]*ur[j]+qr[i];        /*得到pr,qr*/
               }
           pr[i]=pr[i]/hr;
           qr[i]=qr[i]/hr;
           }

           for (i=1;i<m+1;i++)
           {tr=tr+pr[i]*ur[i];}
	       tr=tr/hr;                                  /*得到tr*/

            for (i=1;i<m+1;i++)                       /*得到wr*/
                {wr[i]=qr[i]-tr*ur[i];}

            for (i=1;i<m+1;i++)
                {
                    for (j=1;j<m+1;j++)
                    {
                        a[i][j]=a[i][j]-wr[i]*ur[j]-ur[i]*pr[j];                /*得到新的a*/
                        }
                    }
	}/*小else循环结束*/

	/***变量重新赋为零****/
 dr=0;
 cr=0;
 hr=0;
 tr=0;
zero=0;
for (i=1;i<11;i++)
{ vr[i]=0;
pr[i]=0;
qr[i]=0;
ur[i]=0;
wr[i]=0;}

     }/*小for 循环*/


		  for (i=1;i<11;i++)
		  {
			  for (j=1;j<11;j++)
			  { mk[i][j]=0;}
		  }
		  /*变量重新赋为零*/

	}/* 大else 循环*/ 

	 } /* 大else 循环*/ 




 }/*大FOR循环*/

 /**************************************带双步移位QR分解算法结束**********************************************************/

  if (m==1)
       {solve[k]=a[1][1];        /*得到一个a的特征值*/
        }
	  printf("\n");
	  printf("\n");

	  printf("\n");
	  for (i=1;i<11;i++)
	  {solve3[i]=solve[11-i];
	  solve4[i]=solve2[11-i];
	  printf("%13e,j%13e",solve[i],solve2[i]);   /*输出a的两个特征值*/
	  printf("\n");}
   
printf("\n");

printf("\n");
printf("\n");
for (i=1;i<11;i++)
    {for (j=1;j<11;j++)
	{  if (a[i][j]<=error)
	   a[i][j]=0;
		printf("%13e  ",a[i][j]);}
        printf("\n");
		printf("\n");
        printf("\n");
    }



/******************************************反幂法求特征向量*************************************************/
i=1;

 for (k=1;k<11;k++)
 {
	 if (solve2[k]==0)
	 { solves[i]=solve[k];
	 i++;}

 }
 num=i-1;/*得到实特征值及个数*/
 
 

 for (k=1;k<num+1;k++)
 {
	 for (j=k+1;j<num+1;j++)

 { if (solves[k]>solves[j])
   { middle=solves[j];
   solves[j]=solves[k];
   solves[k]=middle;}

	 }
 }/*实特征值按从小到大排列*/

 printf("\n");
for (i=1;i<7;i++)
{
	printf("%13e",solves[i]);
	printf("\n");}



/*求解特征向量*/


for (ik=1;ik<7;ik++)    /*对于6个实特征值求解特征向量*/
{   
  p=solves[ik]-0.01;    /*选择合适的移位p*/

for (i=1;i<11;i++)
{ if (i==1) u0[i]=0;
else u0[i]=1;}
/*初始化u0*/


	for (i=1;i<11;i++)
{
	for (j=1;j<11;j++)
	{ a3[i][j]=aa[i][j];
	   if (i==j)
		   a3[i][j]=a3[i][j]-p;
	   
	}
	
	}
	     /*得到新的A*/

for (k=1;k>0;k++)                /*用反幂法求解过程*/
{  
	
   eit=0;
	for (i=1;i<11;i++)
	{ eit=eit+u0[i]*u0[i];}
	eit=sqrt(eit);          /*利用u0得到eit*/

	for (i=1;i<11;i++)
	{ y[i]=u0[i]/eit;}       /*利用u0得到y*/

/*********************利用Doolittle分解法反解方程求得新的u0********************************************/
/******分解过程***********/
for (j=1;j<11;j++)
{
	u1[1][j]=a3[1][j];
}

for (i=2;i<11;i++)
{
	l1[i][1]=a3[i][1]/u1[1][1];
}


for (im=1;im<11;im++)
{
	for (j=im;j<11;j++)
	{

		middle=0.0;
		for (tt=1;tt<im;tt++)
		{ middle=middle+l1[im][tt]*u1[tt][j];}
		u1[im][j]=a3[im][j]-middle;
	}

	if (im<10)
	{
		for (i=im+1;i<11;i++)
		{
			middle=0.0;
			for (tt=1;tt<im;tt++)
			{
		    middle=middle+l1[i][tt]*u1[tt][im];}
			l1[i][im]=(a3[i][im]-middle)/u1[im][im];
		}

	}
}
/*************分解过程结束************************************************/
/*****求解过程*****/
l[1]=y[1];
for (i=2;i<11;i++)
{
	middle=0.0;
	for (tt=1;tt<i;tt++)
	{middle=middle+l1[i][tt]*l[tt];}

	l[i]=y[i]-middle;}

u0[10]=l[10]/u1[10][10];

for (i=9;i>0;i--)
{
	middle=0;
	for (tt=i+1;tt<11;tt++)
	{middle=middle+u1[i][tt]*u0[tt];}
	u0[i]=(l[i]-middle)/u1[i][i];
}

/****求解过程结束，得到新的u0************/


beta0=0.0;
for (i=1;i<11;i++)
{   
beta0=beta0+y[i]*u0[i];}         /*得到新的beta2*/

beta2=beta0;
t11=1.0/beta1;
t22=1.0/beta2;


if (k>1)
{
if (abs(t22-t11)/abs(t22)<=10e-13)     /**判断是否停止迭代**********/
{
goto loop;
}
else {beta1=beta2;}
}

else beta1=beta2;

}/*小 for 循环结束*/

loop:t22=1.0/(beta2)+p;

printf("%13e",t22);
printf("\n");
for (i=1;i<11;i++)
{
	x[i]=y[i];
	
	printf("%13e   ",x[i]);

}
/*输出特征向量*/
/************************************************反幂法结束********************************************************/

printf("\n");
printf("\n");




} /*大for 循环*/


     
}/*main*/

/****************************************************************主程序结束*************************************/