我的Gauss-Seidel函数计算结果存在偏差，求问题排查

问题背景 #

编写了一个Gauss-Seidel函数，接收系数矩阵、常数向量数组以及整数n（n-1为矩阵长度），返回变量解数组，但计算结果与预期存在明显偏差。

测试用例 #

系数矩阵 #

double matrix[5][5] = 
{ {-1.601675, 0.700000, 0.000000, 0.000000, 0.000000},
  {0.700000, -1.201675, 0.500000, 0.000000, 0.000000},
  {0.000000, 0.500000, -0.801675, 0.300000, 0.000000},
  {0.000000, 0.000000, 0.300000, -0.401675, 0.100000},
  {0.000000, 0.000000, 0.000000, 1.000000, -1.008375},
};

常数向量 #

double sol[5] = 
{
  -180.041874,
  -0.041874,
  -0.041874,
  -0.041874,
  -0.209372,
};

预期结果 #

通过Wolfram Alpha计算，预期未知量结果：

double T_expected[5] = 
{
  198.298282,
  196.616925,
  194.965693,
  193.373744,
  192.046373,
};

实际运行结果 #

函数返回的结果：

double T[5] = 
{
  199.674915,
  199.257945,
  198.676137,
  197.711841,
  196.277411,
};

函数代码 #

#define EPSILON 1e-06
double *runGaussSeidel(double **matrix, double *sol, int n)
{
    int i, j, converged, m = n - 1;
    double sum1 = 0.0, sum2 = 0.0, temp = 0.0;
    double *T = malloc(sizeof(double) * m);

    // Initial guess.
    for (i = 0; i < m; i++)
        T[i] = 0.0;

    do {   
        converged = 1; 
        for (i = 0; i < m; i++)
        {
            sum1 = sum2 = 0.0;
            for (j = 0; j < i; j++)
            {
                sum1 += matrix[i][j] * T[j];
            }   
            for (j = i + 1; j < m; j++)
            {
                sum2 += matrix[i][j] * T[i];
            }

            double temp = T[i];

            T[i] = (sol[i] - sum1 - sum2) / matrix[i][i];                                                                   
        
            if (T[i] - temp >= EPSILON)
            {
                converged = 0;
            }                    
        }                                                                      
    } while (!converged); 
    
    return T;
}

偏差原因分析 #

1. sum2计算逻辑错误（核心问题） #

Gauss-Seidel迭代公式中，第i个变量更新时，需要计算所有j>i项的和，即sum2 = Σ(matrix[i][j] * T[j])（j从i+1到m-1）。但代码中错误地写成了matrix[i][j] * T[i]，相当于把T[i]乘以所有j>i的matrix[i][j]之和，完全违背迭代公式，直接导致计算值严重偏离正确结果。

修正后的sum2循环代码：

for (j = i + 1; j < m; j++)
{
    sum2 += matrix[i][j] * T[j];
}

2. 收敛判断不严谨 #

当前代码仅判断T[i] - temp >= EPSILON，只考虑解值增大且变化量超阈值的情况，忽略了解值减小且变化量绝对值超阈值的场景（即temp - T[i] >= EPSILON）。这会导致某些迭代步骤中，解的变化虽足够大但方向为负时，程序误判为收敛，提前终止迭代，结果精度不足。

修正后的收敛判断代码：

if (fabs(T[i] - temp) >= EPSILON)
{
    converged = 0;
}

（需确保代码包含<math.h>头文件以使用fabs函数）

3. 矩阵维度易混淆 #

函数中m = n - 1作为矩阵维度，但测试用例是5x5矩阵、5个元素的常数向量和解。若调用时传入n=5，m=4只会计算4个变量，与实际需求矛盾。建议调整函数逻辑，让n直接表示矩阵维度，避免通过n-1推导，减少出错概率。

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