はじめに
scipy integrateのsolve_ivpにより、常微分方程式を解くことができる。ここでは、空気抵抗をうける斜方投射を例としてその使い方を説明する。
コード
解説
モジュールのインポート
バージョン
定数を定義
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の時の結果を見ると、終端速度は一定値に収束することがわかる。
コードをダウンロード(.pyファイル) コードをダウンロード(.ipynbファイル)参考
コメント