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

2022/11/7 1:24:02

本文主要是介绍OpenFOAM 编程 | 求解捕食者与被捕食者模型(predator-prey model)问题(ODEs),对大家解决编程问题具有一定的参考价值,需要的程序猿们随着小编来一起学习吧!

0. 写在前面

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

1. 问题描述

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

2. 解析求解

设任意 tt 时刻山兔与山猫的数量分别是 ϕϕ 和 ψψ ,二者的变化服从下面动力学方程

 

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

 

其中,k1k1,k2k2,μμ 和 νν 都是正常数。

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

  1. k1ϕk1ϕ 表示山兔种群的增长率,与山兔种群数量成正比。
  2. −μϕψ−μϕψ 表示山兔被山猫吃掉而导致的减少率,与乘积 ϕψϕψ (可表示两种动物的相遇概率)成正比。
  3. νϕψνϕψ 表示山猫种群的增长率,由于其数量增长取决于捕食(相遇才有可能),因此 νν 为正值。
  4. −k2ψ−k2ψ 表示山猫种群的死亡率,与其种群数量成正比。

方程组(1)因为含有乘积项,因此是非线性的。现采用线性化的特殊方法求解,即研究种群数量 ϕϕ 和 ψψ 在其稳定值附近的微小涨落。设方程组(1)的稳态解为 ϕ=ϕ0ϕ=ϕ0,ψ=ψ0ψ=ψ0,它们由下面条件决定

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

也就是

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

代数方程(2)的解为

ϕ0ψ0=k2ν=k1μϕ0=k2νψ0=k1μ

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

ϕψ=ϕ0+ξ=ψ0+ηϕ=ϕ0+ξψ=ψ0+η

其中,ξξ 和 ηη 与 ϕ0ϕ0 和 ψ0ψ0 相比都是小量。将上述解带入方程组(1)中可以得到关于变量 ξξ 和 ηη 的方程组

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

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

dξdtdηdt=−k2μνη=k1νμξdξdt=−k2μνηdηdt=k1νμξ

解耦后可得到

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

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

d2ydt2+k2y=0d2ydt2+k2y=0

其通解为

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

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

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

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

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

3. 数值求解

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

dϕdtdψdt=f1(ϕ,ψ)=f2(ϕ,ψ)dϕdt=f1(ϕ,ψ)dψdt=f2(ϕ,ψ)

其中

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

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

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

其中,

 

fi1fi2fi3fi4=fi(ϕk,ψk)=fi(ϕk+Δt2f11,ψk+Δt2f21)=fi(ϕk+Δt2f12,ψk+Δt2f22)=fi(ϕk+Δtf11,ψk+Δtf21)    i=1,2

这篇关于OpenFOAM 编程 | 求解捕食者与被捕食者模型(predator-prey model)问题(ODEs)的文章就介绍到这儿,希望我们推荐的文章对大家有所帮助,也希望大家多多支持为之网!


扫一扫关注最新编程教程