空気抵抗を受けながら落下する小球の速さをPythonとSympyで計算したお話し

#はじめに
　大学の研究室に転がっていた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に渡してやるとグラフも描けてしまいます。さらに極限を求めて漸近線を引くことができます。