일차 미분 방정식의 수치 적분

예를 들어, 다음과 같은 미분 방정식을 생각해 보자.

[eqn: {dx} over {dt} =f(x)]

위의 식은 [eqn: dx=f(x)dt]로 적을 수 있으므로, [eqn: dx = x(t+dt) - x(t)]임을 이용하면

[eqn: x(t+dt)=x(t)+dtf(x(t))]의 형태가 된다. 이 식의 우변은 시간 [eqn: t]에서의 정보에만 관계하므로, 이를 이용해 [eqn: t+dt]의 정보를 얻게 되는 꼴이다. 이 식이 정확히 등호로 성립하는 조건은 당연히 [eqn: dt]가 0으로 수렴하는 경우인데, 컴퓨터 프로그램을 이용해서 [eqn: dt]라는 시간이 경과한 다음의 정보를 얻기 위해서는 충분히 작지만 너무 작지는 않은 [eqn: dt]를 이용할 수밖에 없다. 위의 방식으로 미분방정식을 적분하는 것을 오일러(Euler) 방법이라고 부른다. 즉,

[eqn: x(t+ DELTA t)`=x(t)+ DELTA tf(x(t))]

위의 식에서 좌변을 테일러 전개하면

[eqn: x(t+ DELTA t)=x(t)+ DELTA t dx overdt right|_t + 1 over2 Deltat^2 {d^2 x} over dt^2 right|_t + O(Deltat^3 )]

이므로 식의 우변과 비교하면, 오일러 방법의 오차는 [eqn: O(Deltat^2 )]이다. 즉, 오일러 방법의 정확도는 [eqn: Deltat]에 비례한다. 수치적분 방법의 정확도가 [eqn: DELTA t^n]에 비례하면 [eqn: n]차 방법([eqn: n]-th order)이라 한다. 따라서, 오일러 방법은 1차 적분방법이다.

먼저 살펴볼 간단한 고전역학 문제는 바로, 공기의 저항을 받으며 아래로 떨어지는 물체의 운동이다. 이 물체에 작용하는 힘은 아래를 향한 중력([eqn: mg])과 물체의 속도의 반대 방향의 공기의 저항력으로 구성되는데, 물체의 속도가 아주 빠르지 않을 때에는 저항력을 간단히 [eqn: -bv]로 적을 수 있어서(속도가 클 때는 [eqn: v ^{2}]에 비례하는 저항력이 중요해진다) 전체 힘은 [eqn: F = -bv + mg]가 된다. 즉, 이 경우 물체의 속도는 아래의 미분 방정식을 만족한다.

[eqn: F=ma``` -> ```m dv over dt = -bv + mg]

위 식에 오일러 방법을 적용하면 [eqn: v(t + Delta t) = v(t) + Deltat (g - gamma v)]가 되는데, 여기서 [eqn: gamma = b/m]이다. 자, 이 식을 적분하는 프로그램을 작성해 보자.

import numpy as np, matplotlib.pyplot as plt

g = 9.8; gamma=1.0; dt = 0.05; maxt = 10.0

v = 0.0

t_arr = []

v_arr = []

for t in np.arange(0, maxt, dt):

t_arr.append(t)

v_arr.append(v)

v += dt*(g-gamma*v)

plt.plot(t_arr,v_arr)

plt.show()

프로그램의 한 줄에 여러 코드를 적으려면 g = 9.8; gamma=1.0처럼 세미콜론(;)을 이용하면 된다. np.arange(a,b,d)는 a에서 b까지 d의 간격으로 나열된 숫자들의 배열을 얻는 방법이다(앞 장의 예제 참고). 사실, 위의 미분방정식의 해는 해석적으로도 어렵지 않게 얻을 수 있다. [eqn: dv over {( g - gamma v)} = dt]의 꼴로 적고 양변을 적분하면 [eqn: int {dv} over {v–g/ gamma } =- int gamma dt]이므로 [eqn: ln(v - g/gamma) = -gamma t + C]이고, 따라서 [eqn: v - g/gamma = C’ exp(-gamma t)]인데, 처음 조건인 [eqn: v(t=0) =0]을 대입하면 [eqn: v= {g} over {gamma } (1–e ^{- gamma ```t} )]이 운동방정식의 해가 된다.

== 예제 1 ==

오일러 방법을 이용한 위 프로그램을 이용해서 [eqn: v= {g} over {gamma } (1–e ^{- gamma ```t} )]와 수치적분의 결과의 차이(오차)의 절대값을 시간의 함수로 그리시오. dt를 0.001, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5로 바꾸면서 각각의 dt에서의 오차의 최댓값을 구하고, 이를 dt의 함수로 그리시오.

위에서 살펴본 오일러 방법은 실제 연구나 계산에서는 널리 쓰이지는 않는다. 충분히 정확한 결과를 얻기 위해 사용해야 하는 [eqn: Delta t]의 값이 너무 작아서, 주어진 시간(위 프로그램의 maxt)까지 적분하기 위해 필요한 계산량이 많기 때문이다. 오일러 방법에 비하면 계산량이 두 배 정도에 불과하지만 계산의 정확도는 [eqn: Deltat^2]에 비례해서 훨씬 정확한 2차 방법이 널리 쓰이는데, 가장 대표적인 것으로는 2차 룬제-쿠타 방법(2nd-order Runge-Kutta)이 있다. [eqn: x(t+dt)=x(t)+dtf(x)]를 직접 이용하는 것이 오일러 방법이라면 룬제-쿠타 방법에서는 우변의 [eqn: f(x)]대신에 [eqn: f(x(t+ {dt} over {2} ))]를 이용하는 거다. 물론 [eqn: x(t+ dt over2 )]를 먼저 구해놓아야 이 식을 쓸 수 있기 때문에 룬제-쿠타 방법은 다음과 같은 과정으로 이루어진다. 이중 [eqn: t+ dt over 2]에서의 [eqn: x]의 값을 구하는 첫 번째 식은 위에서 설명한 오일러 방법과 같다.

[eqn: {bar{x}} `` == `x(t+ {Delta t} over {2})=x(t)+ {Delta t} over {2} f(x)# x(t+ Delta t)=x(t)+ Delta tf( {bar{x}} )]

직관적으로 룬제-쿠타 방법을 오일러 방법과 비교해 보자. 오일러 방법은 물체의 속도를 [eqn: t]에서 구하고, 이 속도를 가지고 물체가 [eqn: t+dt]까지 같은 속도로 등속 운동한다고 가정하는 것에 해당한다. 한편, 2차 룬제-쿠타 방법에서는 입자의 속도를 [eqn: t+ dt over 2]에서 구하고 이 속도를 가지고 입자가 [eqn: t]에서 [eqn: t+ dt]까지 등속 운동한다고 가정하는 것에 해당한다. 위에 적은 룬제-쿠타 방법의 수식의 좌변과 우변을 각각 테일러 전개해 비교하면 오차가 [eqn: O( Delta t ^{3} )]임을 쉽게 보일 수 있다(꼭 계산해서 확인해 보기를 추천한다). 즉, 위에서 소개한 룬제-쿠타 방법은 2차 방법(2nd-order method)이다.

룬제-쿠타 방법말고도 자주 이용하는 2차 방법이 있다. 수정-오일러 방법(modified Euler method)이라고 부른다. 수정-오일러 방법에서는 물체의 속도를 [eqn: t+ dt over 2]에서 구하고 이 속도로 물체가 등속운동한다고 가정하는 룬제-쿠타 방법과 달리, [eqn: t]에서의 속도와 [eqn: t+dt]에서의 속도의 평균값을 가지고 물체가 [eqn: t]에서 [eqn: t+dt]까지 등속운동한다고 가정하는 방식이다. 이를 식으로 적으면 아래와 같다.

[eqn: {bar{x}} =x(t)+ Delta tf(x)# x(t+ Deltat)=x(t)+ Deltat {f(x)+f( {bar{x}} )} over {2}]

수정-오일러 방법도 마찬가지로 2차 방법이라는 것을 위 식의 양변을 테일러 전개해 직접 확인해보기를 추천한다. 자, 이제 위에서 오일러 방법을 이용해 풀었던 공기 중에서 저항력을 받으며 낙하하는 물체의 운동을 수정-오일러 방법을 써서 프로그램으로 구현해보자.

import numpy as np, matplotlib.pyplot as plt

g = 9.8; gamma=1.0; dt = 0.05; maxt = 10.0

v = 0.0

t_arr = []

v_arr = []

for t in np.arange(0, maxt, dt):

t_arr.append(t)

v_arr.append(v)

a = g - gamma*v

vtem = v + dt*a

atem = g - gamma*vtem

v += dt*(a + atem)/2.0

plt.plot(t_arr,v_arr)

plt.show()

== 예제 2 ==

수정-오일러 방법을 이용한 위 프로그램을 이용해서 [eqn: v= {g} over {gamma } (1–e ^{- gamma ```t} )]와 수치적분의 결과의 차이(오차)의 절대값을 시간의 함수로 그리시오. dt를 0.001, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5로 바꾸면서 각각의 dt에서의 오차의 최댓값을 구하고, 이를 dt의 함수로 그리시오. 위에서 살펴본 오일러 방법을 이용해 그린 같은 그래프와 비교해 수정-오일러 방법의 장점에 대해 생각해 보시오.

== 예제 3 ==

2차 룬제-쿠타 방법을 적용해 위의 예제를 따라하고 오일러 방법, 2차 룬제-쿠타 방법, 수정-오일러 방법을 비교하시오.

== 과제 1 ==

교재 홈페이지에서 “mers.txt”파일을 다운 받으시오. 이 파일에는 과거 우리나라에서 전염병 MERS가 전파될 때의, 하루 하루의 누적환자수가 들어있다. 최종적으로 [eqn: K]명의 환자가 전염되었다고 하면, 주어진 시간 [eqn: t]에서 병에 걸린 [eqn: N(t)]명의 환자가 하루에 전염시키는 환자의 숫자는 [eqn: N(1-N/K)]에 비례하게 된다. [eqn: N]명의 환자가 아직 병에 걸리지 않은 환자를 만날 확률은 [eqn: (1 - N/K)]에 비례하기 때문이다. 이로부터 하루하루 어떻게 누적환자수가 늘어나는 지를 기술하는 미분방정식은 [eqn: dN over dt = r N left( 1 - N over K right)]로 적을 수 있다. [eqn: N(t=0) = 1]을 초기조건으로 하고,

[eqn: r = 0.25], [eqn: K = 186]으로 해서, 이 미분방정식을 2차 룬제-쿠타 방법을 적용해 수치적분하는 프로그램을 작성해서 그 결과를 mers.txt에 담긴 실제의 데이터와 함께 그리시오.

== 파이선 요약 ==

g = 9.8; gamma=1.0; dt = 0.05  # 여러 줄의 코드를 연결해 한 줄로 적을 때에는 세미콜론(;)을 이용한다.

np.arange(a,b,d) # a에서 b까지 d의 간격으로 나열된 숫자들의 배열을 얻는다.

np.max(np.abs(error)) # 배열 error의 각 요소의 절댓값으로 이루어진 배열을 얻고 그 안에서 가장 큰 값을 돌려준다. 이 장의 예제에서 오차의 절댓값의 최댓값을 얻을 때 유용하게 이용할 수 있다.