Access Current Time Step In Scipy.integrate.odeint Within The Function

Is there a way to access what the current time step is in scipy.integrate.odeint? I am trying to solve a system of ODEs where the form of the ode depends on whether or not a popula

Solution 1:

Of course you can hack something together that might work. You could log t but you have to be aware that the values might not be constantly increasing. This depends on the ODE algorithm and how it works (forward, backward, and central finite differences).

But it will give you an idea where about you are.

logger = [] # visible in odefuncdefodefunc(P,t):

    x = P[0]
    y = P[1]
    z = P[2]

    print(t)
    logger.append(t)

    if logger: # if the list is not emptyif logger[-1] > 2.5: # then read the last value print('hua!')

    if A * x - B * x * y < 0.6:
        dxdt = A/5 * x
        dydt = -C * y + D * x * y
        dzdt = - B * z * y
    else:
        dxdt = A * x - B * x * y
        dydt = -C * y + D * x * y
        dzdt = 0


    dPdt = np.ravel([dxdt, dydt, dzdt])

    return dPdt

print(logger)

Solution 2:

As pointed out in the another answer, time may not be strictly increasing at each call to the ODE function in odeint, especially for stiff problems.

The most robust way to handle this kind of discontinuity in the ode function is to use an event to find the location of the zero of (A * x - B * x * y) - 0.6 in your example. For a discontinuous solution, use a terminal event to stop the computation precisely at the zero, and then change the ode function. In solve_ivp you can do this with the events parameter. See the solve ivp documentation and specifically the examples related to the cannonball trajectories. odeint does not support events, and solve_ivp has an LSODA method available that calls the same Fortran library as odeint.

Here is a short example, but you may want to additionally check that sol1 reached the terminal event before solving for sol2.

from scipy.integrate import solve_ivp

tend = 10defdiscontinuity_zero(t, y):
    return y[0] - 10

discontinuity_zero.terminal = Truedefode_func1(t, y):
    return y

defode_func2 (t, y):
    return -y**2

sol1 = solve_ivp(ode_func1, t_span=[0, tend], y0=[1], events=discontinuity_zero, rtol=1e-8)
t1 = sol1.t[-1]
y1 = [sol1.y[0, -1]]
print(f'time={t1}  y={y1}   discontinuity_zero={discontinuity_zero(t1, y1)}')

sol2 = solve_ivp(ode_func2, t_span=[t1, tend], y0=y1, rtol=1e-8)

plt.plot(sol1.t, sol1.y[0,:])
plt.plot(sol2.t, sol2.y[0,:])
plt.show()

This prints the following, where the time of the discontinuity is accurate to 7 digits.

time=2.302584885712467  y=[10.000000000000002]   discontinuity_zero=1.7763568394002505e-15