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,μμ 和 νν 都是正常数。
在上述方程中有几点需要注意:
- k1ϕk1ϕ 表示山兔种群的净增长率,与山兔种群数量成正比。
- −μϕψ−μϕψ 表示山兔被山猫吃掉而导致的减少率,与乘积 ϕψϕψ (可表示两种动物的相遇概率)成正比。
- νϕψνϕψ 表示山猫种群的增长率,由于其数量增长取决于捕食(相遇才有可能),因此 νν 为正值。
- −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
原创文章,作者:ItWorker,如若转载,请注明出处:https://blog.ytso.com/292645.html