[SciPy] 13. scipy integrateのsolve_ivpを使って常微分方程式を解く その1

matplotlib

はじめに

scipy integrateのsolve_ivpにより、常微分方程式を解くことができる。ここでは、空気抵抗をうける斜方投射を例としてその使い方を説明する。

コード

コードをダウンロード(.pyファイル)

コードをダウンロード(.ipynbファイル)

解説

モジュールのインポート

バージョン

定数を定義

mは質量、kは抗力係数、gは重力加速度とした。

初期値の設定

v0を初期値として設定する。(x,y,vx,vy)として、(0,0)の位置から(10,10)の速度で打ち上げることとした。

関数の定義

微分方程式を定義する。
空気抵抗による力はF = -ku’で表されるので、運動の第2法則から以下の式がたてられる。

\[m \cdot u” = F + mg\]

式を変換すると、以下のようになる。

\[u” = -\frac{k}{m} u’ + g\]

x成分は右辺の1項目のみで、y成分は重量加速度が入るので1項目+2項目となる。

solve_ivpで解く

t_spanは時間の(開始、終了)で、tは実際に計算に用いる配列を設定する。t_evalにtを設定するが、設定しない場合は自動的に適度な間隔で計算が行われる。その他の変数(ここではk)はargs=(k,)のように設定しておく。v0は初期値となる。
solve_ivp(f, t_span, v0, t_eval=t,args=(k,))で微分方程式を解く。
結果は以下のようになる。

結果の表示

ind = solved_data.y[1]>=0でyが正となる値だけを取り出して表示している。
解いたデータはsolved_data.y[0][ind]のようにして得ることができる。

抗力係数k=1の時のxとyの関係は下のようになる。空気抵抗により、x方向の伸びが抑制されていることがわかる。

kを変えた時の変化

for k in np.linspace(0., 5., 5):としてkの値を変化させて、計算を行った結果は以下のようになる。

この時のx方向の速度の経時変化は以下のようになる。
時間はsolved_data.tで得られる。抵抗が大きい場合、速度は0に収束していく様子が見てとれる。

この時のy方向の速度の経時変化は以下のようになる。

k=5の時の結果を見ると、終端速度は一定値に収束することがわかる。

参考

scipy.integrate.solve_ivp — SciPy v1.6.2 Reference Guide
IPython Cookbook - 12.3. Simulating an ordinary differential equation with SciPy
IPython Cookbook,

コメント