微分方程式を「地面に着いたら停止」のような条件で止めたいなら、solve_ivpにイベント関数を渡す。イベントは連続関数が0になる時点として定義し、terminalで停止、directionで通過方向を指定する。到達時刻は通常の出力配列の末尾ではなく、t_eventsから取得する。
最小例で確かめる
以下は説明用に作った小さなデータである。実測データや実行速度の測定結果ではない。コード全体をexample.pyとして保存すれば、入力ファイルを別途用意せずに実行できる。assertは、この例で成り立つべき形や値を確認するために入れてある。
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
g = 9.8
height0 = 10.
def motion(t, state):
height, velocity = state
return [velocity, -g]
def ground(t, state):
return state[0]
ground.terminal = True
ground.direction = -1
result = solve_ivp(motion, (0., 3.), [height0, 0.], events=ground,
t_eval=np.arange(0., 3.01, 0.25), dense_output=True,
max_step=0.1, rtol=1e-9, atol=1e-11)
assert result.success and result.status == 1
assert len(result.t_events[0]) == 1
hit_t = result.t_events[0][0]
hit_state = result.y_events[0][0]
print("status:", result.status)
print("event time:", f"{hit_t:.6f}")
print("last requested output:", f"{result.t[-1]:.6f}")
print("impact velocity:", f"{hit_state[1]:.6f}")
np.testing.assert_allclose(hit_t, np.sqrt(2 * height0 / g), atol=1e-8)
np.testing.assert_allclose(hit_state, [0., -14.], atol=1e-8)
assert result.t[-1] < hit_t
plot_t = np.linspace(0., hit_t, 200)
fig, ax = plt.subplots(figsize=(7, 4), layout="constrained")
ax.plot(plot_t, result.sol(plot_t)[0], color="#a83262", label="ODE solution")
ax.plot(result.t, result.y[0], "o", color="#26778c", label="t_eval output")
ax.plot([hit_t], [0], "*", markersize=13, color="#e38619", label="ground event")
ax.axhline(0, color="0.4", linewidth=1)
ax.set(xlabel="Time [s]", ylabel="Height [m]", title="Stop at downward ground crossing")
ax.legend()
fig.savefig("PY083-event.png", dpi=160)
plt.close(fig)
実行結果
status: 1
event time: 1.428571
last requested output: 1.250000
impact velocity: -14.000000
イベント関数は真偽値を返さない
状態は高さと速度の二つで、微分方程式は高さの変化が速度、速度の変化が-gという空気抵抗なしの落下モデルである。groundは高さそのものを返すので、地面を横切ると正から負へ変わる。「高さが0以下」というTrue/Falseを返す関数にすると、連続な根の探索という前提から外れる。
terminal=Trueはイベント検出時に積分を停止する指定で、direction=-1は正から負への通過だけを対象にする。上向きの通過を無視したい、しきい値を下回ったときだけ停止したい、といった方向条件を表せる。初期状態がちょうど0のときには開始時のイベント判定も関係するので、別途小例で確認する。
イベント時刻と出力時刻は違う
t_evalは保存したい時刻の列であり、ソルバー内部の刻み幅を直接指定するものではない。例では0.25秒ごとの出力を要求したが、地面への到達は約1.428571秒なので、通常の出力の最後は1.25秒になる。result.t[-1]を到達時刻として読むと誤る。
t_eventsにはイベントの時刻、y_eventsにはその時点の状態が入る。イベントが複数種類ある場合は、渡した関数の順に配列が分かれる。イベントが起きなければ対応する配列は空なので、0番目を取り出す前に件数を確かめる。successに加え、statusがイベントによる停止を示す1であることを確認した。
イベントを見落とす条件もある
イベント検出は内部ステップの端での符号変化を探すため、一つのステップ内で複数回0を横切ると見落とすことがある。接するだけで符号が変わらない場合も単純な通過検出では扱いにくい。t_evalを細かくしただけで、こうした問題が必ず解消するわけではない。
max_stepは内部の最大刻み幅を制限する。例では0.1秒としたが、どんな問題でもこの値なら安全という意味ではない。現象の時間スケールに合わせて値を変え、刻み幅や許容誤差を厳しくしたときにイベント時刻が安定するか調べる。連続性や単位、イベント関数の定義も同時に点検する。
解析解や保存則を使って検査する
この落下モデルでは到達時刻sqrt(2h/g)と衝突速度を解析的に求められるため、数値結果を独立に比較できる。図では連続出力による曲線、t_evalで指定した点、イベント位置を区別した。dense_outputで得た補間解は、計算した時間範囲の中で使い、停止後へ安易に延長しない。
実際の問題では解析解がないことも多いが、非負であるべき状態、保存量、既知の極限などは検査に使える。result.messageも保存して終了理由を確認したい。ソルバーのsuccessは数値手続きの終了状態であり、空気抵抗などを省いたモデルが現実を十分に表すことまで保証するものではない。
動作確認環境と参考資料
Linux・CPython 3.12.14、NumPy 2.3.5、pandas 2.2.3、SciPy 1.17.0、Matplotlib 3.10.8の環境で掲載コードを実行した。使用するライブラリはコード冒頭のimportを参照してほしい。公式資料の最新版と、この実行確認版は区別している。数値の末尾や表の表示幅は環境によって変わることがある。