predator-prey model OpenFOAM 編程 | 求解捕食者與被捕食者模型問(wèn)題（ODEs） _生活百科

0. 寫在前面本文問(wèn)題參考自文獻(xiàn) \(^{[1]}\) 第一章例 6，并假設(shè)了一些條件，基于 OpenFOAM-v2206 編寫程序數(shù)值上求解該問(wèn)題。筆者之前也寫過(guò)基于 OpenFOAM 求解偏分方程的帖子，OpenFOAM 編程 | One-Dimensional Transient Heat Conduction 。
1. 問(wèn)題描述假設(shè)一群山貓（捕食者）和一群山兔（被捕食者）生活在同一片區(qū)域，那么我們可以知道，山貓吃了山兔，繁殖力會(huì)增強(qiáng)，山貓的數(shù)量會(huì)增加。這樣一來(lái)，山兔的數(shù)量會(huì)隨之減少。接下來(lái)，山貓由于食物短缺而數(shù)量減少，進(jìn)而導(dǎo)致山兔遇到山貓的機(jī)會(huì)減少（被吃掉的概率降低），結(jié)果山兔的數(shù)量又逐漸增加，這樣山貓得到食物的機(jī)會(huì)也隨之增加，其數(shù)量又再一次增加，而山兔的數(shù)量又會(huì)再一次隨之減少，如此不斷循環(huán) 。
2. 解析求解設(shè)任意 \(t\) 時(shí)刻山兔與山貓的數(shù)量分別是 \(\phi\) 和 \(\psi\) ，二者的變化服從下面動(dòng)力學(xué)方程
\[\begin{aligned}\frac{\mathrmss220yg\phi}{\mathrmemks2kkt} &= k_1 \phi - \mu\phi\psi \\\frac{\mathrmo8q4aky\psi}{\mathrmcusa20qt} &= \nu\phi\psi - k_2 \psi\end{aligned}\tag1\]其中，\(k_1\)，\(k_2\)，\(\mu\) 和 \(\nu\) 都是正常數(shù) 。
在上述方程中有幾點(diǎn)需要注意：

\(k_1\phi\) 表示山兔種群的凈增長(zhǎng)率，與山兔種群數(shù)量成正比。
\(-\mu\phi\psi\) 表示山兔被山貓吃掉而導(dǎo)致的減少率，與乘積 \(\phi\psi\) （可表示兩種動(dòng)物的相遇概率）成正比。
\(\nu\phi\psi\) 表示山貓種群的增長(zhǎng)率，由于其數(shù)量增長(zhǎng)取決于捕食（相遇才有可能），因此 \(\nu\) 為正值。
\(-k_2\psi\) 表示山貓種群的死亡率，與其種群數(shù)量成正比。

方程組（1）因?yàn)楹谐朔e項(xiàng)，因此是非線性的。現(xiàn)采用線性化的特殊方法求解，即研究種群數(shù)量 \(\phi\) 和 \(\psi\) 在其穩(wěn)定值附近的微小漲落。設(shè)方程組（1）的穩(wěn)態(tài)解為 \(\phi=\phi_0\)，\(\psi=\psi_0\)，它們由下面條件決定
\[\begin{aligned}\left . \frac{\mathrmesiiwmm\phi}{\mathrmooemcust} \right |_{\phi=\phi_0,\psi=\psi_0} &= 0 \\\left . \frac{\mathrmewmuka2\psi}{\mathrmy244u4et} \right |_{\phi=\phi_0,\psi=\psi_0} &=0\end{aligned}\]也就是
\[\begin{aligned}k_1 \phi_0 - \mu\phi_0\psi_0 &= 0 \\\nu\phi_0\psi_0 - k_2 \psi_0&=0\end{aligned}\tag2\]代數(shù)方程（2）的解為
\[\begin{aligned}\phi_0 &= \frac{k_2}{\nu} \\\psi_0 &=\frac{k_1}{\mu}\end{aligned}\]現(xiàn)在，將方程組（1）的解寫為下面形式
\[\begin{aligned}\phi &= \phi_0+ \xi \\\psi &= \psi_0 + \eta\end{aligned}\]其中，\(\xi\) 和 \(\eta\) 與 \(\phi_0\) 和 \(\psi_0\) 相比都是小量。將上述解帶入方程組（1）中可以得到關(guān)于變量 \(\xi\) 和 \(\eta\) 的方程組
\[\begin{aligned}\frac{\mathrmogowms2\xi}{\mathrmuc2csi4t} &= k_1\xi-\mu\phi_0\eta-\mu\psi_0\xi-\mu\xi\eta\\\frac{\mathrmewuksiq\eta}{\mathrmuussaq2t} &= \nu\phi_0\eta + \nu\psi_0\xi - k_2\eta+\nu\xi\eta\end{aligned}\tag3\]其中非線性項(xiàng) \(\mu\xi\eta\) 和 \(\nu\xi\eta\) 為二階小量，可以忽略；再將穩(wěn)態(tài)解代入可得線性化的耦合方程組
\[\begin{aligned}\frac{\mathrm02ggmsa\xi}{\mathrmi2m22got} &= -k_2\frac{\mu}{\nu}\eta\\\frac{\mathrmmeemck2\eta}{\mathrmokiiyg0t} &= k_1\frac{\nu}{\mu}\xi\end{aligned}\]解耦后可得到
\[\begin{aligned}\frac{\mathrmmg2ouki^2\xi}{\mathrmsom2gwet^2} +k_1k_2\xi&= 0\\\frac{\mathrmisaaqqm^2\eta}{\mathrm0u6asyot^2} +k_1k_2\eta&= 0\end{aligned}\tag4\]可以知道，式（4）與 L-C 震蕩電路及單擺問(wèn)題同屬于相同的數(shù)學(xué)模型
\[\frac{\mathrmsisyu06^2y}{\mathrmgmqyqu8t^2} + k^2 y = 0\]其通解為
\[y(t) = E\sin(kt+\delta)\ \ \ \ 或\ \ \ \ y(t) = E\cos(kt+\delta)\]其中，\(E\) 和 \(\delta\) 為振幅和初相位，與具體問(wèn)題有關(guān) 。
那么我們也可以得到本問(wèn)題的最終解的形式為
\[\begin{aligned}\phi &= \frac{k_2}{\nu} + E_1 \sin\left(\sqrt{k_1k_2}t+\delta_1\right)\\\psi &= \frac{k_1}{\mu} +E_2 \sin\left(\sqrt{k_1k_2}t+\delta_2\right) \\\end{aligned}\]其中，每個(gè)公式中振幅與初相位取決于各自的初始條件。