Lu分解法的C語言實作

将系數矩陣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[][] ,double U[][], int i, int j ){
    double sU = ;
    for (int k = ; k <= i- ; k++)
    {
        sU += L[i-][k-] * U[k-][j-];
    }

    return sU;
}//計算求和1

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

    return sL;

}//計算求和2

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

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

    double a[][] = { {,,,,},
                    {,,,-,},
                    {,-,-,,},
                    {,,-,}};//将系數存入二維數組
    double L[][] = {};
    double U[][] = {};//初始化部分
    double b[] = {,,,};
    int n = ;//n階
    //輸出[Ab]
    printf("[A]:\n");
    for (int i = ; i <= n; i++)
    {
        for (int j = ; j <= n; j++)
        {
            printf("%f\t", a[i-][j-]);
        }
        printf("\n");
    }

    //計算L,U
    for (int i = ; i <= n; i++)
    {
        L[i-][i-] = ;//對角線元素為1
        for (int j = i; j <= n; j++)
        {
            //由于數組下标從0開始 是以i-1,j-1
            U[i-][j-] = a[i-][j-] - sumU(L,U,i,j);

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

    }





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

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

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

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

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



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

運作結果：

是以

x1=1

x2= 2

x3=-1

x4 = 3