SciPy ODE(상미분 방정식) 풀이

SciPy ODE(상미분 방정식) 풀이

상미분 방정식(Ordinary Differential Equation, ODE)은 미지의 함수와 그 함수의 도함수 간의 관계를 나타내는 수학적 방정식입니다. 과학, 공학, 경제학 등 다양한 분야에서 시스템의 동적 거동을 설명하는 데 사용되며, 미분 방정식을 푸는 과정은 이러한 시스템을 이해하고 예측하는 데 핵심 역할을 합니다.

Python의 SciPy 라이브러리는 scipy.integrate 모듈을 통해 ODE를 효과적으로 풀 수 있는 다양한 방법을 제공합니다. 이 포스팅에서는 SciPy의 solve_ivp 함수를 중심으로 상미분 방정식을 푸는 방법을 단계별로 설명하겠습니다.

---

1. 상미분 방정식(ODE)란 무엇인가?

상미분 방정식(ODE)은 하나의 독립 변수에 대해 한 개 이상의 종속 변수와 그 도함수 간의 관계를 설명하는 방정식입니다. 다음과 같은 형태로 표현됩니다.

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

여기서:

• $ t $는 독립 변수입니다.
• $ y $는 종속 변수입니다.
• $ \frac{dy}{dt} $는 $ y $의 $ t $에 대한 1차 도함수입니다.
• $ f(t, y) $는 주어진 시스템의 변화율을 정의하는 함수입니다.

예: 간단한 1차 ODE

가장 간단한 형태의 ODE는 다음과 같이 표현됩니다.

$$
\frac{dy}{dt} = -ky
$$

여기서 $ k $는 상수입니다. 이 식은 감쇠되는 시스템(예: 라디오 활성 붕괴)을 설명할 때 자주 사용됩니다.

---

2. SciPy를 이용한 ODE 풀이 방법

SciPy는 ODE를 풀기 위해 scipy.integrate.solve_ivp 함수를 제공합니다. 이 함수는 초기값 문제(initial value problem)를 풀 때 사용되며, 다양한 수치적 방법을 선택할 수 있습니다.

solve_ivp 함수의 기본 형식

scipy.integrate.solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, args=())

• fun: 미분 방정식의 우변을 정의하는 함수 $ f(t, y) $.
• t_span: 해를 구할 시간 구간 $[t_0, t_f]$.
• y0: 초기값 $ y(t_0) $.
• method: 수치적 방법(기본값: 'RK45').
• t_eval: 특정 시간에서의 해를 원하는 경우 시간 배열.
• args: 추가 인자 전달.

---

3. 예제 1: 단순 1차 미분 방정식 풀이

문제 설명

다음 1차 미분 방정식을 풀이해보겠습니다.

$$
\frac{dy}{dt} = -2y, \quad y(0) = 5
$$

코드 구현

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# 미분 방정식 정의
def model(t, y):
    return -2 * y

# 초기 조건 및 시간 구간 설정
t_span = (0, 5)  # t=0에서 t=5까지
y0 = [5]         # 초기값 y(0) = 5
t_eval = np.linspace(0, 5, 100)

# 미분 방정식 풀이
sol = solve_ivp(model, t_span, y0, method='RK45', t_eval=t_eval)

# 결과 시각화
plt.figure(figsize=(10, 6))
plt.plot(sol.t, sol.y[0], label='y(t)')
plt.xlabel('Time t')
plt.ylabel('y(t)')
plt.title('1차 ODE: dy/dt = -2y')
plt.legend()
plt.grid(True)
plt.show()

결과 해석

이 방정식의 해는 다음과 같이 주어집니다.

$$
y(t) = y(0) e^{-2t}
$$

시각화된 결과에서, 초기값 5에서 시작하여 시간이 지남에 따라 지수적으로 감소하는 모습을 확인할 수 있습니다.

---

4. 예제 2: 2차 ODE 풀이 (질량-스프링 시스템)

문제 설명

질량-스프링 시스템은 다음과 같은 2차 미분 방정식으로 표현됩니다.

$$
\frac{d^2 x}{dt^2} + 2 \zeta \omega_0 \frac{dx}{dt} + \omega_0^2 x = 0
$$

여기서:

• $ \zeta $: 감쇠 비율
• $ \omega_0 $: 고유 진동수

이를 1차 미분 방정식 두 개로 변환하면 다음과 같이 나타낼 수 있습니다.

$$
\begin{cases}
\frac{dx}{dt} = v \
\frac{dv}{dt} = -2 \zeta \omega_0 v - \omega_0^2 x
\end{cases}
$$

코드 구현

# 질량-스프링 시스템의 미분 방정식 정의
def mass_spring(t, y, zeta, omega0):
    x, v = y
    dxdt = v
    dvdt = -2 * zeta * omega0 * v - omega0**2 * x
    return [dxdt, dvdt]

# 파라미터 설정
zeta = 0.2    # 감쇠 비율
omega0 = 2.0  # 고유 진동수

# 초기 조건 및 시간 구간
t_span = (0, 10)
y0 = [1.0, 0.0]  # x(0)=1, v(0)=0
t_eval = np.linspace(0, 10, 1000)

# 미분 방정식 풀이
sol = solve_ivp(mass_spring, t_span, y0, args=(zeta, omega0), method='RK45', t_eval=t_eval)

# 결과 시각화
plt.figure(figsize=(10, 6))
plt.plot(sol.t, sol.y[0], label='x(t)')
plt.plot(sol.t, sol.y[1], label='v(t)', linestyle='--')
plt.xlabel('Time t')
plt.ylabel('Amplitude')
plt.title('질량-스프링 시스템의 동적 거동')
plt.legend()
plt.grid(True)
plt.show()

결과 해석

위 그래프에서 위치 $ x(t) $는 시간이 지남에 따라 점점 감쇠하는 진동을 보입니다. 이는 감쇠된 질량-스프링 시스템에서 흔히 관찰되는 현상입니다.

---

5. 다양한 수치적 방법

solve_ivp는 다양한 수치적 방법을 지원합니다.

• RK45: 기본값으로 많이 사용되는 4차-5차 Runge-Kutta 방법.
• RK23: 2차-3차 Runge-Kutta 방법.
• BDF: 강성(stiff) 방정식에 적합한 방법.
• LSODA: 자동으로 강성 여부를 판단하여 적절한 방법 선택.

예를 들어, BDF 방법을 사용하는 경우 다음과 같이 설정할 수 있습니다.

sol = solve_ivp(model, t_span, y0, method='BDF', t_eval=t_eval)

---

6. 결론

SciPy의 solve_ivp를 활용하면 상미분 방정식을 효과적으로 풀 수 있습니다. 간단한 1차 방정식부터 복잡한 다변수 시스템까지 다양한 문제를 다룰 수 있으며, 수치적 방법을 적절히 선택하면 정확성과 속도를 모두 만족시킬 수 있습니다.