非线性系统入门与综合建模

上一章，我们把一个平面系统看成一张“水流地图”：矩阵给每个位置配上速度箭头，特征值告诉我们原点附近是向内汇聚、向外散开，还是绕着转。这套办法很好用，但它还有一个尚未回答的问题：真实世界的箭头，为什么一定要由同一个矩阵决定？

猎物增加，会让捕食者更容易找到食物；捕食者增加，又会改变猎物的增长速度。两种数量的作用乘在一起，就出现了 $xy$。一根摆杆偏离竖直方向，重力产生的回复作用取决于 $\sin\theta$，也不可能在所有角度都等于 $\theta$。这些变化没有推翻前面的知识，只是提醒我们：一个线性系统往往是一幅局部放大的图，完整地图还需要保留非线性项。

这一章我们先弄清楚这幅局部图怎么画、能相信到什么程度，再用捕食者与猎物、摆的运动把分析做完整。最后回到一个带环境容量的生态模型，从“想问什么”一直走到“计算结果究竟说明了什么”。这也是本课程最后一次把建模、解析计算、几何分析和数值方法放在同一道问题里。

---

先找系统可能停下来的地方

假设我们研究两个竞争同一类资源的种群。食物充足时，数量较多的种群通常也能产生更多新个体；种群过密时，自己和竞争者都会让增长变慢。为了先看清结构，把种群数量和时间按固定尺度换成无量纲变量，考虑一个教学模型：

$$\begin{cases} x'=x(2-x-y),\\ y'=y(3-x-2y). \end{cases}$$

这里 $x,y\geq0$。括号里的 $2-x-y$ 是第一个种群的单位数量净增长率；$3-x-2y$ 是第二个种群的单位数量净增长率。系数只是为了展示分析过程选定的值，并不是某个真实生态区域的测量结果。

与线性系统相比，$x^2$、$xy$、$y^2$ 让箭头随位置的变化不再是线性的。不过，相平面的读法完全没变。写成一般形式就是

$$\begin{cases} x'=f(x,y),\\ y'=g(x,y). \end{cases}$$

右端不显含时间，称为自治系统。同一个状态 $(x,y)$，无论几点来到这里，都得到同一个速度向量 $(f(x,y),g(x,y))$。这正是我们能把运动画在一张固定地图上的原因。

本章假定向量场在所研究区域内连续可微。这保证每个初值附近有唯一的局部解，不同轨线不能在普通点相交后各走各的。它不保证任意非线性方程的解都能存在到无限远；要作长期判断，还得检查解会不会逃出定义域或者在有限时间内变得无界。

零增长曲线怎样帮助我们读图

先看第一个种群什么时候暂时不变。$x'=0$ 给出

$$x=0\quad\text{或}\quad x+y=2.$$

同样，$y'=0$ 给出

$$y=0\quad\text{或}\quad x+2y=3.$$

这些线叫零增长曲线，也常叫零轨线。第一组上水平速度为零，第二组上竖直速度为零。只有两组同时满足，整个状态才停下来。因此把两组逐一相交，得到四个平衡点：

$$(0,0),\qquad(2,0),\qquad\left(0,\frac32\right),\qquad(1,1).$$

前面三个在边界上，分别表示全部消失、只剩第一种、只剩第二种；最后一个才是共存。求解时如果上来就把两式除以 $x$ 和 $y$，三个边界平衡就会全丢掉。我们在分离变量时遇到过同样的问题：除去一个因子之前，要先安顿好它等于零的情形。

接着看箭头。在正象限内，$x+y<2$ 时 $x'>0$，箭头向右；$x+y>2$ 时箭头向左。 时 ，箭头向上； 时箭头向下。例如在 ，

$$x'=0.8(2-0.8-0.6)=0.48,\qquad y'=0.6(3-0.8-1.2)=0.6.$$

此刻两种数量都增加。在 $(1.2,0.8)$，则有 $x'=0$、$y'=0.16$，箭头竖直向上。它通常会穿过直线 $x$，随后水平速度改变符号。因此，；它更像地图上提示“某个速度分量在这里换号”的分界标记。坐标轴在这个模型中确实是轨线所在的不变边界，因为消失的种群不会凭空重新出现。

---

把平衡点挪到原点，再看最先起作用的项

现在问一个更具体的问题：如果两个种群原本都在 $1$ 附近，稍微多一点或少一点，会不会自己恢复？我们不必先解出全部轨线，可以先只看共存点周围。

令

$$u=x-1,\qquad v=y-1.$$

$u,v$ 是相对平衡的偏差。平衡数量不随时间变化，所以 $u'=x'$、$v'=y'$。把 、 真正代进原方程，而不急着背矩阵公式：

$$\begin{aligned} u'&=(1+u)(-u-v)=-u-v-u^2-uv,\\ v'&=(1+v)(-u-2v)=-u-2v-uv-2v^2. \end{aligned}$$

当偏差都只有原尺度的百分之一时，一次项是百分之一量级，乘积项则是万分之一量级。先保留一次项，就得到

$$\begin{pmatrix}u'\\v'\end{pmatrix} \approx \begin{pmatrix}-1&-1\\-1&-2\end{pmatrix} \begin{pmatrix}u\\v\end{pmatrix}.$$

这就是线性化的来历：我们在平衡点附近，先保留对小偏差起主要作用的线性部分。原方程仍在那里，二次项没有真的消失。

一般地，若平衡点为 $(x_*,y_*)$，定义偏差向量

$$\mathbf{u}=\begin{pmatrix}x-x_*\\y-y_*\end{pmatrix}.$$

多元函数的一阶展开给出

$$\mathbf{u}'=J_*\mathbf{u}+\mathbf{R}(\mathbf{u}), \qquad \frac{\|\mathbf{R}(\mathbf{u})\|}{\|\mathbf{u}\|}\longrightarrow0 \quad(\mathbf{u}\longrightarrow\mathbf{0}),$$

其中

$$J_*=J(x_*,y_*),\qquad J(x,y)= \begin{pmatrix} f_x(x,y)&f_y(x,y)\\ g_x(x,y)&g_y(x,y) \end{pmatrix}.$$

$J$ 叫 Jacobian 矩阵，中文也叫雅可比矩阵。第一列回答“只把 $x$ 增加一点，两种变化率分别怎么动”，第二列回答同样的 $y$ 扰动问题。偏导数因而有明确的模型含义，并不只是计算矩阵的手续。若右端有连续二阶导数，余项通常还能估成 $O(\|\mathbf{u}\|^2)$；只有连续一阶导数时，用上面的更一般表达就够了。

对于竞争模型，

$$J(x,y)= \begin{pmatrix} 2-2x-y&-x\\ -y&3-x-4y \end{pmatrix}.$$

代入 $(1,1)$，正好得到刚才逐项展开的矩阵。它的特征方程是

$$\lambda^2+3\lambda+1=0, \qquad \lambda_{1,2}=\frac{-3\pm\sqrt5}{2}.$$

两根都负。稍后说明的线性化结论允许我们判断：共存平衡局部渐近稳定。这里“局部”修饰初值范围，“渐近”表示随着时间增加逐步回去，并不承诺某个有限时刻突然精确回到平衡。

另外三个点也可以一起算完，避免把模型只研究一半：

两个边界平衡都是鞍点，这在情境里有一个直观读法：只剩一种的状态沿着自身边界可以恢复，但另一种只要以很小的正数量进入，就有增长的可能。讨论实际种群时只取非负方向的扰动；这里的不稳定方向确实进入可行区域，并非来自负种群这种无意义的状态。

---

哪些判断能从线性系统带回来

先把两个容易混的词分开。稳定表示初值足够近，就能一直待在指定的小邻域内；渐近稳定还要求这些近处的解最终趋向平衡。绕着平衡打转而不靠近，可以稳定，却不渐近稳定。距离也不必每时每刻都减小，局部渐近稳定的轨线允许先有一段偏移再回落。

若 $J_*$ 的所有特征值实部都不为零，这个平衡点叫双曲平衡点。对本章假定的连续可微自治系统，线性化能给出可靠的局部稳定性结论：全部实部为负，原平衡局部渐近稳定；存在实部为正的特征值，原平衡不稳定。二维一正一负时是鞍点，附近有分别趋近和远离的方向结构。

我们沿用上一章的迹与行列式，记

$$T=\operatorname{tr}J_*,\qquad D=\det J_*,\qquad \Delta=T^2-4D.$$

由 $\lambda^2-T\lambda+D=0$，可把常见情形放在一起：

最后一行并不是说任何结论都没有：如果另一个根为正，仍可判断不稳定。但若根为 $0$ 和负数，单凭这些根不能判定稳定性。也别把“迹为零”一律当成不能判断；$T=0,D<0$ 明明还是鞍点。

为什么边界情形格外麻烦？因为线性项恰好没有提供向内或向外的主导作用，那些看似很小的非线性项就有机会在长时间里累积起来。

同一个线性中心，可以藏着三种答案

用已经按尺度化简的坐标，考虑下面一族系统，$\sigma$ 可取 $-1,0,1$：

$$\begin{cases} x'=-y+\sigma x(x^2+y^2),\\ y'=x+\sigma y(x^2+y^2). \end{cases}$$

三种情形在原点的 Jacobian 完全相同，都是

$$\begin{pmatrix}0&-1\\1&0\end{pmatrix}, \qquad\lambda=\pm i.$$

但令 $\rho=\sqrt{x^2+y^2}$，在 $\rho>0$ 处直接求导：

$$\rho\rho'=xx'+yy'=\sigma(x^2+y^2)^2, \qquad \rho'=\sigma\rho^3.$$

交叉项 $-xy+xy$ 抵消了，留下的恰恰是高阶项。若 $\sigma=-1$，分离变量得到

$$\rho(t)=\frac{\rho_0}{\sqrt{1+2\rho_0^2t}},$$

半径趋向零，原点渐近稳定。若 $\sigma=0$，半径恒定，得到稳定但不渐近稳定的中心。若 $\sigma=1$，则

$$\rho(t)=\frac{\rho_0}{\sqrt{1-2\rho_0^2t}}, \qquad 0\leq t<\frac1{2\rho_0^2},$$

任意小的正半径最终都会扩大，原点不稳定。这三种情况的角速度都等于 $1$，所以向内螺旋、原地绕圈、向外螺旋，全都能躲在同一个线性中心后面。

零根也一样。例如 $x'=-x^3,\ y'=-y$ 与 $x'=x^3,\ y'=-y$ 在原点都有特征值 ，前者的两个变量都趋向零，后者沿 轴向外离开。线性化测试失去决定性时，应该继续分析原方程，而不是把“暂时判不了”写成“稳定”。

---

猎物先变多，捕食者为什么晚一步跟上

现在让两个种群从竞争者变成捕食者与猎物。先讲清假设：没有捕食者时，猎物的食物暂且视为充足，它按比例增长；没有猎物时，捕食者按比例减少；两者混合得足够均匀，相遇频率近似与两种数量的乘积成正比。暂不加入年龄差异、迁入迁出和捕食饱和。

令 $x$ 为猎物数量，$y$ 为捕食者数量，就得到

$$\begin{cases} x'=\alpha x-\beta xy=x(\alpha-\beta y),\\ y'=\delta xy-\gamma y=y(\delta x-\gamma), \end{cases} \qquad \alpha,\beta,\gamma,\delta>0.$$

$\alpha,\gamma$ 的单位都是时间的倒数。若两种数量分别按“只猎物”“只捕食者”记录，那么 $\beta$ 的单位是“每只捕食者每单位时间”，$\delta$ 的单位是“每只猎物每单位时间”。这样 $\beta y$、$\delta x$ 才能与自然增长率、死亡率相减。

正初值不会在有限时间内穿过坐标轴。例如沿解有

$$x(t)=x_0\exp\!\left(\int_0^t[\alpha-\beta y(s)]\,ds\right),$$

只要解存在且积分有限，右边就严格为正；$y$ 有同样的表达。坐标轴初值则要单独保留：$y=0$ 时 $x=x_0e^{\alpha t}$，$$ 时 。

同时令两种增长率为零，平衡为

$$(0,0),\qquad (x_*,y_*)=\left(\frac\gamma\delta,\frac\alpha\beta\right).$$

正象限内，横线 $y=y_*$ 决定猎物增减，竖线 $x=x_*$ 决定捕食者增减。共存点右下方，猎物多、捕食者少，两种数量都上升；进入右上方后，捕食者已经多到使猎物开始减少，但自己暂时仍在增加。轨线继续绕行，便出现了猎物先达到峰值、捕食者稍后达到峰值的循环次序。

线性化给了什么，又留下了什么

先求一般 Jacobian：

$$J(x,y)= \begin{pmatrix} \alpha-\beta y&-\beta x\\ \delta y&\delta x-\gamma \end{pmatrix}.$$

原点的特征值为 $\alpha,-\gamma$，所以是鞍点。这与坐标轴上的指数增长、衰减正好一致。共存点处则有

$$J_*= \begin{pmatrix} 0&-\beta\gamma/\delta\\ \alpha\delta/\beta&0 \end{pmatrix}, \qquad \lambda^2+\alpha\gamma=0.$$

得到 $\lambda=\pm i\sqrt{\alpha\gamma}$，线性化是中心。经过刚才的反例，我们已经知道：到这里还不能宣布原模型有闭合轨线。

为了进一步分析，把数量除以共存水平：

$$p=\frac{x}{x_*},\qquad q=\frac{y}{y_*}.$$

$p,q$ 都无量纲，系统变成

$$p'=\alpha p(1-q),\qquad q'=\gamma q(p-1).$$

在 $p,q>0$ 时，暂时用 $dq/dp=q'/p'$ 消去时间，可以整理出

$$\alpha\left(\frac1q-1\right)dq =\gamma\left(1-\frac1p\right)dp.$$

两边积分，就提示我们考察

$$H(p,q)=\gamma(p-1-\ln p)+\alpha(q-1-\ln q).$$

我们把对数写在无量纲比值上，避免给“带单位的数量取对数”。刚才消去时间在 $p'=0$ 的位置不能直接用，因此最后一定要回到时间变量核查。链式法则给出

$$\begin{aligned} \frac{dH}{dt} &=\gamma\left(1-\frac1p\right)\alpha p(1-q) +\alpha\left(1-\frac1q\right)\gamma q(p-1)\\ &=\alpha\gamma[(p-1)(1-q)+(q-1)(p-1)]\\ &=0. \end{aligned}$$

这个等式在整个正象限都成立，包括零增长线上。$H$ 就是沿轨线不变的量。

守恒为什么真能推出闭合轨线

令 $\phi(z)=z-1-\ln z$。它在 $z>0$ 上满足

$$\phi'(z)=1-\frac1z,\qquad \phi''(z)=\frac1{z^2}>0.$$

所以它在 $z=1$ 处取得唯一最小值 $0$，并且在 $z\to0^+$ 或 $z\to\infty$ 时都趋向正无穷。于是 $H$ 在 有唯一严格最小值，任意有限的等值线都会被限制在远离坐标轴、也远离无穷远的区域里。

更具体地，$H$ 的 Hessian 是正定对角矩阵，故 $H$ 严格凸。每条非零等值线 $H=h$ 都是围住 $(1,1)$ 的光滑闭曲线：它是紧凸子水平集的边界，边界上梯度不为零。系统速度沿这条曲线切向运动，而且曲线上没有平衡点。光滑紧曲线上处处非零的速度有正的下界，走完有限长度的一圈只需有限时间。这样，每个正的非平衡初值都产生周期运动，就有了方程本身的依据。

共存平衡因此稳定，但不渐近稳定。若初值的 $H>0$，守恒已经排除了它趋向 $H=0$ 的平衡。这里也没有某一条特别吸引周围轨线的周期轨道：每条闭曲线旁边都有别的闭曲线，所以它们不是极限环。极限环首先要求周期轨道在邻近周期轨道中是孤立的，不能把“画成一圈”都叫极限环。

取一组教学参数：时间按月计，$\alpha=0.6$、$\gamma=0.3$，$\beta=0.02$、$\delta=0.005$，单位按前述约定。那么

$$x_*=60,\qquad y_*=30,\qquad \lambda=\pm i\sqrt{0.18}.$$

若初值为 $(90,15)$，则

$$x'(0)=90(0.6-0.02\times15)=27,\qquad y'(0)=15(0.005\times90-0.3)=2.25.$$

两种数量最初都增加，但猎物在捕食者达到 $30$ 时就会停止增加。初值对应 $p_0=1.5,q_0=0.5$，守恒量为

$$H_0=0.3(0.5-\ln1.5)+0.6(-0.5-\ln0.5) \approx0.144249.$$

这条轨线会在这一个等值线上循环。线性化给出的近小振幅周期为 $2\pi/\sqrt{0.18}\approx14.81$ 个月；当前初值并不特别接近平衡，不能把这个近似周期直接当成它的精确周期。

数学上正数量永不归零，不表示真实种群绝无灭绝风险。如果模型算出数量已经远小于一个个体，连续变量的近似就失去直接解释。更不用说食物有限、个体捕食能力有限，这些假设一旦修改，刚才的守恒量和周期结论都要重新核查。

---

摆的能量把局部图接成完整图

把一个刚性细杆的摆稍微推开，它来回摆动；用力推得足够快，它可能越过最高点连续转圈。这两个现象都应当由同一个模型解释。假定杆长 $L$、质量集中在端点，转轴固定，先忽略空气阻力和摩擦。角度 $\theta$ 从竖直向下量起，重力沿切线方向的分量把它拉回去，于是

$$\theta''+\frac gL\sin\theta=0.$$

这里 $t$ 以秒计，$\theta$ 用弧度，$g/L$ 的单位为 $\mathrm{s}^{-2}$。令 $\omega=\theta'$，得到

$$\theta'=\omega,\qquad \omega'=-\frac gL\sin\theta.$$

平衡点是 $(n\pi,0)$。在展开的相平面里它们有无穷多个；若把角度相差 $2\pi$ 的位置视为相同物理位置，相空间可看成一个圆柱。在最低点 $(2n\pi,0)$ 附近，线性化特征值是 $\pm i\sqrt{g/L}$；在最高点 附近，特征值是 ，所以最高点是鞍点、不稳定。

最低点的纯虚根仍需进一步确认。这里最自然的量是机械能除以 $mL^2$：

$$E(\theta,\omega)=\frac12\omega^2+\frac gL(1-\cos\theta).$$

第一项来自动能，第二项来自相对最低点的重力势能；这个经过缩放的 $E$ 单位是 $\mathrm{s}^{-2}$，实际机械能为 $mL^2E$。沿解求导，

$$E'=\omega\omega'+\frac gL\sin\theta\,\theta' =\omega\left(-\frac gL\sin\theta\right) +\frac gL\sin\theta\,\omega=0.$$

最低点是能量的严格局部最小点，附近的小正能量等值线闭合，故它是稳定中心，不渐近稳定。理想摆没有任何消耗能量的项，也就没有理由越摆越小。

从最低点赋予初速度 $\omega_0$，初始能量是 $E_0=\omega_0^2/2$；到达最高点至少需要势能 。所以临界角速度是

$$|\omega_c|=2\sqrt{\frac gL}.$$

当 $0<|\omega_0|<|\omega_c|$ 时，能量不足以越顶，摆动到 $\omega=0$ 后折返，最大偏角满足

$$\theta_{\max}=2\arcsin\!\left(\frac{|\omega_0|}{2\sqrt{g/L}}\right).$$

大于临界速度时，$\omega$ 永不变号，摆连续转圈。恰好等于临界速度时，轨线属于分界轨线：摆无限接近最高点，速度无限趋近零，却不会在有限时间到达平衡后再穿过去。若有限时间到达，唯一性会迫使它一直是最高点的常值解，与先前从最低点运动矛盾。

例如取 $L=1.25\,\mathrm m$、$g=9.8\,\mathrm{m/s^2}$，则 $\sqrt{g/L}=2.8\,\mathrm{s}^{-1}$，临界角速度为 。从最低点以 推出，有 。以 推出则能越顶，经过最高点时

$$|\omega|=\sqrt{6^2-4\times7.84}=\sqrt{4.64}\approx2.154\,\mathrm{rad/s}.$$

小角度公式 $T\approx2\pi\sqrt{L/g}\approx2.244\,\mathrm s$ 只适合小幅摆动。若从 $\theta_0\in(0,\pi)$ 静止释放，能量给出精确周期的积分表达

$$T=4\sqrt{\frac{L}{2g}} \int_0^{\theta_0}\frac{d\theta}{\sqrt{\cos\theta-\cos\theta_0}}.$$

积分上端速度为零，是可积的端点奇性；当振幅趋近 $\pi$ 时，周期则趋向无穷。线性近似丢掉的振幅信息，正是在这里重新出现。

加上阻尼，稳定怎样变成渐近稳定

真实转轴会损失能量。用与角速度成正比、方向相反的阻尼转矩作近似，把除以转动惯量后的阻尼系数记为 $\mu>0$，单位为 $\mathrm{s}^{-1}$：

$$\theta'=\omega,\qquad \omega'=-\mu\omega-\frac gL\sin\theta.$$

平衡位置不变，Jacobian 变成

$$J(\theta,\omega)= \begin{pmatrix}0&1\\-(g/L)\cos\theta&-\mu\end{pmatrix}.$$

最低点的特征方程为 $\lambda^2+\mu\lambda+g/L=0$，两根实部都负。所以无论弱阻尼的螺旋，还是强阻尼的结点，最低点都局部渐近稳定；临界重根的实部也仍为负。最高点满足 $\lambda^2+\mu\lambda-g/L=0$，两根异号，仍是鞍点。

同一个能量函数现在满足

$$E'=-\mu\omega^2\leq0.$$

它把“阻尼消耗能量”写成了可逐步检查的等式。但 $E'\leq0$ 本身还不足以立即断言所有状态都收敛到某个指定最低点：转折瞬间 $\omega=0$ 也有 $E'=0$，而且精确位于最高点的解会永远停在那里。局部收敛已经由特征值得证；若还想研究某个大范围内最终落进哪一个势阱，就要结合初始能量和边界继续分析。

能量还能排除一个常见误画。若阻尼摆有非平衡周期解，经过一个周期 $T$，应有 $E(T)=E(0)$，可是

$$0=E(T)-E(0)=-\mu\int_0^T\omega(t)^2\,dt$$

迫使 $\omega$ 始终为零，矛盾。因此这个无外力的阻尼模型没有非平衡周期运动。若图上画出了永远等幅的闭圈，应先检查是否误用了无阻尼方程。

---

把一个带环境容量的模型真正做完

现在回到生态问题，但把问题说得更具体：猎物资源有限时，少量捕食者能否建立种群？若能，两种数量受到小扰动后，会持续循环还是逐渐恢复？这两个问题分别涉及边界平衡和共存平衡，比笼统地问“生态会怎样”更容易检验。

令 $x(t)$ 为猎物数量，$y(t)$ 为捕食者数量，时间以月计。保留均匀混合、线性捕食相遇项，加入猎物的环境容量限制：

$$\begin{cases} x'=rx\left(1-\frac{x}{K}\right)-axy,\\ y'=baxy-my. \end{cases}$$

$r$ 是猎物低密度净增长率，$K$ 是没有捕食者时的环境容量，$a$ 是捕食相遇系数，$m$ 是捕食者自然死亡率。$b$ 把“损失多少猎物”换成“产生多少捕食者增长”。若两物种数量单位分别记录，$b$ 的单位是“捕食者数量/猎物数量”；把它简单说成无量纲效率，容易忽略两种数量单位之间的转换。

图：先明确问题和数量单位，再用方程、平衡、局部分析与数值核验逐步回答。

先让模型在现实允许的区域里站得住

$r,m$ 的单位是 $\mathrm{月}^{-1}$，$K$ 与 $x$ 同单位，$a$ 为“每只捕食者每月”。于是 $rx(1-x/K)$ 与 都是猎物数量每月， 与 都是捕食者数量每月。初值应满足 ，参数取正。

两个方程都有各自数量作因子，因此非负象限不变。还有 $x'\leq rx(1-x/K)$，说明猎物不能越过 $\max(x_0,K)$ 向上无界增长。为了同时检查捕食者，令

$$Z=bx+y.$$

相加时捕食项正好抵消：

$$Z'=brx(1-x/K)-my =-mZ+bx(r+m-rx/K).$$

后面的二次式在 $x\geq0$ 上有最大值，故

$$Z'\leq-mZ+\frac{bK(r+m)^2}{4r}.$$

用前面的一阶线性比较思路，可知 $Z$ 有上界，进而 $x,y$ 都有上界。右端是光滑多项式，解又被限制在有界非负区域，所以非负初值的解可延续到所有 $t\geq0$。这一步不是多余的：我们现在有资格继续谈长时间行为，而不只是假设计算程序能一直往后跑。

平衡是否可行，比特征值更早一步

由 $y'=y(bax-m)=0$，先分开 $y=0$ 和 $x=m/(ba)$。边界平衡为 、，正的共存平衡为

$$x_*=\frac{m}{ba},\qquad y_*=\frac ra\left(1-\frac{m}{baK}\right).$$

它只有在 $baK>m$ 时才位于正象限。若算出 $y_*<0$，不能把它报告成真实共存数量。$baK=m$ 时共存表达与 $(K,0)$ 重合，是需要另查的临界情形。

一般 Jacobian 是

$$J(x,y)= \begin{pmatrix} r-2rx/K-ay&-ax\\ bay&bax-m \end{pmatrix}.$$

原点特征值 $r,-m$，是鞍点。在 $(K,0)$，特征值为 $-r,baK-m$。所以 $baK<m$ 时，只有猎物的平衡局部渐近稳定； 时它是鞍点。这个阈值也可以不用矩阵读出来：猎物已接近 时，少量捕食者的单位数量净增长率近似为 ，正好决定它能否在稀少时增加。

若正共存平衡存在，在该点用平衡关系化简得

$$J_*= \begin{pmatrix} -rx_*/K&-ax_*\\ bay_*&0 \end{pmatrix}.$$

其迹 $T=-rx_*/K<0$，行列式 $D=ba^2x_*y_*>0$，故共存平衡局部渐近稳定。环境容量项让迹从经典模型中的零变为负，但闭轨线是否消失不能靠“看起来有阻尼”说明；这里的局部恢复由非零负实部给出，远处轨线则要继续研究。

一组完整的数值例子

取教学参数 $r=1$、$K=200$、$a=0.02$、$b=0.25$、$m=0.4$，各自单位遵循上文约定。这组值不代表实测拟合。此时 ，共存平衡为

$$x_*=80,\qquad y_*=30,$$

而

$$J_*= \begin{pmatrix}-0.4&-1.6\\0.15&0\end{pmatrix}, \qquad \lambda=-0.2\pm i\sqrt{0.2}.$$

因此附近会以振荡方式恢复；线性近似的振幅包络按 $e^{-0.2t}$ 缩小，衰减到原来 $1/e$ 的时间尺度约为 $5$ 个月，小扰动振荡周期约为 $2\pi/\sqrt{0.2}\approx14.05$ 个月。这些时间尺度针对靠近平衡的偏差，不能要求任意大扰动都严格遵守。

从 $(x_0,y_0)=(100,20)$ 出发，先做手算检查：

$$x'(0)=100(1-100/200)-0.02\times100\times20=10,$$ $$y'(0)=0.005\times100\times20-0.4\times20=2.$$

程序画出的第一小段应当右上行进。若第一步就猎物下降，首先检查代码中的参数、符号和初值。

用前面学过的经典四阶 Runge–Kutta 法，从 $t=0$ 算到 $20$ 个月，在相同时刻比较步长 $h=0.1$ 与 $h=0.05$ 的结果：

六位小数相同是舍入后的结果，并非两次计算逐位完全一致。再把步长减半到 $0.025$，这些显示位仍一致。这支持所列时刻的数值已经较为稳定。轨线先越过共存水平，再摆回来，与局部螺旋吸引的方向一致。

核验还应包括整个采样过程的非负性、转折位置是否与零增长曲线一致，以及多组初值的表现。前面的经典捕食者模型可以检查 $H$ 是否近似守恒；当前加入环境容量的模型已经改变了方程，不能继续要求旧的 $H$ 不变。阻尼摆则应检查能量衰减，而不是守恒。检验量也必须跟着模型变。

这一次数值比较支持这组参数、这几个采样时刻的计算结果；它没有证明所有初值都收敛，也没有证明模型适用于真实种群。有限时间图像更不能排除极慢的漂移。我们已有的严格结论是可行域保持、有界性，以及各平衡附近的稳定性；如果要把吸引范围扩大成整个正象限，还需要额外的全局论证。

参数变化之后，原答案还剩多少

保持其余参数不变，把 $m$ 从 $0.4$ 提高到 $1.2$，则 $baK=1<1.2$，正共存平衡消失，$(200,0)$ 的特征值变成 $-1,-0.2$，局部渐近稳定。模型表达的是：在仅有猎物的环境容量附近，稀少捕食者的增长不足以补偿死亡。

真实应用中，$r,K$ 可以先从没有捕食者的猎物变化趋势估计，$m$ 需要与捕食者缺乏食物时的减少机制相联系，$a,b$ 则需要相互作用数据。只有一个平衡点的数量，最多给出若干参数组合的约束，通常不能把五个参数分别确定。若只记录猎物、不观察捕食者，也可能有不同参数组合产生相近曲线。此时应承认数据不足，并用参数范围研究结论会不会改变。

最终交给使用模型的人，不应只有一张漂亮相图。应当能明确说出：研究的时间范围是什么，变量怎么量，哪些机制被保留，参数如何获得，阈值是否对参数误差敏感，哪些结论已经证明，哪些只是计算支持。这样，当观察结果与模型不符时，才知道应该检查数值方法、重新估计参数，还是修改捕食饱和、季节变化等机制。

---

自己把局部判断和完整模型接起来

练习：不要漏掉边界平衡

对于 $x'=x(1-x-y)$、$y'=y(2-x-y)$，找出非负象限内所有平衡点并分类，说明少量第二种群进入 后会怎样。

显示答案

两组零增长线分别是 $x=0$ 或 $x+y=1$，$y=0$ 或 $x+y=2$。后两条平行，故没有正共存平衡，只有 。Jacobian 为

练习：纯虚根之后再走一步

分析 $x'=y$、$y'=-x+x^3$ 的三个平衡点，并证明原点附近确有闭合轨线。

显示答案

平衡点为 $(0,0),(1,0),(-1,0)$，Jacobian 为 $\begin{pmatrix}0&1\\-1+3x^2&0\end{pmatrix}$。在 处，特征值为 ，是鞍点；在原点为 ，尚不能直接判中心。

练习：周期运动的两个易混概念

对于经典捕食者模型，若正初值不在共存平衡，能否最后收敛到共存点？其闭轨线能否称为稳定极限环？

显示答案

不能收敛。正非平衡初值的 $H=h>0$，沿解保持不变；若趋向共存点，连续性又要求 $H\to0$，矛盾。共存点的稳定只保证足够小的扰动不会跑远。闭轨线也不是极限环，因为任意一条旁边都有另一能量值对应的闭轨线，缺少周期轨道的孤立性。守恒模型没有把不同能量的初值吸到同一闭轨线上的机制。

练习：阻尼会不会让倒立位置也稳定

对 $\theta''+\theta'+4\sin\theta=0$，判断最低点和最高点的局部稳定性，并计算能量变化率。

显示答案

令 $\omega=\theta'$。最低点的特征方程为 $\lambda^2+\lambda+4=0$，两根 $(-1\pm i\sqrt{15})\mathrm{/}$ 的实部均负，是局部渐近稳定的螺旋型平衡。最高点的方程为 ，两根 异号，仍是鞍点。

练习：零增长线不等于平衡点集合

在一个没有新易感者进入的简化传播模型中，$S'=-\beta SI$、$I'=\beta SI-\gamma I$，其中 、。感染者在哪些位置停止增加？这些位置是否都是平衡点？

显示答案

$I'=I(\beta S-\gamma)=0$ 给出 $I=0$ 或 $S=\gamma/\beta$。当 时，若 ，感染者增加；若 ，感染者减少。

练习：综合结果怎样被检查

在正文带环境容量的数值模型中，若程序从 $(100,20)$ 出发画出第一个时刻两种数量都减少，后来还出现 $y<0$，你会检查什么？若改正后用两个步长算出的图像几乎重合，是否就证明了全局收敛？

显示答案

手算已得 $x'(0)=10,y'(0)=2$，所以足够小的第一步应向右上。先核对 $0.02xy$、$$ 的符号、参数对应及初值，再检查 RK4 各阶段是否使用相应中间状态。精确解保持非负，出现负数要检查步长和方法是否适合，不能直接解释成负种群，也不应只用截断到零掩盖误差。

我们最初解 $y'=2y$ 时，第一次把“未知数”换成了整个函数。走到这里，未知的已经是一组互相牵动的函数，而我们也不再只等着一个显式公式。一个平衡位置、一条不能越过的边界、一个守恒量、一个经核验的数值轨迹，都可能回答模型里真正的问题。

下次遇到一个变化过程，可以从一张纸开始：先写清楚你在追踪什么、每单位时间为什么增加或减少，再给它初值和单位；接着寻找能算出的结构，并检查剩下的计算是否可信。最后把结论翻译回最初的现象，同时带上它成立的条件。能把这个过程独立走通，微分方程才真正从一道道习题，变成了你可以反复使用的建模方法。