Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to use events in scipy.integrate.solve_ivp

Tags:

People also ask

What is the difference between Odeint and solve_ivp?

odeint came first and is uses lsoda from the FORTRAN package odepack to solve ODEs. solve_ivp is a more general solution that lets use decide which integrator to use to solve ODEs. If you define the method param as method='LSODA' then this will use the same integrator as odeint .

How does Scipy integrate?

For n-fold integration, scipy provides the function nquad . The integration bounds are an iterable object: either a list of constant bounds, or a list of functions for the non-constant integration bounds. The order of integration (and therefore the bounds) is from the innermost integral to the outermost one.

What does solve_ivp return?

Now, solve_ivp returns a ndarray named y (in the present case, it will be position. y , of shape (4, 104) ) which according to the scipy. integrate. solve_ivp documentation, gives the "Values of the solution at t".

How do I import an Odeint function in Python?

The Python code first imports the needed Numpy, Scipy, and Matplotlib packages. The model, initial conditions, and time points are defined as inputs to ODEINT to numerically calculate y(t). An optional fourth input is args that allows additional information to be passed into the model function.


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.