捕食者-猎物模型

捕食者-猎物模型（predator–prey model）描述两种群因捕食关系而相互影响的数量变化。经典的洛特卡—沃尔泰拉模型（Lotka–Volterra model）用一对非线性微分方程表达这种关系：猎物为捕食者提供食物，捕食者又减少猎物数量。它给出了种群持续振荡的一种数学机制。

猎物增加以后，捕食者会怎样变化

设 $x(t)$ 是猎物数量，$y(t)$ 是捕食者数量，$t$ 以月计。把数量视为连续量，相当于对个体的出生、死亡和捕食作平均描述。以下选取的参数、曲线与初值均为合成教学模型，不是某个生态系统的实测结果。

先把相互作用暂时拿掉。假设没有捕食者时，猎物的每个个体具有固定净增长率 $a>0$，因此 $x^{\prime }=ax$；假设没有猎物时，捕食者以固定人均死亡率 $c>0$ 减少，因此 $y^{\prime }=-cy$。

现在加入捕食。若双方充分混合，单个捕食者遇到猎物的频率与猎物数量成正比，所有捕食者合计的捕获量就与乘积 $xy$ 成正比。于是猎物损失率写成 $bxy$，其中 $b>0$。再假定由捕食支持的捕食者增长与捕获量成正比，记相应系数为 $d>0$，得到 $$\begin{array}{cc}\frac{dx}{dt}& =ax-bxy=x(a-by),\\ \frac{dy}{dt}& =dxy-cy=y(dx-c).\end{array}$$ 若两种群都以个体数计，$a,c$ 的单位为每月，$b$ 为每捕食者每月，$d$ 为每猎物每月。这样 $by$ 和 $dx$ 都是人均变化率，每一行右边的单位与左边一致。两种群的个体尺度不同，$b$ 与 $d$ 无须相同；比值 $d/b$ 表示模型中的转化关系，而不是“一只猎物必定变成一只捕食者”。

这组方程还假定参数不随季节改变，没有迁入迁出，猎物的资源没有在方程中设置上限，捕食也没有因捕食者吃饱而饱和。它描述的是这些假设共同形成的系统。Jeffrey R. Chasnov，《Mathematical Biology》§1.4，印刷页7–13从种群增长、损失和接触项建立了同类方程。

先算出一个瞬时状态

取一组具体参数：$a=c=0.5$ 每月，$b=0.01$ 每捕食者每月，$d=0.005$ 每猎物每月。初始时有200只猎物、50只捕食者，即 $x(0)=200,y(0)=50$。代入得到 $$x^{\prime }(0)=0.5\times 200-0.01\times 200\times 50=0,$$ $$y^{\prime }(0)=0.005\times 200\times 50-0.5\times 50=25.$$ 猎物当前净增长为零，捕食者却正以每月25只的瞬时速率增加。捕食者一旦超过50，猎物的人均净增长率 $0.5-0.01y$ 就变成负数，于是猎物开始下降。

这里的25是此刻的导数，不能直接断言一个月后捕食者恰好增加25只，因为随后猎物和捕食者都在变。要得到一段时间后的数量，需要把相互影响连续积累起来。

方程也保留了非负性。例如在正解存在期间，第一式可以写成 $$x(t)=x(0)\exp ({\int }_0^t[a-by(s)]ds)>0.$$ 第二式有相同形式。若某种群初始为零，相应坐标轴上的解始终为零；模型不包含从零凭空产生该种群的机制。正初值能否长期存在且保持有界，将由后面的守恒量说明。

用零增长线读出四种变化方向

平衡要求两个变化率同时为零。原点 $(0,0)$ 是一个平衡；对正种群，必须同时满足 $a-by=0$ 和 $dx-c=0$，因此共存平衡为 $$x_{\ast }=\frac{c}{d},y_{\ast }=\frac{a}{b}.$$ 本例得到 $x_{\ast }=100,y_{\ast }=50$。为便于比较，令 $u=x/100$、$v=y/50$，即把各自平衡数量当作单位；再令无量纲时间 $\tau =0.5t$。例如 $u=2$ 表示200只猎物，$v=2$ 表示100只捕食者。

由链式法则 $dx/dt=50du/d\tau$，代入第一式；第二式同理，便得到 $$\frac{du}{d\tau }=u(1-v),\frac{dv}{d\tau }=v(u-1),(u(0),v(0))=(2,1).$$ 以下图中的撇号均指对 $\tau$ 求导。在正象限内，$u>0,v>0$，所以 $u^{\prime }$ 的正负只由 $1-v$ 决定，$v^{\prime }$ 的正负只由 $u-1$ 决定。

水平线 $v=1$ 是猎物的零增长线，竖直线 $u=1$ 是捕食者的零增长线。站在一条零增长线上，只有对应的一种群瞬时不变；两线交点 $(1,1)$ 才是双方同时不变的状态。坐标轴也是原方程的零增长集合的一部分，下图着重画正象限内的两条线。

箭头表示相图中的变化方向，长度经过统一处理，不代表实际变化速率大小。

从右侧的初值 $(2,1)$ 出发，轨迹先向上进入右上区，猎物下降而捕食者增加；越过 $u=1$ 后，捕食者也开始下降；再越过 $v=1$，猎物转而增加。沿四个区域追踪，得到逆时针转动的趋势。不过，方向图尚未说明它究竟闭合、向内收缩，还是向外展开。

沿轨迹保持不变的量

要分清上述三种可能，需要寻找约束轨迹的关系。观察方程，猎物变化含有 $1-v$，捕食者变化含有 $u-1$。选择 $$H(u,v)=u-\log u+v-\log v-2,u>0,v>0,$$ 它的两个偏导数分别是 $1-1/u$ 和 $1-1/v$，与方程里的 $u,v$ 因子相乘时正好约去分母。利用链式法则逐步计算： $$\begin{array}{cc}\frac{dH}{d\tau }& =(1-\frac{1}{u})u^{\prime }+(1-\frac{1}{v})v^{\prime }\\ & =(u-1)(1-v)+(v-1)(u-1)\\ & =(u-1)[(1-v)+(v-1)]=0.\end{array}$$ 因此 $H$ 沿每一条正轨迹保持常数，称为第一积分或守恒量。它是模型中某个组合保持不变，并不是两种群总数量不变；事实上 $(u+v{)}^{\prime }=u-v$，一般不等于零。

若想看这个表达式怎样被找到，在 $u^{\prime }\ne 0$ 的一段轨迹上，可以消去时间： $$\frac{dv}{du}=\frac{v(u-1)}{u(1-v)}.$$ 分离变量为 $(1/v-1)dv=(1-1/u)du$。分别积分，得到 $\log v-v=u-\log u+C$，整理就是 $H$ 为常数。消去时间时暂时排除了 $v=1$ 上的点，而刚才直接求导的验证在整个正象限都成立，补上了这些点。

初值 $(2,1)$ 确定 $$H_0=2-\log 2+1-0-2=1-\log 2\approx \mathrm{0.3068528194.}$$ 以后轨迹只能在这条等值线上运动。相同守恒结构也见佛罗里达大学 PHZ4710 课程，Closed orbits and oscillations；加州大学戴维斯分校 Math 207A，2014年期末解答第2题用消去时间的方法分析了这一系统。

为什么确实是闭轨道，而不是螺旋

令 $\varphi (s)=s-\log s-1$。有 ${\varphi }^{\prime }(s)=1-1/s$：在 $0<s<1$ 上为负，在 $s>1$ 上为正。因此 $\varphi$ 在1取得唯一最小值0，而且当 $s\to 0^{+}$ 或 $s\to \mathrm{\infty }$ 时都趋向无穷。守恒量正是 $H=\varphi (u)+\varphi (v)$。

由此先得到有界性：固定有限的 $H_0$ 时，两个坐标都不能趋向零或无穷，否则某个非负项 $\varphi$ 将超过 $H_0$。轨迹被限制在正象限内部的一个闭有界区域里，既不碰到坐标轴，也不在有限时间发散，所以解能继续延伸。

再看等值线的形状。${\varphi }^{″}(s)=1/s^2>0$，因此从 $(1,1)$ 沿任何直线方向向外走，$H$ 从0严格增加，直到趋向无穷。每条射线与 $H=H_0>0$ 恰交一次，组成包围中心的闭合曲线。除中心外，$H$ 的梯度不为零，等值线光滑；向量场也没有其他正平衡点，所以沿这条曲线的速度处处非零。闭合曲线是紧的，速度大小有正的下界，绕行一周所需时间有限，因而轨迹周期性回到出发状态。

实线是初值(2,1)的轨道；两条虚线分别从(1.2,1)和(1.5,1)出发。不同初值确定不同的守恒值。

这也解释了稳定性的准确含义。共存平衡 $(1,1)$ 是 $H$ 的严格最小点，足够小的守恒值把轨迹限制在足够小的邻域内，因此它是稳定的。更具体地，围绕中心取一条小圆，其上 $H$ 有正的最小值；若初值的 $H$ 更小，就不可能越过这条圆。

但它不是渐近稳定的。非平衡初值的 $H_0>0$ 始终不变，若轨迹趋向中心，就应有 $H\to 0$，与守恒矛盾。一次扰动通常会把系统送到另一条闭轨道，而不是使它回到原来的轨道。这里存在一整族闭轨道，也就不是一条孤立并吸引邻近轨迹的极限环。常说的“中性振荡”指的正是这种不向中心衰减的行为。

原点的情况不同。没有捕食者而引入少量猎物时，$x^{\prime }=ax$ 会让它远离原点；从线性化看，原点的两个特征值为 $a$ 与 $-c$，在整个平面中是鞍点，在生物学使用的非负区域内也不稳定。

同一条轨道上的数量极值与时间

主轨道上的A点为 $(2,1)$。这里 $u^{\prime }=0$，且轨迹从 $v<1$ 一侧进入 $v>1$ 一侧，猎物由增加转为减少，所以A是猎物最大值。捕食者的最大值发生在B点：轨迹从 $u>1$ 进入 $u<1$，故 $v^{\prime }$ 从正变负。

在B处 $u=1$，守恒式化为 $$v-\log v-1=1-\log 2.$$ 其中大于1的解是 $v=2$，所以B恰为 $(1,2)$。小于1的另一个解记为 $q$，满足 $$q-\log q=2-\log 2,q\approx \mathrm{0.4063757400.}$$ 由 $\varphi$ 在0到1严格递减、从无穷降到0，可知小根唯一，可以用二分法求出。例如代入0.4063与0.4065，会把目标值夹在两边，再逐次缩小区间。由于本例的 $H$ 对 $u,v$ 对称，C为 $(q,1)$，D为 $(1,q)$。

换回原数量，猎物在约40.638与200之间振荡，捕食者在约20.319与100之间振荡。它们的最大值都等于各自平衡数量的两倍，但并不同时发生；本例这种相同倍数来自特定参数和初值，不是所有捕食模型都具有的关系。

相图上的一圈，展开成时间轴上的一周期。纵轴为各自平衡数量的倍数，不能直接把两条线的高度当作原始个体数比较。

对该初值进行数值积分，可得以下事件。表中的 $\tau$ 是无量纲时间，实际月份为 $t=2\tau$。

事件	$\tau$（约）	$(u,v)$	含义
A	0	$(2,1)$	猎物最大，捕食者仍在增加
B	1.146456	$(1,2)$	捕食者最大，猎物仍在减少
C	2.763557	$(q,1)$	猎物最小，捕食者仍在减少
D	4.991381	$(1,q)$	捕食者最小，猎物仍在增加
返回A	6.608483	$(2,1)$	完成一周期

所以周期约为13.216965个月，捕食者高峰比猎物高峰晚约2.292912个月。四个阶段用时不相等，不能因为相图大致围成一圈，就把它们当作等长的四分之一周期。

在平衡很近的地方，令 $u=1+\xi ,v=1+\eta$。代入得 ${\xi }^{\prime }=-\eta -\xi \eta$、${\eta }^{\prime }=\xi +\xi \eta$。忽略二次小量 $\xi \eta$ 后，${\xi }^{\prime }=-\eta$、${\eta }^{\prime }=\xi$，再求导得到 ${\xi }^{″}=-\xi$。这个线性近似的无量纲周期是 $2\pi$，在本例为约12.566个月。它与幅度较大的主轨道周期13.217个月不同。一般参数的平衡附近小振幅周期为 $2\pi /\sqrt{ac}$，不能把这个线性化结果当作所有振幅的精确周期；闭轨道的存在由前面的非线性守恒论证保证。

数值算法为什么可能画出向外螺旋

最直接的时间推进是显式欧拉法。选无量纲步长 $h>0$，在同一个旧状态 $(u_n,v_n)$ 上计算两个变化率： $$\begin{array}{cc}{u}_{n+1}& =u_n+h{u}_n(1-v_n),\\ v_{n+1}& =v_n+h{v}_n(u_n-1).\end{array}$$ 例如 $h=0.1$，从 $(2,1)$ 出发，第一步为 $(2,1.1)$；第二步为 $$u_2=2+0.1\times 2(1-1.1)=1.98,v_2=1.1+0.1\times 1.1(2-1)=\mathrm{1.21.}$$ 这里两个更新都使用第一步的旧值。若先算出新的 $u$ 再带入同一步的 $v$，就已经换了算法。

第一步欧拉值虽然接近真解，却已改变守恒量： $$H(2,1.1)-H(2,1)=0.1-\log 1.1\approx 0.00468982>0.$$ 这种偏离会逐步积累。事实上，令 $A=h(1-v_n)$、$B=h(u_n-1)$，则 $u_{n+1}=u_n(1+A)$、$v_{n+1}=v_n(1+B)$。在下一步两个坐标仍为正时，直接相减得 $$\begin{array}{cc}{H}_{n+1}-H_n& =u_{n}A+v_{n}B-\log (1+A)-\log (1+B)\\ & =A+B-\log (1+A)-\log (1+B).\end{array}$$ 第二行用了 $u_{n}A+v_{n}B=h(u_n-v_n)=A+B$。对于 $s>-1$，$\log (1+s)\le s$，等号只在 $s=0$ 成立。因此正象限里的每个非平衡欧拉步都会增大 $H$，把数值状态推向更外侧的等值线。向外螺旋来自算法的误差，不是原连续模型的种群振幅逐渐增大。

上图积分至无量纲时间20；下图用同一初值比较h为0.1、0.05、0.025。理论守恒量误差应为零。

在 $\tau =20$，这三个步长的 $H-H_0$ 分别约为1.605442、0.467689和0.188184。这里展示的是固定算法在有限区间里的计算结果，减半步长并未使它自动成为保持守恒量的方法。步长更大时，还可能因 $1+h(1-v_n)\le 0$ 或 $1+h(u_n-1)\le 0$ 而生成零或负数量。

平衡附近还有一个更简短的解释。线性化的欧拉更新为 ${\xi }_{n+1}={\xi }_n-h{\eta }_n$、${\eta }_{n+1}={\eta }_n+h{\xi }_n$。平方相加，中间的交叉项抵消，得到 $${\xi }_{n+1}^2+{\eta }_{n+1}^2=(1+h^2)({\xi }_n^2+{\eta }_n^2).$$ 任何正步长都令这个线性近似的振幅逐步放大。它与“步长趋于零时，在固定有限时间内可以收敛到真解”并不矛盾：固定步长看无限长时间，与固定时间缩小步长，是不同的极限。相应数值稳定性分析见Driscoll 与 Braun，§11.3，例11.3.7。

本文用于主轨道和时间图的参考计算采用高阶自适应积分，并另用固定小步长四阶Runge–Kutta法对照。检查包括减小步长、守恒量偏差和返回时间的一致性。守恒值接近只能说明轨道层面误差较小，仍需检查沿轨道的时间位置；一个计算点可以落在正确等值线上，却提前或滞后到达。

参数改变与真实生态系统

前面的归一化不限于这组参数。对任意 $a,b,c,d>0$，取 $u=x/(c/d)$、$v=y/(a/b)$、$\tau =at$，得到 $$u^{\prime }=u(1-v),v^{\prime }=\rho v(u-1),\rho =\frac{c}{a}>0.$$ 此时守恒量变为 $H_{\rho }=\rho \varphi (u)+\varphi (v)$。求导后仍有 $${H_{\rho }}^{\prime }=\rho (u-1)(1-v)+\rho (v-1)(u-1)=0.$$ 两个权重都正，所以前面的正性、有界性和闭轨道论证仍成立；轨道形状、不同阶段的用时则会随参数比和初值改变。

若观测中出现振幅衰减、长期漂移或灭绝，就需要检查模型假设。猎物的资源约束可以通过洛吉斯蒂模型中的 $ax(1-x/K)$ 引入；捕食饱和、季节变化、年龄结构和随机个体事件也会改变方程。添加这些机制后，原有守恒量通常不再守恒，因而不应把经典模型的闭轨道结论直接套用。

经典方程的正解不会在有限时间到达零，却不代表很小的真实种群能够避免随机灭绝。拟合到两条起伏曲线，也不足以证明起伏只由捕食关系造成：食物、气候、迁移和观测方式都可能影响数据。相图首先帮助解释给定方程的机制，生态解释还需要独立观测与参数检验。

历史与进一步阅读

阿尔弗雷德·洛特卡在1920年的论文《Analytical Note on Certain Rhythmic Relations in Organic Systems》中研究了生物系统的持续振荡，并在1925年《Elements of Physical Biology》中扩展相关分析。维托·沃尔泰拉独立研究了捕食关系，于1926年发表《Fluctuations in the Abundance of a Species Considered Mathematically》。这些具体发表脉络可见科学史学者 Sharon Kingsland 在 PNAS 发表的研究回顾及洛特卡1920年论文原文。

• Jeffrey R. Chasnov：Mathematical Biology，§1.4，模型建立、平衡和无量纲化。
• BingKan Xue：PHZ4710: Lotka-Volterra System，Closed orbits and oscillations、Limit Cycle两节，守恒量与修改模型后的差别。
• John K. Hunter：Math 207A, Fall 2014, Final Solutions，第2题，印刷页2–3，相轨道的隐式方程。
• Toby A. Driscoll、Richard J. Braun：Fundamentals of Numerical Computation，§11.3，振荡系统的欧拉法稳定性。
• Alfred J. Lotka：上引1920年原文，PNAS 6(7):410–415，DOI：10.1073/pnas.6.7.410。
• Sharon Kingsland：Alfred J. Lotka and the origins of theoretical population ecology，PNAS 112(31):9493–9495（2015），DOI：10.1073/pnas.1512317112。