求解前一时间步的ODE函数（延迟微分方程）

Question

我有这组微分方程：

dy/dt = a*y   - b*x*y  
dx/dt = b*x*y - c*y(t - t_0)

t_0是一个恒定的时间，当t<t_0时，这个词被忽略了。在给定初始条件和所有系数的情况下，如何使用numpy / scipy在python中解决这个问题？

编辑：y(t-t_0)是y时t-t_0的值，而不是y倍t-t_0

Answer 1

在该问题的早期版本中，该问题仅表示一个简单的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()

我们得到了这个漂亮的图：

Answer 2

似乎在全局变量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');

带有示例参数的输出图是：