2014-11-08 5 views
0

파이썬에서 Dirichlet 조건을 가진 2D Laplace에 대한 SOR 솔루션이 있습니다. 오메가가 1.0으로 설정되면 (Jacobi 방법으로), 솔루션이 잘 수렴합니다. 그러나 주어진 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
내 경험에

(하지만 난 정말에 대한 설명을 볼 수있는 시간을 찍어 본적) 융합이 같은 그리드 내에서, 소위 레드 - 블랙 업데이트 방식을 사용하면 더 좋을 것입니다. 즉, 격자를 바둑판 패턴으로 배치 한 것으로보고 처음에는 한 색상을 가진 셀을 업데이트 한 다음 다른 색상의 셀을 업데이트합니다.

코드에서 이와 비슷한 경우, 수렴하는 것처럼 보입니다.

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 

이는 '레드'와 '블랙'그리드 셀을 선택하는 매우 비효율적 인 방법은 아마도,하지만 난 그게 무엇 분명 희망 : 두 번째 그리드는 더 이상 사용되지 않기 때문에 err의 의미는 약간 변경되었습니다 이쪽으로가.

+0

고마워요, 지금은 사실 들쭉날쭉하지 않고 정말 수렴 중입니다! – epicenter