三角分解法求线性方程.c

来自「计算方法的几个经典算法」· C语言 代码 · 共 81 行

#include<stdio.h>
#include<math.h>
 main()
{
   int n,i,j,k,p;
   double x[20],a[20][20],c;
   printf("输入方程组的阶n:");
   scanf("%d",&n);
   printf("输入系数矩阵:\n");


   for(i=1;i<=n;i++)
	   for(j=1;j<=n+1;j++)
		   scanf("%f",&a[i][j]); 


   if(a[1][1]==0)  
   {
	   printf("输出方法失败");
		   return 0;
   } 

   for(i=2;i<=n;i++)
	   a[i][1]=a[i][1]/a[1][1];

   for(k=2;k<=n-1;k++)
   {
	   for(j=k;j<=n+1;j++)
	   {
		   c=0;
		   for(p=1;p<=k-1;p++)
			   c=c+a[k][p]*a[p][j];
		   a[k][j]=a[k][j]-c;
	   }
		if(a[k][k]==0)
		{
		     printf("输出方法失败");
		        return 0;

		}
		c=0;
	   for(i=k+1;i<=n;i++)
	   {
		   c=0;
		   for(p=1;p<=k-1;p++)
			   c=c+a[i][p]*a[p][k];
		   a[i][k]=(a[i][k]-c)/a[k][k];
	   }
   }
   for(j=n;j<=n+1;j++)
   {
	   c=0;
	   for(p=1;p<=n-1;p++)
		   c=c+a[n][p]*a[p][j];
	   a[n][j]=a[n][j]-c;
   }
   if(a[n][n]==0)
   {
	  printf("输出方法失败");
		        return 0;
   }
   x[n]=a[n][n+1]/a[n][n];
   for(i=n-1;i>=1;i--)
   {
	   c=0;
	   for(j=i+1;j<=n;j++)
		   c=c+a[i][j]*x[j];
	   x[i]=(a[i][n+1]-c)/a[i][i];
   }
   for(i=1;i<=n;i++)
   {
	   printf("x[%d]=%lf",i,x[i]);
	   if((i%5)==0) printf("\n");
   }
   return 1;
}