为什么我的高斯-赛德尔函数略有偏差？

Question

我编写了一个 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;
}

Answer 1

我认为你有一个错误 - 你更改 T[i]，然后将其用于下一个变量/迭代。

尽管如此，您的 T_expected[5] 看起来也不正确。你尝试过测试一下吗？

这是我的代码，接近 C 但 C++

auto iteration(const double mtx[5][5],
               const double* rhs,
               const int n, double* res) -> void {
    double* tr = (double*)alloca(n * sizeof(double));
    memcpy(tr, res, n*sizeof(double));
    
    for(int k = 0; k != n; ++k) {
        double s = 0.0;
        for(int i = 0; i != n; ++i) {
            s += mtx[k][i] * tr[i];
        }
        
        res[k] = (rhs[k] - s)/mtx[k][k] + tr[k];
    }
}

auto test(const double mtx[5][5], const double* rhs, const int n, const double* res) -> double {
    double d = 0.0;
    
    double s = 0.0;
    for(int k = 0; k != n; ++k) {
        double s = 0.0;
        for(int i = 0; i != n; ++i) {
            s += mtx[k][i]*res[i];
        }
        
        d += fabs(rhs[k] - s);
    }
    return d;
}

1000次迭代后的解决方案是

198.568 197.141 195.721 194.307 192.901

它计算的绝对差之和（test()函数）等于

7.65291e-14