LU分解原理解释及C语言实现

        本文为笔者在编写spice时，需利用LU求解线性方程组时所写。

1、杜尔里特算法：

        从下至上对矩阵A做初等行变换，将对角线左下方的元素编成零，即得到U矩阵。然后可以证明这些行变换的效果等同于左乘一系列单位下三角矩阵，这一系列单位下三角矩阵的乘积的逆就是L矩阵，它也是一个单位下三角矩阵。即杜尔特里算法证实了矩阵A可以唯一分解成为L和U矩阵。

2、LU分解原理概要：

        将方程组的系数矩阵分解为两个三角矩阵的乘积，以四阶为例进行讨论，即

                        

     可以改写为

                                                                    

    引入中间变量,

                                                                       

                                                                       

    先由一式求出,再带入二式求出即可。

3、LU各矩阵求解图解：

         A、此处展示了求解的过程

       B、求解

      

        C、求解

 

4、算法到代码的思路讲解

                

5、C语言代码

5.1 主函数中进行LU分解

此所求解的方程组为：

T=[1 3 2 5

      1 1  1 1

       2  1  3  4

       1   1   2  1];

B=[33,10,29,13];

解出X=[1,2,3,4]

#include <stdio.h>  
#include <stdlib.h>  
  
int main()  
{  
  int n=4;  
  int i,j,k,i_u,j_u;  
  double T[4][4]={1,3,2,5,1,1,1,1,2,1,3,4,1,1,2,1};  
  double B[4]={33,10,29,13};  
  double X[n],L[n][n],U[n][n],y[n];  
  //1、将L,U矩阵全部赋值为0  
  for(i=0;i<=n-1;i++)  
  {  
      for(j=0;j<=n-1;j++)  
      {  
          L[i][j]=0;  
          U[i][j]=0;  
      }  
  }  
  //2、进行最初的赋值  
  for(i=0;i<=n-1;i++)  
  {  
      U[0][i]=T[0][i];  
      L[i][i]=1;  
  }  
  //3、进行规律性赋值  
  for(j=0;j<=n-2;j++)  
  {  
      for(i=j+1;i<=n-1;i++)  
      {  
          //学会西格玛求和的编程方式  
          L[i][j]=T[i][j];  
          for(k=0;k<=j-1;k++)  
          {  
              L[i][j]=L[i][j]-L[i][k]*U[k][j];  
          }  
          L[i][j]=L[i][j]/U[j][j];  
      }  
      for(i=j+1;i<=n-1;i++)  
      {  
        i_u=j+1;  
        j_u=i;  
        U[i_u][j_u]=T[i_u][j_u];  
        for(k=0;k<=i_u-1;k++)  
        {  
            U[i_u][j_u]=U[i_u][j_u]-L[i_u][k]*U[k][j_u];  
        }  
      }  
  }  
  //4、求y  
  for(i=0;i<=n-1;i++)  
  {  
      y[i]=B[i];  
      for(k=0;k<=i-1;k++)  
      {  
          y[i]=y[i]-L[i][k]*y[k];  
      }  
  }  
  //5、求X  
  for(i=n-1;i>=0;i--)  
  {  
      X[i]=y[i];  
      for(k=i+1;k<=n-1;k++)  
      {  
          X[i]=X[i]-U[i][k]*X[k];  
      }  
      X[i]=X[i]/U[i][i];  
  }  
  return 0;  
}  

代码效果展示：

5.2 编写LU_decompositon子函数，通过调用函数的形式进行LU分解

#include <stdio.h>
#include <stdlib.h>

int main()
{
  int n=4;
  int i,j,k,i_u,j_u;
  double T[4][4]={1,3,2,5,1,1,1,1,2,1,3,4,1,1,2,1};
  double B[4]={33,10,29,13};
  double X[n],L[n][n],U[n][n],y[n];
  void LU_decomposition(int n,double(*T)[n],double *X,double *B);
  LU_decomposition(n,T,X,B);//调用LU分解函数进行LU分解求解方程组
  return 0;
}

void LU_decomposition(int n,double(*T)[n],double *X,double *B)
{
/*************************************************************
 函数名称：LU_decomposition
 函数功能描述：对TX=B进行LU分解求解
 形参说明：n:T方阵的阶数；T：混合系数矩阵；X：未知数矩阵；B：RHS
 作者：Mr.Liu @SDU
 创造时间：2021-08-27 21:44
**************************************************************/
  int i,j,k,i_u,j_u;
  double L[n][n],U[n][n],y[n];
  //1、将L,U矩阵全部赋值为0
  for(i=0;i<=n-1;i++)
  {
      for(j=0;j<=n-1;j++)
      {
          L[i][j]=0;
          U[i][j]=0;
      }
  }
  //2、进行最初的赋值
  for(i=0;i<=n-1;i++)
  {
      U[0][i]=*(*(T+0)+i);
      L[i][i]=1;
  }
  //3、进行规律性赋值
  for(j=0;j<=n-2;j++)
  {
      for(i=j+1;i<=n-1;i++)
      {
          //学会西格玛求和的编程方式
          L[i][j]=*(*(T+i)+j);
          for(k=0;k<=j-1;k++)
          {
              L[i][j]=L[i][j]-L[i][k]*U[k][j];
          }
          L[i][j]=L[i][j]/U[j][j];
      }
      for(i=j+1;i<=n-1;i++)
      {
        i_u=j+1;
        j_u=i;
        U[i_u][j_u]=*(*(T+i_u)+j_u);
        for(k=0;k<=i_u-1;k++)
        {
            U[i_u][j_u]=U[i_u][j_u]-L[i_u][k]*U[k][j_u];
        }
      }
  }
  //4、求y
  for(i=0;i<=n-1;i++)
  {
      y[i]=*(B+i);
      for(k=0;k<=i-1;k++)
      {
          y[i]=y[i]-L[i][k]*y[k];
      }
  }
  //5、求X
  for(i=n-1;i>=0;i--)
  {
      *(X+i)=y[i];
      for(k=i+1;k<=n-1;k++)
      {
          *(X+i)=*(X+i)-U[i][k]*(*(X+k));
      }
      *(X+i)=*(X+i)/U[i][i];
  }
  return 0;
}

代码效果展示：