[SciPy-User] How to create a request to change a function?

[Fabrice Silva]

Le vendredi 31 janvier 2020, Augusto Dufloth a écrit :
> Hello,
> 
> I have a request to change one of scipy’s function.
> 
> I am trying to use odeint to integrate a set of equations of motion.
> My state equations rely on time dependent inputs. And these inputs
> are numpy arrays.
> 
> When I do a flight path reconstruction, the inputs comes from the
> inertial system of the aircraft and I plug it into the state
> equations, so I get the deterministic position of the aircraft and
> than I can compare it to the measured values.
> 
> The issue is that I can’t input an array. Searching stackoverflow I
> saw solutions using for-loops and assigning an arg for each time
> step. This is counter productive and create a huge time calculation
> for long arrays.
> 
> I created my own C-dll to deal with this for loop issue, but it’s
> super counter productive. 
> 
> Can we get an option that the input args are arrays if same size of
> t, and as the state is integrated, the index of the arg follows the
> states?

Hi,
As far as I understand, you want to use time-varying parameters (which
includes right-hand side terms of ODEs) specified as numpy array (i.e.,
at discrete time values) while ODEs are defined in the continuous-time
domain. Specification of an extrapolation method is required for the
well-posedness of the problem: the most trivial is the "Sample-and-
hold" process which maintain values until next sample step.

One simple solution (maybe one that you arleady tested) is to use the
ode class and its integrate method.
Adapting the example given in the docs 
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html#scipy.integrate.ode

from scipy.integrate import ode
def f(t, y, arg1):    return [1j*arg1*y[0] + y[1], -arg1*y[1]**2]def
jac(t, y, arg1):    return [[1j*arg1, 1], [0, -arg1*2*y[1]]]
r = ode(f, jac).set_integrator('zvode', method='bdf')
t = np.arange(0, 10, 1)results = np.empty((len(t), 1 + len(y0)),
dtype=float) + np.nan
r.set_initial_value([1.0j, 2.0], t[0])for ind in range(len(t) -
1):    r.set_f_params(param[ind]).set_jac_params(param[ind])    tmp_res
ult = r.integrate(t[ind+1])    if not
r.successful():        break    results[ind] = tmp_result
This introduced a reasonnable overhead in performance while begin quite
readable.

Best regards

Fabrice
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-user/attachments/20200131/432266d4/attachment.html>