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