改进的平方根法

import numpy as np
from fractions import Fraction


def modified_cholesky_decomposition(A):
    """
    对对称矩阵 A 进行改进平方根分解 A = LDL^T
    参数：
    A : 对称矩阵 (n x n)
    返回：
    L : 下三角矩阵 (对角线为1)
    D : 对角矩阵
    """
    n = A.shape[0]
    L = np.zeros((n, n), dtype=object)  # 使用 object 类型存储 Fraction
    D = np.zeros(n, dtype=object)  # 对角矩阵 D

    # 初始化 L 的对角线为 1
    for i in range(n):
        L[i, i] = Fraction(1)

    for i in range(n):
        # 计算 D[i]
        sum_d = sum(D[k] * L[i, k] * L[i, k] for k in range(i))
        D[i] = Fraction(A[i, i]) - sum_d

        if abs(float(D[i])) < 1e-10:  # 检查是否接近零
            raise ValueError("矩阵不可分解或需要进一步调整")

        # 计算 L 的第 i 列
        for j in range(i + 1, n):
            sum_l = sum(D[k] * L[j, k] * L[i, k] for k in range(i))
            L[j, i] = (Fraction(A[j, i]) - sum_l) / D[i]

    return L, D


def solve_with_modified_cholesky(A, b):
    """
    使用改进平方根分解求解 Ax = b
    参数：
    A : 对称矩阵 (n x n)
    b : 常数向量 (n x 1)
    返回：
    x : 解向量
    """
    L, D = modified_cholesky_decomposition(A)

    n = len(b)
    # 前向代入：Ly = b
    y = np.zeros(n, dtype=object)
    for i in range(n):
        y[i] = Fraction(b[i]) - sum(L[i, j] * y[j] for j in range(i))

    # 中间步骤：Dz = y
    z = np.zeros(n, dtype=object)
    for i in range(n):
        z[i] = y[i] / D[i]

    # 回向代入：L^T x = z
    x = np.zeros(n, dtype=object)
    for i in range(n - 1, -1, -1):
        x[i] = z[i] - sum(L[j, i] * x[j] for j in range(i + 1, n))

    return x



if __name__ == "__main__":
    # A = np.array([[1, 2, 1],
    #               [2, 5, 0],
    #               [1, 0, 14]], dtype=float)
    #
    # b = np.array([4, 7, 15], dtype=float)
    A = np.array([[2, -1, 1],
                  [-1, -2, 3],
                  [1, 3, 1]], dtype=float)

    b = np.array([4, 5, 6], dtype=float)

    try:
        # 改进平方根分解
        L, D = modified_cholesky_decomposition(A)
        print("L 矩阵（分数形式）:")
        for row in L:
            print([str(f) for f in row])
        print("\nD 矩阵（分数形式）:")
        print([str(d) for d in D])

        # 验证分解
        L_float = L.astype(float)
        D_float = np.diag(D.astype(float))
        print("\n验证 L * D * L^T:")
        print(np.dot(np.dot(L_float, D_float), L_float.T))
        print("原始 A:")
        print(A)

        # 求解
        x = solve_with_modified_cholesky(A, b)

        # 输出结果
        print("\n解（分数形式）:")
        print(f"x = {x[0]}")
        print(f"y = {x[1]}")
        print(f"z = {x[2]}")

        print("\n解（浮点数形式）:")
        print(f"x = {float(x[0]):.4f}")
        print(f"y = {float(x[1]):.4f}")
        print(f"z = {float(x[2]):.4f}")

        # 验证
        print("\n验证结果:")
        print("Ax =", np.dot(A, x.astype(float)))
        print("b =", b)

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