平方根法矩阵分解以及直接解法.txt

数值方法课程中的程序

#include<iostream.h>
#include<math.h>
void Cholesky(double A[][4],double L[][4],int n)
{  int i,j,k;
   double s=0,t=0;
   for(j=0;j<n;j++)
   {  for(i=0;i<n;i++)
      {  if(i==j)
         {  s=0;
            for(k=0;k<=j-1;k++)
            {  s=s+L[j][k]*L[j][k];}
            L[j][j]=sqrt(A[j][j]-s);  
         }


         if(i>j)
         {   t=0;
             for(k=0;k<=j-1;k++)
             {  t=t+L[i][k]*L[j][k];}
             L[i][j]=(A[i][j]-t)/L[j][j];     
         }
         
         else if(i<j)  L[i][j]=0;
 
      }
   } 

}

void zhijiejiefa(double L[][4],double U[][4],double b[],double x[],int n)
{  int i,k;
   double s,t;
   double y[4];
   
   y[0]=b[0]/L[0][0];
   for(i=1;i<n;i++)
   {  s=0.0;
      for(k=0;k<i;k++)
      {  s=s+L[i][k]*y[k];}
      y[i]=(b[i]-s)/L[i][i];
   }

   x[n-1]=y[n-1]/U[n-1][n-1];
   for(i=n-2;i>=0;i--)
   {  t=0.0;
      for(k=i+1;k<=n-1;k++)
      {  t=t+U[i][k]*x[k];}
      x[i]=(y[i]-t)/U[i][i];
   }
}

void main()
{  double A[4][4]={{5,-4,1,0},{-4,6,-4,1},{1,-4,6,-4},{0,1,-4,5}};
   double L[4][4],U[4][4];
   double x[4],b[4]={2,-1,-1,2};
   
   int n=4,r,i;
   Cholesky(A,L,4);
   
   for(r=0;r<n;r++)
   {  for(i=0;i<n;i++)
      {  U[r][i]=L[i][r];
      }
   }
   
   cout<<"L矩阵为:"<<endl;
   for(r=0;r<n;r++)
   {  for(i=0;i<n;i++)
      {  cout<<L[r][i]<<" "<<" "<<" ";
      }
      cout<<endl;
   }
   
   
  
   cout<<endl<<"U矩阵(L的为转置矩阵）:"<<endl;
   for(r=0;r<n;r++)
   {  for(i=0;i<n;i++)
      {
         cout<<U[r][i]<<" "<<" "<<" ";
      }
      cout<<endl;
   }
  
  zhijiejiefa(L,U,b,x,4);
  cout<<endl<<"方程组的解为:"<<endl;
  for(i=0;i<n;i++)
  cout<<x[i]<<endl;

}