이번엔 왜 6번이냐고 물으신다면 몇번까지 했는지 기억이 나지 않기 때문...이라고.

# Elementary Numerical analysis 6
# based on Python
# (C) 2013. Keehwan Nam, Dept. of physics, KAIST.
# snowall@gmail.com / snowall@kaist.ac.kr

# Relaxation method
# http://snowall.tistory.com/2561 이 글의 개선된 버전이 되겠다.

# 라플라스 방정식이나 푸아송 방정식을 풀기 위해서 사용하는 가장 대표적인 방법은 뭐니뭐니해도 가우스 법칙이다. 가우스 법칙은 표면에서 필드를 적분하면 내부에 존재하는 전하량에 비례한다는 사실을 말해준다. 여기서 완화법(relaxation method)이 왜 작동하는지 설명할 수 있다.

# 라플라스 방정식의 해가 조화 함수(Harmonic function)이라고도 불리운다는 사실을 기억하자. 조화함수의 가장 중요한 특성 중 하나는 평균값 정리인데, 이 정리는 조화 함수가 잘 정의되는 임의의 점에서 성립한다. 주어진 임의의 점에서 거리가 일정한 점의 집합, 즉 어떤 '구면'(spherical surface)을 정의할 수 있다. 이 구면에서(3차원이나 2차원 표면이 아니라도) 조화함수의 함수값의 평균값은 구면의 중심점에서의 함수값과 같다.

# 그런데, 우리가 해결하고자 하는 문제의 공간은 격자로 잘게 쪼개져 있다. 따라서,  '구면'은 바로 옆에 있는 점들로 정의된다. 자, 이제 2차원 구면을 생각해 보자.

# 함수 f(x, y)가 잘 정의되고, x0, y0에서 h만큼 떨어져 있는 '구면'은 다음과 같이 정의된다.

# 구면 = {(x0+h, y0), (x0-h, y0). (x0, y0+h), (x0, y0-h)}
# 2차원이고 네모난 격자이므로 이렇게 점 4개가 '구면'을 정의한다. 이 점 위에서의 함수값들의 평균값은 다음과 같이 정의된다.

# f_mean(x0, y0) = (f(x0+h, y0)+f(x0-h, y0)+f(x0, y0+h)+f(x0, y0-h))/4
# 적분은 덧셈으로 바뀌었고, 4개 있으니까 4로 나눈 것이다. 아주 쉽다.
# 조화 함수의 특성으로부터, f_mean(x0, y0) = f(x0, y0) 이어야 한다.

# 하지만, f(x,y)를 알고 있다면 당연히 이게 성립하는걸 확인할 수 있겠지만 모르는 마당에 어떻게 확인할 수 있을까?
# 라플라스 방정식이나 푸아송 방정식을 공부해 본 사람이라면 누구나 알고 있겠지만, 이 방정식들의 해는 만약 존재한다면 유일하게 존재한다. 즉, 수단과 방법을 가리지 않고 그런 해를 찾아냈다면 그 해는 우리가 찾아 헤메이던 바로 그것이다.

# 만약 함수값 f(x0, y0)을 위에서 계산한 f_mean(x0, y0)으로 대신한다면, 그래도 아무것도 모르는 것 보다는 답에 조금 더 가까울 것이다. f_mean(x0, y0)을 그렇게 두고, f_mean(x0+h, y0)을 계산하고, f_mean(x0, y0+h)를 계산하고, 그렇게 계산할 수 있다. 이렇게 h씩 움직이면서 모든 점을 한바퀴 돌고 나면, 아마 원래 알고 있던 함수보다는 조금 더 진짜 함수에 가까워진 모양이 될 것이다. 이 짓을 더이상 함수값이 변하지 않을 때 까지 반복하자.

# 이 방법은 왜 정답에 수렴할까? 왜냐하면, 경계조건이 있기 때문이다. 경계에서 주어진 함수값은 정해져 있고, 따라서 경계 근처에서 평균을 내면 어떻게든 경계의 함수값을 따라갈 수밖에 없다. 경계 근처에서 라플라스, 푸아송 방정식을 만족한다면, 그 옆에서도, 그 옆에서도, 어떤 위치에서도 다 만족할 수 밖에 없다. 충분히 오래 반복한다면 답은 반드시 나올 것이다.

# 이제, 예를 들어서 답을 찾아보자.

import numpy as np
f = []
for i in range(52):
    f.append([])
    for j in range(52):
        f[-1].append(0)
        
# 0으로 가득 찬 52 x 52 짜리 행렬을 만들었다.

for i in range(52):
    f[0,i] = sin(float(i))
    f[-1,i] = sin(float(51-i))
    f[i,0] = sin(float(51-i))
    f[i,-1] = sin(i)
    
# 네모난 경계에 경계조건을 주었다.

desire = 1.e-10 # 내가 원하는 오차의 한계
error = 1.
while error>desire:
    error = 0.
    for i in range(1, 51):
        for j in range(1, 51):
            tmp_error= f[i, j]
            f[i, j] = (f[i-1, j]+f[i+1, j]+f[i, j+1]+f[i, j-1])/4. # 문제를 해결하는 핵심 코드.
            error+= abs(tmp_error-f[i, j])

result = open("solution.txt", "w")
result.write(f)
result.close()
# 이게 루틴의 끝이다.
# 심지어 문제를 푸는 루틴이 경계조건을 부여하는 루틴보다 간단하다!

# 푸아송 방정식이라면?
while error>desire:
    error = 0.
    for i in range(1, 51):
        for j in range(1, 51):
            tmp_error= f[i, j]
            f[i, j] = -h*h*rho(x[i], y[j])+(f[i-1, j]+f[i+1, j]+f[i, j+1]+f[i, j-1])/4. # 문제를 해결하는 핵심 코드.
            error+= abs(tmp_error-f[i, j])


# -h*h*rho(x[i], y[i])이 추가된 것 외에 똑같다. rho는 푸아송 방정식에서 inhomogenius인 부분을 나타내는 함수이다. 저게 0이면 라플라스 방정식과 똑같다. x[i]와 y[j]는 i, j에서의 실제로 주어진 x, y값이다. 이것도 속도를 빠르게 하고 싶으면 미리 rho(x, y)를 rho[i, j]로 만들어 놓고 시작하면 된다.
            
# 완화법의 문제는 느리다는 점이다. 이것을 해결하기 위해서, 해가 대충 알려져 있는 경우 f(x, y)를 위에서 쓴 것처럼 0으로 배치하지 말고 대충 알려진 근사함수로 미리 넣고 시작하면 보다 빠르게 수렴할 것이다. 사실 미리 넣어주는 함수를 아무 함수나 집어넣어도 함수가 무한대로 뻗치지 않는 한 무조건 수렴한다. 그러나 빠르게 값을 얻고 싶다면 미리미리 비슷해 보이는 함수에서 시작하는 것이 조금 더 빠른 답을 얻는데 도움이 된다.

# 미리 넣고 시작하는 근사 함수를 얻기 위해서, 격자를 우선 성기게 잘라서 대충 답을 구하고, 이 답을 잘게 쪼갠 격자에 대입해서 시작하면 보다 빠르게 답을 얻을 수 있다. 이런 방법을 다중격자(multigrid) 방법이라고 한다.


# 이 방법을 확장한 Successive Over-relaxation (SOR) method 라는 방법이 있다. 이것은 주변의 평균값을 더할 때 가중치를 주어서 더 빠르게 수렴하도록 하는 방법이다.


# 라플라스 방정식과 푸아송 방정식 등, 타원꼴 편미분방정식은 이 방법으로 언제나 풀 수 있다. 고차원에서도 성립하므로, 2차원보다 높은 차원에 대해서도 손쉽게 확장할 수 있다. 문제는 격자를 쪼개는 것이 사각형이 아니면 매우 곤란해진다는 점이다. 이 문제를 해결하기 위해 유한요소법(Finite element method, FEM)이 도입되었다. 그러나 다음 시간에는 유한차 시간영역 법(Finite difference time domain, FDTD)을 다룰 것이다. FEM은 완화법을 보다 일반적으로 확장한 것이다. 어쨌든, 다음 시간에.


by snowall 2013. 9. 18. 02:53