ODE with time-varying coefficients in scipy

user_1506850 0
Viewed 130

I am evaluating a set of ODEs with time varying coefficients

def deriv(y, t, N, coefficients):
    S, I, R = y
    dSdt = coefficients['beta'](t) * S * I / N * -1
    dIdt = coefficients['beta'](t) * S * I / N - coefficients['gamma']* I
    dRdt = coefficients['gamma'] * I
    return dSdt, dIdt, dRdt

In particular, I have 'beta' values in a pre-calculated array, of size equal to int(max(t)).

coefficients = {'beta' : beta_f,'gamma':0.1}
def beta_f(t):
     return  mybetas.iloc[int(t)]

# Initial conditions vector
y0 = (S0, I0, R0)
# Integrate the SIR equations over the time grid, t.
ret = odeint(deriv,y0,t,args=(N,coefficients))

When I run odeint, it is evaluate also for value beyond max(t), raising an index out of bound error in beta_f. How to limit the evaluation span for odeint?

ode scipy numpy python
1 Answers

Since len(mybetas) == int(max(t)), you can get an out-of-bounds error even for values of t which are not beyond max(t). For example, mybetas.iloc[int(max(t))] will give you the out-of-bounds error, even though int(max(t)) <= max(t) for positive values of t.

But to your point, odeint does indeed check some values outside of the domain of integration. I had to deal with a problem similar to yours just a few weeks ago, and the following two discussions on stackoverflow were really helpful:

integrate.ode sets t0 values outside of my data range

Solve ODEs with discontinuous input/forcing data

The second link explains why it might be computationally faster to solve the ODE with odeint over each individual integer time step one after the other in a for loop, instead of letting odeint deal with the discontinuities in your derivative caused by jumps in the values of your betas.

Otherwise, if this is appropriate for your study case, you can interpolate your betas, and let the function beta_f return interpolated values of beta. Of course, you will have to extend the interpolation domain slightly beyond your integration domain, since odeint might want to evaluate the derivative for some t larger than max(t).

user_159647770
Related