I am not sure if the event handling in scipy.integrate.solve_ivp is working correctly. In the example below, I integrated a derivative which should result in a cubic polynomial with three roots at x=-6, x=-2 and x=2. I set up an event function that returns y, which will be zero at those x-values. I expected to see three entries in the t_events attribute of the solution, but I only see one, even though it is clear the solution crosses the x-axis three times.
Am I doing something wrong here?
Reproducing code example:
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])
The code above results in:
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. ]])
```
The problem is you can see from the sol.y array there should be three zeros (there are three sign changes), but only one event is recorded.
Scipy/Numpy/Python version information:
1.0.0 1.13.3 sys.version_info(major=3, minor=6, micro=0, releaselevel='final', serial=0)
[UPDATE]: if I use the optional max_step argument to solve_ivp, and make it small enough, then I see all three roots. It seems that the event functions are not called at the t_eval steps, but rather only on internal solver steps, which are far fewer than the t_eval steps, and which end up skipping some roots. That doesn't seem very useful, as you have to know how to set max_steps to avoid missing roots.
