LU分解法

LU分解的原理是将一个矩阵分解为一个下三角矩阵L和一个上三角矩阵U的乘积。

下面是LU分解的原理：
定一个n * n的矩阵A，我们希望将它分解为下三角矩阵L和上三角矩阵U的乘积，即A = LU。 我们可以在矩阵A上进行列主元高斯消元，消元的过程中会用到矩阵的初等变换。我们把消元过程中记录下来的消元主元放到一个n * n的矩阵U的对角线上，将消元过程中得到的消元因子放到一个n * n的矩阵L的对角线以下，得到LU分解。

import numpy as np

# 定义函数进行LU分解
def LU_decomposition(A):
    n = len(A[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:
            # 求U的第i行
            for j in range(i, n):
                sum = 0
                for k in range(0, i):
                    sum += L[i][k] * U[k][j]
                U[i][j] = A[i][j]-sum
            # 求L的第i列
            for j in range(i+1, n):
                sum = 0
                for k in range(0, i):
                    sum += L[j][k] * U[k][i]
                L[j][i] = (A[j][i] - sum)/U[i][i]
    return L, U

# 定义函数进行前/回代
def solve(L, U, b):
    n = len(b)
    y = np.zeros(n)
    x = np.zeros(n)
    # 前代
    for i in range(n):
        s = sum(L[i][j] * y[j] for j in range(i))
        y[i] = b[i] - s

    # 回代
    for i in range(n - 1, -1, -1):
        s = sum(U[i][j] * x[j] for j in range(i, n))
        x[i] = (y[i] - s) / U[i][i]

    return x


# 测试
A = np.array([[1, 2, -1], [2, 1, -2], [-3, 1, 1]])
b = np.array([3, 3, -6])

# 进行LU分解
L, U = LU_decomposition(A)
# 求解线性方程组
x = solve(L, U, b)

# 输出结果
print('系数矩阵A:\n', A)
print('向量b:', b)
print('解向量x:', x)
print('库函数解的结果', np.linalg.solve(A, b))

"""
#L*U=A
for i in range(1,n+1):
    for j in range(1,n+1):
        sum = 0
        for k in range(1,n+1):
            sum+=L[i][k]*U[k][j]
        A[i][j]=sum

step1:用矩阵L的第一行,去乘矩阵U的第i列,得到矩阵U的第1行
        U[1][i]=A[1][i]         (矩阵L的第1行L[1][1]==1,L[1][j]==0)
      用矩阵L的第j行,去乘矩阵U的第1列,得到矩阵L的第1列
        L[j][1]=A[j][1]/U[1][1] (矩阵U的第1列U[1][1]=A[1][1],U[j][1]==0)
next_step:用矩阵L的第k行去乘矩阵U的第i列(i>=k),得到矩阵U的第k行
        U[k][i]=A[k][i]-sum(L[k][m]*U[m][k] for k in range(1,k))
      用矩阵L的第j行,去乘矩阵U的第k列,得到矩阵L的第k列
      L[j][k]=(A[j][k]-sum(L[j][m]*U[m][k] for m in range(1,k)))/U[k][k]
goto next_step

LU分解法先求出矩阵U的一行元素,再求出矩阵L的一列元素,就这样逐行逐列交替求得矩阵中的所有元素。
"""