소개

대학의 연구실에 구르고 있던 Mathematica에서 놀았던 기억이 있습니다. 그리고 수십 년의 세월이 흐르는 동안 연구생활에서 몇번이나 미분방정식이 눈앞에 나타나, 수식처리소프트가 있다고 생각하는 일이 있었지만, 만일이 되면 수치적분 굉장히 좋다고 타락한 채 환력 눈앞까지 되어 버렸던 노 기술자입니다. 지금 연도 잘 부탁드립니다.
그런데 최근 변태 동료의 다나카 씨라는 수학자가 Sympy를 소개해 주고, 이것을 사용하고 있는 동안 젊은 날의 추억이 되살아나고, 히다 신지에 꺼낸, 아니 이 수식 처리 라이브러리를 사용해 고교 물리 의 문제에서도 풀어 보려고 생각한 대로입니다.

환경

Windows10 pro
Anaconda + Python 3.7.6 (64bit) + Sympy 1.5.1
Visual Studio Code 1.44.2

기본식

고등학교 시절, 노기술자가 처음으로 미분방정식의 위력을 알게 된 것이 공기 저항 중의 낙하 문제.

질량 $ m $ (kg)의 작은 구가 수직 아래쪽으로 떨어집니다. 소구는 시간 $0$(s)에 초속 $0$에서 낙하를 개시하고, 낙하의 도중, 소구에는 속도 $v$(m/s)에 비례하는 공기 저항 $kv$(N)이 속도 의 방향과 역방향에 작용하는 것으로서, 시각 $t$(s)에 있어서의 소구의 속도를 구해라.

옛날, 스루다이에는 사카마사라고 하는 위대한 물리의 선생님이 계시고, 미분 방정식을 온 글자로 살짝 써 주신다. 고등학교 1학년 때 오전부에 숨어 있던 저는 순식간에 흩날리는 겁니다. 하지만 어떻게든 노트에 써서, 라는 날들이 계속 맑아져 현역으로 이과 일류에, 라고 이끼는 자랑도 넣으면서, 우선은 운동 방정식을 세웁시다. 시각 $t$(s)의 소구의 속도를 $v(t)$(m/s)로서,

운동 방정식

$$ m\frac{dv\left(t\right)}{d t} {} = mg-kv{\left(t\right)}$$
이것을 풀어
$$v{\left(t\right)} =\frac{mg}{k} +\frac{e^{\frac{C_{1} k}{m}} }{k}e^{-\frac{k t}{m}}$$
상수를
$$A=\frac{e^{\frac{C_{1} k}{m}} }{k}$$
정리하면
$$v{\left(t\right)} =\frac{mg}{k} + A e^{-\frac{k t}{m}}$$
라는 일반해가 됩니다. 초기 조건 $t=0(s) 에서 v{\left(t\right)}=0$ 을 넣어 $A$를 구하면,
$$A=-\frac{mg}{k}$$
되고,
$$v{\left(t\right)} =\frac{mg}{k} (1- e^{-\frac{k t}{m}})$$
얻을 수 있습니다. 속도 $v{\left(t\right)} $의 $t→\infty$의 극한, 즉 종단 속도는
$$\lim_{t\to\infty} v(t)=\frac{mg}{k}$$
됩니다.

이것을 Sympy로 풀어 그래프도 그린다

#%%
import sympy as sym
from IPython.display import Math, display
from sympy.solvers import ode
from sympy import oo
import matplotlib.pyplot as plt

# 変数を定義
t, v, k, m, g = sym.symbols("t v k m g")
v = sym.Function('v')

# 運動方程式（常微分方程式）
eq = sym.Eq(m * v(t).diff(t,1) - m * g + k * v(t), 0)

# 一般解
ans = sym.expand(sym.dsolve(eq))
display(Math(f"一般解: {sym.latex(ans)}"))

# 特殊解(文字式)
ans2 = sym.expand(sym.dsolve(eq, ics={v(0):0}))
display(Math(f"特殊解(文字式）: {sym.latex(ans2)}"))

# 特殊解（数値）
m1 = 0.4
k1 = 0.2
g1 = 9.8
eq2 = eq.subs([(m, m1), (k, k1), (g, g1)])
ans3 = sym.dsolve(eq2, ics={v(0):0})
p1 = sym.plot(ans3.rhs, (t, 0, 10), show=False) # rhs は右辺のこと
display(Math(f"特殊解(m={m1}, k={k1}, g={g1}): {sym.latex(ans3)}"))

# 終端速度
voo = sym.limit(ans3.rhs, t, oo)
display(Math(r"終端速度：\lim_{t \to \infty} v(t)" \
    f"={voo:.2f}"))
p2 = sym.plot(voo, (t, 0, 10), show=False)
p1.extend(p2)
p1[1].line_color = "0.8"
p1.show()
fig = p1._backend.fig
fig.legend([f"f(t)={ans3.rhs}", f"f(t)={voo:.1f}"], loc='center')
fig

이렇게하면,



특수해를 요구하는 것이 Sympy의 좋은 곳입니다.
게다가 그 특수해의 식을 sympy.plot에 건네주면 그래프도 그려 버립니다. 또한 극한을 구하여 점근선을 그릴 수 있습니다.