利用矩阵的LU分解求解线性方程组

线性方程组可以表示为Ax=b的形式，其中A是一个n×n的方阵，x和b是n维向量。LU分解是将矩阵A分解为一个下三角矩阵L和一个上三角矩阵U的乘积，即A=LU。LU分解的算法可以使用高斯-约旦消元法或Doolittle分解法等方法实现。

假设已经完成了LU分解，那么可以通过以下步骤求解线性方程组：

1. 求解Ly=b，其中y是一个n维向量。这可以通过前代法（forward substitution）实现，即从上到下依次求解。

2. 求解Ux=y，其中x是一个n维向量。这可以通过回代法（back substitution）实现，即从下到上依次求解。

下面给出一个Python实现代码：

import numpy as np

def lu_solve(A, b):
    """
    利用LU分解求解线性方程组Ax=b
    输入：
        A: n×n的矩阵
        b: n维向量
    输出：
        x: n维向量，满足Ax=b
    """
    n = A.shape[0]
    # LU分解
    L = np.zeros((n, n))
    U = np.zeros((n, n))
    for i in range(n):
        L[i, i] = 1.0   # L的对角线元素为1
        for j in range(i, n):
            U[i, j] = A[i, j] - np.dot(L[i, :i], U[:i, j])
        for j in range(i+1, n):
            L[j, i] = (A[j, i] - np.dot(L[j, :i], U[:i, i])) / U[i, i]
    # 前代求解Ly=b
    y = np.zeros(n)
    for i in range(n):
        y[i] = b[i] - np.dot(L[i, :i], y[:i])
    # 回代求解Ux=y
    x = np.zeros(n)
    for i in range(n-1, -1, -1):
        x[i] = (y[i] - np.dot(U[i, i+1:], x[i+1:])) / U[i, i]
    return x

使用这个函数可以很方便地求解线性方程组，例如：

A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 10]])
b = np.array([6, 15, 29])
x = lu_solve(A, b)
print(x)

输出结果为：

[-1.  2.  1.]

这正是线性方程组的解。