LU分解C++实现

• 1、待定系数法LU分解
• 2、总结LU分解公式
• 3、方程求解

设 A = L U A=LU A=LU,即
[ a 11 a 12 a 13 ⋯ a 1 n a 21 a 22 a 23 ⋯ a 2 n a 31 a 32 a 33 ⋯ a 3 n ⋮ ⋱ a n 1 a n 2 a n 3 ⋯ a n n ] begin{bmatrix} 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} \ vdots&& ddots&\a_{n1}&a_{n2}&a_{n3}& cdots&a_{nn} \ end{bmatrix} ⎣ ⎡​a11​a21​a31​⋮an1​​a12​a22​a32​an2​​a13​a23​a33​⋱an3​​⋯⋯⋯⋯​a1n​a2n​a3n​ann​​⎦ ⎤​= [ 1 l 21 1 l 31 l 32 1 ⋮ ⋱ ⋮ l n 1 l n 2 ⋯ l n , n − 1 1 ] begin{bmatrix} 1&&&& \ l_{21}&1&& & \ l_{31}&l_{32}&1& & \ vdots&& ddots&vdots\l_{n1}&l_{n2}& cdots&l_{n,n-1}&1 \ end{bmatrix} ⎣ ⎡​1l21​l31​⋮ln1​​1l32​ln2​​1⋱⋯​⋮ln,n−1​​1​⎦ ⎤​ [ u 11 u 12 u 13 ⋯ u 1 n u 22 u 23 ⋯ u 2 n u 33 ⋯ u 3 n ⋮ ⋱ u n n ] begin{bmatrix} u_{11}&u_{12}&u_{13}&cdots&u_{1n} \ &u_{22}&u_{23}&cdots & u_{2n}\ &&u_{33}& cdots&u_{3n} \&&&vdots&ddots\&&&&u_{nn} \ end{bmatrix} ⎣ ⎡​u11​​u12​u22​​u13​u23​u33​​⋯⋯⋯⋮​u1n​u2n​u3n​⋱unn​​⎦ ⎤​

1、比较等式两边的第一行，可得
u 1 j = a 1 j , j = 1 , 2 , . . . . , n . u_{1j}=a_{1j}, j=1,2,....,n. u1j​=a1j​,j=1,2,....,n.
再比较等式两边的第一列，可得
a i 1 = l i 1 u 11 ⇒ l i 1 = a i 1 / u 11 , i = 2 , 3 , . . . , n a_{i1}=l_{i1}u_{11} Rightarrow l_{i1}=a_{i1}/u_{11}, i=2,3,...,n ai1​=li1​u11​⇒li1​=ai1​/u11​,i=2,3,...,n
2、比较等式两边的第二行，可得
a 2 j = l 21 u 1 j + u 2 j ⇒ u 2 j = a 2 j − l 21 u 1 j , j = 2 , 3 , . . . , n a_{2j}=l_{21}u_{1j}+u_{2j}Rightarrow u_{2j}=a_{2j}-l_{21}u_{1j}, j=2,3,...,n a2j​=l21​u1j​+u2j​⇒u2j​=a2j​−l21​u1j​,j=2,3,...,n
再比较等式两边的第二列，可得
a i 2 = l i 1 u 12 + l i 2 + u 22 ⇒ l i 1 = ( a i 2 − l i 1 u 12 ) / u 22 , i = 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 ai2​=li1​u12​+li2​+u22​⇒li1​=(ai2​−li1​u12​)/u22​,i=2,3,...,n
3、依次类推比较等式两边的第k行，可得
u k j j = a k j − ( l k 1 u 1 j + ⋯ + l k , k − 1 u k − 1 , j ) , j = k , k + 1 , . . . , n u_{kjj}=a_{kj}-(l_{k1}u_{1j}+cdots+l_{k,k-1}u_{k-1,j}), j=k,k+1,...,n ukjj​=akj​−(lk1​u1j​+⋯+lk,k−1​uk−1,j​),j=k,k+1,...,n
依次类推比较等式两边的第k列，可得
l i k = ( a i k − l i 1 u 1 k − ⋯ − l i , k − 1 u k − 1 , k ) / u k k , i = k , k + 1 , . . . , n l_{ik}=(a_{ik}-l_{i1}u_{1k}-cdots-l_{i,k-1}u_{k-1,k})/u_{kk}, i=k,k+1,...,n lik​=(aik​−li1​u1k​−⋯−li,k−1​uk−1,k​)/ukk​,i=k,k+1,...,n
直到第n步算出L和U的所有元素。

u 1 i = a 1 i ( i = 1 , 2 , 3 , ⋯ , n ) u_{1i}=a_{1i}(i=1,2,3,⋯,n) u1i​=a1i​(i=1,2,3,⋯,n)
l i 1 = a i 1 / u 11 ( i = 2 , 3 , ⋯ , n ) l_{i1}=a_{i1}/u_{11}(i=2,3,⋯,n) li1​=ai1​/u11​(i=2,3,⋯,n)
u r i = a r i − ∑ k = 1 r − 1 l r k u k i ( i = r , r + 1 , r + 2 , ⋯ , n ) u_{ri}=a_{ri}−sum_{k=1} ^{r−1}l_{rk}u_{ki}(i=r,r+1,r+2,⋯,n) uri​=ari​−k=1∑r−1​lrk​uki​(i=r,r+1,r+2,⋯,n)
l i r = ( a i r − ∑ k = 1 r − 1 l i k u k r ) / u r r ( i = r + 1 , ⋯ , n ; 且 r ≠ n ) l_{ir}=(a_{ir}−sum_{k=1}^{r−1}l_{ik}u_{kr})/u_{rr}(i=r+1,⋯,n;且r≠n) lir​=(air​−k=1∑r−1​lik​ukr​)/urr​(i=r+1,⋯,n;且r=n)
C++实现方面，首先LU分解的函数传入两个参数，方阵的一阶数组和方阵的阶数（方阵用一维数组的行优先表示）。

计算步骤：

1. 初始化LU矩阵，L矩阵上三角为0，对角线为1，U矩阵下三角为0.

2. 计算U矩阵的第一行和L矩阵的第一列

3. 循环计算U矩阵和L矩阵，第一层循环表示U计算到第几行，同时也表示L计算到第几列，先计算U后计算L。第二层循环分别表示U矩阵改行的第几列元素。第三层循环就是公式中的求和符号部分。
下面展示一些 内联代码片。

#include 
// 参数：一个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);
}

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

矩阵A可做如下分解 A = L U A=LU A=LU
则 A x = b ⇔ L U x = b Ax=bLeftrightarrow LUx=b Ax=b⇔LUx=b
令
U x = y Ux=y Ux=y
则
A x = b ⇔ L y = b ⇒ y = L − 1 b U x = y ⇒ x = U − 1 y Ax=bLeftrightarrow Ly=bRightarrow y=L^{-1}b Ux=yRightarrow x=U^{-1}y Ax=b⇔Ly=b⇒y=L−1bUx=y⇒x=U−1y

算法总共分为两步：

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