2014-11-08 136 views
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 

回答

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 

這可能是在選擇「紅色」和「黑色」網格單元一個非常低效的方式,但我希望這是清楚什麼繼續這樣。

+0

謝謝,它現在真的在收斂,沒有任何鋸齒邊緣! – epicenter 2014-11-09 11:07:23