Gauss-Zeidel迭代法代码返回错误结果，求问题排查方案

我编写了一段用Gauss-Seidel法求解线性方程组的Python代码，但运行结果不正确，无法定位问题根源，代码如下：

import numpy as np

A = np.matrix([[10, 3, 0], 
               [3, 15, 1], 
               [0, 1, 7]]).astype(float)
f = np.matrix([2, 12, 5]).astype(float)
f = np.reshape(f, (3, 1))
x_new = np.matrix([0, 0, 0]).astype(float)
x_new = np.reshape(x_new, (3, 1))
N = 100
D = np.matrix(np.zeros(A.shape))
for n in range(A.shape[0]):
    D[n, n] = A[n, n]
A1 = np.tril(A, -1)
A2 = np.triu(A, 1)
D_inv = np.linalg.inv(D)
x = np.matrix([0, 0, 0]).astype(float)
x = np.reshape(x, (3, 1))
eps = 1e-5
dim = len(A)

def gauss_zeidel(x):
    i = 0
    t = 0
    s1 = 0
    s2 = 0
    while t <= N:
        x_new = np.copy(x)
        for i in range(dim):
            print("--new cycle--")
            print("===============================")
            for j in range (i-1):
                s1 += A[i,j] * x_new.item(j)
            for j  in range(i+1, dim):
                s2 += A[i,j] * x.item(j)
            print("iteration result")
            x_new.itemset(i, (f.item(i) - s1 - s2)/A[i,i])
            print((x.item(i)), i)
        print(t)
        print("===============================")
        x = x_new
        t+=1
    print(x)

问题分析及修正点 #

• 累加变量未重置：s1和s2在函数开头定义后，每次循环i时没有清零，导致累加了之前所有迭代的结果，完全偏离了Gauss-Seidel法的计算逻辑。必须将s1、s2的初始化放到每个i的循环内部。

• 循环范围错误：计算s1时，range(i-1)会漏掉j=i-1这个已更新的变量，正确的范围应该是range(i)，确保包含从0到i-1的所有已更新变量。

• 缺少收敛判断：代码仅执行固定N次迭代，没有利用eps判断是否收敛，既浪费计算资源，也无法确定结果是否稳定。

• 冗余变量冗余：函数外定义的x_new、D、A1、A2、D_inv等变量在函数内完全未使用，属于无效代码，可直接删除。

修正后的代码 #

import numpy as np

A = np.matrix([[10, 3, 0], 
               [3, 15, 1], 
               [0, 1, 7]]).astype(float)
f = np.matrix([2, 12, 5]).astype(float).reshape(3, 1)
x = np.matrix([0, 0, 0]).astype(float).reshape(3, 1)
N = 100
eps = 1e-5
dim = len(A)

def gauss_zeidel(x):
    t = 0
    while t <= N:
        x_new = np.copy(x)
        max_diff = 0
        for i in range(dim):
            # 每次计算前重置累加变量
            s1 = 0
            s2 = 0
            # 累加已更新的前i个变量
            for j in range(i):
                s1 += A[i, j] * x_new.item(j)
            # 累加未更新的后续变量
            for j in range(i+1, dim):
                s2 += A[i, j] * x.item(j)
            # 更新当前变量
            new_val = (f.item(i) - s1 - s2) / A[i, i]
            x_new.itemset(i, new_val)
            # 记录最大更新误差
            diff = abs(new_val - x.item(i))
            if diff > max_diff:
                max_diff = diff
        # 判断是否达到收敛精度
        if max_diff < eps:
            print(f"迭代{t}次后收敛")
            break
        x = x_new
        t += 1
    print("最终求解结果：")
    print(x)
    return x

# 调用迭代函数
gauss_zeidel(x)

内容的提问来源于stack exchange，提问作者LevGor