LU分解的C++實現

　　又是一次數值科學與計算方法的實驗題目，LU分解的推導就不贅述，其核心公式如下：

$u_{1i}=a_{1i}   (i=1,2,3,\cdots ,n) $

$l_{i1}=a_{i1}/u_{11} ( i=2,3,\cdots ,n)$

$u_{ri}=a_{ri}-\sum_{k=1}^{r-1}l_{rk}u_{ki}    (i=r,r+1,r+2,\cdots ,n)$

$l_{ir}=(a_{ir}-\sum_{k=1}^{r-1}l_{ik}u_{kr}) /u_{rr}    (i=r+1,\cdots ,n;且r\ne n)$

　　C++實現方面，首先LU分解的函數傳入兩個參數，方陣的一階數組和方陣的階數（方陣用一維數組的行優先表示）。

　　計算步驟：

　　1. 初始化LU矩陣，L矩陣上三角為0，對角線為1，U矩陣下三角為0.

　　2. 計算U矩陣的第一行和L矩陣的第一列

　　3. 循環計算U矩陣和L矩陣，第一層循環表示U計算到第幾行，同時也表示L計算到第幾列，先計算U后計算L。第二層循環分別表示U矩陣改行的第幾列元素。第三層循環就是公式中的求和符號部分。

程序如下：

#include <iostream>

// 參數：一個order階矩陣，和矩陣的階數
void lowerUpperFactor(double *matrix, int order) {
    printf("--------原矩陣：--------\n");
    printMatrix(matrix, order,order);
    
    // 結果變量  L矩陣和U矩陣都是order階矩陣
    double *L = new double[order*order];
    double *U = new double[order*order];

    // 初始化全為0
    for (int i = 0; i < order; i++) {
        // 初始化U下三角為0
        for (int j = 0; j < i; j++) {
            *(U + i * order + j) = 0;
        }
        //初始化L對角線為1，上三角為0
        *(L + i * order + i) = 1;
        for (int j = i + 1; j < order; j++) {
            *(L + i * order + j) = 0;
        }
    }
    // 計算U的第一行和L的第一列
    int i = 0;
    for (i = 0; i < order; i++) {
        *(U + i) = *(matrix + i);
    }
    for (i = 1; i < order; i++) {
        *(L + i * order) = *(matrix + i * order) / *U;
    }
    // 計算其余行列
    int temp;
    for (int i = 1; i < order; i++) {
        // 計算矩陣U
        for (int j = i; j < order; j++) {
            temp = 0;
            for (int k = 0; k < i; k++) {
                temp+= (*(U + k * order + j) * (*(L + i * order + k)));
            }
            *(U + i * order + j) = *(matrix + i * order + j) - temp;
        }
        // 計算矩陣L
        for (int j = i+1; j < order; j++) {
            temp = 0;
            for (int k = 0; k < i; k++) {
                temp += *(U + k * order + i) * (*(L + j * order + k));
            }
            *(L + j * order + i) = (*(matrix +j * order + i) - temp) / (*(U+i* order + i));
        }
    }

    printf("------矩陣U------\n");
    printMatrix(U, order,order);
    printf("------矩陣L------\n");
    printMatrix(L, order, order);

    if (L) {
        delete[] L;
    }
    if (U) {
        delete[] U;
    }
}

void printMatrix(double *matrix, int row, int column) {
　　for (int i = 0; i < row; i++) {
　　　　for (int j = 0; j < column; j++) {
　　　　　　printf("%6.2lf ", *(matrix + i * column + j));
　　　　}
　　　　printf("\n");
　　}
}

int main() {

    //double matrix[] = { 1,2,3,2,5,2,3,1,5 };
    //int order = 3;
    double matrix1[] = { 1,1,-1,2,1,2,0,2,-1,-1,2,0,0,0,-1,1 };
    int order1 = 4;

    lowerUpperFactor(matrix1, order1);
}

提供兩個算例驗證其正確性和通用性：

1. 這個取自《數值分析》第四版孫家廣中的例題

通過程序計算如下：

2. 這個取自《Numerical Analysis-Burden Faires》中的例題

程序計算如下：