SOR 方法不收敛

Question

我有一个在 Python 中具有狄利克雷条件的二维拉普拉斯的 SOR 解决方案。 如果 omega 设置为 1.0（使其成为 Jacobi 方法），则解收敛得很好。 但是使用给定的欧米茄，无法达到目标误差，因为解决方案在某个时候变得疯狂，无法收敛。为什么它不与给定的欧米茄公式收敛？ 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

这可能是选择“红色”和“黑色”网格单元格的一种非常低效的方法，但我希望清楚这种方式发生了什么。