SymPy를 이용하여 미분방정식 해결하기

SymPy 란?

SymPy는 심볼릭 수학을 위한 Python 라이브러리입니다. 가능한 한 이해하기 쉽고 확장하기 편하게 컴퓨터 대수 시스템(CAS)을 지원합니다. 라이브러리는 Python으로 작성되었습니다. (공식 홈페이지 소개 번역)

가령 인수분해, 부정적분을 구하는 문제 등은 수치해석 라이브러리인 SciPy 등으로 해결할 수 없는 문제입니다. SymPy는 이러한 문제에 '마치 손으로 푸는 것과 같은' 솔루션을 제시합니다.

SymPy 설치하기

본 글에서는 ODE만을 다룹니다.

SymPy는 내장 라이브러리가 아닙니다. 따라서 이를 설치해야 합니다.

• pip 환경에서

pip install sympy

• conda 환경에서

conda install sympy

SymPy로 미분방정식 해결하기

1. 모듈 로드

SymPy 모듈을 로드합니다.

import sympy as sp

2. 가독성 개선

보기 편하게 하기 위해 use_latex를 사용하도록 하겠습니다.

sp.init_printing(use_latex=True)

use_unicode를 사용해도 되고, 사용하지 않아도 가독성이 감소하는 것을 제외하면 실행하지 않아도 되는 코드입니다.

3. 함수 설정

미분방정식에 사용할 함수를 만들도록 하겠습니다.

함수는 매개변수와 종속변수간의 관계를 나타냅니다. 본 글에서는 한 개의 매개변수에 대한 함수만을 다룹니다. 일반적으로 통용되는 미지수 x와 y를 각각 매개변수와 종속변수로, f를 함수로 사용하겠습니다.

SymPy에서는 변수를 나타내기 위해 일련의 과정을 거쳐야 합니다. 변수는 상수, 변수, 함수 등 모든 것을 포괄하는 개념입니다. symbols 메서드를 이용해 변수를 정의할 수 있습니다.

x = sp.symbols('x')
f = sp.symbols('f', cls=sp.Function) # 함수임을 나타냄

이때 종속변수 y는 함수 f에 x를 대입했을 때의 값이므로 다음과 같이 코드를 입력합니다.

y = f(x)

4. 미분방정식 풀이

이제 미분방정식을 풀어보도록 하겠습니다.

다음 선형 미분방정식을 해결하기 위해, 먼저 방정식 자체를 정의하도록 하겠습니다. Equation 메서드는 첫 번째 인자에 좌변의 식을, 두 번째 인자에 우변의 식을 넣어 등식을 정의할 수 있습니다.

eq = sp.Eq(x * sp.diff(y, x) + y, 2 * x)

이 때 sp.diff(y, x)는 f의 x에 대한 도함수입니다. 보다 직관적인 표기를 위해 sp.diff(f(x), x)를 사용할 수도 있으나, 가독성의 문제로 이전에 y를 정의하였습니다.

이제 등식(미분방정식)이 정의되었으므로, 이를 해결합니다. dsolve 함수는 미분방정식의 해를 제공합니다.

sp.dsolve(eq, y)
# >> f(x)=C1/x+x

이 또한 마찬가지로 미분방정식 eq를 y, 즉 f(x)에 관해 풀라는 의미입니다.

이번엔 고계 미분방정식을 풀어보도록 하겠습니다.

eq = sp.Eq(sp.diff(y, x, 2) + y, sp.sec(x))
sp.dsolve(eq, y)
# >> f(x)=(C1+x)sin(x)+(C2+log(cos(x))cos(x))

sp.diff에 세 개의 인자가 들어가면 첫 번째 인자(y)를 두 번째 인자(x)에 대해 세 번째 인자(2)에 해당하는 횟수의 미분, 즉 이계 도함수를 의미합니다. sec함수는 SymPy에서 제시하는 함수로, 어떤 SymPy 변수를 인자로 받는 함수입니다. 즉, sp.sec(x)는 SymPy의 식을 나타냅니다.

5. 초기조건 설정

때론 미분방정식의 초기 조건이 제시되기도 합니다.

그 경우 dsolve에 초기 조건에 대한 식을 전달하여 해를 명확히 할 수도 있습니다.

첫 번째 미분방정식에 초기조건 $y(1)=2$을 적용해 봅시다.

sp.dsolve(eq, y, ics={f(1): 2})
# >> f(x)=x+1/x

dsolve에 ics 인자에 좌변을 key로, 우변을 value로 한 dict 개체를 전달할 수 있습니다. 여기서 $y(1)$은 f에 1을 대입한 값이므로 y(1)이 아닌 f(1)로 입력해야 함이 마땅합니다.

두 번째 미분방정식은 이계 미분방정식이므로 두 개의 초기조건 $y(0)=0$, $y'(0)=0$을 적용해 봅시다.

sp.dsolve(eq, y, ics={f(0): 0, sp.diff(y, x).subs(x, 0): 0})
# >> f(x)=xsin(x)+log(cos(x))cos(x)

subs 메서드는 식의 해당하는 미지수(첫 번째 인자; x)에 값(두 번째 인자; 0)를 대입한 식을 반환합니다. 따라서 sp.diff(y, x).subs(x, 0)은 $y'\big\vert_{x=0}=y'(0)$을 의미합니다.

Note

symbols 메서드는 여러개의 변수를 한 번에 정의할 수 있습니다. 가령,

x, y = sp.symbols('x y')

과 같이 사용 가능합니다.