Python无第三方库实现Doolittle LU分解出现除零错误如何解决

错误直接原因 #

你运行时代码除以0的根本原因是分解到第二行（索引i=1）时，上三角矩阵U的对角元upper[1][1]计算结果为0，后续计算下三角矩阵L的对应列时需要除以该对角元，从而触发报错。

问题详细分析 #

1 理论层面的原生Doolittle分解限制 #

不带选主元的原生Doolittle LU分解有严格的适用前提：待分解方阵的所有顺序主子式都不为0，否则分解过程中必然会出现对角主元为0的情况，无法继续计算。
你使用的测试矩阵二阶顺序主子式为[[4,2],[8,4]]，行列式为4*4 - 2*8 = 0，完全不满足原生Doolittle分解的前提条件，所以必然会出现零主元报错，这是核心问题。

2 代码实现的额外错误 #

你在计算lower[k][i]时强制增加了int()转换，这一操作完全不符合算法要求：

• 如果计算结果不是整数，强制转int会丢失精度，导致后续U矩阵计算出现误差，甚至本来不为0的主元被错误计算为0
• 就算输入矩阵元素全为整数，LU分解的结果也不一定是整数，强制类型转换没有任何合理性

修复方案 #

方案1：仅修复代码逻辑错误，保留原生Doolittle逻辑 #

把int()转换去掉，改成浮点计算即可，该方案仍仅支持顺序主子式全为非0的矩阵，你当前的测试矩阵仍会触发零主元报错：

# 将lower[k][i]计算行修改为
lower[k][i] = (mat[k][i] - soma) / upper[i][i]

方案2：增加部分选主元逻辑，兼容一般非奇异方阵 #

如果需要处理存在零主元的普通非奇异方阵，需要在每轮分解前增加选主元、交换行的逻辑，避免出现零主元，修改后的完整代码如下：

def luDecomposition(mat, n):
    # 复制原矩阵避免修改输入
    mat = [row.copy() for row in mat]
    lower = [[0.0 for _ in range(n)] for _ in range(n)]
    upper = [[0.0 for _ in range(n)] for _ in range(n)]
    # 记录行交换的置换数组
    perm = list(range(n))

    for i in range(n):
        # 部分选主元：找i列下方绝对值最大的行
        max_row = i
        for k in range(i, n):
            if abs(mat[k][i]) > abs(mat[max_row][i]):
                max_row = k
        # 交换当前行和最大主元行
        mat[i], mat[max_row] = mat[max_row], mat[i]
        perm[i], perm[max_row] = perm[max_row], perm[i]
        # 同步交换L矩阵前i列的对应行
        for j in range(i):
            lower[i][j], lower[max_row][j] = lower[max_row][j], lower[i][j]

        # 计算U的第i行
        for k in range(i, n):
            soma = 0.0
            for j in range(i):
                soma += lower[i][j] * upper[j][k]
            upper[i][k] = mat[i][k] - soma

        # 计算L的第i列
        lower[i][i] = 1.0
        for k in range(i+1, n):
            soma = 0.0
            for j in range(i):
                soma += lower[k][j] * upper[j][i]
            lower[k][i] = (mat[k][i] - soma) / upper[i][i]

    print("     L    \t\t          U    ")
    for i in range(n):
        for j in range(n):
            print(f"{lower[i][j]:.4f}", end="\t")
        print("", end="\t")
        for j in range(n):
            print(f"{upper[i][j]:.4f}", end="\t")
        print("")
    print("置换顺序（原矩阵行索引）：", perm)

# 驱动测试代码
mat = [[4, 2, -1, 1],
       [8, 4, 2, -2],
       [1, 4, 2, -2],
       [2, 1, -2, 6]]
luDecomposition(mat, 4)

内容的提问来源于stack exchange，提问作者Bryan L. Andrade