我最近在研究一个描述生物系统的耦合非线性常微分方程（ODE）系统，具体形式如下：
$$
\begin{cases}
\dfrac{dp}{dt} = -\gamma p f,\
\
\dfrac{df}{dt} = \gamma p f,\
\
\dfrac{dT}{dt} = \left( 1 - \dfrac{k}{T} \right) \left{ b (1 - T) \dfrac{r^n}{fn + r^n} - m \dfrac{f^n}{fn + r^n} \right},
\end{cases}
$$
其中参数取值为：$\gamma = 0.5432$，$k = 0.0026$，$b = 0.0885$，$n = 14.9832$，$r = 0.8265$，$m = 0.9780$。

为了分析系统的稳定性，我首先发现前两个方程可以推出 $p + f$ 是一个非负常数，记为$w$；结合我的参数和初始条件，这里$w$约为0.801。通过消去$p$，原系统可以简化为：
$$
\begin{cases}
\dfrac{df}{dt} = \gamma (w - f) f,\
\
\dfrac{dT}{dt} = \left( 1 - \dfrac{k}{T} \right) \left{ b (1 - T) \dfrac{r^n}{fn + r^n} - m \dfrac{f^n}{fn + r^n} \right}.
\end{cases}
$$

这个简化后的系统有四个$(f, T)$形式的平衡点：
$$(0, k), (0, 1), (w, k), \left( w, 1- \frac{m}{b} (w / r)^n \right).$$

我通过计算雅可比矩阵的特征值得出：前两个平衡点是不稳定的，后两个则是稳定平衡点。

但实际模拟时却出现了矛盾：我用Mathematica对原系统进行模拟，无论设置什么初始条件，解似乎永远不会收敛到第四个平衡点 $\left( w, 1- \frac{m}{b} (w / r)^n \right)$，只会收敛到$(w,k)$。我使用的模拟代码如下：

s = NDSolve[{p'[t] == -gamma p[t] f[t], f'[t] == gamma p[t] f[t], 
T'[t] == b (1 - T[t]) (r^n/(f[t]^n + r^n)) (1 - 0.0026/T[t]) - m (f[t]^n/(f[t]^n + r^n)) (1 - 0.0026/T[t]), 
p[0] == p0, f[0] == 0.001, T[0] == T0}/.{gamma -> 0.5432, b -> 0.0885, m -> 0.9780, r -> 0.8265, n -> 14.9832, p0 -> 0.8, T0 -> 0.8}, {p, f, T}, {t, 0, 720}];

Plot[Evaluate[{f[t], T[t]}/.s], {t, 0, 720}, PlotStyle -> {Blue, Red}, AxesOrigin -> {-1, 0}, PlotRange -> All]

模拟结果的图像显示：代表$f(t)$的蓝色曲线最终趋近于$w≈0.8$，代表$T(t)$的红色曲线趋近于$k=0.0026$，也就是解长期收敛于平衡点$(w,k)$。

*我的核心问题是：*为什么系统永远无法收敛到那个理论上的稳定平衡点？一个系统有可能存在永远无法到达的稳定平衡点吗？

补充说明（原系统转化的原因） #

直接对原系统做稳定性分析会遇到一个问题：原系统的平衡点是$(f \to 0, T \to k)$和$(f \to 0, T \to 1)$，但$p$可以取任意值，这让雅可比矩阵的计算变得模糊不清。而通过转化系统消去$p$后，这个问题就得到了解决。

备注：内容来源于stack exchange，提问作者user1303236