0

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

1 回答 1

0

以我的经验（但我从未花时间真正寻找解释）在同一网格内使用所谓的红黑更新方案时收敛似乎更好。这意味着您将网格视为棋盘图案，首先更新具有一种颜色的单元格，然后更新另一种颜色的单元格。

如果我对您的代码进行这样的处理，它似乎确实会收敛。的含义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

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

于 2014-11-08T22:13:56.140 回答