如何用JiTCDDE求解多维延迟微分方程？含具体系统实现疑问

你选JiTCDDE完全没问题——它就是专门处理这类包含多个状态变量延迟的微分方程的工具，你的困惑主要是没掌握多变量延迟的调用方式。咱们一步步来解决这个问题：

第一步：明确方程结构 #

先把你的延迟微分方程整理清楚（修正了原方程里的笔误）：
$$
\left{
\begin{array}{l}
\dot{v}(t) = y(t) \
\dot{y}(t) = \frac{a_1\alpha}{\omega_1} \cdot y(t-\tau) \cdot \left(1-\tanh^2[v(t-\tau)]\right) - v(t) - \frac{1}{Q_1} \cdot y(t)
\end{array}
\right.
$$
参数：$a_1=70,\ Q_1=50,\ \omega_1=2260,\ \alpha=10,\ \tau \in [0,8\times10^{-3}]$

在JiTCDDE中，我们可以把状态变量定义为一个向量：$x[0] = v(t)$，$x[1] = y(t)$，这样就能统一处理两个变量的当前值和延迟值。

第二步：用JiTCDDE定义方程组 #

核心是用.past(τ)方法调用对应状态变量的延迟值，不管是$y(t-\tau)$还是$v(t-\tau)$，都可以通过这个方法获取。下面是完整的代码框架：

from jitcdde import jitcdde, y
import numpy as np
from matplotlib import pyplot as plt

# 定义参数
a1 = 70
Q1 = 50
omega1 = 2260
alpha = 10

# 定义微分方程组
def system():
    # x[0] = v(t), x[1] = y(t)
    yield y(1)  # 对应v的导数：dv/dt = y(t)
    yield (a1*alpha/omega1) * y(1, tau) * (1 - np.tanh(y(0, tau))**2) - y(0) - y(1)/Q1  # 对应y的导数

# 准备绘制分岔图的数据
tau_values = np.linspace(0, 8e-3, 50)  # 取50个τ值，可根据精度调整
v_bifurcation = []

for tau in tau_values:
    # 初始化DDE求解器，传入方程组和延迟值τ
    DDE = jitcdde(system, delays=[tau])
    
    # 设置初始历史：假设初始时刻v和y都为0（可根据实际需求调整）
    DDE.constant_past([0.0, 0.0])
    
    # 先积分一段时间让系统达到稳态/周期态
    DDE.adjust_diff()
    DDE.integrate_blindly(0.1)  # 先积分0.1秒，足够系统稳定
    
    # 收集稳态下v的极值（或均值，根据需求）
    v_values = []
    t_start = DDE.t
    for _ in range(100):
        t = DDE.t + 0.001
        DDE.integrate(t)
        v_values.append(DDE.y[0])
    
    # 提取极值作为分岔点数据
    v_max = max(v_values)
    v_min = min(v_values)
    v_bifurcation.append( (v_min, v_max) )

# 绘制分岔图
plt.figure(figsize=(10,6))
for tau, (v_min, v_max) in zip(tau_values, v_bifurcation):
    plt.scatter(tau, v_min, color="blue", s=5)
    plt.scatter(tau, v_max, color="blue", s=5)

plt.xlabel(r"$\tau$ (s)")
plt.ylabel(r"$v(t)$")
plt.title("Bifurcation Diagram of the DDE System")
plt.grid(True)
plt.show()

关键细节解释 #

• 多变量延迟的调用：y(0, tau)就是$v(t-\tau)$，y(1, tau)就是$y(t-\tau)$，JiTCDDE会自动处理不同状态变量的延迟历史。
• 初始历史设置：constant_past是最简单的初始条件，如果你需要更复杂的初始历史，可以用past_from_function自定义。
• 稳态判断：integrate_blindly先让系统运行一段时间，跳过暂态过程，之后收集的数据才是稳态/周期解的信息。
• 分岔数据收集：对于每个τ，我们收集v的最大值和最小值，这样能展示分岔时的周期倍化或混沌行为。

常见问题排查 #

如果运行时出现错误，先检查：

• 延迟值τ的单位是否和时间积分的单位一致（这里都是秒，没问题）
• 初始历史是否合理，比如如果系统有固有稳态，初始值可以设为该稳态值
• 积分时间是否足够长，确保系统已经脱离暂态

内容的提问来源于stack exchange，提问作者A.MdP