predator-prey model OpenFOAM 编程 | 求解捕食者与被捕食者模型问题（ODEs） _生活百科

0. 写在前面本文问题参考自文献 \(^{[1]}\) 第一章例 6，并假设了一些条件，基于 OpenFOAM-v2206 编写程序数值上求解该问题。笔者之前也写过基于 OpenFOAM 求解偏分方程的帖子，OpenFOAM 编程 | One-Dimensional Transient Heat Conduction 。
1. 问题描述假设一群山猫（捕食者）和一群山兔（被捕食者）生活在同一片区域，那么我们可以知道，山猫吃了山兔，繁殖力会增强，山猫的数量会增加。这样一来，山兔的数量会随之减少。接下来，山猫由于食物短缺而数量减少，进而导致山兔遇到山猫的机会减少（被吃掉的概率降低），结果山兔的数量又逐渐增加，这样山猫得到食物的机会也随之增加，其数量又再一次增加，而山兔的数量又会再一次随之减少，如此不断循环。
2. 解析求解设任意 \(t\) 时刻山兔与山猫的数量分别是 \(\phi\) 和 \(\psi\) ，二者的变化服从下面动力学方程
\[\begin{aligned}\frac{\mathrm{d}\phi}{\mathrm{d}t} &= k_1 \phi - \mu\phi\psi \\\frac{\mathrm{d}\psi}{\mathrm{d}t} &= \nu\phi\psi - k_2 \psi\end{aligned}\tag1\]其中，\(k_1\)，\(k_2\)，\(\mu\) 和 \(\nu\) 都是正常数。
在上述方程中有几点需要注意：

\(k_1\phi\) 表示山兔种群的净增长率，与山兔种群数量成正比。
\(-\mu\phi\psi\) 表示山兔被山猫吃掉而导致的减少率，与乘积 \(\phi\psi\) （可表示两种动物的相遇概率）成正比。
\(\nu\phi\psi\) 表示山猫种群的增长率，由于其数量增长取决于捕食（相遇才有可能），因此 \(\nu\) 为正值。
\(-k_2\psi\) 表示山猫种群的死亡率，与其种群数量成正比。

方程组（1）因为含有乘积项，因此是非线性的。现采用线性化的特殊方法求解，即研究种群数量 \(\phi\) 和 \(\psi\) 在其稳定值附近的微小涨落。设方程组（1）的稳态解为 \(\phi=\phi_0\)，\(\psi=\psi_0\)，它们由下面条件决定
\[\begin{aligned}\left . \frac{\mathrm{d}\phi}{\mathrm{d}t} \right |_{\phi=\phi_0,\psi=\psi_0} &= 0 \\\left . \frac{\mathrm{d}\psi}{\mathrm{d}t} \right |_{\phi=\phi_0,\psi=\psi_0} &=0\end{aligned}\]也就是
\[\begin{aligned}k_1 \phi_0 - \mu\phi_0\psi_0 &= 0 \\\nu\phi_0\psi_0 - k_2 \psi_0&=0\end{aligned}\tag2\]代数方程（2）的解为
\[\begin{aligned}\phi_0 &= \frac{k_2}{\nu} \\\psi_0 &=\frac{k_1}{\mu}\end{aligned}\]现在，将方程组（1）的解写为下面形式
\[\begin{aligned}\phi &= \phi_0+ \xi \\\psi &= \psi_0 + \eta\end{aligned}\]【predator-prey model OpenFOAM 编程 | 求解捕食者与被捕食者模型问题（ODEs）】其中， \(\xi\) 和 \(\eta\) 与 \(\phi_0\) 和 \(\psi_0\) 相比都是小量。将上述解带入方程组（1）中可以得到关于变量 \(\xi\) 和 \(\eta\) 的方程组
\[\begin{aligned}\frac{\mathrm{d}\xi}{\mathrm{d}t} &= k_1\xi-\mu\phi_0\eta-\mu\psi_0\xi-\mu\xi\eta\\\frac{\mathrm{d}\eta}{\mathrm{d}t} &= \nu\phi_0\eta + \nu\psi_0\xi - k_2\eta+\nu\xi\eta\end{aligned}\tag3\]其中非线性项 \(\mu\xi\eta\) 和 \(\nu\xi\eta\) 为二阶小量，可以忽略；再将稳态解代入可得线性化的耦合方程组
\[\begin{aligned}\frac{\mathrm{d}\xi}{\mathrm{d}t} &= -k_2\frac{\mu}{\nu}\eta\\\frac{\mathrm{d}\eta}{\mathrm{d}t} &= k_1\frac{\nu}{\mu}\xi\end{aligned}\]解耦后可得到
\[\begin{aligned}\frac{\mathrm{d}^2\xi}{\mathrm{d}t^2} +k_1k_2\xi&= 0\\\frac{\mathrm{d}^2\eta}{\mathrm{d}t^2} +k_1k_2\eta&= 0\end{aligned}\tag4\]可以知道，式（4）与 L-C 震荡电路及单摆问题同属于相同的数学模型
\[\frac{\mathrm{d}^2y}{\mathrm{d}t^2} + k^2 y = 0\]其通解为
\[y(t) = E\sin(kt+\delta)\ \ \ \ 或\ \ \ \ y(t) = E\cos(kt+\delta)\]其中，\(E\) 和 \(\delta\) 为振幅和初相位，与具体问题有关。
那么我们也可以得到本问题的最终解的形式为
\[\begin{aligned}\phi &= \frac{k_2}{\nu} + E_1 \sin\left(\sqrt{k_1k_2}t+\delta_1\right)\\\psi &= \frac{k_1}{\mu} +E_2 \sin\left(\sqrt{k_1k_2}t+\delta_2\right) \\\end{aligned}\]其中，每个公式中振幅与初相位取决于各自的初始条件。
3. 数值求解从上一节可知，我们需要数值求解一个耦合的常微分方程组，可以用RungeKutta法\(^{[2]}\) 。简单推导过程如下：
\[\begin{aligned}\frac{\mathrm{d}\phi}{\mathrm{d}t} &= f_1\left( \phi,\psi \right) \\\frac{\mathrm{d}\psi}{\mathrm{d}t} &= f_2\left( \phi,\psi \right) \\\end{aligned}\]其中，
\[\begin{aligned}f_1\left( \phi,\psi \right) &=k_1 \phi - \mu\phi\psi \\f_2\left( \phi,\psi \right) &= \nu\phi\psi - k_2 \psi \\\end{aligned}\]四阶Runge-Kutta方法可以表示为：
\[\begin{aligned}\phi^{k+1} &= \phi^{k} + \frac{\Delta t}{6} \left(f_{11} + 2f_{12} + 2f_{13} + f_{14} \right) \\\psi^{k+1} &= \psi^{k} + \frac{\Delta t}{6} \left(f_{21} + 2f_{22} + 2f_{23} + f_{24} \right) \\\end{aligned}\]其中，