2.lu分解法.cpp

LU分解法

#define N 3
#include <stdio.h>
 void main()
{   
	 double a[N][N+1],b,X[N],sum1=0,sum2=0,s=0,sum=0,y[N];
     int i,j,k,t;
     printf("please input the data:");
     for(i=0;i<=N-1;i++)
     for(j=0;j<=N;j++)
	{  scanf("%lf",&b);
       a[i][j]=b;
	}
     for(i=0;i<=N-1;i++)
     for(j=0;j<=N;j++)
	{  printf("%lf ",a[i][j]);
       if((j+1)%(N+1)==0)
       printf("\n");
	}
printf("\n");
     for(k=0;k<=N-1;k++)
	{
       for(j=k;j<=N-1;j++)
	   {  
		   for(t=0;t<=k-1;t++)
		   {   sum1+=a[k][t]*a[t][j];
               sum2+=a[j+1][t]*a[t][k];
		   }
          a[k][j]=a[k][j]-sum1;
          a[j+1][k]=(a[j+1][k]-sum2)/a[k][k];
          sum1=0;sum2=0;
	   }
	 }
      for(i=0;i<=N-1;i++)
      for(j=0;j<=N;j++)
	 { 
		  printf("%lf ",a[i][j]);
          if((j+1)%(N+1)==0)
          printf("\n");
	 }
printf("\n");
y[0]=a[0][N];
      for(i=1;i<=N-1;i++)
	 {
		 {  for(t=0;t<=i-1;t++)
            s+=a[i][t]*y[t];
		 }
          y[i]=a[i][N]-s;
          s=0;
	 }
      for(i=0;i<=N-1;i++)
	 { 
		  printf("y[%d]=%lf ",i+1,y[i]);
	 }
      printf("\n");
      X[N-1]=y[N-1]/a[N-1][N-1];
      for(k=N-2;k>=0;k--)
      for(j=k+1;j<=N-1;j++)
	 {  
		  sum+=a[k][j]*X[j];
          X[k]=(y[k]-sum)/a[k][k];
	 }
for(i=0;i<=N-1;i++)
printf("X[%d]=%lf ",i+1,X[i]);
printf("\n");
}