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)
Author And Source
この問題について(Scipyのsolve_ivpで条件を満たすまでシミュレーションを行う), 我々は、より多くの情報をここで見つけました https://qiita.com/sirococoa/items/df57f4cdb4891aa6f055著者帰属:元の著者の情報は、元のURLに含まれています。著作権は原作者に属する。
Content is automatically searched and collected through network algorithms . If there is a violation . Please contact us . We will adjust (correct author information ,or delete content ) as soon as possible .