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

技术分享 3年前 (2022-11-06) 0 999+

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 - muphipsi \ frac{mathrm{d}psi}{mathrm{d}t} &= nuphipsi - k_2 psi end{aligned} tag1 ]

其中，(k_1)，(k_2)，(mu) 和 (nu) 都是正常数。

在上述方程中有几点需要注意：

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

方程组（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 - muphi_0psi_0 &= 0 \ nuphi_0psi_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} ]

其中，(xi) 和 (eta) 与 (phi_0) 和 (psi_0) 相比都是小量。将上述解带入方程组（1）中可以得到关于变量 (xi) 和 (eta) 的方程组

[begin{aligned} frac{mathrm{d}xi}{mathrm{d}t} &= k_1xi-muphi_0eta-mupsi_0xi-muxieta\ frac{mathrm{d}eta}{mathrm{d}t} &= nuphi_0eta + nupsi_0xi - k_2eta+nuxieta end{aligned} tag3 ]

其中非线性项 (muxieta) 和 (nuxieta) 为二阶小量，可以忽略；再将稳态解代入可得线性化的耦合方程组

[begin{aligned} frac{mathrm{d}xi}{mathrm{d}t} &= -k_2frac{mu}{nu}eta\ frac{mathrm{d}eta}{mathrm{d}t} &= k_1frac{nu}{mu}xi end{aligned} ]

解耦后可得到

[begin{aligned} frac{mathrm{d}^2xi}{mathrm{d}t^2} +k_1k_2xi&= 0\ frac{mathrm{d}^2eta}{mathrm{d}t^2} +k_1k_2eta&= 0 end{aligned} tag4 ]

可以知道，式（4）与 L-C 震荡电路及单摆问题同属于相同的数学模型

[frac{mathrm{d}^2y}{mathrm{d}t^2} + k^2 y = 0 ]

其通解为

[y(t) = Esin(kt+delta) 或 y(t) = Ecos(kt+delta) ]

其中，(E) 和 (delta) 为振幅和初相位，与具体问题有关。

那么我们也可以得到本问题的最终解的形式为

[begin{aligned} phi &= frac{k_2}{nu} + E_1 sinleft(sqrt{k_1k_2}t+delta_1right)\ psi &= frac{k_1}{mu} +E_2 sinleft(sqrt{k_1k_2}t+delta_2right) \ end{aligned} ]

其中，每个公式中振幅与初相位取决于各自的初始条件。

3. 数值求解

从上一节可知，我们需要数值求解一个耦合的常微分方程组，可以用RungeKutta法(^{[2]})。简单推导过程如下：

[begin{aligned} frac{mathrm{d}phi}{mathrm{d}t} &= f_1left( phi,psi right) \ frac{mathrm{d}psi}{mathrm{d}t} &= f_2left( phi,psi right) \ end{aligned} ]

其中，

[begin{aligned} f_1left( phi,psi right) &= k_1 phi - muphipsi \ f_2left( phi,psi right) &= nuphipsi - 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} ]

其中，

[begin{aligned} f_{i1} &= f_i left( phi_k, psi_k right) \ f_{i2} &= f_i left( phi_k+frac{Delta t}{2}f_{11}, psi_k+frac{Delta t}{2}f_{21} right) \ f_{i3} &= f_i left( phi_k+frac{Delta t}{2}f_{12}, psi_k+frac{Delta t}{2}f_{22} right) \ f_{i4} &= f_i left( phi_k+{Delta t}f_{11}, psi_k+{Delta t}f_{21} right) \ end{aligned} i=1,2 ]

求解代码采用 Python 编写，如下所示

#!/usr/bin/python3 # -*- coding:utf-8 -*-  import numpy as np  k1 = 0.7 k2 = 0.5 mu = 0.1 nu = 0.02  def f1(phi,psi):     return k1*phi-mu*phi*psi  def f2(phi,psi):     return nu*phi*psi-k2*psi  tStart = 0 tEnd   = 100.0 n      = 100000 deltaT = tEnd / n halfDeltaT = deltaT / 2.0 Solution = np.ndarray([n+1,2]) Solution[0] = [30,20]   for i in range(n):     f11 = f1(Solution[i][0], Solution[i][1])     f21 = f2(Solution[i][0], Solution[i][1])      f12 = f1(Solution[i][0] + halfDeltaT * f11, Solution[i][1] + halfDeltaT * f21)     f22 = f2(Solution[i][0] + halfDeltaT * f11, Solution[i][1] + halfDeltaT * f21)      f13 = f1(Solution[i][0] + halfDeltaT * f12, Solution[i][1] + halfDeltaT * f22)     f23 = f2(Solution[i][0] + halfDeltaT * f12, Solution[i][1] + halfDeltaT * f22)      f14 = f1(Solution[i][0] + deltaT * f11, Solution[i][1] + deltaT * f21)     f24 = f2(Solution[i][0] + deltaT * f11, Solution[i][1] + deltaT * f21)      Solution[i+1][0] = Solution[i][0] + deltaT / 6.0 * (f11 + 2*f12 + 2*f13 + f14)     Solution[i+1][1] = Solution[i][1] + deltaT / 6.0 * (f21 + 2*f22 + 2*f23 + f24)     print((i+1)*deltaT,Solution[i+1][0],Solution[i+1][1])

4. OpenFOAM 求解

使用OpenFOAM 数值求解常微分方程（组）主要用到 ODESystem.H（构造微分方程系统）和 ODESolver.H（求解器）；此外，在 OpenFOAM 中需要对常微分方程（组）进行整理(^{[3]})，进而方便编写代码进行求解。

对于任意阶常微分方程可以转化为一系列一阶常微分方程，这个过程称为降阶，一阶常微分方程的个数与原方程的阶数相等（对于耦合常微分方程组，其阶数等于所有方程阶数之和）。对于某个 (n) 阶常微分方程，可按下面形式降阶

[y^{(n)}(x) = f left( x, y^{(0)}, y^{(1)},ldots,y^{(n-1)} right) ]

其中，(n) 为阶数，(y^{(0)}=y) 。

进一步，引入符号 (mathrm{D}) 对各阶导数重新定义，此过程称为转换

[mathrm{D}_j = y^{(j-1)} j=1,2,ldots,n-1 ]

最终，使用新符号重新表达原系统，此过程称为诱导

[begin{aligned} mathrm{D}'_j &= mathrm{D}_{j+1} \ mathrm{D}'_n = y^{(n)} &= fleft( x, mathrm{D}_1, mathrm{D}_2,ldots,mathrm{D}_n right) end{aligned} ]

在 OpenFOAM 中，存在另外一个过程，该过程仅与刚性系统求解器相关，这类求解器需要雅可比矩阵和对自变量的偏导数，即

[J = begin{bmatrix} frac{partial mathrm{D}'_1}{partial mathrm{D}_1} & frac{partial mathrm{D}'_1}{partial mathrm{D}_2} & cdots & frac{partial mathrm{D}'_1}{partial mathrm{D}_n}\ frac{partial mathrm{D}'_2}{partial mathrm{D}_1} & frac{partial mathrm{D}'_2}{partial mathrm{D}_2} & cdots & frac{partial mathrm{D}'_2}{partial mathrm{D}_n}\ vdots & vdots & ddots & vdots \ frac{partial mathrm{D}'_n}{partial mathrm{D}_1} & frac{partial mathrm{D}'_n}{partial mathrm{D}_2} & cdots & frac{partial mathrm{D}'_n}{partial mathrm{D}_n}\ end{bmatrix} 和 frac{partial mathrm{D}'_1}{partial x},frac{partial mathrm{D}'_2}{partial x},,ldots, frac{partial mathrm{D}'_n}{partial x} ]

接下来，我们看一下如何实现相关求解代码。首先看一下如何构造方程系统。系统代码需要继承 Foam::ODESystem 抽象类，并且需要全部实现三个方法nEqns()、 derivatives() 和 jacobian()，其中 jacobian() 方法对于非刚性求解器可以将实现置空（空函数体）。

让我们重新回顾一下公式（1），可知 nEqns() 应该返回 2；此外，定义 (Y=[phi,psi]^{mathrm{T}}) ，公式（1）可整理成如下向量形式

[frac{mathrm{d}Y}{mathrm{d}t} = begin{bmatrix} k_1 & -muphi \ nupsi & -k_2 \ end{bmatrix} Y ]

因此，导数可按照公式（1）编写即可，只不过需要注意是向量形式。最后，对应之前的描述的降阶过程，可以知道

[Y' = fleft( t, Yright) ]

进而可以知道， (D_1 = Y, D'_1=Y')，可得到雅可比矩阵和对自变量的偏导数分别为

[frac{partial mathrm{D}'_1}{partial mathrm{D}_1} = frac{partial Y'}{partial Y} = begin{bmatrix} k_1 & -muphi \ nupsi & -k_2 \ end{bmatrix}, frac{partial mathrm{D}'_1}{partial t} = 0 ]

需要注意的是，雅可比矩阵只有一个元素 (frac{partial mathrm{D}'_1}{partial mathrm{D}_1})，只不过这个元素是一个块的形式。

具体代码实现如下所示

#include "ODESystem.H"  class ODEs : public Foam::ODESystem { public:     ODEs() {}     ~ODEs() {}     // 初始化参数     ODEs(const Foam::scalar k1, const Foam::scalar mu, const Foam::scalar k2,          const Foam::scalar nu)     {         k1_ = k1;         mu_ = mu;         k2_ = k2;         nu_ = nu;     }     // 方程个数     Foam::label nEqns() const override { return 2; }     // 求导     void derivatives(const Foam::scalar x, const Foam::scalarField& y,                      Foam::scalarField& dydx) const override     { // 两个未知量存成向量，y[0] -> phi, y[1] -> psi         dydx[0] = k1_ * y[0] - mu_ * y[0] * y[1];         dydx[1] = nu_ * y[0] * y[1] - k2_ * y[1];     }     // 计算符号的雅可比矩阵和关于自变量的导数     void jacobian(const Foam::scalar x, const Foam::scalarField& y, Foam::scalarField& dfdx,                   Foam::scalarSquareMatrix& dfdy) const override     {         dfdx[0] = 0;         dfdx[1] = 0;          dfdy[0][0] = k1_;         dfdy[0][1] = -mu_ * y[0];          dfdy[1][0] = nu_ * y[1];         dfdy[1][1] = -k2_;     }  private:     Foam::scalar k1_;     Foam::scalar mu_;     Foam::scalar k2_;     Foam::scalar nu_; };

对应的，我们实现下主函数

#include <iostream> #include <memory>  #include "ODESystem.H" #include "ODESolver.H"  class ODEs : public Foam::ODESystem {     // 这里的代码在上边已经介绍，此处省略 };  int main(int argc, char* argv[]) {     const Foam::scalar startTime = 0.0;         // 开始时间     const Foam::scalar endTime   = 100.0;       // 结束时间     const Foam::scalar phi0      = 30;          // 山兔初始值     const Foam::scalar psi0      = 20;          // 山猫初始值     const Foam::label n          = 100000;      //     const Foam::scalar deltaT    = endTime / n; // 步长     // 系数，参考自文献[4]     const Foam::scalar k1 = 0.7;     const Foam::scalar mu = 0.1;     const Foam::scalar k2 = 0.5;     const Foam::scalar nu = 0.02;     // 构造对象     ODEs odes(k1, mu, k2, nu);      // 构造求解器，具体使用的算法通过参数传递     Foam::dictionary dict;     dict.add("solver", argv[1]);     Foam::autoPtr<Foam::ODESolver> solver = Foam::ODESolver::New(odes, dict);      // 初始化一些变量     Foam::scalar tStart = startTime;     Foam::scalarField PhiPsi(odes.nEqns()); // 因变量     PhiPsi[0] = phi0;     PhiPsi[1] = psi0;     Foam::scalarField ddt(odes.nEqns()); // 保存导数值      // 计算过程     for (Foam::label i = 0; i < n; ++i)     {         Foam::scalar dtEst = deltaT / 2;         Foam::scalar tEnd  = tStart + deltaT;         //         odes.derivatives(tStart, PhiPsi, ddt);         solver->solve(tStart, tEnd, PhiPsi, dtEst);         //         tStart = tEnd;         //         Foam::Info << tStart << "," << PhiPsi[0] << "," << PhiPsi[1] << Foam::endl;     }      return 0; }

此外，CMakeLists.txt 文件可参考笔者之前的随笔，如 OpenFOAM编程 | Hello OpenFOAM 和 OpenFOAM 编程 | One-Dimensional Transient Heat Conduction，此处不再赘述。

5. 数据分析

笔者通过命令行参数分别采用RKCK45 算法和 seulex 算法（需要用到雅可比矩阵）对该问题进行求解，从下图可见二者求解得到的结果是一致的。

同时运行笔者之前提到的 Python 代码后得到的数值结果与 OpenFOAM 计算结果绘制在同一张图中，二者高度重合。

同时，解析解法（线性化的特殊解法）得到的结论是二者均按照 (sqrt{k_1k_2}) 圆频率震荡，那么对应的周期为 $T = 2pi / sqrt{k_1k_2} = 2 pi / sqrt{0.7*0.5} approx 10.62 $，而数值解中得到的周期为 12.425，笔者认为在本文的条件假设下，其中的差距来自于线性解法中没有考虑非线性，但这个解法仍然具有实际意义；读者可以尝试改用绝对值较小的系数来降低其非线性程度。

另外，感兴趣的读者可以尝试使用 Matlab 或 GNU Octave 求解该问题。