LU分解C++实现

1、待定系数法LU分解

设$A=LU$,即
$\left[\begin{array}{ccccc}{a}_{11}& a_{12}& a_{13}& \cdots & a_{1n}\\ a_{21}& a_{22}& a_{23}& \cdots & a_{2n}\\ a_{31}& a_{32}& a_{33}& \cdots & a_{3n}\\ ⋮& & \ddots & & \\ a_{n1}& a_{n2}& a_{n3}& \cdots & a_{nn}\end{array}\right]$=$\left[\begin{array}{ccccc}1& & & & \\ l_{21}& 1& & & \\ l_{31}& l_{32}& 1& & \\ ⋮& & \ddots & ⋮& \\ l_{n1}& l_{n2}& \cdots & l_{n,n-1}& 1\end{array}\right]$ $\left[\begin{array}{ccccc}{u}_{11}& u_{12}& u_{13}& \cdots & u_{1n}\\ & u_{22}& u_{23}& \cdots & u_{2n}\\ & & u_{33}& \cdots & u_{3n}\\ & & & ⋮& \ddots \\ & & & & u_{nn}\end{array}\right]$

1、比较等式两边的第一行，可得
$$u_{1j}=a_{1j},j=1,2,....,n.$$
再比较等式两边的第一列，可得
$$a_{i1}=l_{i1}{u}_{11}\Rightarrow l_{i1}=a_{i1}/u_{11},i=2,3,...,n$$
2、比较等式两边的第二行，可得
$$a_{2j}=l_{21}{u}_{1j}+u_{2j}\Rightarrow u_{2j}=a_{2j}-l_{21}{u}_{1j},j=2,3,...,n$$
再比较等式两边的第二列，可得
$$a_{i2}=l_{i1}{u}_{12}+l_{i2}+u_{22}\Rightarrow l_{i1}=(a_{i2}-l_{i1}{u}_{12})/u_{22},i=2,3,...,n$$
3、依次类推比较等式两边的第k行，可得
$$u_{kjj}=a_{kj}-(l_{k1}{u}_{1j}+\cdots +l_{k,k-1}{u}_{k-1,j}),j=k,k+1,...,n$$
依次类推比较等式两边的第k列，可得
$$l_{ik}=(a_{ik}-l_{i1}{u}_{1k}-\cdots -l_{i,k-1}{u}_{k-1,k})/u_{kk},i=k,k+1,...,n$$
直到第n步算出L和U的所有元素。

2、总结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);
}

3、方程求解

现在已经知道矩阵 $Ax=b$的LU分解了，就要去求解方程的解了，已知方程$Ax=b$

矩阵A可做如下分解 $A=LU$
则$$Ax=b\iff LUx=b$$
令
$$Ux=y$$
则
$$Ax=b\iff Ly=b\Rightarrow y=L^{-1}bUx=y\Rightarrow x=U^{-1}y$$

算法总共分为两步：

向前回代求解$LUx=b$
向后回代求解 $Ux=y$
如果矩阵 A 不进行分解直接通过$A^{-1}b$ 来求解方程运算量会比较大，因为通常来说 是不容易求的。而将 A 分解为 LU 的形式，单位下三角矩阵L 的逆和上三角矩阵 U 的逆是容易求的，因此很容易可以求出 y和x ，运算量将非常小。