常微分方程（十五）：种群动力学

为什么猞猁和雪兔的数量会以诡异的规律每 10 年循环一次？ 为什么引入一个新物种有时会导致整个生态系统崩溃？为什么相似的竞争者有时能共存，有时却会相互灭绝？答案并不在物种本身，而在于描述这些物种之间关系的方程。本章将带你梳理数学生态学中的经典模型：从单一种群的 Logistic 模型和 Allee 效应，到多物种间的竞争、捕食-被捕食振荡，再到年龄结构与空间扩散。

本章要点#

• 单种群增长的三大模型：Malthus（指数增长）、Logistic（饱和增长）和 Allee 效应（存在灭绝阈值）
• Lotka-Volterra 捕食-被捕食模型：闭合周期轨道、守恒量，以及“富营养化悖论”
• Holling 功能反应：如何打破严格的周期性
• 两种群竞争模型的四种典型结局：稳定共存、单方胜出（两种情形）和双稳态（奠基者效应）
• 年龄结构模型：Leslie 矩阵、主特征值作为长期增长率，以及稳定年龄分布
• 集合种群模型（Levins）：斑块占用动态与区域灭绝阈值
• 空间扩散模型：Fisher-KPP 行波，其最小传播速度为 $c_{\min} = 2\sqrt{Dr}$

先修知识：第 7 章 的相平面分析、第 8 章 的非线性稳定性理论，以及第 13 章 的偏微分方程引论（用于空间扩散部分）。

单种群增长#

设 $N(t)$ 表示种群数量。

Malthus 模型（1798 年）：$\dot N = r N$ ，其解为 $N(t) = N_0 e^{rt}$ 。该模型数学形式简洁，但在生物学上仅适用于少数几代——长期来看并不现实。

种群增长速率在 $N = K/2$ 处达到最大值 $rK/4$ 。

此时，$N=0$ 和 $N=K$ 都是稳定平衡点，而 $N=A$ 是一个不稳定的阈值。若初始种群低于 $A$ ，则走向灭绝；若高于 $A$ ，则趋向承载力 $K$ 。这种现象称为双稳态——同一个方程拥有两个吸引域。

Allee 效应解释了为何一些濒危物种即使在栖息地恢复后仍难以复苏，也说明了为何小规模的野化放归常常失败。其数学形式与力学中的三次势函数完全相同，动力学行为相当于在双势阱景观中进行梯度下降。

左上图：人均增长率 $\dot N / N$ 随种群密度的变化。Logistic 模型（蓝色）在 $N < K$ 时始终为正；强 Allee 模型（红色）在 $N < A$ 时为负。右上图：$\dot N$ 与 $N$ 的关系——强 Allee 模型存在三个平衡点。左下图：不同初始种群的演化轨迹——起始值低于 $A$ 的种群最终灭绝，高于 $A$ 的则增长至 $K$ 。右下图：对应的“势能”函数 $V(N)$ ，满足 $\dot N = -V'(N)$ 。两个稳定状态（$N=0$ 和 $N=K$ ）对应势阱，而阈值 $A$ 则是分隔两者的势垒。若系统受到随机扰动，就可能在两个吸引域之间“隧穿”，从而为生态学中的“状态跃迁”提供清晰的数学模型。

捕食-被捕食模型：Lotka-Volterra#

• $\alpha$ ：无捕食者时被捕食者的内禀增长率
• $\gamma$ ：无食物时捕食者的死亡率
• $\beta, \delta$ ：相遇与转化效率参数

令导数为零，可得两个平衡点：$(0, 0)$ （灭绝状态，为鞍点）和共存平衡点 $(\gamma/\delta,\ \alpha/\beta)$ （一个中性中心）。

一个守恒量#

直接计算可得沿任意解轨有 $\dot H = 0$ 。因此，所有 Lotka-Volterra 轨道都位于 $H$ 的某个等值线上。这些等值线是围绕共存点的闭合曲线，这正是解呈现周期性的原因——该中心确实是中性的，而非渐近稳定的。

但这种守恒性极为脆弱：对方程稍作扰动就会破坏它。真实生态系统并不严格满足 Lotka-Volterra 假设，因此自然界中的种群振荡通常是极限环（如下文 Holling 响应所示），而非 LV 模型中连续的一族闭轨。

雪兔-猞猁数据#

哈德逊湾公司 1845–1930 年的毛皮收购记录清晰地展示了雪兔与加拿大猞猁之间约 10 年的周期性波动。Lotka-Volterra 模型在定性上能很好地复现这一现象。不过，现代生态学家普遍认为，这一周期实际上主要由雪兔与植被的相互作用驱动，猞猁只是被动跟随。尽管如此，历史上模型与数据的高度吻合仍使 Lotka-Volterra 名声大噪。

左上图：经典的异相振荡——被捕食者先达到峰值，捕食者滞后约四分之一周期。右上图：相平面上五条不同初值的闭合轨道，均围绕共存点 $(\gamma/\delta,\ \alpha/\beta)$ ，向量场指示旋转方向。左下图：沿一条轨道计算的 $H$ 值——仅有数值漂移，验证了其守恒性。右下图：当采用 Holling II 型功能响应并加入 Logistic 被捕食者增长时，提高承载力 $K$ 会导致平衡点失稳，产生一个不断扩大的极限环——这就是著名的“富营养化悖论”。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

def lv(t, s, a=1.0, b=0.1, d=0.05, g=0.5):
    x, y = s
    return [a*x - b*x*y, d*x*y - g*y]

sol = solve_ivp(lv, (0, 60), [10, 2], t_eval=np.linspace(0, 60, 4000))
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].plot(sol.y[0], sol.y[1]); ax[0].set_xlabel('被食者'); ax[0].set_ylabel('捕食者')
ax[1].plot(sol.t, sol.y[0], label='被食者')
ax[1].plot(sol.t, sol.y[1], label='捕食者'); ax[1].legend()
plt.tight_layout(); plt.show()

功能反应（Holling）#

现实中，捕食者存在摄食饱和现象——一只狼不可能一天吃掉无限多只兔子。Holling 将捕食者的功能反应分为三类：

• I 型：线性响应 $g(x) = a x$ （即 Lotka-Volterra 默认形式，生物学上罕见）
• II 型：$g(x) = \frac{a x}{1 + a h x}$ ，摄食率趋于饱和值 $1/h$ ，最为常见
• III 型：$g(x) = \frac{a x^2}{1 + a h x^2}$ ，呈 S 形，在低密度时允许捕食者切换猎物

将 Lotka-Volterra 中的线性相互作用替换为 Holling-II 型响应，便得到 Rosenzweig-MacArthur 模型。该模型的共存平衡点会随着被捕食者承载力 $K$ 的增大，通过 Hopf 分岔失去稳定性。更高的 $K$ 意味着更多食物 → 更多捕食者 → 更剧烈的被捕食者崩溃 → 更大的极限环振幅 → 最终可能导致灭绝。换句话说，食物越多反而越危险。

富营养化悖论#

在固定其他参数的情况下，Hopf 分岔发生在 $K$ 超过某一临界值 $K_c$ 之时，此后出现的极限环振幅大致按 $\sqrt{K - K_c}$ 增长。当振幅足够大时，种群数量会逼近坐标轴，此时即使微小的随机扰动也可能导致某一物种灭绝。这并非纯粹的理论反直觉——湖泊富营养化实验中已观察到类似现象。

两种群竞争#

无量纲的竞争系数 $\alpha_{ij}$ 衡量的是：一个物种 $j$ 的个体对物种 $i$ 的人均增长率的抑制程度，相对于一个同种个体的抑制效果。

该系统最多存在四个平衡点：

• $(0, 0)$ —— 当两物种均可独立增长时，此点总是不稳定的
• $(K_1, 0)$ —— 物种 1 胜出
• $(0, K_2)$ —— 物种 2 胜出
• 内部平衡点 $(\hat N_1, \hat N_2)$ （若存在）—— 共存

这些平衡点的稳定性取决于：当另一物种处于其承载力时，当前物种是否能够成功入侵。由此可推导出四种典型结局：

竞争排斥原理（Gause，1934）是上述结论的一个推论：若两个物种的生态位完全相同（即 $\alpha_{12} = \alpha_{21} = 1$ 且 $K_1 = K_2$ ），则无法稳定共存。共存的前提是生态位分化——种内竞争必须强于种间竞争。

每幅图均包含零增长线（蓝色/红色）、向量场（灰色箭头）、五条不同初值的轨迹（彩色），以及平衡点（黑点或星号）。左上：共存情形——零增长线以正确方式相交，内部星号为稳定点。右上与左下：排斥情形——内部平衡点不存在，所有轨迹最终汇入某一条坐标轴。右下：双稳态情形——内部平衡点为鞍点*，其稳定流形构成吸引域的分界线；初始条件落在哪一侧，哪一方就获胜。*

年龄结构模型#

• $F_i$ ：年龄为 $i$ 的个体在单位时间内的生育率（后代数量）
• $P_i$ ：从年龄 $i$ 存活至 $i+1$ 的概率

该模型本质上是线性迭代 $n_{t+1} = L\,n_t$ ，其长期行为由 Leslie 矩阵 $L$ 的主特征值 $\lambda$ 决定：

• $\lambda > 1$ ：种群以几何速率 $\lambda$ 增长
• $\lambda = 1$ ：种群规模稳定
• $\lambda < 1$ ：种群逐渐衰退至灭绝

对应的右特征向量即为稳定年龄分布——即系统达到动态平衡后，各年龄组所占的比例。无论初始年龄结构如何，种群最终都会收敛到这一分布（对于非负不可约的 Leslie 矩阵，Perron-Frobenius 定理保证主特征值为实数且唯一主导，排除了持续振荡的可能性）。

为何 $\lambda$ 是关键指标#

净生殖率定义为 $R_0 = \sum_i \ell_i F_i$ ，其中 $\ell_i = \prod_{j < i} P_j$ 是存活至年龄 $i$ 的概率。可以证明 $R_0 = 1$ 当且仅当 $\lambda = 1$ 。但种群增长不仅取决于生育总数，还取决于生育时机——越早出生的后代对增长的贡献越大，因为它们能更早开始繁殖。这种“时间价值”效应已被特征值 $\lambda$ 自动纳入考量。

左上：Leslie 矩阵 $L$ 的热力图——首行表示生育率 $F_i$ ，次对角线表示存活率 $P_i$ 。右上：特征值在复平面上的分布；主特征值（红星）位于单位圆外，表明种群增长。左下：从 1000 名新生儿出发的种群轨迹——约 10–15 年后进入稳定的几何增长阶段。右下：稳定年龄分布（红色柱状图）与 $t=30$ 时刻的实际年龄结构（黑点）高度吻合。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

import numpy as np

F = np.array([0.0, 0.5, 2.5, 1.5, 0.2])           # 生育率
P = np.array([0.6, 0.7, 0.85, 0.5])               # 存活率
L = np.zeros((5, 5)); L[0, :] = F
for i in range(4): L[i + 1, i] = P[i]

lam, vec = np.linalg.eig(L)
i = np.argmax(np.abs(lam))
print('长期增长率:', lam[i].real)
print('稳定年龄分布:', np.abs(vec[:, i].real / vec[:, i].real.sum()))

集合种群：多斑块系统#

其中 $p$ 为被占据斑块的比例，$c$ 为定殖率（每个被占斑块向空斑块扩散的速率），$e$ 为局部灭绝率。系统平衡点为 $p^* = 1 - e/c$ ，仅当 $c > e$ 时为正。

该模型在结构上与 SIS 传染病模型完全一致：斑块相当于“个体”，定殖过程相当于“传播”。集合种群得以维持的条件是 $c > e$ ，这构成了一个区域灭绝阈值。

这对保护生物学的启示是：即使每个栖息地斑块本身健康完好，若斑块间的连通性不足（即定殖率 $c$ 过低），仍可能导致整个区域的物种灭绝。道路和围栏虽未破坏单个斑块，却通过切断扩散廊道降低了 $c$ ，从而加剧了灭绝风险。

空间扩散：Fisher-KPP 模型#

其证明仅需两步：在前沿区域（$N \ll K$ ）对模型线性化；代入行波试探解 $N = e^{-\lambda(x - ct)}$ ，可得色散关系 $c = D\lambda + r/\lambda$ ，对其关于 $\lambda$ 求最小值得到 $c_{\min} = 2\sqrt{Dr}$ 。而实际被选择的波速恰好就是这个线性最小值——这便是著名的 KPP 选择原理。

该公式在生态学中应用广泛：可用于预测气候变化下植物分布区的扩张速度、解释 1920 年代麝鼠入侵欧洲的实测波速（约为 $\sqrt{Dr}$ ），以及描述遗传学中有利等位基因在种群中固定的传播速度。

左上：不同灭绝率 $e$ 下 Levins 模型的 $p(t)$ 曲线。$e$ 越大，平衡占比 $p^*$ 越低；当 $e > c$ 时，集合种群发生区域性灭绝。右上：$p^*$ 随 $e/c$ 的变化——在 $e = c$ 处存在明显的灭绝阈值。左下：Fisher-KPP 行波在 $t = 0, 20, 40, 60, 80$ 时刻的快照，波前保持恒定形状，以速度 $c_{\min} = 2\sqrt{Dr}$ 匀速推进。右下：数值追踪的波前位置；其线性拟合斜率与理论值 $2\sqrt{Dr}$ 的偏差小于 2%。

超越 Fisher：图灵斑图#

当两个物种（如激活子与抑制子，或捕食者与被捕食者）以不同速率扩散时，原本均匀的空间分布可能因微小扰动而失稳，自发形成空间斑图。这就是 Turing 不稳定性（1952）。该机制被假定为动物皮毛花纹、半干旱地区植被条带，以及某些贻贝床空间格局的成因，也是数学形态发生学的核心支柱。

总结#

数学生态学的魅力在于：方程数量虽少，却能涌现出极其丰富的行为。仅靠几条组合规则——扩散加反应、单种群加年龄结构、均质假设加异质性——便足以构建整个领域的理论大厦。

练习题#

概念题

1. 从强 Allee 方程出发，显式推导势函数 $V(N)$ ，并指出两个吸引域及势垒的位置。
2. 为何 Lotka-Volterra 模型中的中心是中性的，而非渐近稳定的？请利用守恒量 $H(x,y)$ 进行论证。
3. 两个竞争物种满足 $\alpha_{12} = 0.9$ 、$\alpha_{21} = 0.9$ ，$K_1 = K_2 = 100$ ，$r_1 = r_2 = 0.5$ 。计算其内部平衡点，并判断其稳定性。

计算题

4. 通过精细网格数值模拟 Fisher-KPP 行波，追踪 50% 密度等值面，验证理论波速 $c_{\min} = 2\sqrt{Dr}$ 。
5. 在 Rosenzweig-MacArthur 模型中改变 $K$ ，数值求解 Hopf 分岔点，并与解析预测结果对比。
6. 构造一个 Leslie 矩阵，其生育率 $F = (0, 0, 1.2, 1.0, 0.4)$ ，存活率 $P = (0.7, 0.85, 0.85, 0.5)$ 。求其长期增长率和稳定年龄分布。

编程题

7. 制作双稳态竞争系统的相平面动画：让初始条件连续跨越鞍点的稳定流形，观察吸引域的切换过程。
8. 在 1000 个斑块上实现随机 Levins 模型，比较其平均占用率与确定性解 $p^* = 1 - e/c$ ，特别关注 $e/c \to 1$ 时的行为。
9. 编写二维 Fisher-KPP 模拟器，从局部初始条件出发观察波的传播，并与一维理论波速进行比较。
10. 实现一个生育率 $F_i(t)$ 随时间周期性变化（模拟季节效应）的 Leslie 模型，观察种群年龄结构如何响应外部强迫。

下一步#

下一章进入控制理论——一个把 ODE 当工具去主动设计系统行为的领域。PID 控制、状态反馈、最优控制、卡尔曼滤波——它们都是在问’怎么调节输入让 ODE 产生我想要的输出’。这一视角把 ODE 从描述工具升级成设计工具。

参考文献#

• Murray, Mathematical Biology I & II, Springer (2002, 2003)
• Edelstein-Keshet, Mathematical Models in Biology, SIAM Classic (2005)
• Kot, Elements of Mathematical Ecology, Cambridge (2001)
• Hastings, Population Biology: Concepts and Models, Springer (1997)
• Lotka, Elements of Physical Biology, Williams & Wilkins (1925)
• Volterra, “Variazioni e fluttuazioni del numero d’individui in specie animali conviventi,” Mem. R. Acad. Naz. dei Lincei (1926)
• Fisher, “The wave of advance of advantageous genes,” Ann. Eugenics 7 (1937)
• Hanski, Metapopulation Ecology, Oxford (1999)
• Turing, “The chemical basis of morphogenesis,” Phil. Trans. Roy. Soc. B 237 (1952)