Best way to check if SOR solution has converged for 2d matrix

Question

I have written a SOR algorithm to solve the Laplace equation on a 2d grid. The outside of the grid is fixed at 0 and the central square is fixed at 10.

I can obtain the fully converged solution for some values of the relaxation parameter $\beta$ but not others. In the other cases the convergence condition is never met and thus iterates to infinity. I suspect it is related to floating point/precision error but I don't know how to avoid it.

The condition is as follows: $$r = \sum_{i}\sum_{j}\left|V_{i,j}^{n+1}-V_{i,j}^{n}\right| \\ \frac{r}{n}<c$$ Where $n$ is the number of non-fixed points and $c$ is a small finite value. Basically, it should converge when the average pointwise change in the solution is smaller than $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;
}

The grid is $N\times N$ (defined elsewhere in the program). Ideally, c=DBL_EPSILON. It converges well when $\beta=1$ (Gauss-Seidel) and some other values but not universally. It converges more consistently when c is made larger but it still fails in certain areas.

It has a tendency to blow up around $\beta\geq1.5$. For larger grid sizes this means the optimal parameter is not found since convergence is never reached. In the attached image $c$ is machine epsilon.

UPDATE: I suspect $\beta=1.5$ is significant to the problem. When the numbers become very small, multiplying by anything larger than 1.5 will cause numbers to be rounded up significantly. Effectively, a $\beta$ of 1.5 will start to behave like a $\beta=2$, for which SOR will never converge...

The question is: How do I avoid this?

Answer 1

I think the primary issue you're seeing is related to how you're checking for convergence.

The chosen value for your tolerance is likely too tight.

The tolerance should at a minimum allow 1 double precision (DP) epsilon up or down of floating point rounding error.

However, 1 DP epsilon is relative to the magnitude of the floating point value.

In this case, we happen to know that the maximum value is at the central node (assuming 0 boundary conditions, though those might end up being the max value).

So tol can be chosen as $\epsilon = 2^{-52} C$, where $C$ is the max (central) value.

The reason this is $2^{-52}$ and not $2^{-53}$ is because we want to allow 1 DP epsilon above $C$, not just 1 DP epsilon below $C$.

This value is approximately double tol = 2.23e-16 * C; in decimal. In practice, I usually pick something larger than this to allow for a few DP epsilon up or down, especially when there are more floating point operations involved per degree of freedom. double tol = 1e-15 * C; is a very reasonable choice, though I do often relax this tolerance even further when dealing with iterative solvers for complex systems.

A second thing you can do other than checking that v has converged is to check how well the error metric itself has converged/stagnated (e.g. compare r_prev with r). If you're not making any progress one iteration to the next, then the method has effectively converged to the best possible value of v it can resolve. Whether that converged solution is any good or not is up to you (methods can converge to invalid/undesirable solutions).

Here's my test code if you want to try and compare results:

#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);
}