低松弛法求解线性代数方程组的代码问题求助

问题根源分析 #

你的代码在处理4阶及以上矩阵时失效，核心问题有三个：

1. 迭代逻辑错误：误用同步更新而非异步更新 #

当前代码采用同步更新（类似Jacobi迭代），所有phi[i]的更新依赖上一轮的完整数组。但低松弛法是Gauss-Seidel迭代的变体，需要异步更新——计算第i个变量时，要使用已更新的phi[0]到phi[i-1]，以及未更新的phi[i+1]到phi[n-1]。同步更新会大幅降低收敛速度，甚至在非对角占优矩阵下直接发散。

2. 测试矩阵不满足迭代法收敛条件 #

迭代法（包括低松弛、Gauss-Seidel）收敛的关键前提是矩阵严格对角占优（每行对角线元素的绝对值大于该行其他元素绝对值之和）。你的测试矩阵A：

• 第一行：|5|=5 < |8|+|-3|+|21|=32，不满足
• 第三行：|3|=3 < |-3.6|+|0.34|+|-76|=79.94，不满足
• 第四行：|9|=9 < |-5.7|+|2.26|+|-5.6|=13.56，不满足

这种非对角占优的矩阵，即便迭代逻辑正确也大概率无法收敛。

3. 松弛因子与初始猜测的选择问题 #

欠松弛法的omega必须在(0,1)区间内，取值不当（如接近0或1）会导致收敛极慢；初始猜测若离真实解过远，也会加剧发散风险。

修复后的代码 #

针对上述问题，修改迭代逻辑为异步更新，并添加收敛条件检查提示：

import numpy as np

def relax_solver(A, b, omega, initial_guess, convergence_criteria, max_iter, print_iter):
    """
    用低松弛法（欠松弛，omega ∈ (0,1)）求解线性方程组，基于Gauss-Seidel异步迭代
    输入：
        A: n阶numpy矩阵
        b: n维numpy向量
        omega: 松弛因子，需满足 0 < omega < 1
        initial_guess: 初始解猜测
        max_iter: 最大迭代次数
        print_iter: 是否打印每轮迭代信息
    返回：
        收敛后的解向量，迭代次数
    """
    n = A.shape[0]
    phi = np.copy(initial_guess)  # 避免修改原始输入
    residual = np.linalg.norm(A @ phi - b)
    iteration = 0

    # 检查矩阵是否严格对角占优，给出提示
    diag_dominant = True
    for i in range(n):
        diag_abs = abs(A[i,i])
        row_sum = sum(abs(A[i,j]) for j in range(n) if j != i)
        if diag_abs <= row_sum:
            diag_dominant = False
            print(f"警告：第{i+1}行不满足严格对角占优，迭代可能不收敛")
    
    while residual > convergence_criteria and iteration < max_iter:
        # 异步更新：逐个修改phi，后续计算使用已更新的值
        for i in range(n):
            sigma = 0
            # 前半部分：已更新的phi[0..i-1]
            for j in range(i):
                sigma += A[i,j] * phi[j]
            # 后半部分：未更新的phi[i+1..n-1]
            for j in range(i+1, n):
                sigma += A[i,j] * phi[j]
            # 低松弛更新公式
            phi[i] = (1 - omega) * phi[i] + (omega / A[i,i]) * (b[i] - sigma)
        
        residual = np.linalg.norm(A @ phi - b)
        iteration += 1
        if print_iter:
            print(f'迭代 {iteration}, 残差: {residual:.6f}, 解: {", ".join(f"{val:.6f}" for val in phi)}')
    
    if iteration >= max_iter:
        print("警告：已达最大迭代次数，未收敛")
    print("\n最终解:", ", ".join(f"{val:.6f}" for val in phi))
    print(f"迭代次数: {iteration}")
    return phi, iteration

测试与使用建议 #

1. 验证收敛性的测试用例 #

先使用严格对角占优的矩阵验证代码有效性：

# 严格对角占优的4阶矩阵，真实解为[1,1,1,1]
A_test = np.array([[10, 1, 2, 3],
                   [1, 11, 1, 2],
                   [2, 1, 12, 1],
                   [3, 2, 1, 13]])
b_test = np.array([14, 15, 16, 17], dtype=float)

# 调用函数：omega取0.8，初始猜测全0，收敛阈值1e-6，最大迭代100次
relax_solver(A_test, b_test, omega=0.8, initial_guess=np.zeros(4), 
             convergence_criteria=1e-6, max_iter=100, print_iter=True)

该用例会快速收敛到真实解。

2. 原测试矩阵的处理建议 #

如果必须求解原矩阵A：

• 先对矩阵做行变换（交换行、缩放行），使其满足严格对角占优
• 用直接解法np.linalg.solve(A,B)获取真实解，再调整初始猜测和omega值（可尝试0.1~0.3区间）
• 增加迭代次数，观察残差是否逐步下降

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