scipy/solve_ivp における初期条件のイベント検知回避

問題点

solve_ivpのイベント検知関数 find_active_events はデフォルトで以下のように定義されている．（コメントは省略した）

ivp.py

def find_active_events(g, g_new, direction):
    g, g_new = np.asarray(g), np.asarray(g_new)
    up = (g <= 0) & (g_new >= 0)
    down = (g >= 0) & (g_new <= 0)
    either = up | down
    mask = (up & (direction > 0) |
            down & (direction < 0) |
            either & (direction == 0))

    return np.nonzero(mask)[0]

上記定義に従うと，初期条件 $(t_0, y_0)$ が g=event(t0,y0)=0 を満たす場合でもイベントを検知してしまうため，例えば Poincaré 断面 $\Pi_0 = \{ \boldsymbol{x} \in R^n ~|~ x = 0\}$ 上の点を初期条件として再び $\Pi_0$ 上へ戻る軌道を計算したい場合などには不便である．

解決策

このイベントを回避するためには，ivp.py を直接，以下のように変更すれば良い．

ivp.pyの変更

# up = (g <= 0) & (g_new >= 0)
up = (g <= -1e-12) & (g_new >= 0)

# down = (g >= 0) & (g_new <= 0)
down = (g >= 1e-12) & (g_new <= 0)

原理

g の値は g_new の古い値をストアしているに過ぎないため，必ず g == 0 よりも（時系列的な意味で）先に g_new == 0 が検知される．
言い換えれば，g == 0 によって検知されるイベントはすべて「初期条件により生じるイベント」といえる．

そこで，abs(g) < 1e-12 の場合にはイベントを検知しないよう，上記の変更を施せば，初期条件由来のイベント検知を回避でき，かつ元々の機能も損なわない．

補足

ivp.py の格納場所は inspect モジュールの getfile 関数によって確認できる．

>>> from scipy.integrate import solve_ivp
>>> import inspect
>>> inspect.getfile(solve_ivp)