计算科学 - 检查 SOR 解决方案是否已针对 2d 矩阵收敛的最佳方法 - 吾爱随笔录

检查 SOR 解决方案是否已针对 2d 矩阵收敛的最佳方法

计算科学迭代法收敛 C 精确

2021-12-06 02:33:35

我编写了一个 SOR 算法来求解二维网格上的拉普拉斯方程。网格外侧固定为 0，中心方格固定为 10。

的某些值，我可以获得完全收敛的解，但不能获得其他值。在其他情况下，永远不会满足收敛条件，因此迭代到无穷大。我怀疑它与浮点/精度错误有关，但我不知道如何避免它。 $\beta$

条件如下：其中是非固定点的数量，是一个小的有限值。基本上，当解中的平均逐点变化小于时，它应该收敛。

r = \sum_{i} \sum_{j} | V_{i, j}^{n + 1} - V_{i, j}^{n} | \frac{r}{n} < c

$r = \sum_{i}\sum_{j}\left|V_{i,j}^{n+1}-V_{i,j}^{n}\right| \\ \frac{r}{n}<c$

n

$n$

c

$c$

c

$c$

//Successive over-relaxation.  argument b is the relaxation parameter.
unsigned int sor(double *v,double b) {
  short i,j,t=1;
  unsigned int s=0;
  double l,r,n=((N-2)*(N-2)-1);
  while (t) {//t is true when non-converged
    r = 0.0;
    for (i=1;i<N-1;i++) {
      for (j=1;j<N-1;j++) {
        if (i!=N/2 || j!=N/2) {
          l = (*(v+N*(i+1)+j) + *(v+N*(i-1)+j) + *(v+N*i+j+1) + *(v+N*i+j-1))*0.25;
          r += fabs(b*(l-*(v+N*i+j))); //add to residual
          *(v+N*i+j) += b*(l-*(v+N*i+j)); //update new value
        }
      }
    }
    if (r/n<c) {t = 0;} //check convergence
    s++;
    if (s==USHRT_MAX) {t = 0;}  //iteration limit is reached for some values of b.  It will go higher than this, but USHRT_MAX is well beyond what is expected
  }
return s;
}

网格是（在程序的其他地方定义）。理想情况下，. (Gauss-Seidel) 和其他一些值时，它收敛得很好，但不是普遍的。当变大时，它会更一致地收敛，但在某些区域仍然会失败。 $N\times N$ c=DBL_EPSILON $\beta=1$ c

它有在附近爆发的趋势。对于较大的网格尺寸，这意味着找不到最佳参数，因为永远不会达到收敛。在附图中，是机器 epsilon。 $\beta\geq1.5$ $c$

更新：我怀疑对问题很重要。当数字变得非常小时，乘以任何大于 1.5 的值将导致数字显着向上取整。实际上，一个将开始表现得像一个，SOR 永远不会收敛...... $\beta=1.5$ $\beta$ $\beta=2$

问题是：我该如何避免这种情况？

1个回答

我认为您看到的主要问题与您如何检查收敛有关。

为您的容差选择的值可能太紧了。

容差应至少允许浮点舍入误差向上或向下 1 个双精度 (DP) epsilon。

但是，1 DP epsilon 与浮点值的大小有关。

在这种情况下，我们碰巧知道最大值在中心节点处（假设边界条件为 0，尽管最终可能是最大值）。

所以 tol 可以选择为，其中是最大（中心）值。 $\epsilon = 2^{-52} C$ $C$

这是而不是的原因是因为我们希望允许上方 1 DP epsilon ，而不仅仅是下方 1 DP epsilon 。 $2^{-52}$ $2^{-53}$ $C$ $C$

该值大约double tol = 2.23e-16 * C;为十进制。在实践中，我通常选择比这更大的东西来允许向上或向下一些 DP epsilon，特别是当每个自由度涉及更多浮点运算时。double tol = 1e-15 * C;是一个非常合理的选择，尽管在处理复杂系统的迭代求解器时，我经常会进一步放宽这种容忍度。

除了检查是否v已经收敛之外，您可以做的第二件事是检查错误度量本身收敛/停滞的程度（例如与比较r_prev）r。如果您在一次迭代到下一次迭代中没有取得任何进展，那么该方法已经有效地收敛到v它可以解决的最佳可能值。该融合解决方案是否有用取决于您（方法可能会收敛到无效/不受欢迎的解决方案）。

如果您想尝试比较结果，这是我的测试代码：

#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>

size_t sor(double *v, double beta, double central_value, size_t width) {
  size_t iter = 1;

  // 2^-52 to allow 1 epsilon above or below central_value
  // in practice, you might want to allow a few epsilon
  double tol = central_value * 2.220446049250313e-16;

  for (;;) {
    double err = 0;
    // perform 1 iteration of SOR
    for (size_t row = 1; row < width - 1; ++row) {
      for (size_t col = 1; col < width - 1; ++col) {
        double vnp1 = (1 - beta) * v[row * width + col];

        if (row == width / 2 && col == width / 2) {
          // replace center cell's equation with identity
          vnp1 += beta * central_value;
        } else {
          vnp1 += beta *
                  (v[row * width + col - 1] + v[row * width + col + 1] +
                   v[(row - 1) * width + col] + v[(row + 1) * width + col]) /
                  4;
        }
        // accumulate L1 error metric of v
        err += fabs(vnp1 - v[row * width + col]);
        v[row*width+col] = vnp1;
      }
    }
    // check convergence
    if (err / ((width - 2) * (width - 2) - 1) <= tol) {
      break;
    }
    // max iterations
    if (iter >= 1000000) {
      break;
    }
    ++iter;
  }

  return iter;
}

int main(int argc, char **argv) {
  // test driver
  size_t width = 15;
  double beta = 1.5;

  double *v = (double *)malloc(sizeof(double) * width * width);
  for (size_t i = 0; i < width * width; ++i) {
    v[i] = 0;
  }
  size_t iter = sor(v, beta, 10, width);
  printf("%zu\n", iter);
  free(v);
}

其它你可能感兴趣的问题

上一篇从有限元模型计算等效热阻和热容下一篇CVX能否解决凸问题？