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

Posted 2023-03-12

【中文标题】如何使用 scipy.integrate.odeint 求解具有时间相关变量的 ODE 系统

【英文标题】：How to solve a system of ODEs with scipy.integrate.odeint with a time-dependent variable

【发布时间】：2017-08-16 08:30:51

【问题描述】：

我正在使用 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 个数字。

【问题讨论】：

对于我正在做的现实生活中的问题，我需要用真实数据列表替换系统中的 1 个微分方程。我也尝试在 Y 函数定义之外定义 Si，但得到相同的错误 codenumbers = [345, 299, 933, 444, 265, 322] codefor t in [0, 5]: @ 987654329@ Si = 数字

不要谈论用数字列表替换微分方程（这对我来说没有意义），而是尝试从一开始就根据您拥有的信息来解释问题。例如，类似“我有一个数字列表。这些数字的含义是 [...]。我想建模一个过程，其中这些数字用于 [...]”

你是对的，我用一种特别愚蠢的方式表达了我的问题，我深表歉意。我确实理解你的观点和 python 错误，我不能将我的数字乘以一个列表。我试图了解如何在每个时间步 t 向 f0 和 f1 方程提供列表中的项目。感谢您抽出宝贵时间提供帮助。

我没有清楚地表达我想要做什么。我是一名生物学家，我有一个简单的 ODE 系统，它模拟蛋白质如何被激活，然后与其他几种蛋白质相互作用。我现在有关于给定时间跨度内激活蛋白质水平的真实数据，我想用真实的时间序列数据替换代表模型中该蛋白质激活的方程。在僵尸启示录模型的类比中，我想用不同时间点的人类流行列表替换描述人口出生率和死亡率的方程。

我已编辑问题标题和问题。感谢您的宝贵时间。

【参考方案1】：

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

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

【讨论】：

感谢您的回答，它清晰易懂。不幸的是 interporlate 对我不起作用，我的数据非常嘈杂，它基本上是一片云。是否有可以使用固定步长的 odeint 替代方法？如果不需要，我是否需要使用循环来明确指定步长为 1？

@BobbyM：噪音并不重要；您仍然可以尝试拟合您的数据（我在下面添加了代码）。

【参考方案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()

【讨论】：

非常感谢您抽出宝贵时间做出回应。我会研究它，我相信我会学到很多东西。我应该更清楚我的犹豫（我正在学习编码，甚至谈论编码都需要非常精确）。您的绿点和蓝点都导致类似的拟合方程的事实是问题所在。我的目标是找出生物过程中每个步骤对输出噪声的相对贡献。通过将真实数据强制输入模型，我希望看到它在最终输出中产生了多少噪音。我可能需要重新考虑我的方法

@BobbyM：绿点没有用于合身。它们代表“真实”数据，即没有任何测量误差的数据。在现实生活中，这是不现实的，你总会有测量误差，生物变异（基因表达是随机的）。为了模仿这一点，我在由蓝点表示的真实数据中添加了噪声。红线是通过仅考虑蓝点（即测量数据）生成的。可以看出，尽管有噪音，但可以很好地再现真实数据。这就是我们在生物学中一直在做的事情：处理噪音:)

@BobbyM：研究进展如何？对上面的代码有任何疑问吗？它解决了你的问题吗？

嗨，对不起，我没有早点回复，上周我没有解决这个问题（很多期中评分）。我现在了解您的解决方案并且代码运行良好。你对测量误差的评论让我思考了很多。虽然我在查看实验数据时经常考虑测量误差，但我从未想过在建模过程中也需要考虑它。我仍在努力寻找最好的方法。非常感谢您的帮助。

【参考方案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)