计算方法：矩阵三角分解法解线性方程组

#include<stdio.h>
#define N 100
int print(double a[N][N],int n)
{
   int i,j;
   for(i=1;i<=n;i++)
   {
    for(j=1;j<=n;j++)
     printf("%.4f\t",a[i][j]);
   printf("\n");
  }
  return 0;
}  
      
int main()
{
  double a[N][N+1]={0},l[N][N]={0},u[N][N]={0},y[N]={0},x[N]={0},sum;
   int i,j,k,n;
   printf("请输入方程的阶数:\n");
  scanf("%d",&n);
  printf("请输入方程增广矩阵:\n");  
  for(i=1;i<=n;i++)//输入方程增广矩阵
   for(j=1;j<=n+1;j++)
    scanf("%lf",&a[i][j]);
    for(i=1;i<=n;i++)//初始化l数组的对角线上的元素，全部为1
          l[i][i]=1;    
  for(k=1;k<=n;k++)
  {
             for(i=k;i<=n;i++)//用L的第K行去乘U的每一列，得到U的第K行，i是第K行的第i列
             {   
      for(j=1,sum=0;j<=k-1;j++)//j是用来计算第i列之前的数的和
          sum=sum+l[k][j]*u[j][i];
      u[k][i]=a[k][i]-sum;
       } 
       for(i=k+1;i<=n;i++)//用L的每一行去乘U的第K列，得到L的第K列，其中i是第K列的第i行
       {   
      for(j=1,sum=0;j<=k-1;j++)//j是指第i列前面的数的和
          sum=sum+l[i][j]*u[j][k];
      l[i][k]=(a[i][k]-sum)/u[k][k];
       }
   }
   printf("L阵为:\n");
   print(l,n);
   printf("U阵为:\n");
   print(u,n);
  for(i=1;i<=n;i++)
  {
        for(j=1,sum=0;j<=i-1;j++)
       sum=l[i][j]*y[j]+sum;//注意乘的是y[j]而不是y[i]，
        y[i]=a[i][n+1]-sum; 
  }printf("中间向量y为:\n");
  for(i=1;i<=n;i++)
   printf("y[%d]=%f\n",i,y[i]);
  for(i=n;i>=1;i--)
  {
       for(j=n,sum=0;j>i;j--)
           sum=u[i][j]*x[j]+sum;
      x[i]=(y[i]-sum)/u[i][i];
  }
  printf("方程的根为:\n");
  for(i=1;i<=n;i++)
      printf("x%d=%f\n",i,x[i]);
  system("pause");
  return 0;
}           
                          
                    
 /*3
11 -3 -2 3
-23 11 1 0
1 -2 2 -1
x3=1.001,x2=2.001,x1=1.000

3
1 4 -5 3
1 3 -2 2
6 -1 18 2
x1=1.333333.x2=0,x3=-0.33333

3
1 -1 1 -4
5 -4 3 -12
2 1 1 11
x1=3,x2=6,x3=-1

 

3
2 3 5 5
3 4 7 6
1 3 3 5
x1=-4,x2=1,x3=2(书上答案有误)

4
1.1161 0.1254 0.1397 0.1490 1.5471
0.1582 1.1675 0.1768 0.1871 1.6471
0.1968 0.2071 1.2168 0.2271 1.7471
0.2368 0.2471 0.2568 1.2671 1.8471
x1=1.0406,x2=0.9870,x3=0.9351,x4=0.8813

3
2 1 2 6
4 5 4 18
6 -3 5 5
x1=1,x2=2,x3=1
*/