如何在scipy.integrate.solve_ivp中使用事件

Question

我不确定scipy.integrate.solve_ivp中的事件处理是否正常工作.在下面的例子中,我整合了一个导数,它应该产生一个三次多项式,在x = -6,x = -2和x = 2时有三个根.我设置了一个返回y的事件函数,该函数在这些x值处为零.我希望在解决方案的t_events属性中看到三个条目,但我只看到一个,即使很明显解决方案穿过x轴三次.

我在这里做错了吗？

再现代码示例:

def fprime(x, y):
    return 3 * x**2 + 12 * x - 4

def event(x, y):
    return y

import numpy as np
from scipy.integrate import solve_ivp
sol = solve_ivp(fprime, (-8, 4), np.array([-120]), t_eval=np.linspace(-8, 4, 10), events=[event])

上面的代码导致:

message: 'The solver successfully reached the interval end.'
      nfev: 26
      njev: 0
       nlu: 0
       sol: None
    status: 0
   success: True
         t: array([-8.        , -6.66666667, -5.33333333, -4.        , -2.66666667,
        -1.33333333,  0.        ,  1.33333333,  2.66666667,  4.        ])
  t_events: [array([-6.])]
         y: array([[-120.        ,  -26.96296296,   16.2962963 ,   24.        ,
           10.37037037,  -10.37037037,  -24.        ,  -16.2962963 ,
           26.96296296,  120.        ]])

```

问题是你可以从sol.y数组看到应该有三个零(有三个符号变化),但只记录了一个事件.

Scipy/Numpy/Python版本信息:

1.0.0 1.13.3 sys.version_info(major=3, minor=6, micro=0, releaselevel='final', serial=0)

[更新]:如果我将可选的max_step参数用于solve_ivp,并使其足够小,那么我会看到所有三个根.看起来事件函数不是在t_eval步骤中调用的,而是仅在内部求解器步骤上调用,这些步骤远远少于t_eval步骤,最终会跳过一些根.这似乎不是很有用,因为您必须知道如何设置max_steps以避免错过根.

Answer 1

事件函数不会在t_eval步骤中调用，而仅在内部求解器步骤中调用

这是对的。找到新的一(t, y)对后，将根据它们计算事件并将其传递给事件函数，将事件函数的符号与上一步进行比较。find_active_events

t_eval根本不被事件计算使用，稍后处理。因此，事件看到的符号变化就是您在带有 的输出中看到的符号变化t_eval=None，其中只有一个。

y: array([[-120.        , -110.49687882,  -35.93785936,   94.46893375,
     120.        ]]) 

似乎值得在 SciPy 跟踪器上提出问题；有一些关于 的悬而未决的问题solve_ivp，这仍然是相当新的。