如何求解具有 scipy.integrate.odeint 且具有随时间变化的变量的 ODE 系统
How to solve a system of ODEs with scipy.integrate.odeint with a time-dependent variable
我正在使用 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 个数字。
您无法先验地知道数值积分器在哪些点计算 ODE 函数。积分器（odeint 和其他未明确 "fixed step-size" 的积分器）动态生成一个内部点列表，该列表可能比给定的采样点列表具有更小或有时更大的步长。输出值是从内部列表中插入的。
如果您想用函数替换 ODE 的一部分，则必须将示例数据转换为函数。这可以通过插值来完成。使用 scipy.interpolate.interp1 函数生成函数对象，然后您可以像使用任何其他标量函数一样使用这些对象。
据我从你问题下方的评论中了解到，你试图合并可能有噪音的测量数据。您可以使用这些数据来适应您的时间进程，而不是直接插入数据。在这里，我显示变量 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()
要具体执行我在问题中提出的问题，即使用值代替 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)
我正在使用 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 个数字。
您无法先验地知道数值积分器在哪些点计算 ODE 函数。积分器（odeint 和其他未明确 "fixed step-size" 的积分器）动态生成一个内部点列表，该列表可能比给定的采样点列表具有更小或有时更大的步长。输出值是从内部列表中插入的。
如果您想用函数替换 ODE 的一部分，则必须将示例数据转换为函数。这可以通过插值来完成。使用 scipy.interpolate.interp1 函数生成函数对象，然后您可以像使用任何其他标量函数一样使用这些对象。
据我从你问题下方的评论中了解到，你试图合并可能有噪音的测量数据。您可以使用这些数据来适应您的时间进程，而不是直接插入数据。在这里，我显示变量 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()
要具体执行我在问题中提出的问题，即使用值代替 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)