2번이 왜 이제 올라오느냐도 역시 의미 없는 질문이다.


--

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

# Differential equation
# 추천 참고자료: http://www.swarthmore.edu/NatSci/echeeve1/Ref/NumericInt/FrameNumInt.html

# 미분방정식은 미분 변수가 1개인 상미분방정식(Ordinary differential equation:ODE)과 2개 이상인 편미분방정식(Partial differential equation:PDE)으로 나눠진다.
# 일단 1개인 미방을 어떻게 푸는지 알아보자. 연립 미방은 아직 다루지 않고, 나중에 다룬다.

# ODE는 선형 미방과 비선형 미방으로 나눠진다. 선형미방은 미방의 해가 두개 이상 있을 때, 두 해의 선형결합도 해가 되는 경우이다.
# 즉, f1과 f2가 미방의 연산자 L에 대해서 L(f1)=L(f2)=0일 때 L(f1+f2)=0이 만족되는 경우이다. 물론 그 역을 선형 미방이라고 하지는 않는다.
# 비선형 미방은 풀기 골치아프므로 일단 선형 미방을 어떻게 푸는지 알아보자.

# 선형 미방은 다시 최대 차수인 도함수에 따라 n차 미방으로 구분된다. 즉, 가령 2차 도함수가 최대차수인 도함수라면 2차 미분방정식이 된다.

# 1차 이상의 미분을 포함하는 고차 미분방정식은 여러개의 중간 변수를 도입하여 1차 연립 미분 방정식으로 간주할 수 있다. 따라서, 1차 미분방정식을 잘 풀 수 있다면 고차 미분방정식은 손쉽게 풀린다. 그러므로 1차 미방에 집중하자.
# 1차 미분방정식을 이산화(discretize)하게 되면 차분방정식(difference dquation)이 된다. 차분방정식은 고등학교 때 배우는 계차수열의 한 형태이다.

# 그럼, 이제 차분방정식을 풀어보자. 가장 간단하게, 오일러 풀이법이 있다.

# 오일러 풀이법에서는 dy/dt + a*y(t) = f(t) 라는 문제를 풀게 된다. a는 값을 알고 있는 상수이고, f(t)는 이미 알고 있는 함수이다.

# 오일러 풀이법은 정말 간단하다. 주어진 시간 t로부터 h만큼의 시간이 지난 후의 값인 y(t+h)를 다음과 같이 근사한다.
# y(t+h)=y(t)+(dy/dt)h
# 아주 간단한 1차 근사가 되겠다. 그럼, 이제 y(t+h)를 알았으니, 같은 방법으로 한번 더 계산하면 y(t+2h)도 알아낼 수 있다. 이 작업을 반복하면, y(t=t0)를 알고 있을 때 우리는 아주 많은 문제를 풀 수 있게 된다.

# y(t+h) = y(t)+(dy/dt)h = y(t) + (- a*y(t) + f(t))*h
# 잘 생각해보면, 이것은 아주 간단한 반복작업이 된다. 이를 위한 파이썬 코드를 쓴다면, 다음과 같을 것이다.

import numpy as np
a=2.
f = lambda x, 3.*np.exp(-4.*x)
y0 = 0. # t=t0일때의 값이라고 하자.
h = 0.1 # time slice가 된다. 이 값은 작을수록 더 정확한 답을 얻게 되지만, 계산이 느려진다.
y = [y0] # 일단 여기서 시작해 보자.
t = 0.
tmax = 10. # 얻고싶은 시간의 마지막 순간. 즉, 이 시간 이후에는 관심이 없고 그 전까지만 관심이 있다는 뜻이다.
while t<tmax:
    y.append(y[-1]+(f(t)-a*y[-1])*h)
    t+=h

# while 안에 있는 코드가 가장 중요하다. 잘 보면, 별거 없다는 사실을 알 수 있을 것이다.

# 이제 2차 미방을 1차 미방 2개로 쪼개는 기법을 아주 간단히 이해하고 넘어가자.
# 주어진 2차 미방이 d(dy/dt)/dt + a*dy/dt + b*y(t) = f(t) 라고 하면, 일단 dy/dt = z 라고 하자. 그리고 z를 대입하면
# dz/dt + a*z + b*y = f(t)
# dy/dt = z

# 이렇게 2개의 미분방정식으로 쪼개졌다. 한번 풀어보도록 하자.

# 일단 y0와 z0=dy/dt(t=0)는 알려져 있다고 하자.
y = [y0]
z = z0.
# z는 우리가 문제를 풀기를 원하는 값이 아니라 중간 저장 변수이므로 그냥 이렇게 두면 된다. 만약 2차 미분방정식인데 1차 미분도 같이 구해야 하는 경우(뉴턴 방정식에서 속도와 위치를 둘 다 구하는 경우, 회로 방정식에서 전하와 전류를 둘 다 구하는 경우)  등등에서는 z를 리스트로 선언해서 써도 된다.

while t<tmax:
    y.append(y[-1]+z*h)
    z+=(-a*z-b*y+f(t))*h
    t+=h
   
# 2차 미방이지만 위와 같이 딱 1줄 더 늘어난 아름다운 코드로 만들 수 있다.

# 오일러 풀이법은 가장 간단하지만 h가 상당히 작아야 정확한 답을 얻을 수 있다. 미분이나 적분을 할 때 함수값의 대칭성을 이용해서 같은 구간 h에 대해서도 보다 정확한 계산을 할 수 있었는데, 여기서도 마찬가지 테크닉을 사용할 수 있다.

# 룽게-쿠타(Runge-Kutta) 풀이법은 한발 앞서 간 도함수와 평균을 내서 다시 구한다. 즉, 기울기 계산을 두번한다. 말로는 알아듣기 어려울 것이고, 직접 코드를 보는 것이 보다 이해하기 쉬울 것이다.
# y(t+h) = y(t)+(dy/dt)h = y(t) + (- a*y(t) + f(t))*h
# dy/dt + a*y(t) = f(t)
y = [y0]
while t<tmax:
    y1= -a*y[-1]+f(t) # 일단 기울기 계산
    y.append(y1*h+y[-1]) # 일단 값 집어넣고
    y2 = -a*y[-1]+f(t+h) # 한칸 간 곳에서의 기울기 계산
    y[-1] = y[-2]+0.5*(y1+y2)*h # 실제 값은 두 기울기의 평균만큼만 증가하는 것으로.
    t+=h
   
# 계산 시간을 아주 조금이라도 더 절약하고 싶으면 0.5*h를 미리 h2같은걸로 정의해놓고 0.5*h대신 h2를 곱하는 것으로 바꾸면 조금 빨라진다.
# y2를 계산하는 부분은 그 다음줄로 한번에 집어넣을 수도 있지만, 읽기가 힘들어지므로 굳이 그렇게 할 필요는 없다.

# 2차 룽게-쿠타 방법 뿐만 아니라 더 높은 차수의 룽게-쿠타 방법도 가능하다. 위와 같이 기울기 가중치를 어떻게 주느냐를 잘 생각하면 되는데, 아이디어는 비슷하므로 크게 어려울 것 없다.

# 조금 더 일반화하면 predictor-corrector 방법이 된다.
# 미분방정식이 dy/dt+f(y(t), t)=0으로 주어진 경우에 사용하는데, 위의 미분방정식은 f(y, t) = a*y-f(t)로 주어진 하나의 사례에 해당한다.
# Predictor-corrector 방법은 그럼 RK방법을 그대로 갖고 오면 된다.
y = [y0]
while t<tmax:
    y1= f(y[-1], t) # 일단 기울기 계산
    y.append(y1*h+y[-1]) # 일단 값 집어넣고
    y2 = f(y[-1], t+h) # 한칸 간 곳에서의 기울기 계산
    y[-1] = y[-2]+0.5*(y1+y2)*h # 실제 값은 두 기울기의 평균만큼만 증가하는 것으로.
    t+=h

# 물리학에서 풀어야 하는 수많은 미분방정식들은 대체로 운동방정식이다. 운동방정식은 2차 미분방정식이므로, 2차 도함수를 알고 있다는 사실로부터 뭔가를 직접 해결할 수 있다.
# 그래서 나온 방법이 벌레(Verlet) 방법이다. Verlet은 사람 이름이고, 프랑스 사람이라 t가 묵음이다.
# 벌레 알고리즘에 대해서는 이 글을 참고해 볼 수 있다.
# http://www.physics.udel.edu/~bnikolic/teaching/phys660/numerical_ode/node5.html
# 일단, 테일러 전개로부터 다음과 같은 사실을 알 수 있다.
# y(t-h) = y(t)-v(t)*h+a(t)*h*h/2
# 여기서 v=dy/dt, a=dv/dt 이다.

# -h가 아니라 +h에 대해서 다시 쓰면
# y(t+h) = y(t)+v(t)*h+a(t)*h*h/2
# 이렇게 될 것이다.

# 두개를 더해보자.
# y(t-h)+y(t+h) = 2y(t)+a(t)*h*h
# v가 사라졌다!

# 이제 다음과 같이 y(t+h)를 추정할 수 있다.
# y(t+h) = 2y(t)-y(t-h)+a(t)*h*h
# 물리 문제를 해결할 때, 이 공식은 매우 기가막힌 방법이다. 2차 미분방정식을 풀어야 하는데, 속력을 구하지 않고 힘으로부터 곧바로 위치를 얻어내기 때문이다.

# 또한, 두개를 빼면 속력을 얻을 수 있다.
# v(t) = (y(t+h)-y(t-h))/2h
# 이 공식은 미분에서 공부했던 3점 공식과 같다.

# 이렇게 계산하는 방법을 original Verlet 알고리즘이라고 한다. 그런데, 오리지널 벌레는 문제가 있다. y(t=0)일때를 알고 있어도, y(t+h)를 구하기 위해서는 y(t-h)를 알아야 한다. 즉, 최초 2점을 알아야 한다는 문제가 있다. 그래서 이 방법을 개선하기 위해서 velocity Verlet 알고리즘이 등장한다. 일단 공식부터 보고 시작하자.

# y(t+h) = y(t)+v(t)*h+0.5*a(t)*h*h
# v(t+h) = v(t)+0.5*(a(t+h)+a(t))*h

# 왠지 앞에서 보았던 오일러 풀이법으로 2차 미분방정식을 푸는 것과 비슷해 보이지만, 그보다 더 정확하다. 그리고 이 경우 알지도 못하는 과거인 y(t-h)를 몰라도 문제를 해결할 수 있으므로 고민이 한층 덜어졌다.

# 이 문제에서 a(t) 가 알려져 있다고 하면, 코드를 다음과 같이 짜볼 수 있을 것이다.
y=[y0]
v=v0 # 2차 미방이므로 y0와 v0는 알고 있어야 한다.
while t<tmax:
    y.append(y[-1]+v*h+0.5*a(t)*h*h)
    v+=0.5*(a(t+h)+a(t))*h
    t+=h
# 이외에도 많은 알고리즘이 있지만 대체로 이런 아이디어들의 변형과 개선을 도모한 방법들이다. 실제로 자신이 풀려고 하는 문제에 적합한 알고리즘을 선정하고, 그 문제에 적합하게 알고리즘을 개조/개선해서 적용하는 것이 좋다.

# 수치해석적인 방법의 장점은 선형 미분방정식 뿐만 아니라 비선형 미분방정식의 경우에도 적용할 수 있다는 것이다. 사실, 선형 제차 미분방정식의 해는 매우 쉽게 해석적으로 구할 수 있다. 그러나 비제차 미분방정식이나 비선형 미분방정식의 경우 해가 알려진 것이 매우 적으나, 실용적 관점에서 수치해석을 이용하여 문제를 해결할 수 있다.
# 수치해석으로 비선형 미분방정식을 해결할 때 주의해야 하는 것은 오차의 전파이다. 비선형 효과에 의해 오차가 예상보다 커질 수 있으므로 오차를 줄이기 위해서 기존의 공식을 적절히 변형하는 것이 필요할 수 있다.

by snowall 2013. 9. 10. 02:19