如何使用 scipy.integrate.odeint 求解具有时间相关变量的 ODE 系统

Question

我正在使用 scipy 食谱中的 Zombie Apocalypse example 来学习如何在 python 中求解 ODE 系统。

在此模型中，有一个方程可以根据出生率、死亡率和初始人口提供每天的人口数量。然后根据人口计算创建和杀死多少僵尸。

我有兴趣将人口微分方程替换为一个数据列表，该数据列表告诉我们每个时间步长的人口数量。我收到以下错误：

TypeError: can't multiply sequence by non-int of type 'float'

正如人们所指出的那样，这是因为将单个数字乘以列表是没有意义的。我不确定如何在每次 T 时将列表中的数字提供给微分方程。

这是两次尝试的代码

# solve the system dy/dt = f(y, t)
def f(y, t):
    Si = [345, 299, 933, 444, 265, 322] # replaced an equation with list
    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]

我也试过了

numbers = [345, 299, 933, 444, 265, 322]
for t in [0, 5]:
    Si = numbers

# 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]

这两种尝试都有相同的问题，即向 f0 和 f1 提供整个列表，而不是从列表中迭代地提供 1 个数字。

Answer 1

您无法先验地知道数值积分器在什么点计算 ODE 函数。积分器（odeint 和其他没有明确“固定步长”的）动态生成一个内部点列表，其步长可能比给定的采样点列表更小，有时也更大。输出的值是从内部列表中插入的。

如果您想用函数替换 ODE 的一部分，那么您必须将样本数据转换为函数。这可以通过插值来完成。使用 scipy.interpolate.interp1 函数生成函数对象，然后您可以像使用任何其他标量函数一样使用它们。

Answer 2

据我从您的问题下方的 cmets 中了解到，您尝试合并可能有噪声的测量数据。您可以使用这些数据来适应您的时间进程，而不是直接插入数据。在这里，我显示了您的变量 S 的结果：

green dots 是从您提供的 ODE 系统的解决方案中抽取的。为了模拟测量误差，我在这些数据中添加了一些噪声 (blue dots)。然后，您可以让您的 ODE 系统尽可能好地重现这些数据 (red line)。

对于这些任务，您可以使用lmfit。重现情节的代码如下所示（一些解释可以在 inline cmets 中找到）：

# 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()

Answer 3

要具体执行我在问题中提出的问题，即使用值代替 hte ODES 之一，您需要使用一个循环，在该循环中使用 odesolver 求解系统 1 秒，然后取输出作为循环下一次迭代的初始条件。这种方法的代码如下。然而，正如许多人指出的那样，在大多数情况下，最好使用 Cleb 和其他人所描述的插值

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)