Scipyのsolve_ivpで条件を満たすまでシミュレーションを行う

概要

scipyライブラリのsolve_ivpのevent機能について試すために自由落下のシミュレーションを行います。
今回は地面に到着した瞬間までシミュレーションを行います。

式

自由落下の式は次のようになります。

\begin{align}
\frac{dy}{dt} &= v \\
\frac{dv}{dt} &= -mg
\end{align}

ここで、物体の高さを$y$、速度を$v$とします。
$m$と$g$は定数です。

$y=0$のとき地面と衝突したと考えます。
solve_ivpのeventはfloatを返す関数で定義され、eventの発生時に0を返します。
今回はそのまま物体の高さ$y$を返します。

プログラム

上記の式を用いてシミュレーションプログラムを作成します。

from scipy.integrate import solve_ivp
import numpy
from matplotlib import pyplot

# 物理定数
m = 1
g = 1


def free_fall(t, y):
    '''
    自由落下の方程式
    '''
    return numpy.asarray([y[1], -m * g])

def ground(t, y):
    '''
    地面との衝突の判定式
    '''
    return y[0]
ground.terminal = True  # 条件が満たされたらシミュレーションを終了する


if __name__ == '__main__':
    # 初期値
    y0 = numpy.asarray([10, 0])
    # 0から10までの時間を100分割してシミュレーションを行う
    sol = solve_ivp(free_fall, [0, 10], y0, t_eval=numpy.linspace(0, 10, 100), events=ground)
    print(sol.t)
    print(sol.y)
    # グラフ出力
    pyplot.plot(sol.t, sol.y[0])
    pyplot.ylabel('y')
    pyplot.xlabel('t')
    pyplot.show()

結果

プログラムを実行すると次の出力とグラフが得られます。

[0.         0.1010101  0.2020202  0.3030303  0.4040404  0.50505051
 0.60606061 0.70707071 0.80808081 0.90909091 1.01010101 1.11111111
 1.21212121 1.31313131 1.41414141 1.51515152 1.61616162 1.71717172
 1.81818182 1.91919192 2.02020202 2.12121212 2.22222222 2.32323232
 2.42424242 2.52525253 2.62626263 2.72727273 2.82828283 2.92929293
 3.03030303 3.13131313 3.23232323 3.33333333 3.43434343 3.53535354
 3.63636364 3.73737374 3.83838384 3.93939394 4.04040404 4.14141414
 4.24242424 4.34343434 4.44444444]
[[10.          9.99489848  9.97959392  9.95408632  9.91837568  9.87246199
   9.81634527  9.75002551  9.6735027   9.58677686  9.48984797  9.38271605
   9.26538108  9.13784308  9.00010203  8.85215794  8.69401082  8.52566065
   8.34710744  8.15835119  7.9593919   7.75022957  7.5308642   7.30129579
   7.06152433  6.81154984  6.55137231  6.28099174  6.00040812  5.70962147
   5.40863177  5.09743904  4.77604326  4.44444444  4.10264259  3.75063769
   3.38842975  3.01601877  2.63340475  2.2405877   1.8375676   1.42434445
   1.00091827  0.56728905  0.12345679]
 [ 0.         -0.1010101  -0.2020202  -0.3030303  -0.4040404  -0.50505051
  -0.60606061 -0.70707071 -0.80808081 -0.90909091 -1.01010101 -1.11111111
  -1.21212121 -1.31313131 -1.41414141 -1.51515152 -1.61616162 -1.71717172
  -1.81818182 -1.91919192 -2.02020202 -2.12121212 -2.22222222 -2.32323232
  -2.42424242 -2.52525253 -2.62626263 -2.72727273 -2.82828283 -2.92929293
  -3.03030303 -3.13131313 -3.23232323 -3.33333333 -3.43434343 -3.53535354
  -3.63636364 -3.73737374 -3.83838384 -3.93939394 -4.04040404 -4.14141414
  -4.24242424 -4.34343434 -4.44444444]]

地面に衝突したタイミングでシミュレーションが終了していることがわかります。
ground.terminal = Trueの行をコメントアウトしたり、値をFalseに設定するとt=10までシミュレーションが継続します。

参考文献

scipy.integrate.solve_ivp — SciPy v1.6.3 Reference Guide (https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html)