파이썬 Sympy 기호수학 - 응용 1. 상미분방정식

SymPy 사용법¶

SymPy Tutorial

sympy 라이브러리를 불러오고, 사용할 기호 변수를 선언한다. 출력을 LaTex 수식으로 나타내고, 맷플롯립 모듈을 불러온다.

In [1]:

from sympy import *

x, y, z, t = symbols('x y z t')
f, g, h = symbols('f, g, h', cls=Function)

init_printing()
%matplotlib inline 

상미분 방정식¶

미분 방정식의 해 구하기 $$y'(t) = -2y , \qquad y(0)=0.5$$

y 를 t의 함수 y(t) 로 사용하려면, 함수로 선언하여야 한다.

In [2]:

y = symbols('y', cls=Function)
y(t)

Out[2]:

$\displaystyle y{\left(t \right)}$

In [3]:

deq = Eq( y(t).diff(t), -2*y(t) )
deq

Out[3]:

$\displaystyle \frac{d}{d t} y{\left(t \right)} = - 2 y{\left(t \right)}$

미분방정식을 dsolve() 함수로 풀면 일반해가 얻어진다.

In [4]:

dsolve( deq, y(t) )

Out[4]:

$\displaystyle y{\left(t \right)} = C_{1} e^{- 2 t}$

초기조건을 ics= {...} 으로 전달하면, 초기조건이 적용된 특수해(particular solution)를 얻는다.

In [5]:

psol = dsolve( deq, ics= {y(0):0.5} )
psol

Out[5]:

$\displaystyle y{\left(t \right)} = 0.5 e^{- 2 t}$

윗 식의 오른쪽 부분인 .rhs 가 해이며, 이를 그래프로 그린다.

In [6]:

plot( psol.rhs, (t,0,2))

Out[6]:

<sympy.plotting.plot.Plot at 0x24cee24ba90>

[예제] 비제차 미분방정식

$$y' + y \tan x = \sin 2x , \qquad y(0)= 1$$

y 를 x의 함수 y(x) 로 선언하고, 미분 방정식을 기술한다.

In [7]:

y(x)
deq = Eq( y(x).diff(x) + y(x)*tan(x), sin(2*x) )
deq

Out[7]:

$\displaystyle y{\left(x \right)} \tan{\left(x \right)} + \frac{d}{d x} y{\left(x \right)} = \sin{\left(2 x \right)}$

dsolve() 함수로 일반해를 얻는다.

In [8]:

dsolve( deq, y(x) )

Out[8]:

$\displaystyle y{\left(x \right)} = \left(C_{1} - 2 \cos{\left(x \right)}\right) \cos{\left(x \right)}$

초기조건을 ics= {...} 로 설정하면, 초기조건이 적용된 특정해를 얻는다.

In [9]:

psol = dsolve( deq, ics= {y(0):1} )
psol

Out[9]:

$\displaystyle y{\left(x \right)} = \left(3 - 2 \cos{\left(x \right)}\right) \cos{\left(x \right)}$

전개하여 보기 좋게 만든다.

In [10]:

expand( psol.rhs )

Out[10]:

$\displaystyle - 2 \cos^{2}{\left(x \right)} + 3 \cos{\left(x \right)}$

2계 미분방정식¶

$$ 10 y'' + 10 y' + 90 y = 0 , \qquad y(0)= 0.16, \quad y'(0)=0 $$

함수 y(t) 를 선언하고, 미분 방정식을 기술한다.

In [11]:

y(t)
deq = Eq( 10*y(t).diff(t,2) + 10*y(t).diff(t) + 90*y(t), 0 )
deq

Out[11]:

$\displaystyle 90 y{\left(t \right)} + 10 \frac{d}{d t} y{\left(t \right)} + 10 \frac{d^{2}}{d t^{2}} y{\left(t \right)} = 0$

dsolve() 함수의 인자에 초기조건을 ics= {...} 로 전달하여, 초기조건이 적용된 특정해를 얻는다.

In [12]:

sol = dsolve( deq, ics= { y(0):0.16, y(t).diff(t).subs(t,0):0 } )
sol

Out[12]:

$\displaystyle y{\left(t \right)} = \left(0.0270449361513125 \sin{\left(\frac{\sqrt{35} t}{2} \right)} + 0.16 \cos{\left(\frac{\sqrt{35} t}{2} \right)}\right) e^{- \frac{t}{2}}$

In [13]:

plot( sol.rhs, xlim=(0,8), ylim=(-0.2,0.2) )

Out[13]:

<sympy.plotting.plot.Plot at 0x24cefdfa3a0>

2계 비제차 미분방정식¶

$$ y'' + 2 y' + 0.75 y = 2 \cos x - 0.25 \sin x + 0.09 x , \qquad y(0)= 2.78, \quad y'(0)=-0.43 $$

함수 y(x) 를 선언하고, 미분 방정식을 기술한다.

In [14]:

y(x)
deq = Eq( y(x).diff(x,2) + 2*y(x).diff(x) + 0.75*y(x), 2*cos(x)-0.25*sin(x)+0.09*x )
deq

Out[14]:

$\displaystyle 0.75 y{\left(x \right)} + 2 \frac{d}{d x} y{\left(x \right)} + \frac{d^{2}}{d x^{2}} y{\left(x \right)} = 0.09 x - 0.25 \sin{\left(x \right)} + 2 \cos{\left(x \right)}$

dsolve() 함수의 인자에 초기조건을 ics= {...} 로 전달하여, 초기조건이 적용된 특정해를 얻는다.

In [15]:

sol = dsolve( deq, ics= { y(0):2.78, y(x).diff(x).subs(x,0):-0.43 } )
sol

Out[15]:

$\displaystyle y{\left(x \right)} = 0.12 x + 1.0 \sin{\left(x \right)} - 0.32 + 3.1 e^{- 0.5 x}$

In [16]:

plot( sol.rhs, (x,0,20) )

Out[16]:

<sympy.plotting.plot.Plot at 0x24cefb0ae50>

연립 미분 방정식¶

\begin{align} {y_1} ' &= -0.02 y_1 + 0.02 y_2 \\ {y_2} ' &= 0.02 y_1 - 0.02 y_2 \end{align}$$ \qquad \qquad \qquad \qquad \qquad 초기조건 \qquad y_1 (0) = 150, \quad y_2 (0) = 0 $$

함수 $y_1(t), y_2 (t)$ 를 선언한다.

In [17]:

y1, y2 = symbols('y1 y2', cls=Function)
y1(t)
y2(t)

Out[17]:

$\displaystyle \operatorname{y_{2}}{\left(t \right)}$

미분 방정식을 기술한다.

In [18]:

eq1 = Eq( y1(t).diff(t), -0.02*y1(t)+0.02*y2(t) )
eq1

Out[18]:

$\displaystyle \frac{d}{d t} \operatorname{y_{1}}{\left(t \right)} = - 0.02 \operatorname{y_{1}}{\left(t \right)} + 0.02 \operatorname{y_{2}}{\left(t \right)}$

In [19]:

eq2 = Eq( y2(t).diff(t), 0.02*y1(t)-0.02*y2(t) )
eq2

Out[19]:

$\displaystyle \frac{d}{d t} \operatorname{y_{2}}{\left(t \right)} = 0.02 \operatorname{y_{1}}{\left(t \right)} - 0.02 \operatorname{y_{2}}{\left(t \right)}$

dsolve() 함수로 일반해를 얻는다.

In [20]:

dsolve( [ eq1, eq2 ] )

Out[20]:

$\displaystyle \left[ \operatorname{y_{1}}{\left(t \right)} = 1.0 C_{1} - 1.0 C_{2} e^{- 0.04 t}, \ \operatorname{y_{2}}{\left(t \right)} = 1.0 C_{1} + 1.0 C_{2} e^{- 0.04 t}\right]$

dsolve() 함수에 초기조건을 ics= {...} 으로 전달하여, 특정해를 얻는다.

In [21]:

sol = dsolve( [ eq1, eq2 ], ics= { y1(0):150, y2(0):0 } )
sol

Out[21]:

$\displaystyle \left[ \operatorname{y_{1}}{\left(t \right)} = 75.0 + 75.0 e^{- 0.04 t}, \ \operatorname{y_{2}}{\left(t \right)} = 75.0 - 75.0 e^{- 0.04 t}\right]$

그래프로 그린다.

In [22]:

plot( sol[0].rhs, sol[1].rhs , (t,0,100) )

Out[22]:

<sympy.plotting.plot.Plot at 0x24cefeae220>

비제차 연립 미분 방정식¶

\begin{align} {y_1} ' &= -3 y_1 + y_2 - 6 e ^ {-2t} \\ {y_2} ' &= \quad y_1 -3 y_2 + 2 e ^ {-2t} \end{align}

함수 $y_1(t), y_2 (t)$ 를 선언한다.

In [23]:

y1, y2 = symbols('y1 y2', cls=Function)
y1(t)
y2(t)

Out[23]:

$\displaystyle \operatorname{y_{2}}{\left(t \right)}$

미분 방정식을 기술한다.

In [24]:

eq1 = Eq( y1(t).diff(t), -3*y1(t)+y2(t)-6*exp(-2*t) )
eq1

Out[24]:

$\displaystyle \frac{d}{d t} \operatorname{y_{1}}{\left(t \right)} = - 3 \operatorname{y_{1}}{\left(t \right)} + \operatorname{y_{2}}{\left(t \right)} - 6 e^{- 2 t}$

In [25]:

eq2 = Eq( y2(t).diff(t), y1(t)-3*y2(t)+2*exp(-2*t) )
eq2

Out[25]:

$\displaystyle \frac{d}{d t} \operatorname{y_{2}}{\left(t \right)} = \operatorname{y_{1}}{\left(t \right)} - 3 \operatorname{y_{2}}{\left(t \right)} + 2 e^{- 2 t}$

dsolve() 함수로 일반해를 얻는다.

In [26]:

dsolve( [ eq1, eq2 ] )         

Out[26]:

$\displaystyle \left[ \operatorname{y_{1}}{\left(t \right)} = - C_{2} e^{- 4 t} - 2 t e^{- 2 t} + \left(C_{1} - 2\right) e^{- 2 t}, \ \operatorname{y_{2}}{\left(t \right)} = C_{2} e^{- 4 t} - 2 t e^{- 2 t} + \left(C_{1} + 2\right) e^{- 2 t}\right]$

last updated in 2022.3