Метод SOR не сходится

Question

У меня есть решение SOR для 2D Лапласа с условием Дирихле в Python. Если для омеги установлено значение 1,0 (что делает метод Якоби), решение сходится нормально. Но используя данные омеги, целевая ошибка не может быть достигнута, потому что решение в какой-то момент просто сходит с ума, не сходясь. Почему бы не сходиться с данной формулой омега? живой пример на 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

По моему опыту (но я никогда не тратил время на поиски объяснения), кажется, что конвергенция лучше при использовании так называемой красно-черной схемы обновления внутри той же сетки. Это означает, что вы смотрите на сетку, расположенную на шахматном узоре, и сначала вы обновляете ячейки одного цвета, а затем обновляете ячейки другого цвета.

Если мне что-то подобное с вашим кодом, то оно, похоже, сходится. Значение err был немного изменен, потому что вторая сетка больше не используется:

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

Вероятно, это очень неэффективный способ выбора "красных" и "черных" ячеек сетки, но я надеюсь, что понятно, что происходит таким образом.