时滞微分方程

时滞微分方程（delay differential equation）把过去的状态放进当前的变化规律。它适合描述信号传输、运输、成熟期等需要时间的过程。与普通微分方程相比，确定未来通常需要一段已知历史。

温控器收到的消息晚了两分钟

设一个理想化温控系统的目标温度固定。令 $x (t)$ 表示实际温度减目标温度，单位为摄氏度，时间 $t$ 用分钟计。正值表示偏热，负值表示偏冷。为单独研究反馈延迟，假设调节产生的净温度变化率正比于检测到的温差，方向相反；这个线性模型只描述目标附近的小偏差。

如果测量立即到达，反馈规律是 $x^{'} = - a x (t)$ ，其中 $a > 0$ 的单位为每分钟。从 $x (0) = x_{0}$ 出发，解为 $x_{0} e^{- a t}$ ，温差单调衰减。

若检测与执行总共落后 $τ$ 分钟，控制器使用的是 $x (t - τ)$ ，于是 $x^{'} (t) = - a x (t - τ) .$ 取教学参数 $a = 0.4$ 每分钟、 $τ = 2$ 分钟。在时刻0，如果控制器收到两分钟前的温差为0.8度，则当前变化率是 $- 0.32$ 度每分钟。即使当前已经不热，它仍可能继续按旧信息降温。

图中向右的箭头表示时间经过，反馈箭头把过去的测量送到当前。 $τ$ 是模型中的延迟，不是数值求解时随意选择的步长。

时间轴上过去t减tau的测量经过固定延迟来到当前t，当前温差变化率由过去温差乘负a确定 — 反馈作用在现在，依据却来自过去。

同一个当前值，为何会有不同未来

在最初两分钟内，方程不断读取区间 $[- 2, 0]$ 上的温差。因此除了 $x (0)$ ，还要指定 $x (s) = ϕ (s), - τ \leq s \leq 0.$ 函数 $ϕ$ 称为初始历史。比较两种历史： $ϕ_{1} (s) = 0.8, ϕ_{2} (s) = 0.8 (1 + s / 2), - 2 \leq s \leq 0.$ 两者在当前都等于0.8，但两分钟前分别是0.8与0。故 ${x_{1}}^{'} (0 +) = - 0.32$ ，而 ${x_{2}}^{'} (0 +) = 0$ 。当前值相同，当前变化率已经不同。

对第一种历史， $0 \leq t \leq 2$ 时有 $x_{1} (t) = 0.8 - 0.32 t$ 。对第二种历史， $ϕ_{2} (t - 2) = 0.4 t$ ，所以 ${x_{2}}^{'} = - 0.16 t$ ，积分得 $x_{2} (t) = 0.8 - 0.08 t^{2}$ 。两分钟后分别为0.16与0.48度。这个差别由历史决定。

两段历史在零时刻同为零点八，常值历史之后得到直线下降，斜坡历史之后得到抛物线下降，两分钟后分别为零点一六和零点四八 — 虚线表示事先给定的历史，实线表示方程决定的未来。

把过去算好的部分接着使用

对于常值历史，前两分钟已经得到 $x (t) = 0.8 - 0.32 t$ 。在 $2 \leq t \leq 4$ 内，延迟项读取的正是这段解： $x^{'} (t) = - 0.4 [0.8 - 0.32 (t - 2)] = - 0.32 + 0.128 (t - 2) .$ 从 $x (2) = 0.16$ 积分，得到 $x (t) = 0.16 - 0.32 (t - 2) + 0.064 (t - 2)^{2}, 2 \leq t \leq 4.$ 于是 $x (4) = - 0.224$ 。温差已经越过零，这和没有延迟的单调回归不同。继续计算 $[4, 6]$ 时，再将刚求出的二次多项式作为已知延迟输入。

这种分步法以一个延迟长度为区间，把已经知道的过去用于下一段。在更一般的方程 $x^{'} (t) = F (t, x (t), x (t - τ))$ 中，下一段仍可能需要求一个关于当前状态的普通微分方程；已知的是延迟项，并非整个右端。Marc R. Roussel 的时滞方程讲义讨论了历史、分步法和线性化。

如果 $F$ 连续并对状态变量局部 Lipschitz 连续，连续历史可以在这些条件成立的范围内逐段确定唯一解。全局存在还需排除有限时间发散。历史在0处的左导数不必符合方程给出的右导数，因此求解器还需要处理由接点传播的光滑性变化。

真正控制响应的是反应速度乘以延迟

令 $u = t / τ$ 、 $z (u) = x (τ u) / 0.8$ ，则链式法则给出 $z^{'} (u) = - q z (u - 1), q = a τ .$ $q$ 没有单位。它衡量控制器在“等待一个延迟”的时间内会作出多强的修正。增大反应强度或增加延迟，都可能把同一类反馈推向振荡。

常值历史 $z (u) = 1$ （ $- 1 \leq u \leq 0$ ）下，可比较 $q = 0.4, 1.2, 1.8$ 。前两组的偏差最终衰减，但衰减过程形状不同；第三组出现增长的振荡。下图中的水平时间单位是一个延迟长度。

三个上下排列的图用相同初始历史比较q为零点四、一点二和一点八，前两组衰减，第三组振荡放大 — 振荡可以来自模型自身的延迟；判断是否由数值误差造成，还需缩小计算步长。

临界值为什么与圆周率有关

寻找指数形式的解 $z (u) = e^{λ u}$ ，代入得 $λ + q e^{- λ} = 0.$ 这称为特征方程。指数项使它通常具有无穷多个复根。根的实部决定指数模式增长还是衰减，虚部决定振荡。

先证明 $0 < q < π / 2$ 时没有实部非负的根。设 $λ = α + i β$ 且 $α \geq 0$ ，分离实虚部： $α = - q e^{- α} \cos β, β = q e^{- α} \sin β .$ 第二式给出 $| β | \leq q < π / 2$ ，故 $\cos β > 0$ 。第一式却要求 $α < 0$ ，矛盾。对这个常系数滞后型方程，特征根全在左半平面的谱判据因而给出零解的渐近稳定性。

临界时令 $λ = i ω$ 、 $ω > 0$ ，得到 $q \cos ω = 0, ω = q \sin ω .$ 第一对虚轴根出现在 $ω = q = π / 2$ 。此时有不衰减的正弦模式，因此零解不再渐近稳定。沿根分支对 $λ e^{λ} = - q$ 求导， $\frac{d λ}{d q} = \frac{λ}{q (1 + λ)} .$ 在 $λ = i π / 2$ 处，实部为正；随 $q$ 增大，这对根进入右半平面。其余虚轴穿越发生在 $q = π / 2 + 2 k π$ （非负整数 $k$ ），穿越方向也相同。因此本方程的渐近稳定范围恰为 $0 < a τ < π / 2$ ，超过它则不稳定。这是所写线性延迟方程的结论，不是所有控制系统的统一阈值。

如何判断计算出的振荡可信

求解 $x (t - τ)$ 通常需要保存过去的计算点，并对网格间的值插值。简单地把它换成“上一个点”，会把延迟错误地改成求解步长。

本例的头两个区间有明确多项式，可以先用它们检查程序；更长时间则比较步长减半后的曲线。还可以把 $τ = 0$ 的无延迟模型单独求解，检查是否得到指数衰减。有限区间的数值对照支持计算精度，长期稳定性仍由前面的理论判断。

温控系统真实的散热、功率上限和执行器饱和会改变方程。当计算温差增长到不再“小”时，应回到建模假设。延迟还可出现在捕食者-猎物模型的成熟期、库存补货和神经信号中；它们需要各自的机制方程。

来源与继续阅读

Marc R. Roussel，Delay-differential equations：历史函数、分步法、特征方程与根的位置。
先修：一阶线性微分方程、平衡点与稳定性、欧拉公式；后续：分岔、微分方程数值解、动态系统建模。