1. 线性方程组直接法对比：高斯消元与LU分解

问题设定与两种直接法的实现

本节针对同一个稠密线性方程组 Ax = b，分别用高斯消元（Gauss Elimination, GE）与 LU 分解（Doolittle，带部分主元）求解，并对比计算量、数值稳定性与主元策略的影响，最后用残差 ‖b − Ax̂‖∞ 验证解的正确性。

为便于复现，取一个 4×4 的示例矩阵（教学示例数据，非实测）：

A = [[ 4, -2,  1,  3],
     [ 8, -3,  2,  7],
     [-4,  1, -3,  1],
     [ 6, -5,  4, 12]]
b = [ 7, 17, -9, 21]

两种方法的核心差异在于：高斯消元把右端项 b 与系数矩阵同步消元，一次性得到解；LU 分解只对 A 做分解得到 L、U 与置换 P，再通过前代/回代处理 b。因此同一 A 面对多个右端项时，LU 分解只需分解一次。

关键实现（Python，可直接运行）

import numpy as np

def gauss_elimination(A, b):
    A = A.astype(float).copy()
    b = b.astype(float).copy()
    n = len(b)
    swaps = 0
    for k in range(n - 1):
        # 部分主元：选列主元绝对值最大者
        p = k + np.argmax(np.abs(A[k:, k]))
        if p != k:
            A[[k, p]] = A[[p, k]]
            b[[k, p]] = b[[p, k]]
            swaps += 1
        if abs(A[k, k]) < 1e-14:
            raise ValueError("矩阵奇异或接近奇异")
        for i in range(k + 1, n):
            m = A[i, k] / A[k, k]
            A[i, k:] -= m * A[k, k:]
            b[i] -= m * b[k]
    # 回代
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (b[i] - A[i, i+1:] @ x[i+1:]) / A[i, i]
    return x, swaps

def lu_decompose(A):
    A = A.astype(float).copy()
    n = A.shape[0]
    L = np.eye(n)
    U = A.copy()
    P = np.eye(n)
    swaps = 0
    for k in range(n - 1):
        p = k + np.argmax(np.abs(U[k:, k]))
        if p != k:
            U[[k, p]] = U[[p, k]]
            P[[k, p]] = P[[p, k]]
            L[[k, p], :k] = L[[p, k], :k]   # 只交换已生成的乘子
            swaps += 1
        if abs(U[k, k]) < 1e-14:
            raise ValueError("矩阵奇异或接近奇异")
        for i in range(k + 1, n):
            L[i, k] = U[i, k] / U[k, k]
            U[i, k:] -= L[i, k] * U[k, k:]
    return P, L, U, swaps

def lu_solve(P, L, U, b):
    b = P @ b.astype(float)
    n = len(b)
    y = np.zeros(n)
    for i in range(n):
        y[i] = b[i] - L[i, :i] @ y[:i]
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - U[i, i+1:] @ x[i+1:]) / U[i, i]
    return x

A = np.array([[ 4, -2,  1,  3],
              [ 8, -3,  2,  7],
              [-4,  1, -3,  1],
              [ 6, -5,  4, 12]], dtype=float)
b = np.array([7, 17, -9, 21], dtype=float)

x_ge, sw_ge = gauss_elimination(A, b)
P, L, U, sw_lu = lu_decompose(A)
x_lu = lu_solve(P, L, U, b)

print("GE 解:", x_ge, " 交换次数:", sw_ge)
print("LU 解:", x_lu, " 交换次数:", sw_lu)
print("GE 残差:", np.linalg.norm(b - A @ x_ge, np.inf))
print("LU 残差:", np.linalg.norm(b - A @ x_lu, np.inf))

预期输出（示例数据可直接推出）：两种方法给出同一解 x ≈ [1, 2, -1, 1]，残差均在 1e-14 量级（浮点舍入范围内）。若残差明显偏大，说明主元选择或回代实现有误。

计算量与主元策略对比

结果判读与选型建议

• 解一致性：GE 与 LU 得到的 x̂ 应在浮点舍入范围内一致；若差异明显，优先检查 LU 中交换行时是否同步交换了 L 已生成的乘子。
• 残差验证：用 ‖b − A x̂‖∞ 判断，示例数据下应接近机器精度；残差大通常意味着主元过小导致放大误差，或回代索引写错。
• 计算量：单右端项时两者同阶（≈ n³/3）；多右端项时 LU 分解一次即可复用，优势随右端项数量线性放大。
• 稳定性：部分主元使乘子 |l_ik| ≤ 1，避免小主元放大误差；对病态矩阵仍需结合条件数评估。