CODEWITHSUNDEEP
pythonscipy
John Kitchin
John Kitchin

Reputation: 2433

How to use events in scipy.integrate.solve_ivp

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.

Upvotes: 16

Views: 16770

Answers (1)

user6655984
user6655984

Reputation:

the event functions are not called at the t_eval steps, but rather only on internal solver steps

This is correct. After a new (t, y) pair is found, the events are computed from them and passed to find_active_events that will compare the signs of event-functions to the previous step.

t_eval is not used by the event computation at all, it is dealt with later. So the sign changes that events see are those that you see in the output with t_eval=None, of which there is just one.

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

Seems worth raising an issue on SciPy tracker; there are a few open issues regarding solve_ivp, which is still pretty new.

Upvotes: 5

Related Questions