SOR method does not converge

Question

I have a SOR solution for 2D Laplace with Dirichlet condition in Python. If omega is set to 1.0 (making it a Jacobi method), solution converges fine. But using given omegas, target error cannot be reached because solution just goes wild at some point, failing to converge. Why wouldn't it converge with given omega formula? live example on repl.it

from math import sin, exp, pi, sqrt

m = 16

m1 = m + 1
m2 = m + 2
grid = [[0.0]*m2 for i in xrange(m2)]
newGrid = [[0.0]*m2 for i in xrange(m2)]

for x in xrange(m2):
    grid[x][0] = sin(pi * x / m1)
    grid[x][m1] = sin(pi * x / m1)*exp(-x/m1)

omega = 2 #initial value, iter = 0
ro = 1 - pi*pi / (4.0 * m * m) #spectral radius
print "ro", ro
print "omega limit", 2 / (ro*ro) - 2/ro*sqrt(1/ro/ro - 1)


def next_omega(prev_omega):
    return 1.0 / (1 - ro * ro * prev_omega / 4.0)

for iteration in xrange(50):
    print "iter", iteration,
    omega = next_omega(omega)
    print "omega", omega,
    for x in range(1, m1):
        for y in range(1, m1):
            newGrid[x][y] = grid[x][y] + 0.25 * omega * \
                (grid[x - 1][y] + \
                 grid[x + 1][y] + \
                 grid[x][y - 1] + \
                 grid[x][y + 1] - 4.0 * grid[x][y])
    err = sum([abs(newGrid[x][y] - grid[x][y]) \
        for x in range(1, m1) \
        for y in range(1, m1)])
    print err,
    for x in range(1, m1):
        for y in range(1, m1):
            grid[x][y] = newGrid[x][y]
    print

Answer 1

In my experience (but I never took the time to really look for an explanation) convergence seems to be better when using a so-called red-black update scheme, within the same grid. This means that you look at the grid as being laid out on a checkerboard pattern, and first you update the cells that have one color, and afterwards you update the cells of the other color.

If I something like this with your code, it does seem to converge. The meaning of err was changed slightly because the second grid isn't used anymore:

for iteration in xrange(50):
    print "iter", iteration,
    omega = next_omega(omega)
    err = 0
    print "omega", omega,
    for x in range(1, m1):
        for y in range(1, m1):
            if (x%2+y)%2 == 0: # Only update the 'red' grid cells
                diff = 0.25 * omega * \
                    (grid[x - 1][y] + \
                     grid[x + 1][y] + \
                     grid[x][y - 1] + \
                     grid[x][y + 1] - 4.0 * grid[x][y])

                grid[x][y] = grid[x][y] + diff
                err += diff

    for x in range(1, m1):
        for y in range(1, m1):
            if (x%2+y)%2 == 1: # Only update the 'black' grid cells
                diff = 0.25 * omega * \
                    (grid[x - 1][y] + \
                     grid[x + 1][y] + \
                     grid[x][y - 1] + \
                     grid[x][y + 1] - 4.0 * grid[x][y])

                grid[x][y] = grid[x][y] + diff
                err += diff

    print err

This is probably a very inefficient way in selecting the 'red' and 'black' grid cells, but I hope it's clear what's going on this way.