追赶法+thomas算法

import numpy as np
from fractions import Fraction

def thomas_algorithm(a, b, c, d):
    """
    使用追赶法求解三对角线性方程组 Ax = d
    参数：
    a : 下对角线元素 (长度 n-1)
    b : 主对角线元素 (长度 n)
    c : 上对角线元素 (长度 n-1)
    d : 右侧常数向量 (长度 n)
    返回：
    x : 解向量
    """
    n = len(b)
    # 使用 object 类型存储 Fraction
    c_prime = np.zeros(n - 1, dtype=object)  # 修改后的上对角线
    d_prime = np.zeros(n, dtype=object)  # 修改后的右侧向量
    x = np.zeros(n, dtype=object)  # 解向量

    # 将输入转换为 Fraction
    a = [Fraction(ai) for ai in a]
    b = [Fraction(bi) for bi in b]
    c = [Fraction(ci) for ci in c]
    d = [Fraction(di) for di in d]

    # 前向追赶
    c_prime[0] = c[0] / b[0]
    d_prime[0] = d[0] / b[0]

    for i in range(1, n - 1):
        denom = b[i] - a[i - 1] * c_prime[i - 1]
        if abs(float(denom)) < 1e-10:
            raise ValueError("矩阵不可分解或需要调整")
        c_prime[i] = c[i] / denom
        d_prime[i] = (d[i] - a[i - 1] * d_prime[i - 1]) / denom

    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[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

def lu_decomposition(A):
    """
    对矩阵 A 进行 LU 分解
    参数：
    A : 矩阵 (n x n)
    返回：
    L : 下三角矩阵
    U : 上三角矩阵
    """
    n = A.shape[0]
    L = np.zeros((n, n), dtype=object)
    U = np.zeros((n, n), dtype=object)

    for i in range(n):
        # 上三角矩阵 U
        for j in range(i, n):
            sum0 = sum(L[i][k] * U[k][j] for k in range(i))
            U[i][j] = A[i][j] - sum0

        # 下三角矩阵 L
        for j in range(i, n):
            if i == j:
                L[i][i] = Fraction(1)
            else:
                sum0 = sum(L[j][k] * U[k][i] for k in range(i))
                L[j][i] = (A[j][i] - sum0) / U[i][i]

    return L, U

if __name__ == "__main__":
    a = [-1, -1, -1, -1]  # 下对角线
    b = [2,2,2,2,2]  # 主对角线
    c = [-1,-1,-1,-1]  # 上对角线
    d = [1,0,0,0,0]  # 右侧常数

    try:
        # 求解
        x = thomas_algorithm(a, b, c, d)

        # 输出结果
        print("解（分数形式）:")
        print(f"x1 = {x[0]}")
        print(f"x2 = {x[1]}")
        print(f"x3 = {x[2]}")
        print(f"x4 = {x[3]}")
        print(f"x5 = {x[4]}")

        print("\n解（浮点数形式）:")
        print(f"x1 = {float(x[0]):.4f}")
        print(f"x2 = {float(x[1]):.4f}")
        print(f"x3 = {float(x[2]):.4f}")
        print(f"x4 = {float(x[3]):.4f}")
        print(f"x5 = {float(x[4]):.4f}")

        # 构建完整的矩阵 A
        n = len(b)
        A = np.zeros((n, n), dtype=object)
        for i in range(n):
            A[i][i] = b[i]
            if i > 0:
                A[i][i - 1] = a[i - 1]
            if i < n - 1:
                A[i][i + 1] = c[i]

        # 进行 LU 分解
        L, U = lu_decomposition(A)

        # 输出 L 和 U 矩阵
        print("\nL 矩阵（分数形式）:")
        for row in L:
            print([str(f) for f in row])

        print("\nU 矩阵（分数形式）:")
        for row in U:
            print([str(f) for f in row])

        # 验证 LU 分解
        L_float = L.astype(float)
        U_float = U.astype(float)
        A_float = np.dot(L_float, U_float)
        print("\n验证 L * U:")
        print(A_float)
        print("原始 A:")
        print(A.astype(float))

        # 验证解
        d_array = np.array(d, dtype=float)
        x_float = x.astype(float)
        print("\n验证结果:")
        print("Ax =", np.dot(A_float, x_float))
        print("d =", d_array)

    except ValueError as e:
        print(f"错误: {e}")