解线性方程组的直接解法+代码实现

import numpy as np
"""
Gauss消元法: 基本的矩阵消元方法，可能对病态矩阵不稳定。
列主元消元法: 通过选择最大主元提高数值稳定性。
Doolittle分解: LU分解的一种形式，L对角线元素为1。
Crout分解: LU分解的另一种形式，U对角线元素为1
追赶法: 专为三对角矩阵设计的高效算法。
Cholesky分解: 用于对称正定矩阵，得到A=LL^T。
改进Cholesky分解: LDL^T形式，避免开平方运算
"""

# 1. Gauss消元法（基本高斯消元）
def gauss_elimination(A, b):
    n = len(b)
    # 增广矩阵
    Ab = np.hstack([A, b.reshape(-1, 1)])

    # 消元过程
    for i in range(n - 1):
        for j in range(i + 1, n):
            factor = Ab[j, i] / Ab[i, i]
            Ab[j, i:] = Ab[j, i:] - factor * Ab[i, i:]

    # 回代过程
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (Ab[i, -1] - np.sum(Ab[i, i + 1:n] * x[i + 1:])) / Ab[i, i]

    return x


# 2. 列主元消元法
def pivot_gauss_elimination(A, b):
    n = len(b)
    Ab = np.hstack([A, b.reshape(-1, 1)])

    for i in range(n - 1):
        # 找主元
        max_idx = np.argmax(abs(Ab[i:, i])) + i
        if max_idx != i:
            Ab[[i, max_idx]] = Ab[[max_idx, i]]

        # 消元
        for j in range(i + 1, n):
            factor = Ab[j, i] / Ab[i, i]
            Ab[j, i:] = Ab[j, i:] - factor * Ab[i, i:]

    # 回代
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (Ab[i, -1] - np.sum(Ab[i, i + 1:n] * x[i + 1:])) / Ab[i, i]

    return x


# 3. Doolittle分解 (LU分解的一种形式)
def doolittle_decomposition(A):
    n = A.shape[0]
    L = np.zeros((n, n))
    U = np.zeros((n, n))

    for i in range(n):
        L[i, i] = 1
        # 计算U的上三角部分
        for j in range(i, n):
            U[i, j] = A[i, j] - np.sum(L[i, :i] * U[:i, j])
        # 计算L的下三角部分
        for j in range(i + 1, n):
            L[j, i] = (A[j, i] - np.sum(L[j, :i] * U[:i, i])) / U[i, i]

    return L, U


# 4. Crout分解
def crout_decomposition(A):
    n = A.shape[0]
    L = np.zeros((n, n))
    U = np.zeros((n, n))

    for j in range(n):
        U[j, j] = 1
        # 计算L的下三角部分
        for i in range(j, n):
            L[i, j] = A[i, j] - np.sum(L[i, :j] * U[:j, j])
        # 计算U的上三角部分
        for i in range(j + 1, n):
            U[j, i] = (A[j, i] - np.sum(L[j, :j] * U[:j, i])) / L[j, j]

    return L, U


# 5. 追赶法（用于三对角矩阵）
def thomas_algorithm(a, b, c, d):
    # a: 下对角线
    # b: 主对角线
    # c: 上对角线
    # d: 右端项
    n = len(b)
    c_prime = np.zeros(n - 1)
    d_prime = np.zeros(n)

    # 前向消元
    c_prime[0] = c[0] / b[0]
    d_prime[0] = d[0] / b[0]

    for i in range(1, n - 1):
        temp = b[i] - a[i - 1] * c_prime[i - 1]
        c_prime[i] = c[i] / temp
        d_prime[i] = (d[i] - a[i - 1] * d_prime[i - 1]) / temp

    d_prime[n - 1] = (d[n - 1] - a[n - 2] * d_prime[n - 2]) / (b[n - 1] - a[n - 2] * c_prime[n - 2])

    # 回代
    x = np.zeros(n)
    x[n - 1] = d_prime[n - 1]
    for i in range(n - 2, -1, -1):
        x[i] = d_prime[i] - c_prime[i] * x[i + 1]

    return x


# 6. Cholesky分解（用于对称正定矩阵）
def cholesky_decomposition(A):
    n = A.shape[0]
    L = np.zeros((n, n))

    for i in range(n):
        # 对角元素
        L[i, i] = np.sqrt(A[i, i] - np.sum(L[i, :i] ** 2))
        # 下三角部分
        for j in range(i + 1, n):
            L[j, i] = (A[j, i] - np.sum(L[j, :i] * L[i, :i])) / L[i, i]

    return L
# Cholesky分解 解方程组
def cholesky_solve01(A, b):
    n = A.shape[0]
    L = cholesky_decomposition(A)
    y = np.linalg.solve(L, b)
    x = np.linalg.solve(L.T, y)
    return x


# 使用Cholesky分解求解Ax=b
def solve_cholesky02(A, b):
    # 确保A是对称矩阵
    if not np.allclose(A, A.T):
        raise ValueError("矩阵A必须是对称的")

    # 确保A是正定的（简单检查）
    if np.any(np.linalg.eigvals(A) <= 0):
        raise ValueError("矩阵A必须是正定的")

    # 进行Cholesky分解 A = LL^T
    L = cholesky_decomposition(A)

    n = A.shape[0]
    # 第一步：求解Ly = b（前向代入）
    y = np.zeros(n)
    for i in range(n):
        y[i] = (b[i] - np.sum(L[i, :i] * y[:i])) / L[i, i]

    # 第二步：求解L^Tx = y（回向代入）
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - np.sum(L[i + 1:, i] * x[i + 1:])) / L[i, i]

    return x


# 7. 改进Cholesky分解（LDL^T分解）
def modified_cholesky_decomposition(A):
    n = A.shape[0]
    L = np.zeros((n, n))
    D = np.zeros(n)

    for i in range(n):
        L[i, i] = 1
        # 计算D的对角元素
        D[i] = A[i, i] - np.sum(L[i, :i] ** 2 * D[:i])
        # 计算L的下三角部分
        for j in range(i + 1, n):
            L[j, i] = (A[j, i] - np.sum(L[j, :i] * L[i, :i] * D[:i])) / D[i]

    return L, D


if __name__ == "__main__":
    # # 测试矩阵
    # A = np.array([[4, 12, -16],
    #               [12, 37, -43],
    #               [-16, -43, 98]], dtype=float)
    # b = np.array([1, 2, 3])
    #
    # # Gauss消元
    # print("Gauss消元结果:", gauss_elimination(A.copy(), b.copy()))
    #
    # # 列主元
    # print("列主元消元结果:", pivot_gauss_elimination(A.copy(), b.copy()))
    #
    # # Doolittle
    # L, U = doolittle_decomposition(A.copy())
    # print("Doolittle分解:\nL=\n", L, "\nU=\n", U)
    #
    # # Crout
    # L, U = crout_decomposition(A.copy())
    # print("Crout分解:\nL=\n", L, "\nU=\n", U)
    #






    # # Cholesky
    # Cholesky_A = np.array([[3.1, 1.5,1.0],
    #               [1.5, 2.5,0.5],
    #               [1.0, 0.5, 4.2]], dtype=float)
    # Cholesky_b=np.array([10.83,9.2,17.1])
    # L = cholesky_decomposition(Cholesky_A.copy())
    # x=solve_cholesky02(Cholesky_A.copy(),Cholesky_b.copy())
    #
    # print("Cholesky_A分解:\nL=\n", L)
    # print("解x:\nx=\n", x)
    # # # 验证
    # print("验证Ax:", Cholesky_A @ x)
    # print("原b:", Cholesky_b)
    #
    # # 与NumPy内置函数比较
    # x_numpy = np.linalg.solve(Cholesky_A, Cholesky_b)
    # print("NumPy解:", x_numpy)







    #
    # # 改进Cholesky
    # L, D = modified_cholesky_decomposition(A.copy())
    # print("改进Cholesky分解:\nL=\n", L, "\nD=\n", D)

    # 三对角矩阵
    # a: 下对角线
    # b: 主对角线
    # c: 上对角线
    # d: 右端项
    a = np.array([1,1,1,1])
    b = np.array([2, 2, 2,2,2])
    c = np.array([1, 1, 1,1])
    d = np.array([1, 0, 0,0,1])
    print("追赶法结果:", thomas_algorithm(a, b, c, d))