ICode9

精准搜索请尝试: 精确搜索
首页 > 编程语言> 文章详细

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

2022-11-06 13:11:58  阅读:276  来源: 互联网

标签:java 函数 学习 系统 语言 平台 方法 安装 QML c++ 数据


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

标签:java,函数,学习,系统,语言,平台,方法,安装,QML,c++,数据
来源:

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有