Python에서 ode45() 함수 모방

MATLAB의 ode45() 함수는 적응형 룽게-쿠타 방법으로 상미분방정식(ODE)의 초깃값 문제를 풉니다. Python에는 문자 그대로 ode45()라는 이름의 함수는 없지만, SciPy의 solve_ivp()는 같은 작업 흐름을 다룹니다. 방정식을 정의하고, 솔버에 시간 구간과 시작 상태를 지정한 뒤, 반환된 객체에서 해를 읽으면 됩니다.

MATLAB의 ode45()를 대신해 SciPy의 solve_ivp() 사용

scipy.integrate의 solve_ivp() 함수는 ode45()에 가장 가까운 실용적인 대응 함수입니다. 기본 방법인 RK45는 자동 스텝 크기 제어 기능을 갖춘 5(4)차 명시적 룽게-쿠타 공식으로, MATLAB에서 ode45()가 사용하는 알고리즘 계열과 같습니다. 따라서 해가 천천히 변하는 구간에서는 큰 스텝을 사용하고, 빠르게 변하는 구간에서는 작은 스텝을 사용합니다.

그러나 그대로 교체해 쓸 수 있는 함수는 아닙니다. 인수 이름이 다르고, 결과도 t와 y 배열 한 쌍이 아니라 하나의 객체로 반환됩니다. 그래도 도함수 함수, 시간 범위, 초기 상태라는 전체적인 패턴은 두 언어 간에 그대로 대응합니다. SciPy가 아직 설치되어 있지 않다면 pip install scipy로 설치하세요.

ODE를 정의하고 solve_ivp()을 호출하기

핵심 호출 형태는 solve_ivp(fun, t_span, y0, method="RK45", t_eval=None, args=None, ...)입니다. 주요 입력은 다음과 같습니다.

• fun: 현재 시간 t와 상태 y에 대한 도함수 dy/dt를 반환하는 콜러블 fun(t, y)입니다. 단일 방정식이면 하나의 값을 반환하고, 연립 방정식이면 각 상태 성분에 해당하는 도함수 하나씩을 포함한 배열을 반환합니다.
• t_span: 통합 구간의 시작과 끝을 지정하는 (t_start, t_end) 쌍입니다.
• y0: t_start에서의 상태로, 이 배열의 길이가 시스템의 방정식 개수를 정합니다.
• method (선택 사항): 통합 방식이며 기본값은 'RK45'입니다.
• args (선택 사항): fun(t, y, *args) 형태로 함수에 전달할 추가 값입니다.
• rtol 및 atol (선택 사항): 적응형 스텝 크기를 제어하는 상대 허용 오차와 절대 허용 오차입니다.

다음 예제는 y(0) = 1에서 시작해 구간 [0, 2]에서 dy/dt = 2y를 풉니다. 이 호출은 결과 객체를 반환하며, t와 y 속성을 출력하면 솔버가 선택한 시점과 각 시점의 상태를 확인할 수 있습니다. 이 섹션의 코드와 출력에 표시되는 t: 및 y: 레이블은 영어 그대로 유지됩니다.

이 섹션의 코드와 표시되는 레이블 또는 출력은 영어로 유지됩니다.

from scipy.integrate import solve_ivp


def growth(t, y):
    return 2 * y


result = solve_ivp(growth, (0, 2), [1])
print("t:", result.t)
print("y:", result.y[0])

t: [0.         0.07580098 0.56590702 1.18757739 1.85907272 2.        ]
y: [ 1.          1.16369695  3.10133569 10.75272756 41.18553629 54.59497095]

보고된 시점의 간격이 불균일한 이유는 solve_ivp()이 고정된 증가분으로 진행하는 대신 해의 국소적인 변화에 맞춰 각 스텝을 조정하기 때문입니다.

솔버 결과 읽기

solve_ivp()은 scipy.integrate.OdeResult 객체를 반환합니다. 가장 자주 사용하는 속성은 다음과 같습니다.

• t: 보고된 시점을 담은 NumPy 배열입니다. t_eval를 지정하지 않으면 솔버의 내부 적응형 스텝이 들어 있습니다.
• y: 모양이 (n, len(t))인 2차원 NumPy 배열입니다. i번째 행에는 보고된 각 시점에서의 상태 성분 i가 들어 있으므로, y[:, k]는 t[k]에서의 상태입니다. 단일 방정식에서는 y[0]가 해 곡선입니다.
• sol: 구간 내 임의의 시점에서 해를 평가하기 위해 result.sol(t)처럼 호출할 수 있는 연속 OdeSolution 보간 객체입니다. dense_output=True를 전달한 경우에만 채워지며, 그렇지 않으면 None입니다.
• success와 message: 통합이 완료되었는지 나타내는 부울 값과 종료 방식을 사람이 읽을 수 있게 설명하는 메시지입니다.
• status: 구간이 완료되면 0, 종료 이벤트가 발생하면 1, 솔버가 진행하지 못하면 -1입니다.

y를 사용하기 전에 success(또는 status)를 확인하는 것이 좋습니다. 실패한 실행도 부분적인 배열을 반환할 수 있기 때문입니다.

t_eval로 출력 시점 선택하기

적응형 스텝 시점은 보고하고자 하는 시점과 거의 일치하지 않습니다. 정렬된 시간 배열을 t_eval에 전달하면 solve_ivp()은 내부 스텝 대신 정확히 그 시점들에서의 해를 반환합니다. 이 섹션의 코드와 출력에 표시되는 t: 및 y: 레이블은 영어 그대로 유지됩니다.

이 섹션의 코드와 표시되는 레이블 또는 출력은 영어로 유지됩니다.

from scipy.integrate import solve_ivp


def growth(t, y):
    return 2 * y


times = [0.5, 1.0, 1.5, 2.0]
result = solve_ivp(growth, (0, 2), [1], t_eval=times)
print("t:", result.t)
print("y:", result.y[0])

t: [0.5 1.  1.5 2. ]
y: [ 2.71823524  7.38683161 20.07910026 54.59497095]

이제 t에는 요청한 시점만 들어 있고, y에는 각 시점에서 보간된 상태가 들어 있습니다. t_eval의 모든 값은 t_span 안에 있어야 합니다. 통합 구간 밖의 시점을 요청하면 값이 t_span 내에 있어야 한다는 오류와 함께 호출이 중단됩니다. 끝점은 허용되므로 t_eval에는 t_start와 t_end를 모두 포함할 수 있습니다.