Jacobi与Gauss-Seidel迭代求解器溢出错误排查求助

你的迭代求解器出现溢出错误，核心原因是给定的矩阵不满足Jacobi/Gauss-Seidel迭代的收敛条件，导致迭代过程中数值持续暴涨，最终触发溢出。同时代码存在一些实现细节问题，加剧了错误。

1. 收敛性根本问题 #

Jacobi和Gauss-Seidel迭代收敛的充分条件是矩阵满足严格对角占优（即每行对角元素的绝对值大于该行其他元素绝对值之和），或者矩阵对称正定。

以你的矩阵A1为例：

• 第一行：对角元素绝对值|2|=2，其他元素绝对值和4+2+2=8 → 2 < 8，不满足
• 第二行：对角元素绝对值|2|=2，其他元素绝对值和1+4+3=8 → 2 < 8，不满足
• 第三行：对角元素绝对值|8|=8，其他元素绝对值和3+3+2=8 → 仅满足弱对角占优，不满足严格条件
• 第四行：对角元素绝对值|-3|=3，其他元素绝对值和1+1+6=8 → 3 < 8，不满足

由于矩阵不满足收敛条件，迭代过程中解向量的元素会无限增大，最终导致数值溢出。

2. 代码实现细节修复 #

除了收敛性问题，代码中的一些细节也需要调整：

Jacobi迭代修复 #

• 初始时直接赋值x0=x0_in会导致引用原输入数组，需改为np.copy避免副作用
• 手动计算距离可以替换为np.linalg.norm，更高效且减少溢出风险

修复后的代码：

import numpy as np

def Jacobi(A_in, x0_in, b_in, tol, step_max):
    N = len(A_in)
    x0 = np.copy(x0_in)  # 修复：复制输入向量，避免修改原数据
    b = np.copy(b_in)
    A = A_in.copy()
    it = 0
    x1 = np.copy(x0)
    while True:
        for i in range(N):
            S = 0.0
            for j in range(N):
                if i != j:
                    S += A[i][j] * x0[j]
            x1[i] = (b[i] - S) / A[i][i]
        
        # 用numpy内置范数计算距离，替代手动求和
        dist_x0_x1 = np.linalg.norm(x0 - x1)
        print(f"Iteration {it}: {x1}")
        
        if dist_x0_x1 < tol:
            return x1
        x0 = np.copy(x1)
        it += 1
        if it > step_max:
            print("Diverging")
            break

# 测试：先求精确解对比
A1 = np.array([[2,4,-2,-2],[1,2,4,-3],[-3,-3,8,-2],[-1,1,6,-3]], dtype=np.float64)
b1 = np.array([4,-3,5,1], dtype=np.float64)
exact_sol = np.linalg.solve(A1, b1)
print(f"Exact solution: {exact_sol}")

# 调用Jacobi（注意：由于矩阵不收敛，仍会发散，但代码无语法/复制问题）
Jacobi(A1, np.zeros(4), b1, 1e-5, 20)

Gauss-Seidel迭代修复 #

• 替换手动距离计算为np.linalg.norm，避免手动求和的溢出风险
• 使用更高精度的float64数据类型，避免float32的精度溢出

修复后的代码：

import numpy as np

def GaussSeidel(A_in, x0_in, b_in, tol=1E-5, max_step=100):
    b = np.copy(b_in)
    AC = A_in.copy()
    N = len(A_in)
    it = 0
    x_g = np.copy(x0_in)
    x0 = np.copy(x_g)  # 初始化当前迭代向量
    while True:
        for i in range(N):
            S = 0.0
            for j in range(N):
                if i == j:
                    continue
                S += AC[i][j] * x0[j]
            x0[i] = (b[i] - S) / AC[i][i]
        
        # 用numpy范数计算距离
        dist = np.linalg.norm(x0 - x_g)
        
        if dist < tol:            
            return x0
        x_g = np.copy(x0)
        it += 1
        if it > max_step:
            print("Divergent")
            break

# 测试
A = np.array([[2.0,4.0,-2.0,-2.0],[1.0,2.0,4.0,-3.0],[-3.0,-3.0,8.0,-2.0],[-1.0,1.0,6.0,-3.0]], dtype=np.float64)
b1 = np.array([4.0,-3.0,5.0,1.0], dtype=np.float64)
exact_sol = np.linalg.solve(A, b1)
print(f"Exact solution: {exact_sol}")

# 调用Gauss-Seidel（同样因矩阵不收敛会发散）
GaussSeidel(A, np.zeros(4), b1)

3. 解决收敛问题的建议 #

如果必须使用迭代法求解该方程组，可尝试以下方案：

• 改用收敛的迭代方法：比如GMRES等适用于非对称矩阵的迭代法
• 调整矩阵形式：尝试对原方程组进行行变换，使其满足严格对角占优（比如交换行顺序，或缩放行）
• 使用直接法验证：先用np.linalg.solve求解精确解，确认方程组有唯一解，再尝试迭代法逼近

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