计算机数值方法——LU分解法（C++\Python代码实现）

算法流程

首先，LU分解法就是在高斯消元法的基础上，把矩阵$A$分解为一个上三角矩阵$U$与一个单位下三角矩阵$L$的乘积。

懒得敲LaTeX公式了，书上由具体的推导过程，这里我们重点介绍代码吧

主要就是在高斯消元的过程中标记单位下三角矩阵$L$，算法复杂度$O(N^3)$,没有变化。

C++代码

#include <bits/stdc++.h>
using namespace std;
typedef pair<int, int> PII;
#define int long long
const int N = 1e2 + 10;
int a[N][N], b[N][N];
int n;
void LU_Factorization()
{
    for (int i = 0; i < n; i++)
        b[i][i] = 1;

    for (int k = 0; k < n; k++)
        for (int i = k + 1; i < n; i++)
        {
            double t = a[i][k] / a[k][k];
            b[i][k] = t;
            for (int j = 0; j < n; j++)
            {
                a[i][j] -= a[k][j] * t;
            }
        }
    printf("L矩阵为：\n");
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
            printf("%d ", b[i][j]);
        puts("");
    }
    printf("U矩阵为：\n");
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
            printf("%d ", a[i][j]);
        puts("");
    }
}
signed main()
{
    cin >> n;
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            cin >> a[i][j];
    LU_Factorization();
    return 0;
}

python代码

老师同意调库了，，试一下numpy（虽然本质上没改变什么，但是在学习中成长）

# LU分解
import numpy as np
def LU_Factorization(A):
    n = A.shape[0]
    L = np.zeros([n, n])
    U = np.zeros([n, n])
    for i in range(n):
        L[i][i] = 1
        if i == 0:
            U[0][0] = A[0][0]
            for j in range(1, n):
                U[0][j] = A[0][j]
                L[j][0] = A[j][0] / U[0][0]
        else:
            for j in range(i, n):
                temp = 0
                for k in range(0, i):
                    temp = temp + L[i][k] * U[k][j]
                U[i][j] = A[i][j] - temp
            for j in range(i + 1, n):
                temp = 0
                for k in range(0, i):
                    temp = temp + L[j][k] * U[k][i]
                L[j][i] = (A[j][i] - temp) / U[i][i]
    return L, U
def main():
    A = np.array([[1, 1, 1], [0, 4, -1], [2, -2, -1]])
    L, U = LU_Factorization(A)
    print('L矩阵为：\n',L, '\nU矩阵为：\n', U)
main()