はじめに
scipy integrateのsolve_ivpにより、常微分方程式を解くことができる。ここでは、空気抵抗と回転によるマグヌス効果をうける斜方投射を例として、その使い方を説明する。
マグヌス効果とは回転して進行する球に力が生じる現象のことである。
コード
解説
モジュールのインポート
バージョン
初期値の設定
v0を初期値として設定する。(x,y,vx,vy)として、(0,10)の位置から(10,10)の速度で打ち上げることとした。
定数を定義
mは質量、kは抗力係数、gは重力加速度とした。sはマグヌス効果の係数で、wは回転ベクトルで (0,0,1)とした。
関数の定義
微分方程式を定義する。
運動の第2法則から以下の式がたてられる。
\[m \cdot u” = F + mg\]
空気抵抗による力はF = -ku’で表され、マグヌス効果による力は、ω=(0,0,1)とすると、s(−ωvy, ωvx, 0)となるので式を変換すると、以下のようになる。
\[u_x” = -\frac{k}{m} u_x’ -\frac{s}{m} u_y’ \]
\[u_y” = -\frac{k}{m} u_y’ +\frac{s}{m} u_x’+ g\]
solve_ivpで解く
t_spanは時間の（開始、終了）で、tは実際に計算に用いる配列を設定する。t_evalにtを設定するが、設定しない場合は自動的に適度な間隔で計算が行われる。その他の変数はargs=(k,s,w)のように設定しておく。v0は初期値となる。
solve_ivp(f, t_span, v0, t_eval=t,args=(k,s,w))で微分方程式を解く。
解けたかどうかは、messageを確認することでわかる。
結果の表示
ind = solved_data.y[1]>=0でyが正となる値だけを取り出して表示している。
解いたデータはsolved_data.y[0][ind]のようにして得ることができる。
マグヌス係数s=1の時のxとyの関係は下のようになる。
sを変えた時の変化
for s in [0,0.5,1,1.5]:としてsの値を変化させて、計算を行った結果は以下のようになる。
参考
コメント