글
이번엔 왜 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은 완화법을 보다 일반적으로 확장한 것이다. 어쨌든, 다음 시간에.
댓글