(연립) 미분 방정식의 수치적 해법 (2) : 룽게-쿠타 방법(Runge-Kutta Method)과 로트카-볼테라 방정식(Lotka-Volterra Equation)

이 글에서는 미분 방정식의 수치적 해법 중 룽게-쿠타 방법(Runge-Kutta Method)에 대해 설명하겠습니다.

---

룽게-쿠타 방법(Runge-Kutta Method)

룽게-쿠타 방법은 미분 방정식의 수치적 해법 중 하나로서, 2차 룽게-쿠타 방법, 4차 룽게-쿠타 방법 등이 있습니다. 여기서는 주로 많이 쓰이는 4차 룽게-쿠타 방법(RK4)에 대해 다루도록 하겠습니다.

 

룽게-쿠타 방법은 다변수 함수의 테일러 전개로부터 유도되는 방법입니다. 하지만, 이 글에서는 룽게-쿠타 방법의 유도는 다루지 않고, 간단히 4차 룽게-쿠타 방법의 알고리즘만 소개하도록 하겠습니다.

 

우리가 풀어야 할 미분 방정식을 아래와 같이 설정합니다.

$$\frac{dy}{dt}=f(t, y),\quad y(t_0)=y_0$$

$t_{i+1}=t_i+h$로 설정하고, 아래와 같이 값을 정의할 수 있습니다.

$$k_1=f(t_i, y_i)$$

$$k_2=f(t_i+\frac{1}{2}h, y_i+\frac{1}{2}hk_1)$$

$$k_3=f(t_i+\frac{1}{2}h, y_i+\frac{1}{2}hk_2)$$

$$k_4=f(t_i+h, y_i+hk_3)$$

위와 같이 $k_1, k_2, k_3, k_4$를 정의하면, 아래와 같은 근사를 사용할 수 있습니다.

$$y_{i+1}=y_i+\frac{1}{6}h(k_1+2k_2+2k_3+k_4)$$

이와 같은 근사를 사용하여 함수값을 찾는 방법이 4차 룽게-쿠타 방법(RK4)입니다.

---

연립 미분 방정식의 풀이

미지수가 여러개인 미분 방정식을 여러개 모아놓은 system을 연립 미분 방정식이라고 합니다. 연립 미분 방정식의 수치해법은 미지수가 1개인 일반적인 미분 방정식의 수치해법과 크게 다르지 않습니다. 

 

아래와 같은 $m$개의 미분 방정식을 생각해봅시다.

$$y^{'}(x)_{1}=f_{1}(x, y_1, ... , y_m), \quad y_1(x_0)=\alpha_1$$

$$y^{'}(x)_{2}=f_{2}(x, y_1, ... , y_m), \quad y_2(x_0)=\alpha_2$$

$$\vdots$$

$$y^{'}(x)_{m}=f_{m}(x, y_1, ... , y_m), \quad y_m(x_0)=\alpha_m$$ 

이를 벡터로 표기한다면,

$$\mathbf{y^{'}}(x)=\mathbf{f}(x, \mathbf{y}(x)), \quad \mathbf{y}(x)=\mathbf{\alpha}$$

와 같이 표기할 수 있습니다. 이때

$$\mathbf{y}(x)=\left \{\begin{array}{cc} y_{1}(x) \\ y_{2}(x) \\ \vdots \\ y_{m}(x) \end{array} \right \}, \; \mathbf{f}=\left \{\begin{array}{cc} f_{1}(x, y_1, ... , y_m) \\ f_{2}(x, y_1, ... , y_m) \\ \vdots \\ f_{m}(x, y_1, ... , y_m) \end{array} \right \}, \; \mathbf{\alpha}=\left \{\begin{array}{cc} \alpha_1 \\ \alpha_2 \\ \vdots \\ \alpha_m \end{array} \right \}$$

입니다. 이렇게 된 식은 벡터가 들어간 것을 제외하면 미지수가 1개인 미분방정식과 같은 형태이므로, 벡터의 형태로 원래 우리가 알고 있는 수치적 방법을 적용하면 됩니다. 파이썬의 경우에는, float형 변수 자리에 float형 array를 사용하면 되는 것이죠.

---

로트카-볼테라 방정식(Lotka-Volterra Equation)

어느 생태계에서 포식자와 피식자 사이 관계를 생각해봅시다. 피식자의 수가 줄어들면, 포식자는 먹이가 줄어들어 포식자 개체수 또한 감소하게 됩니다. 그러면 피식자의 개체수가 다시 늘고, 이에 따라 포식자의 수 또한 다시 증가하고, 또 피식자의 수가 줄어들고..... 이렇게 포식자와 피식자의 개체수는 주기적으로 반복되지 않을지 의문을 가질 수 있습니다. 이에 대한, 포식자와 피식자의 개체수에 관한 미분 방정식이 바로 로트카-볼테라 방정식(Lotka-Volterra Equation)입니다. 

 

피식자의 개체수를 $x$, 포식자의 개체수를 $y$라 합니다. 또한 피식자와 포식자 사이 상호작용으로 인한 개체수 변동에 대한 상수를 각각 $b, d$로, 피식자와 포식자 스스로에 의한 개체수 변동에 대한 상수를 각각 $a, c$로 놓습니다. 피식자와 포식자 스스로에 의한 개체수 변동은 피식자와 포식자의 개체수 $x, y$에 비례할 것이고, 피식자와 포식자 사이 상호작용으로 인한 개체수 변동은 포식자와 피식자가 마주치는 횟수, 즉 $x$와 $y$의 곱에 비례할 것임을 알 수 있습니다. 이러한 사실을 바탕으로, 미분 방정식을 기술해보면 아래와 같은 결과를 얻을 수 있습니다.

$$\frac{dx}{dt}=ax-bxy$$

$$\frac{dy}{dt}=-cx+dxy$$

이 연립 비선형 상미분 방정식이 바로 로트카-볼테라 방정식(Lotka-Volterra Equation)입니다.

---

코드 작성

이제 4차 룽게-쿠타 방법을 이용하여 로트카-볼테라 방정식을 풀어보도록 하겠습니다.

우선 순서도를 그려보았습니다.

위 순서도를 바탕으로 아래 코드를 작성해보았습니다.

import numpy as np
import matplotlib.pyplot as plt

a = 3
b = 2
c = 1
d = 2

def lveq(t, y):
    dydt0 = a*y[0] - b*y[0]*y[1]
    dydt1 = -c*y[1] + d*y[0]*y[1]
    return np.array([dydt0, dydt1])

def RK4(f, h, t0, tmax, y0):
    n = int((tmax - t0)/h)
    t = t0
    y = y0
    m = len(y0)
    tsol = np.zeros(n+1)
    ysol = np.zeros([n+1, m])
    tsol[0] = t
    ysol[0, :] = y
    for i in range(n):
        k1 = f(t, y)
        k2 = f(t+0.5*h, y+0.5*h*k1)
        k3 = f(t+0.5*h, y+0.5*h*k2)
        k4 = f(t+h, y+h*k3)
        y = y + h*(k1 + 2*k2 + 2*k3 + k4)/6
        t = t + h
        tsol[i+1] = t
        ysol[i+1, :] = y
    return tsol, ysol

h = 0.01
t0, tmax = 0, 20
y0 = np.array([1.9, 1.7])
t, y = RK4(lveq, h, t0, tmax, y0)
plt.plot(t, y[:,0], label = r'$x$'' (prey)')
plt.plot(t, y[:,1], label = r'$y$'' (predator)')
plt.xlabel(r'$t$')
plt.ylabel(r'$x, y$'' (population)')
plt.legend()
plt.title("Lotka-Volterra Equation")
plt.show()

아래는 코드에 대한 간단한 설명입니다.

풀어야 할 방정식(이 코드의 경우는 로트카-볼테라 방정식)을 입력합니다. 함수의 범위를 t0부터 tmax까지로 설정하고, 초깃값 y0를 설정해줍니다. 설정한 범위를 길이가 h인 구간 n개로 나눕니다. 위에서 설명한 4차 룽게-쿠타 방법의 알고리즘을 바탕으로 k1, k2, k3, k4를 계산하고 이를 통해 y와 t값을 계산한 후 이를 tsol, ysol에 추가해주는 작업을 n번 반복합니다. 이렇게 얻은 tsol과 ysol을 바탕으로 그래프를 그리면 코드는 끝이 납니다.

위 코드를 실행하면 아래와 같은 그래프를 얻을 수 있습니다.

파란색이 피식자, 주황색이 포식자이다.

위 로트카-볼테라 방정식의 설명과 같이, 포식자 수가 증가하여 피식자 수가 감소하고, 이에 따라 포식자 수가 감소하고 피식자 수가 다시 증가하고, 또 다시 포식자 수가 증가하는 변화가 주기적으로 반복됨을 볼 수 있습니다. 

로트카-볼테라 방정식의 예시로는 캐나다에서 눈신토끼와 스라소니의 관계를 예로 들 수 있습니다. 아래 1840년대부터 1930년대까지의 데이터를 보면, 피식자 눈신토끼와 포식자 스라소니의 관계는 로트카-볼테라 방정식을 따름을 알 수 있습니다.

https://www.ealt.ca/blog/hare-vs-lynx-cycle

그런데, 우리의 생활과 바로 밀접한 예시가 하나 더 있습니다. 눈치 채셨나요?

 

바로 코로나 19 확진자 수입니다.

로트카-볼테라 방정식은 감염병에도 적용할 수 있습니다. 피식자 대신 비감염자, 포식자 대신 감염자(전파자)를 대입하여 방정식을 적용할 수 있습니다. 그래서 실제 코로나 19 확진자 수가 계속해서 감소하거나 증가하지 않고, 아래 그래프와 같이 진동하는 것입니다.

https://www.yna.co.kr/view/GYH20220706000300044

 

이렇게 4차 룽게-쿠타 방법을 이용하여, 포식자와 피식자의 개체수에 관련된 비선형 연립 미분 방정식 로트카-볼테라 미분 방정식을 풀어보았습니다.