Как решить систему ODE с scipy.integrate.odeint с переменной, зависящей от времени

Насколько я понял из комментариев под вашим вопросом, вы пытаетесь включить данные измерений, которые могут быть зашумленными. Вместо того, чтобы подключать данные напрямую, вы можете использовать эти данные для соответствия своим временным курсам. Здесь я показываю результат для вашей переменной S:

green dots взяты из решения предоставленной вами системы ODE. Чтобы имитировать ошибки измерения, я добавил к этим данным немного шума (blue dots). Затем вы можете настроить свою систему ODE для максимально точного воспроизведения этих данных (red line).

Для этих задач вы можете использовать lmfit. Код, воспроизводящий сюжет, выглядит следующим образом (некоторые пояснения можно найти во встроенных комментариях):

# zombie apocalypse modeling
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from lmfit import minimize, Parameters, Parameter, report_fit
from scipy.integrate import odeint


# solve the system dy/dt = f(y, t)
def f(y, t, paras):

    Si = y[0]
    Zi = y[1]
    Ri = y[2]

    try:
        P = paras['P'].value
        d = paras['d'].value
        B = paras['B'].value
        G = paras['G'].value
        A = paras['A'].value

    except:
        P, d, B, G, A = paras
    # the model equations (see Munz et al. 2009)
    f0 = P - B * Si * Zi - d * Si
    f1 = B * Si * Zi + G * Ri - A * Si * Zi
    f2 = d * Si + A * Si * Zi - G * Ri
    return [f0, f1, f2]


def g(t, x0, paras):
    """
    Solution to the ODE x'(t) = f(t,x,p) with initial condition x(0) = x0
    """
    x = odeint(f, x0, t, args=(paras,))
    return x


def residual(paras, t, data):
    x0 = paras['S0'].value, paras['Z0'].value, paras['R0'].value
    model = g(t, x0, paras)
    s_model = model[:, 0]
    return (s_model - data).ravel()

# just for reproducibility reasons
np.random.seed(1)

# initial conditions
S0 = 500.              # initial population
Z0 = 0                 # initial zombie population
R0 = 0                 # initial death population
y0 = [S0, Z0, R0]     # initial condition vector
t = np.linspace(0, 5., 100)         # time grid

P = 12      # birth rate
d = 0.0001  # natural death percent (per day)
B = 0.0095  # transmission percent  (per day)
G = 0.0001  # resurect percent (per day)
A = 0.0001  # destroy percent  (per day)

# solve the DEs
soln = odeint(f, y0, t, args=((P, d, B, G, A), ))
S = soln[:, 0]
Z = soln[:, 1]
R = soln[:, 2]

# plot results
plt.figure()
plt.plot(t, S, label='Living')
plt.plot(t, Z, label='Zombies')
plt.xlabel('Days from outbreak')
plt.ylabel('Population')
plt.title('Zombie Apocalypse - No Init. Dead Pop.; No New Births.')
plt.legend(loc=0)
plt.show()

# generate fake data
S_real = S[0::8]
S_measured = S_real + np.random.randn(len(S_real)) * 100
t_measured = t[0::8]

plt.figure()
plt.plot(t_measured, S_real, 'o', color='g', label='real data')

# add some noise to your data to mimic measurement erros
plt.plot(t_measured, S_measured, 'o', color='b', label='noisy data')

# set parameters including bounds; you can also fix parameters (use vary=False)
params = Parameters()
params.add('S0', value=S0, min=490., max=510.)
params.add('Z0', value=Z0, vary=False)
params.add('R0', value=R0, vary=False)
params.add('P', value=10, min=8., max=12.)
params.add('d', value=0.0005, min=0.00001, max=0.005)
params.add('B', value=0.01, min=0.00001, max=0.01)
params.add('G', value=G, vary=False)
params.add('A', value=0.0005, min=0.00001, max=0.001)

# fit model
result = minimize(residual, params, args=(t_measured, S_measured), method='leastsq')  # leastsq nelder
# check results of the fit
data_fitted = g(t, y0, result.params)

plt.plot(t, data_fitted[:, 0], '-', linewidth=2, color='red', label='fitted data')
plt.legend()
# display fitted statistics
report_fit(result)

plt.show()

Вы не можете заранее знать, в каких точках числовой интегратор вычисляет функцию ОДУ. Интегратор (odeint и другие, которые явно не являются «фиксированным размером шага») динамически генерирует внутренний список точек, который может иметь меньший, а иногда и больший размер шага, чем данный список точек выборки. Значения для вывода интерполируются из внутреннего списка.

Если вы хотите заменить часть ODE функцией, вам необходимо преобразовать данные образца в функцию. Это можно сделать с помощью интерполяции. Используйте scipy.interpolate. interp1 для создания функциональных объектов, которые затем можно использовать, как любую другую скалярную функцию.

чтобы конкретно сделать то, что я задавал в вопросе, а именно использовать значения вместо одного из hte ODES, вам нужно будет использовать цикл, в котором вы используете odesolver для решения вашей системы в течение 1 секунды, а затем принимаете выходные данные как начальные условия для следующей итерации цикла. Код для этого подхода приведен ниже. Однако, как отмечали многие, в большинстве ситуаций будет лучше использовать интерполяцию, как описано Клебом и другими.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

Si = [345, 299, 933, 444, 265, 322] # replaced an equation with list

#Parameters
P = 0           # birth rate
d = 0.0001  # natural death percent (per day)
B = 0.0095  # transmission percent  (per day)
G = 0.0001  # resurect percent (per day)
A = 0.0001  # destroy percent  (per day)

# solve the system dy/dt = f(y, t)
def f(y, t):
    Zi = y[0]
    Ri = y[1]
    # the model equations (see Munz et al. 2009)
    f0 = B*Si*Zi + G*Ri - A*Si*Zi
    f1 = d*Si + A*Si*Zi - G*Ri
    return [f0, f1]

# initial conditions
Z0 = 0                      # initial zombie population
R0 = 0                      # initial death population
y0 = [Z0, R0]   # initial condition vector
# a timestep of 1 forces the odesolve to use your inputs at the beginning and provide outputs at the end of the timestep. 
# In this way the problem that LutzL has described is avoided.
t  = np.linspace(0, 1, 2)       
Si =np.array(Si).T

#create a space for storing your outputdata
dataZ =[]
dataR =[]
#use a for loop to use your custom inputs for Si
for Si in Si:
    y0 = [Z0, R0]
    soln = odeint(f, y0, t)
    Z = soln[:, 0]
    R = soln[:, 1]
    #define your outputs as the initial conditions for the next iteration of the loop
    Z_0 = Z[1]
    R_0 = R[1]
    #store your outputs 
    dataZ.append(Z[1])
    dataR.append(R[1])


print (dataZ)
print (dataR)