將係數矩陣A轉變成等價兩個矩陣L和U的乘積
，其中L和U分別是單位下三角矩陣和上三角矩陣。當A的所有順序主子式都不為0時，矩陣A可以分解為A=LU（所有順序主子式不為0，矩陣不一定不可以進行LU分解）。其中L是下三角矩陣，U是上三角矩陣。

LU分解在本質上是高斯消元法的一種表達形式。實質上是將A通過初等行變換變成一個上三角矩陣，其變換矩陣就是一個單位下三角矩陣。這正是所謂的杜爾裡特演算法（Doolittle algorithm）：從下至上地對矩陣A做初等行變換，將對角線左下方的元素變成零，然後再證明這些行變換的效果等同於左乘一系列單位下三角矩陣，這一系列單位下三角矩陣的乘積的逆就是L矩陣，它也是一個單位下三角矩陣。

–

構造矩陣A，b

                {1,2,3,4, }
    A=[         {1,4,2,-8,}
                {1,-1,4,1}   ]
                {1,3,5,2}

b = [14,-17,2,8]^T

程式碼實現：

#include <stdio.h>
#include <stdlib.h>
//LU分解法實現解線性方程組
//copyright @ Mryang

double sumU(double L[4][4] ,double U[4][4], int i, int
 
 j ){
    double sU = 0.0;
    for (int k = 1; k <= i-1 ; k++)
    {
        sU += L[i-1][k-1] * U[k-1][j-1];
    }

    return sU;
}//計算求和1

double sumL(double L[4][4] ,double U[4][4], int i, int j ){
    double sL = 0.0;
    for (int k = 0; k <= j-1; k++)
    {
        sL += L[i-1][k-1] * U[k-1][j-1
 
];
    }

    return sL;

}//計算求和2

double sumY(double L[4][4] ,double y[4],int i){
    double sY=0.0;
    for (int k = 1; k <= i - 1; k++)
    {
        sY += L[i-1][k-1] * y[k-1];
    }
    return sY;
}//計算求和3

double sumX(double U[4][4] ,double x[4],int i ,int m){
    double sX = 0.0;
    for (int k = i+1; k <= m; k++)
    {
        sX += U[i-1][k-1] * x[k-1];
    }
    return sX;
}//計算求和4
int main(){

    double a[4][4] = { {1,2,3,1,},
                    {1,4,1,-1,},
                    {1,-1,-2,3,},
                    {1,3,-1,2}};//將係數存入二維陣列
    double L[4][4] = {0};
    double U[4][4] = {0};//初始化部分
    double b[4] = {8,8,12,19};
    int n = 4;//n階
    //輸出[Ab]
    printf("[A]:\n");
    for (int i = 1; i <= n; i++)
    {
        for (int j = 1; j <= n; j++)
        {
            printf("%f\t", a[i-1][j-1]);
        }
        printf("\n");
    }

    //計算L,U
    for (int i = 1; i <= n; i++)
    {
        L[i-1][i-1] = 1;//對角線元素為1
        for (int j = i; j <= n; j++)
        {
            //由於陣列下標從0開始 所以i-1,j-1
            U[i-1][j-1] = a[i-1][j-1] - sumU(L,U,i,j);

            if(j+1 <= n) L[j][i-1] = (a[j][i-1] - sumL(L,U,j+1,i))/U[i-1][i-1];//i變j+1，j變i 
        }

    }





    //輸出U
    printf("U:\n");
    for (int i = 1; i <= n; i++)
    {
        for (int j = 1; j <= n; j++)
        {
            printf("%f\t",U[i-1][j-1]);
        }
        printf("\n");
    }
    //輸出L
    printf("L:\n");
    for (int i = 1; i <= n; i++)
    {
        for (int j = 1; j <= n; j++)
        {
            printf("%f\t",L[i-1][j-1]);
        }
        printf("\n");
    }

    //由Ly=b 求y
    double y[4] = {0.0};
    y[0] = b[0];//y(1) = b(1);

    for (int i = 2; i <= n; i++)
    {
        y[i-1] = b[i-1] - sumY(L,y,i);
    }

    //由 Ux=y 求x
    double x[4] = {0.0};

    for (int i = n; i >= 1; i--)
    {
        x[i-1] =( y[i-1] - sumX(U,x,i,n))/ U[i-1][i-1];
    }



    //輸出y
    printf("y:\n");
    for (int i = 0; i < n; i++)
    {
        printf("%f\n",y[i]);
    }
    printf("\n");
    //輸出x
    printf("x:\n");
    for (int i = 0; i < n; i++)
    {
        printf("%f\n",x[i]);
    }
    printf("\n");
    system("pause");
    return 0;
}

執行結果：

所以
x1=1
x2= 2
x3=-1
x4 = 3