如何在Scipy solve_ivp中为ODE系统的某一组件引入延迟？

问题背景 #

现有一个基于Scipy solve_ivp 的ODE系统，希望为其中ff_s的f_e输入项添加时间延迟（即使用f_e(t-delay)替代当前的f_e）。原模型代码如下：

def model_step(t,
           y,
           w_0, w_1, w_2, w_3, w_4, w_5, w_6, w_7, w_8, w_9, w_10,
           tau_0, tau_1, tau_2, tau_3,
           threshold,
           power,
           q,
           i_0, i_1, i_2, i_3
           ):
    """
    Builds differential equation model
    """
    
    f_e, f_p, f_s, f_v = y

    ff_e = (min(max((i_0 + w_0 * f_e + w_1 * f1(t) + w_2 * f_p + w_3 * f_s), threshold), 25))
    ff_p = (min(max((i_1 + w_4 * f_e + f2(t)), threshold), 25))
    ff_s = (min(max((i_2 + w_5 * f_e + w_6 * f3(t) - w_7 * f_v), threshold), 25))
    ff_v = (min(max((i_3 + w_8 * f_e - w_9 * f_s +  w_10 * f4(t)), threshold), 25))

    df_edt = ((q * ff_e ** power) - f_e) / tau_0
    df_pdt = ((q * ff_p ** power) - f_p) / tau_1
    df_sdt = ((q * ff_s ** power) - f_s) / tau_2
    df_vdt = ((q * ff_v ** power) - f_v) / tau_3

    dydt = [df_edt, df_pdt, df_sdt, df_vdt]

    return dydt

目标修改：

ff_s = (min(max((i_2 + w_5 * f_e(t-delay) + w_6 * f3(t) - w_7 * f_v), threshold), 25))

由于solve_ivp仅依赖当前时间点的状态y，需要找到合适的方法引入延迟项。

解决方案 #

方法一：手动记录历史状态并插值（适配solve_ivp） #

solve_ivp本身不支持延迟微分方程（DDE），但可以通过手动维护历史时间和状态数据，在计算延迟项时用插值得到f_e(t-delay)的值。步骤如下：

• 初始化历史缓存：求解开始前，预定义延迟时间段内的初始状态（比如t < 0时的f_e值，通常用初始条件或常数）。
• 用类封装历史数据：避免使用全局变量，用类来记录所有已求解的时间点t和对应的f_e值。
• 计算延迟时刻的状态：当t >= delay时，从历史数据中插值获取t-delay时刻的f_e；若t < delay，直接使用初始预设的延迟状态。
• 插值保证精度：推荐用线性插值或样条插值（如scipy.interpolate.interp1d）估算延迟时刻的f_e，减少离散点误差。

示例代码：

import numpy as np
from scipy.interpolate import interp1d
from scipy.integrate import solve_ivp

class DelayedModel:
    def __init__(self, delay, initial_y, f1, f2, f3, f4):
        self.delay = delay
        # 初始化历史：t=0时的状态
        self.t_history = [0.0]
        self.fe_history = [initial_y[0]]
        self.f1 = f1
        self.f2 = f2
        self.f3 = f3
        self.f4 = f4

    def get_delayed_fe(self, t):
        t_delayed = t - self.delay
        if t_delayed <= 0:
            # t小于延迟时间时返回初始状态的f_e
            return self.fe_history[0]
        # 线性插值获取延迟时刻的f_e
        interp_func = interp1d(self.t_history, self.fe_history, kind='linear', fill_value="extrapolate")
        return interp_func(t_delayed)

    def model_step(self, t, y, *params):
        w_0, w_1, w_2, w_3, w_4, w_5, w_6, w_7, w_8, w_9, w_10, tau_0, tau_1, tau_2, tau_3, threshold, power, q, i_0, i_1, i_2, i_3 = params
        
        f_e, f_p, f_s, f_v = y
        # 更新历史数据
        self.t_history.append(t)
        self.fe_history.append(f_e)

        # 获取延迟后的f_e
        fe_delayed = self.get_delayed_fe(t)

        ff_e = (min(max((i_0 + w_0 * f_e + w_1 * self.f1(t) + w_2 * f_p + w_3 * f_s), threshold), 25))
        ff_p = (min(max((i_1 + w_4 * f_e + self.f2(t)), threshold), 25))
        ff_s = (min(max((i_2 + w_5 * fe_delayed + w_6 * self.f3(t) - w_7 * f_v), threshold), 25))
        ff_v = (min(max((i_3 + w_8 * f_e - w_9 * f_s +  w_10 * self.f4(t)), threshold), 25))

        df_edt = ((q * ff_e ** power) - f_e) / tau_0
        df_pdt = ((q * ff_p ** power) - f_p) / tau_1
        df_sdt = ((q * ff_s ** power) - f_s) / tau_2
        df_vdt = ((q * ff_v ** power) - f_v) / tau_3

        return [df_edt, df_pdt, df_sdt, df_vdt]

# 使用示例
delay = 1.0
initial_y = [1.0, 0.0, 0.0, 0.0]
delayed_model = DelayedModel(delay, initial_y, f1, f2, f3, f4)

# 按顺序传入所有参数
params = (w_0, w_1, w_2, w_3, w_4, w_5, w_6, w_7, w_8, w_9, w_10, tau_0, tau_1, tau_2, tau_3, threshold, power, q, i_0, i_1, i_2, i_3)

# 调用solve_ivp，可设置max_step限制步长提升精度
sol = solve_ivp(delayed_model.model_step, [0, 10], initial_y, args=params, method='RK45', max_step=0.1)

方法二：使用专门的延迟微分方程求解器 #

更可靠的方式是用专门针对DDE的求解器，无需手动维护历史，求解器会自动处理延迟项，精度和稳定性更有保障。比如Scipy自带的ddeint：

from scipy.integrate import ddeint
import numpy as np

def dde_model(y, t, delay, params, f1, f2, f3, f4):
    w_0, w_1, w_2, w_3, w_4, w_5, w_6, w_7, w_8, w_9, w_10, tau_0, tau_1, tau_2, tau_3, threshold, power, q, i_0, i_1, i_2, i_3 = params
    
    # y(t)是当前状态，y(t-delay)是延迟状态
    f_e, f_p, f_s, f_v = y(t)
    fe_delayed = y(t - delay)[0] if t >= delay else y(0)[0]

    ff_e = (min(max((i_0 + w_0 * f_e + w_1 * f1(t) + w_2 * f_p + w_3 * f_s), threshold), 25))
    ff_p = (min(max((i_1 + w_4 * f_e + f2(t)), threshold), 25))
    ff_s = (min(max((i_2 + w_5 * fe_delayed + w_6 * f3(t) - w_7 * f_v), threshold), 25))
    ff_v = (min(max((i_3 + w_8 * f_e - w_9 * f_s +  w_10 * f4(t)), threshold), 25))

    df_edt = ((q * ff_e ** power) - f_e) / tau_0
    df_pdt = ((q * ff_p ** power) - f_p) / tau_1
    df_sdt = ((q * ff_s ** power) - f_s) / tau_2
    df_vdt = ((q * ff_v ** power) - f_v) / tau_3

    return [df_edt, df_pdt, df_sdt, df_vdt]

# 定义初始历史函数：t <= 0时的状态
def history(t):
    return initial_y  # initial_y是初始状态列表

# 使用ddeint求解
t = np.linspace(0, 10, 100)
sol = ddeint(dde_model, history, t, args=(delay, params, f1, f2, f3, f4))

内容的提问来源于stack exchange，提问作者Yehor K.