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
    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.

enter image description here

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?

  • $\begingroup$ 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? $\endgroup$ Mar 5 '20 at 3:58
  • $\begingroup$ @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. $\endgroup$ Mar 5 '20 at 19:18
  • $\begingroup$ 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. $\endgroup$ Mar 5 '20 at 19:34
  • $\begingroup$ @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. $\endgroup$ Mar 5 '20 at 19:39
  • $\begingroup$ What is the rationale for if (i!=N/2 || j!=N/2) {? $\endgroup$ Mar 6 '20 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]) /
        // 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) {
    // max iterations
    if (iter >= 1000000) {

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.