この記事は

python(正確にはscipy)のodeintというライブラリを使って運動方程式を解きます。anacondaで入れたpython3.6.3を使っています。

運動方程式(二階微分方程式)を解いていますが、一階の微分方程式など他の微分方程式にも応用可能です。

Pythonで運動方程式を解きたい

人生とは

簡単な運動方程式を解きたい場合が人生の中では出てくる。
特に、解いて運動の様子をplotしたい場合がある。
その時人はどうするか

今回はpythonのライブラリ"odeint"を用いて楽してプロットすることを目的にする。

odeint

おでんitに見える。
ODE integralのことだと思う。
ODE=Ordinary Differential Equation (常微分方程式)

公式reference
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html

odeというライブラリもある
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html

調べるとどうやらもっと新しいpackage"solve_ivp"があるらしく、こちらが推奨されている。いつかこれについても記事を書きたい。
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp

ちなみにodeintと、他二つ(ode, solve_ivp)では引数の順番が違うというエグい仕様になっている。混同しないように気をつけよう。

今回解きたい運動方程式

運動方程式なら何でも良いのだが、今回は簡単な投げ上げを解く。

m\ddot y=-rv+mg

odeintで二階微分方程式を特には、二つの一階微分に直す必要がある。例えば今回の場合は

\dot y=v \\
m\dot v=-rv+mg

という式を実際にはodeintで解く。

コード

全体の流れ

コードの流れは以下の通り
1. 必要なライブラリを読み込む
2. 解きたい微分方程式を記述する
3. 初期条件等を決めて方程式をodeintに解いてもらう
4. 結果を描く

準備

まず、色々ライブラリを入れる。

import numpy as np #numpy
from scipy.integrate import odeint #odeint
import matplotlib.pyplot as plt #to draw graphs

微分方程式を書く

そして、微分方程式を定義する。
空気抵抗r,質点の質量mは次のステップで定義される。
重力定数gは変更しないので今回はここで定義している。

def func(s, t, r, m):

    y, v = s #sは変数yとvの組
    g=9.80665 #m/s^2
    dsdt = [v, (-r*v-m*g)/m]
    return dsdt 

初期条件を決めて、微分方程式を解いてもらう

空気抵抗r,質量mを決め、初期条件を決めて微分方程式を解く。
sol = odeint(func, y0, t, args=(r,m))の最後の引数について、一点注意。
args=(r,m)のようにリストを渡すので、引数が一つの場合はargs=(r,)のようにしなければならない。今回は困らないが、他の運動方程式を解く場合には注意。

r=12
m=100
y0 = [0,10]#位置0から初速10で投げ上げ
t = np.linspace(0, 7, 201)#時刻0から7まで、201step刻みで計算する

sol = odeint(func, y0, t, args=(r,m))

結果を図に描く

結果を可視化してみよう。

plt.plot(t, sol[:, 0], 'b', label='y')#yについてplot
plt.plot(t, sol[:, 1], 'g', label='v')#vについてplot
plt.legend(loc='best')#レジェンドを付ける
plt.xlabel('t')
plt.grid()#格子を付ける
plt.show()

そうするとこんな図が得られるはず。

この次は

最初はtを0から100までなどと変更して終端速度を確認したり、あるいは斜めに投げ上げた場合の問題を解いたりするのがコードを理解する良い練習問題になると思う。