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 .
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.
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".
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?
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.
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.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With