欢迎您的访问
专注架构,Java,数据结构算法,Python技术分享

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

0. 写在前面

本文问题参考自文献

[1]第一章例 6,并假设了一些条件,基于 OpenFOAM-v2206 编写程序数值上求解该问题。笔者之前也写过基于 OpenFOAM 求解偏分方程的帖子,OpenFOAM 编程 | One-Dimensional Transient Heat Conduction

1. 问题描述

假设一群山猫(捕食者)和一群山兔(被捕食者)生活在同一片区域,那么我们可以知道,山猫吃了山兔,繁殖力会增强,山猫的数量会增加。这样一来,山兔的数量会随之减少。接下来,山猫由于食物短缺而数量减少,进而导致山兔遇到山猫的机会减少(被吃掉的概率降低),结果山兔的数量又逐渐增加,这样山猫得到食物的机会也随之增加,其数量又再一次增加,而山兔的数量又会再一次随之减少,如此不断循环。

2. 解析求解

设任意

t时刻山兔与山猫的数量分别是

ϕ和

ψ,二者的变化服从下面动力学方程

 

(1)dϕdt=k1ϕ−μϕψdψdt=νϕψ−k2ψ

其中,

k1,

k2,

μ和

ν都是正常数。

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

  1. k1ϕ表示山兔种群的增长率,与山兔种群数量成正比。
  2. −μϕψ表示山兔被山猫吃掉而导致的减少率,与乘积 ϕψ

    (可表示两种动物的相遇概率)成正比。

  3. νϕψ表示山猫种群的增长率,由于其数量增长取决于捕食(相遇才有可能),因此 ν

    为正值。

  4. −k2ψ表示山猫种群的死亡率,与其种群数量成正比。

方程组(1)因为含有乘积项,因此是非线性的。现采用线性化的特殊方法求解,即研究种群数量

ϕ和

ψ在其稳定值附近的微小涨落。设方程组(1)的稳态解为

ϕ=ϕ0,

ψ=ψ0,它们由下面条件决定

 

dϕdt|ϕ=ϕ0,ψ=ψ0=0dψdt|ϕ=ϕ0,ψ=ψ0=0

也就是

 

(2)k1ϕ0−μϕ0ψ0=0νϕ0ψ0−k2ψ0=0

代数方程(2)的解为

现在,将方程组(1)的解写为下面形式

其中,

ξ和

η与

ϕ0和

ψ0相比都是小量。将上述解带入方程组(1)中可以得到关于变量

ξ和

η的方程组

 

(3)dξdt=k1ξ−μϕ0η−μψ0ξ−μξηdηdt=νϕ0η+νψ0ξ−k2η+νξη

其中非线性项

μξη和

νξη为二阶小量,可以忽略;再将稳态解代入可得线性化的耦合方程组

解耦后可得到

 

(4)d2ξdt2+k1k2ξ=0d2ηdt2+k1k2η=0

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

其通解为

 

y(t)=Esin⁡(kt+δ)    或    y(t)=Ecos⁡(kt+δ)

其中,

E和

δ为振幅和初相位,与具体问题有关。

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

 

ϕ=k2ν+E1sin⁡(k1k2t+δ1)ψ=k1μ+E2sin⁡(k1k2t+δ2)

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

3. 数值求解

从上一节可知,我们需要数值求解一个耦合的常微分方程组,可以用RungeKutta法

[2]。简单推导过程如下:

其中,

 

f1(ϕ,ψ)=k1ϕ−μϕψf2(ϕ,ψ)=νϕψ−k2ψ

四阶Runge-Kutta方法可以表示为:

 

ϕk+1=ϕk+Δt6(f11+2f12+2f13+f14)ψk+1=ψk+Δt6(f21+2f22+2f23+f24)

其中,

 

fi1=fi(ϕk,ψk)fi2=fi(ϕk+Δt2f11,ψk+Δt2f21)fi3=fi(ϕk+Δt2f12,ψk+Δt2f22)fi4=fi(ϕk+Δtf11,ψk+Δtf21)    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(x,y(0),y(1),…,y(n−1))

其中,

n为阶数,

y(0)=y。

进一步,引入符号

D对各阶导数重新定义,此过程称为转换

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

 

Dj′=Dj+1Dn′=y(n)=f(x,D1,D2,…,Dn)

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

 

J=[∂D1′∂D1∂D1′∂D2⋯∂D1′∂Dn∂D2′∂D1∂D2′∂D2⋯∂D2′∂Dn⋮⋮⋱⋮∂Dn′∂D1∂Dn′∂D2⋯∂Dn′∂Dn]    和    ∂D1′∂x,∂D2′∂x,,…,∂Dn′∂x

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

让我们重新回顾一下公式(1),可知 nEqns() 应该返回 2;此外, 定义

Y=[ϕ,ψ]T,公式(1)可整理成如下向量形式

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

进而可以知道,

D1=Y,D1′=Y′,可得到雅可比矩阵和对自变量的偏导数分别为

 

∂D1′∂D1=∂Y′∂Y=[k1−μϕνψ−k2],    ∂D1′∂t=0

需要注意的是,雅可比矩阵只有一个元素

∂D1′∂D1,只不过这个元素是一个块的形式。

具体代码实现如下所示

#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 OpenFOAMOpenFOAM 编程 | One-Dimensional Transient Heat Conduction,此处不再赘述。

5. 数据分析

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

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

同时,解析解法(线性化的特殊解法)得到的结论是二者均按照

k1k2圆频率震荡,那么对应的周期为

T=2π/k1k2=2π/0.7∗0.5≈10.62,而数值解中得到的周期为 12.425,笔者认为在本文的条件假设下,其中的差距来自于线性解法中没有考虑非线性,但这个解法仍然具有实际意义;读者可以尝试改用绝对值较小的系数来降低其非线性程度。

另外,感兴趣的读者可以尝试使用 MatlabGNU Octave 求解该问题。

参考文献

[1] 顾樵. 数学物理方法[M]. 北京:科学出版社, 2012.
[2] Chenglin LI.数值计算(四十七)RungeKutta求解常微分方程组
[3] Hassan Kassem. How to solve ODE in OpenFOAM
[4] 捕食者与被捕食者模型——logistic-volterra


防止迷路,请关注笔者博客 博客园@Fitanium
喜欢的朋友还请点赞、收藏、转发,您的支持将是笔者创作的最大动力。

赞(0) 打赏
版权归原创作者所有,任何形式转载请联系作者;码农code之路 » OpenFOAM 编程 | 求解捕食者与被捕食者模型(predator-prey model)问题(ODEs)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

微信扫一扫打赏