sympyで超基本の微分方程式を解いてみよう

こんにちは森の生活者、ガマ仙人です。
現在この記事は限定公開にしてます。
qiita初めてなので、怖くて怖くて公開できません。
で、Sympyを昨日から勉強してます。
Maximaでできることはだいたいできるようですね。
今日は簡単な微分方程式を解いてみたいと思います。
超簡単な例題でやってみましたが
あら、不思議、あっけなく解けてしましました。
やるなあSympy!
で、とりあえずメモしておきます。

#自由落下運動における高さと加速度の関係

test1.py

h = Symbol('h') #高さ
t = Symbol('t') #時間
g = Symbol('g') #重力加速度

#微分方程式を定義するーー時間ごとの高さを2階微分したものはマイナスの重力加速度である
before =  Eq(h(t).diff(t,2), -g)

#微分方程式を解いてみよう
after =dsolve(before,h(t))

#結果を表示してみる
display(before,after)

結果

\frac{d^{2}}{d t^{2}}  h{\left (t \right )} = - g

h{\left (t \right )} = C_{1} + C_{2} t - \frac{g t^{2}}{2}

というわけで解けてますね

#単振動の微分方程式

test2.py

x = Symbol('x') #振動の位置
w = Symbol('ω') #角速度
t = Symbol('t') #時間

#微分方程式を定義するーーなんと表現していいかわからんので式をそのまま書くよ
before = Eq(x(t).diff(t,2), -w**2*x(t))

#微分方程式を解いてみよう
after =dsolve(before,x(t))

#結果を表示してみる
display(before,after)

結果

\frac{d^{2}}{d t^{2}}  x{\left (t \right )} = - ω^{2} x{\left (t \right )}

x{\left (t \right )} = C_{1} e^{- i t ω} + C_{2} e^{i t ω}

これも解けましたね

#y′ = ay の基本形

test3.py

x=Symbol('x')           
y=Symbol('y') 
a=Symbol('a') 

eq = Eq(y(x).diff(x), a*y(x))
display(eq, dsolve(eq, y(x)))

\frac{d}{d x} y{\left (x \right )} = a y{\left (x \right )}

y{\left (x \right )} = C_{1} e^{a x}

#ロジスティック方程式の基本形

test4.py

x=Symbol('x')           
y=Symbol('y') 
a=Symbol('a') 
b=Symbol('b') 

eq = Eq(y(x).diff(x), (a-b*y(x))*y(x))
display(eq, dsolve(eq, y(x)))

\frac{d}{d x} y{\left (x \right )} = \left(a - b y{\left (x \right )}\right) y{\left (x \right )}

y{\left (x \right )} = \frac{a e^{a \left(C_{1} + x\right)}}{b \left(e^{a \left(C_{1} + x\right)} - 1\right)}

#補足
ちなみに、言う必要もないかもしれまんが
上記のサンプル動かす前に
次のことやっておいてくださいね。
数式が綺麗に出力されますから、気持ちがいいです。

from sympy import *
from sympy.interactive import printing
printing.init_printing(use_latex=True)