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

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} 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?

• But you know from theory that for SOR to converge, $\beta$ needs to be sufficiently small. Isn't what you see what you should expect from theory? – Wolfgang Bangerth Mar 5 at 3:58
• @WolfgangBangerth. No, you are mistaken. For SOR convergence follows when $0<\beta<2$. The optimal $\beta$ lies between 1 and 2 but it depends on the problem. For larger grid sizes the optimal value tends to increase, this trend can be seen to start in the graph. From some testing I determined that the computed solution in those cases does not explode, it just simply fails to trip the convergence condition. – user8384493 Mar 5 at 19:18
• Is this same thing as you see here: en.wikipedia.org/wiki/Successive_over-relaxation#/media/… I think your plot looks reasonable so I'm not sure what's the problem here. – Alone Programmer Mar 5 at 19:34
• @AloneProgrammer. No, that plots the spectral radius versus the relaxation parameter. I am interested in the actual number of iterations required to acheive convergence for different parameters. – user8384493 Mar 5 at 19:39
• What is the rationale for if (i!=N/2 || j!=N/2) {? – Abdullah Ali Sivas Mar 6 at 5:26

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