第五讲 微分方程模型 —— 种群模型

数学建模 讲义 NUDT 2023SP

---

5.1 微分方程模型

• 微分方程（Differential Equation）是一种数学方程，用来描述某一类函数与其导数之间的关系.
• 如果微分方程（组）中的自变量是时间（一般记为 $t$），则对应的系统常被称为 时变系统 或 动态（力）系统 (Dynamical System).
  
  • 如果时变系统的模型中不显示地出现自变量 $t$，这样的系统常常也被称为 定常系统 或 自治系统 (Autonomous System).
• 除了少数形式简单或特殊的微分方程（组），大多数的微分方程（组）都无法得出精确的解析解，对其解的计算和分析常常只能诉诸数值方法或其他的量化分析工具（如微分方程定性理论）.

---

---

种群模型 （Population Ecology Model）

• 种群（生态）模型一般特指用于研究特定物种的个体总数随时间变化的模型.
• 依据研究背景假设的不同，较为经典的种群模型有
  
  • Malthus 模型：单一种群在资源无限情况下的变化
  • Logistic 模型：单一种群在资源有限情况下的变化
  • 可开发的单种群模型：引入固定比例（或数量）的种群损失（开发利用）
  • 多种群模型：多个种群间的相互作用（合作、竞争、弱肉强食）

---

5.2 Malthus 模型

• 得名于 Thomas Robert Malthus，他于 1798 年写下了最早的人口学著作之一 《人口论》（An Essay on the Principle of Population）.
• 预言："人口增长超越食物供应，会导致人均占有食物的减少."—— Malthus 危机（或 人口陷阱）

---

• 设 $p(t)$ 为给定的物种在时刻 $t$ 的种群规模.
• $r(t,p)$ 为该物种在时刻 $t$ 出生率与死亡率之差，称为自然增长率.
• Malthus 模型：假设 $r(t,p)=r$ 为常数，则种群的增长规律可表示为
  
  • $\bbox[yellow]{\dfrac{\mathrm{d} p}{\mathrm{d} t}=r\cdot p}$
• 若设初值为 $p(t_0)=p_0$，则以上方程的解为 $p(t)=p_0\mathrm{e}^{r(t-t_0)}$.
• 由于其增长形式为指数形式，故该模型又称为指数增长模型.

---

模型的特点

• 贡献：
  
  • 首次提出了利用微分方程模型来研究人口的增长规律.
  • 假设人口的增长率与当前人口规模之间存在定量的联系.
  • 在一定的时间范围内，可以认为人口的自然增长率是恒定的，有其合理性.
• 不足：
  
  • 默认假设种群生存依赖的资源是无限的，没有考虑资源受限情况下种群内部竞争对增长率的影响.
  • 更合理的模型形式应该是 $r(t,p)$ 不为常数，且是随着 $p$ 的增长而单调递减.

---

5.3 Logistic 模型

该模型由 Pierre François Verhulst 于 1838年-1847 年间陆续发表在三篇论文中，他在 Adolphe Quetelet 的指导下，通过引入修正项修改了 Malthus 的指数增长模型，使之更加符合现实中资源受限的情况，该模型的解所对应的一类函数被其统称为 Logistics 函数.

---

Logistics 函数

$$\bbox[yellow]{f(t)=\dfrac{L}{1+\mathrm{e}^{-k(t-t_0)}}}$$

• 常被用于描述生物种群的自限性增长（Self-limiting Growth）：
  
  • 增长的初始阶段近似于指数增长；
  • 随着增长逐渐饱和，曲线放缓至接近线性；
  • 在成熟阶段，增长停止.

---

Logistics 模型

$$\dfrac{\mathrm{d} p}{\mathrm{d} t}=rp\bbox[yellow]{-bp^2}.$$

• 给定初值 $p(t_0)=p_0$，可以解得 $p(t)=\dfrac{rp_0}{bp_0+(r-bp_0)\mathrm{e}^{-r(t-t_0)}}$
• Logistic 模型中的 $-bp^2$ 被称为 阻滞项，因此 Logistic 模型也被称为阻滞增长模型.

---

模型的解释

$$p(t)=\dfrac{rp_0}{bp_0+(r-bp_0)\mathrm{e}^{-r(t-t_0)}}$$

• 不难验证 $\lim\limits_{t\to\infty}p(t)=\dfrac{r}{b}$.
• 不论初始的种群规模如何，随着时间的增长，群体规模总是小于并且趋于极限值 $r/b$.
  
  • $r/b$ 称为该种群的 最大容纳量.

---

Logistic 模型的应用

记 $N=r/b$，则 Logistics 模型中的微分方程可以写为

$$\dfrac{\mathrm{d} p}{\mathrm{d} t}=rp\left(1-\frac{p}{N}\right)$$

• 其中固有增长率 $r$ 为常数，$N$ 表示环境资源对该种群的最大容量.
• 1911年 A. G. McKendrick 应用改模型研究了肉汤中细菌的生长情况.
• 有人曾用该模型对 1790～1950 年美国人口数量作过预测，与实际数据相当吻合，误差不超过 $2.5\%$.

---

传染病的传播

• 在人群中未被免疫的新型传染性病原体，通常会在早期呈指数级传播，有大量易感个体尚未被感染.
• 2020 年初，导致 2019 冠状病毒病的 SARS-CoV-2 病毒在多国的感染过程中呈现出指数级增长.
• 此后，易感宿主减少（持续感染直到超过群体免疫阈值）或通过社交距离措施减少潜在宿主的被传染概率等因素，可能使呈指数增长的传染曲线首先线性化，然后趋缓，达到最大值.

---

5.4 可开发的单种群模型

例： 考察一个渔场，我们要建立一个在有捕捞条件下鱼的总量所满足的方程，并且在稳定的前提下讨论如何控制捕捞使持续产量最大.

---

按比例捕捞的渔场鱼量模型

• 设 $t$ 时刻渔场中鱼的总量为 $p(t)$, $r$ 为固有增长率，$N$ 为环境资源允许的最大鱼量.
• 在无捕捞条件下，$p(t)$ 服从 Logistic 模型：$\dfrac{\mathrm{d} p}{\mathrm{d} t}=rp\left(1-\dfrac{p}{N}\right)$
• 假设：单位时间的捕捞量 $h$ 与渔场鱼量成正比，比例系数为 $k$，表示单位时间捕捞率，也即 $h(p)=k\cdot p$.
• 有捕捞条件下渔场的鱼量模型：$\dfrac{\mathrm{d} p}{\mathrm{d} t}=rp\left(1-\dfrac{p}{N}\right)\bbox[yellow]{-kp}$.

---

模型分析

• 记 $F(p)=rp\left(1-\frac{p}{N}\right)-kp$，以下使用微分方程的定性理论，对模型的解进行讨论.
• 平衡点 / 不动点：满足 $F(p)=0$ 的点称为方程的平衡点.
  
  • 解得 $\displaystyle p_1=N\left(1-\frac{k}{r}\right), p_2=0$.
  • 进一步地 $F'(p_1)=k-r, F'(p_2)=r-k$.
• 称平衡点 $p^*$ 是稳定的，是指：对方程的任一个解 $p=p(t)$，恒有 $\bbox[yellow]{\lim\limits_{t\to \infty}p(t)=p^*}$.

---

平衡解稳定性的判断

• 判断平衡点 $p^*$ 是否稳定，一般需要解方程，很多时候是行不通的.
• 根据一阶近似方程判断，改写方程为：$\dfrac{\mathrm{d} p}{\mathrm{d} t}=F(p)= F'(p^*)(p-p^*)$
• 近似方程的一般解为：$p(t)=C\cdot \mathrm{e}^{F'(p^*)t}+p^*$.
  
  • 若 $F'(p^*)<0$，则 $p^*$ 是稳定平衡点；
  • 若 $F'(p^*)>0$，则 $p^*$ 不是稳定平衡点.
  • 若 $F'(p^*)=0$，则不能直接作出判断（可能稳定，也可能不稳定）.

---

结论

本例中，$F'(p_1)=k-r$，$F'(p_2)=r-k$，所以

• 当 $k<r$ 时，$F'(p_1)<0,\ F'(p_2)>0$，故 $p_1$ 是稳定平衡点，$p_2$ 不是；
• 当 $k>r$ 时，$F'(p_1)>0,\ F'(p_2)<0$，故 $p_1$ 不是稳定平衡点，$p_2$ 是.
• 结论：
  
  • 捕捞适度 ($k<r$)情况下，可使渔场鱼量持续稳定在 $p_1=N\left(1-\dfrac{k}{r}\right)$ 上下，并获得持续产量 $h(p_1)=kp_1$.
  • 若发生了捕捞过度（$k>r$），渔场鱼量将减至 $p_2=0$，这意味着开发将是不可持续的，最终对渔场产生不可逆转的破坏.

---

进一步的讨论

问题：如何控制捕捞强度 $k$ 使得持续产量 $h(p_1)=kp_1$ 最大？

• $\displaystyle h(p_1)=kp_1=N\left(1-\dfrac{k}{r}\right)\cdot k\triangleq h(k)$
• $\displaystyle\dfrac{\mathrm{d} h}{\mathrm{d} k}=N\left(1-\frac{2k}{r}\right)=0\ \Rightarrow k_{\max}=\frac{r}{2}$
• 对应的 $\displaystyle h_{\max}=\frac{r}{2}N\left(1-\frac{1}{2}\right) =\frac{rN}{4},\quad p_{\max}=\frac{N}{2}$.
• 结论：控制捕捞强度 $k=\displaystyle\frac{r}{2}$，使渔场鱼量保持在最大鱼量 $N$ 的一半时，可以获得最大的持续产量 $h_m=rN/4$.

---

肿瘤生长模型

• 在医学上，以 $x(t)$ 表示肿瘤的大小，其变化动态可描述为 $x'=r\left(1-{\frac {x}{K}}\right)x$.
• 如果采用化疗等手段，则模型可修改为 $x'=r\left(1-{\frac {x}{K}}\right)x-c(t)x$，其中 $c(t)$ 为治疗引起的肿瘤死亡率.
• 在理想的长期的治疗下，$c(t)$ 可表示为周期为 $T$ 的周期函数或常数函数.
• 结论：若 $\displaystyle {\frac {1}{T}}\int _{0}^{T}c(t)\,dt>r$，则 $\lim _{t\to +\infty }x(t)=0$，即：如果平均治疗引起的肿瘤死亡率大于基线增殖率，则疾病能被根除.

---

5.5 弱肉强食模型

• 弱肉强食模型也称为 Lotka-Volterra 方程，或 捕食者-食饵模型 (Predator Prey Model).
• 经常用来描述生物系统中，掠食者与猎物进行互动时的动态模型，也就是两者族群规模的消长.
• 此方程分别在 1925 年与 1926 年，由 Alfred James Lotka 与 Vito Volterra 独立发表.

---

Lotka-Volterra 方程

$$\begin{cases} \dfrac{\mathrm{d}x}{\mathrm{d}t}=\alpha x-\beta xy & \text{被捕食者（猎物）}\\ \dfrac{\mathrm{d}y}{\mathrm{d}t}=-\gamma y+\delta xy & \text{捕食者（掠食者）} \end{cases}$$

• 假设：
  
  • 猎物所接受的食物供给已经达到最大极限，且除非遭遇掠食者的捕食，否则繁殖数量的增加以指数方式成长.

---

模型的解

---

微分方程求解器

---

种群的相互依存

• 当猎物数量充足的时候，掠食者的族群也会兴旺起来.
• 不过掠食者的族群最后仍然会因为超过猎物所能供给的数量而开始衰减.
• 当掠食者的族群族群缩减，则猎物族群将会再次增大.
• 两者的族群大小便以周期性的成长与衰减进行循环.

---

族群规模的平衡

• 令 $\begin{cases} P(x,y)=\alpha x-\beta xy=0\\ Q(x,y)=-\gamma y+\delta xy=0 \end{cases}$，解得两个平衡点
  
  1. $(x,y)=(0,0)$：两个物种完全灭绝.
  2. $(x,y)=\left(\dfrac{\gamma}{\delta},\dfrac{\alpha}{\beta}\right)$：两个族群能够维持一个不为零的数量，并且能够永久持续.

---

线性动力系统的平衡点

• 对于形如 $\left\{\begin{array}{l}\dfrac{\mathrm{d}x}{\mathrm{d}t}=P(x,y)\\ \dfrac{\mathrm{d}y}{\mathrm{d}t}=Q(x,y)\end{array}\right.$ 的微分动力系统，
• 在其平衡点附近，通常可将其近似为线性动力系统 $\left\{\begin{array}{l}\dfrac{\mathrm{d}x}{\mathrm{d}t}=P'_xx+P'_yy\\ \dfrac{\mathrm{d}y}{\mathrm{d}t}=Q'_xx+Q'_yy\end{array}\right.$，
• 通过对 $\boldsymbol{A}=\boldsymbol{J}(x,y)=\dfrac{\partial(P,Q)}{\partial (x,y)}=\begin{bmatrix}P'_x & Q'_x \\ P'_y & Q'_y\end{bmatrix}$ 的分析，可得到该线性动力系统的解的性态特征，进而对原微分动力系统的性态有所了解.

---

线性自治系统平衡点的分类

---

弱肉强食模型的平衡点

$$\boldsymbol{A}=\boldsymbol{J}(x,y)=\begin{bmatrix}\alpha-\beta y & -\beta x \\ \delta y & \delta x-\gamma\end{bmatrix}$$

• $J(0,0)=\begin{bmatrix}\alpha & 0 \\ 0 & \gamma\end{bmatrix}$
  
  • 特征值 $\lambda_1=\alpha$，$\lambda_2=-\gamma$.
  • $\mathrm{det}\boldsymbol{A}<0$，因此该平衡点是一个鞍点 (Saddle).
• $J\left(\dfrac{\gamma}{\delta},\dfrac{\alpha}{\beta}\right)=\begin{bmatrix}0 & -\frac{\beta\gamma}{\delta} \\ \frac{\alpha\delta}{\beta} & 0\end{bmatrix}$
  
  • 特征值 $\lambda_{1,2}=\pm\mathrm{i}\sqrt{\alpha\gamma}$.
  • $\mathrm{tr}(\boldsymbol{A})=0$ 且 $\mathrm{det}\boldsymbol{A}>0$，因此该平衡点为一个 焦点 (Center).

---

Lotka-Volterra 方程组的相图

为了对弱肉强食模型进行讨论，可以将其转换为一个自治系统，形如：

$$\dfrac{\mathrm{d}y}{\mathrm{d} x}=\dfrac{-y(\gamma-\delta x)}{x(\alpha-\beta y)}.$$

• 在 $xOy$ 平面上画出其解的图像，不同的初值可能对应于不同的封闭曲线（轨线）：$(y^{\alpha}\mathrm{e}^{-\beta y})(x^{\gamma}\mathrm{e}^{-\delta x})=C\;(C\in\mathbb{R})$
• 随着时间的增长，所有轨线上点都围绕焦点逆时针运动.

---

5.6 多种群模型

• 假设：
  
  • 资源受限.
  • 种群间存在相互影响.
• 根据种群间相互影响的不同，典型的多种群模型可分为
  
  • 相互竞争的模型
  • 相互依存模型
  • 弱肉强食模型（捕食者-食饵模型）

---

一般形式的多种群模型

$$\left\{ \begin{array}{l} \displaystyle\frac{\mathrm{d} x_1}{\mathrm{d} t}= r_1 x_1\left(1-\frac{x_1+\alpha x_2}{N_1}\right) \\ \displaystyle\frac{\mathrm{d} x_2}{\mathrm{d} t}= r_2 x_2\left(1-\frac{x_2+\beta x_1}{N_2}\right) \\ \end{array} \right.$$

• $N_1,N_2$ 分别表示两个物种的环境最大承载量.
• $\alpha$，$\beta$ 的不同取值决定了这两个种群的不同关系. 一般来说，
  
  • $\alpha, \beta>0$，对应种群间 相互竞争的模型.
  • $\alpha, \beta<0$，对应种群间 相互依存的模型.
  • $\alpha\cdot \beta<0$，对应 弱肉强食模型.

---

更多不同种群的情况

$$\dfrac{\mathrm{d}x_i}{\mathrm{d}t}=r_ix_i\left(1-\sum\limits_{j=1}^n\dfrac{a_{ij}x_j}{N_i}\right),\; i=1,2,...,n.$$

• $N_i$ 为每个物种的环境最大承载量.
• 矩阵 $\boldsymbol{A}=[a_{ij}]_{n\times n}$ 反映了物种之间的相互影响.

---

Lorenz 方程与奇异吸引子

$$\begin{cases}{\dfrac {\mathrm {d} x}{\mathrm {d} t}}&=\sigma (y-x),\\[6pt]{\dfrac {\mathrm {d} y}{\mathrm {d} t}}&=x(\rho -z)-y,\\[6pt]{\dfrac {\mathrm {d} z}{\mathrm {d} t}}&=xy-\beta z.\end{cases}$$

• 右图参数：$\rho=28$，$\sigma=10$，$\beta=8/3$.

---

5.7 综合案例：鱼群的可持续捕捞

• 为了保护人类赖以生存的自然环境，可再生资源(如渔业、林业资源)的开发必须适度.
• 一种合理简化的策略是，在实现可持续收获的前提下，追求最大的产量或最佳效益.
• 考虑对某种鱼的最优捕捞策略: 假设这种鱼分 4 个年龄组，称 1 龄鱼，...，4 龄鱼。各年龄组每条鱼的平均重量分别为 5.07，11.55，17.86，22.99(克) ，各年龄组鱼的自然死亡率为 0.8(1/年).

---

产卵繁殖

• 这种鱼为季节性集中产卵繁殖，平均每条 4 龄鱼的产卵量为 $1.109\times 10^5$，3 龄鱼的产卵量为这个数的一半，2 龄鱼和 1 龄鱼不产卵.
• 产卵和孵化期为每年的最后 4 个月，卵孵化并成活为 1 龄鱼，成活率( 孵化出的 1 龄鱼条数与产卵总量 $n$ 之比)为：$\dfrac{1.22\times 10^{11}}{1.22\times 10^{11}+n}$.

---

捕捞规定

• 渔业管理部门规定，每年只允许在产卵孵化期前的 8 个月内进行捕捞作业.
• 如果每年投入的捕捞能力(如渔船数、下网次数等)固定不变，这时单位时间捕捞量将与各年龄组鱼群条数成正比，比例系数称为捕捞强度系数.
• 捕鱼通常使用 13mm 网眼的拉网，这种网只能捕捞 3 龄鱼和 4 龄鱼，两种鱼的捕捞强度系数之比为 0.42:1，渔业上称这种方式为固定网力量捕捞.

---

待解决的问题

建立数学模型分析如何实现可持续捕获 (即每年开始捕捞时渔场中各年龄组鱼群条数基本不变), 并且在此前提条件下得到最高的年收获量(捕捞总重量).

---

问题假设

• 假设鱼群每年在 8 月最后一天产卵，卵在 12 月最后一天孵化.
• 假设各龄鱼群数量随时间呈周期性连续变化，周期为一年，因此仅考虑一年内各龄鱼群数量变化，同时将 1 年中的 12 个月映射到 [0,1] 连续区间.
• 仅考虑该鱼群自身繁殖和死亡以及捕捞带来的数量变化.

---

模型参数

• 设 $x_i(t)$ 为 $t$ 时刻 $i(i=1,2,3,4)$ 龄鱼的数量，$k$为 4 龄鱼的捕捞强度系数，$n$ 为每年产卵量.
• 根据假设 2 和假设 3，不同年龄的鱼群在一年内只有自然死亡或被捕捞，不存在鱼群数量增长情况，也不存在不同年龄鱼群之间的数量变化.

---

不同年龄鱼群的种群变化模型

$$\left\{\begin{array}{ll} \dfrac{\mathrm{d}x_i(t)}{\mathrm{d}t}=-0.8x_i(t), & t\in[0,1],\;i=1,2.\\ \dfrac{\mathrm{d}x_i(t)}{\mathrm{d}t}=-0.8x_i(t), & t\in(2/3,1],\;i=3,4.\\ \dfrac{\mathrm{d}x_3(t)}{\mathrm{d}t}=-(0.8+0.42k)x_3(t), & t\in[0,2/3].\\ \dfrac{\mathrm{d}x_4(t)}{\mathrm{d}t}=-(0.8+k)x_4(t), & t\in[0,2/3]. \end{array}\right.$$

---

模型的初值

• 只有 3 龄鱼和 4 龄鱼在 8 月底产卵，因此，产卵量可表达为：
  
  • $n=1.109\times 10^5\times \left[0.5x_3(2/3)+x_4(2/3)\right]$.
• 根据鱼卵的成活率可知，年度开始时1龄鱼数量：
  
  • $x_1(0)=n\times\dfrac{1.22\times 10^{11}}{1.22\times 10^{11}+n}$.
• 2 龄鱼、3 龄鱼和 4 龄鱼在年初时数量满足：
  
  • $x_2(0)=x_1(1)$，$x_3(0)=x_2(1)$，$x_4(0)=x_3(1)$.
• 六个独立微分方程，通过初值之间的关联构成了一个微分方程组.

---

最终，$k$ 和 $x_1(0)$ 满足如下模型：

$$\begin{array}{l|l} \left\{\begin{array}{ll} \dfrac{\mathrm{d}x_i(t)}{\mathrm{d}t}=-0.8x_i(t), & t\in[0,1],\;i=1,2.\\ \dfrac{\mathrm{d}x_i(t)}{\mathrm{d}t}=-0.8x_i(t), & t\in(2/3,1],\;i=3,4.\\ \dfrac{\mathrm{d}x_3(t)}{\mathrm{d}t}=-(0.8+0.42k)x_3(t), & t\in[0,2/3].\\ \dfrac{\mathrm{d}x_4(t)}{\mathrm{d}t}=-(0.8+k)x_4(t), & t\in[0,2/3].\\ n=1.109\times 10^5\times \left[0.5x_3(2/3)+x_4(2/3)\right] \end{array}\right. & \begin{array}{l} \text{初始条件：}\\ x_1(0)=n\times\dfrac{1.22\times 10^{11}}{1.22\times 10^{11}+n}\\ x_2(0)=x_1(1)\\ x_3(0)=x_2(1)\\ x_4(0)=x_3(1) \end{array} \end{array}$$

---

目标函数

• 寻找最佳的捕捞强度，使得年捕捞量达到最大，即确定 $k$ 使得以下目标函数达到最大：
  
  • $\max f(k)=\max\left(17.86\int_{0}^{\frac23}0.42kx_3(t)\mathrm{d}t+22.99\int_0^{\frac23}kx_4(t)\mathrm{d}t\right)$.
• 约束条件由前述微分方程确定，即对于给定的 $k$ 根据前述微分方程组模型确定 $x_1(0)$ 和 $x_3(t)$, $x_4(𝑡)$, $t∈[0, 2/3]$，进而决定了捕捞量.

---

---

结论

• 当 $k=17.36$ 时，总获得的最大捕鱼量为 $3.88\times 10^{11}$g.
• 各年龄段鱼在不同时间点的数量如右表所示

$$\boxed{\begin{array}{c|c} \text{鱼的阶段} & \text{数量}\\ \hline \text{1月初1龄鱼} & 1.196\times 10^{11}\\ \text{1月初2龄鱼} & 5.374\times 10^{10}\\ \text{1月初3龄鱼} & 2.415\times 10^{10}\\ \text{1月初4龄鱼} & 8.396\times 10^{7}\\ \text{8月底3龄鱼} & 1.067\times 10^{8}\\ \text{8月底4龄鱼} & 436 \end{array}}$$

---

课后思考题

1. ($\star\star\star$) 驱逐舰在黑夜中搜索潜艇，据雷达测量发现潜艇在正北三公里处，但潜艇立即下潜. 已知潜艇下潜后立即朝某个未知方向直线行驶，且驱逐舰的速度是潜艇速度的2倍。问驱逐舰应采取何种路线追逐，才能保证它能开到潜艇的正上方？
2. ($\star\star$) 在可开发的单种群模型中，如果渔场鱼量的自然增长仍服从 Logistic 规律，$r$ 为固有增长率，$N$ 为环境资源允许的最大鱼量，而单位时间捕捞量为常数 $q$. 请分别就 $q>rN/4$, $q<rN/4$, $q=rN/4$ 这 3 种情况讨论渔场鱼量方程的平衡点及其稳定性.