0
我有一個SOR解決方案,用於2D拉普拉斯和Python中的Dirichlet條件。 如果歐米茄設置爲1.0(使其成爲雅可比方法),解決方案收斂很好。 但是,使用給定的omegas,目標錯誤無法達成,因爲解決方案在某個時刻只是瘋狂,未能收斂。爲什麼它不會與給定的歐米茄公式一致? live example on repl.itSOR方法不收斂
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
謝謝,它現在真的在收斂,沒有任何鋸齒邊緣! – epicenter 2014-11-09 11:07:23