Cholesky分解法-python

import math

def SolverEqCholesky(matrix_A, vector_b, n):
    # 创建一个n*1的零向量作为结果
    x = [[0.] for _ in range(n)]

    # Cholesky分解
    L = [[0.] * n for _ in range(n)]
    for i in range(n):
        for j in range(i+1):
            if i == j:
                sum_val = 0.
                for k in range(j):
                    sum_val += L[j][k] ** 2
                L[j][j] = math.sqrt(matrix_A[j][j] - sum_val)
            else:
                sum_val = 0.
                for k in range(j):
                    sum_val += L[i][k] * L[j][k]
                L[i][j] = (matrix_A[i][j] - sum_val) / L[j][j]

    # 前向替换 Ly = b
    y = [[0.] for _ in range(n)]
    for i in range(n):
        y[i] = vector_b[i]
        for j in range(i):
            y[i] -= L[i][j] * y[j]
        y[i] /= L[i][i]

    # 后向替换 L^T x = y
    for i in range(n-1, -1, -1):
        x[i] = y[i]
        for j in range(i+1, n):
            x[i] -= L[j][i] * x[j]
        x[i] /= L[i][i]

    return x