我有这组微分方程：
dy/dt = a*y - b*x*y
dx/dt = b*x*y - c*y(t - t_0)
t_0 是一个常数时间，当t<t_0 时该术语被忽略。在给定初始条件和所有系数的情况下，如何在 python 中使用 numpy/scipy 解决这个问题？
编辑：y(t-t_0) 是当时y 的值t-t_0，而不是y 乘以t-t_0
【问题讨论】：
t<t_0，忽略哪个词？ x' 是指相对于t 的派生词吗？你的参数在什么范围内，你想在什么时候评估这个？
在问题的早期版本中，问题仅说明了一个简单的 ODE 系统。然后改为延迟微分方程，下面的答案不再有效。我把它留作将来参考。
要解决系统延迟问题，必须使用额外的 python 包。例如包JiTCDDE 允许求解这种方程。 这里提出了一个相关问题：Solve ODE in Python with a time-delay
scipy 函数 ode 可能就是你要找的东西：
让我们首先定义两个方程组。一种用于t<t0，另一种用于t>t0。我们称这些函数为f 和f2。此外，我们还计算雅可比矩阵，供积分器使用。
def f(t,y,a,b,c,t_0):
return [b*y[0]*y[1]-c*(t-t_0), a*y[1]-b*y[0]*y[1]]
def f2(t,y,a,b,c):
return [b*y[0]*y[1], a*y[1]-b*y[0]*y[1]]
def jac_f(t,y,a,b):
return [[b*y[1],b*y[0]],[-b*y[1],a-b*y[1]]]
然后我们导入ode，调用积分器两次。第一次，我们从起始值（我将其设置为 t=0）积分直到达到t0，然后使用对t>t0 有效的方程系统开始第二次积分。我们将最后计算的值作为初始条件传递给积分器并继续积分，直到达到t=4（任意选择）。
from scipy.integrate import ode
y_res = []
t_list = []
a,b,c,t0=1,1,1,1
y0=[1,2]
t_start=0
t_fin=t0
dt=0.01
r=ode(f2,jac_f).set_integrator("vode", method="adams", with_jacobian=True)
r.set_initial_value(y0, t_start).set_f_params(a,b).set_jac_params(a,b)
while r.successful() and r.t < t_fin:
r.integrate(r.t+dt)
y_res.append(r.y)
t_list.append(r.t)
y0=y_res[-1]
t_start=t0
t_fin=4
dt=0.01
r=ode(f,jac_f).set_integrator("vode", method="adams", with_jacobian=True)
r.set_initial_value(y0, t_start).set_f_params(a,b,c,t0).set_jac_params(a,b)
while r.successful() and r.t < t_fin:
r.integrate(r.t+dt)
y_res.append(r.y)
t_list.append(r.t)
我们现在可以绘制结果曲线：
import matplotlib.pyplot as plt
yy=np.stack(y_res)
plt.plot(t_list, yy[:,0], label="x(t)")
plt.plot(t_list, yy[:,1], label="y(t)")
plt.legend()
plt.show()
【讨论】：
c*max(0,t-t0)。由于它是连续的，求解器在该点上积分应该没有问题，即使它必须将步长调整到 t0 左右。
似乎对全局变量 sol_y 执行插值也可以：
import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
def dudt(t, u, params):
x, y = u
a, b, c, t0 = params
dydt = a*y - b*x*y
if t <= t0:
dxdt = b*x*y
else:
dxdt = b*x*y - c*get_old_y(t-t0)
return [dxdt, dydt]
def get_old_y(old_t):
return np.interp(old_t, sol_t, sol_y)
def jac_dudt(t, u, params):
x, y = u
a, b, c, t0 = params
jac = [[ b*y, b*x-c],
[-b*y, a-b*y]]
return jac
# parameters
t0 = 1
params = 1, 1, 2, t0
u0 = [1, 2]
t_end = 3*t0
dt = 0.05
# integration
r = ode(dudt, jac_dudt).set_integrator("vode",
method="adams",
with_jacobian=True)
r.set_initial_value(u0, 0).set_f_params(params).set_jac_params(params)
sol_t, sol_x, sol_y = [], [], []
while r.successful() and r.t < t_end:
r.integrate(r.t + dt)
sol_x.append(r.y[0])
sol_y.append(r.y[1])
sol_t.append(r.t)
# graph
plt.plot(sol_t, sol_x, '-|', label='x(t)')
plt.plot(sol_t, sol_y, '-|', label='y(t)')
plt.legend(); plt.xlabel('time'); plt.ylabel('solution');
带有示例参数的输出图是：
【讨论】：