【備忘録】Pyhtonで数式計算_常微分方程式

この記事の続きです

yyhhyy.hatenablog.com

 このsympyというライブラリの便利さは異常で、微分方程式も解けます。

例題*1

豚の体重wの成長率は豚の体重に比例すると仮定し、
初日（t=0）での豚の体重を200ポンド、
初日（t=0）での豚の体重の成長率は、5ポンド／日
とする

⇒ 成長率(dw(t) dt)が現時点での体重に比例するので、
dwdt = c * w(t) とおける。
w(0) = 200 かつ dwdtが t=0 のとき 5 になるので
dwdt = (5/200)*w(t)  = (1/40)*w(t)　となる
この微分方程式を解きたい

解きたい式

dwdt = (1/40)*w(t)

コード

# sympyをインポート
from sympy import *
# 変数 t を宣言
t= symbols("t")
# wは式として使うので式(Function)として宣言。
w = symbols("w",cls=Function)
dwdt = w(t).diff(t)
eq1 = Eq ( (dwdt - 1*w(t)/40),0)
# 微分方程式を解く時は、dsolveを使う。
# 変数分離法を使えるので hint="separable"のオプションを加える
ans = dsolve(eq1,w(t),hint="separable")
print("ans",ans)
# 積分定数C1が残るので変数宣言をしておく
C1=symbols("C1")
# 結果を見て改めてwを定義
w = C1*exp(t/40)
# t=0の時 200になる筈なので方程式を作り、C1について解く
eq2 = Eq ( (w.subs(t,0) - 200),0)
C1ans = solve([eq2],[C1])
print("C1",C1ans)
# wが求まったので改めて定める
w = 200*exp(t/40)
#実際に微分してt=0を代入した時に5になるか確認してみる
print(w.diff(t))
print(w.diff(t).subs(t,0))

結果

('ans', w(t) == C1*exp(t/40))
('C1', {C1: 200})
5*exp(t/40)
5

ちょっとかしこすぎないだろうか？sympy

-----------------------------------
微分方程式の解き方について

3.2. Sympy : Python での代数計算 — Scipy lecture notes

tokeigaku.blog.jp

積分定数を宣言しちゃえ！ということに気がついたのはこちらの記事から

note.chiebukuro.yahoo.co.jp

今回(も)引用した例題については、こちらの本のp57より。
この本、感度の計算についても触れてあって、統計の本の中でもマーケティング寄り。
まだ読書中ですが良書の予感がします。

数理モデリング入門 ファイブ・ステップ法 [ マーク・M．ミールスキャーアト ]

*1:例によって前回の記事と同じ本からの引用です