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

Question

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

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

Если я что-то вроде этого с вашим кодом, похоже, сходится. Смысл 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 

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