OpenFOAM 程式設計 | 求解捕食者與被捕食者模型（predator-prey model）問題（ODEs）

0. 寫在前面

本文問題參考自文獻 $^{[1]}$ 第一章例 6，並假設了一些條件，基於 OpenFOAM-v2206 編寫程式數值上求解該問題。筆者之前也寫過基於 OpenFOAM 求解偏分方程的貼文，OpenFOAM 程式設計 | One-Dimensional Transient Heat Conduction。

1. 問題描述

假設一群山貓（捕食者）和一群山兔（被捕食者）生活在同一片區域，那麼我們可以知道，山貓吃了山兔，繁殖力會增強，山貓的數量會增加。這樣一來，山兔的數量會隨之減少。接下來，山貓由於食物短缺而數量減少，進而導致山兔遇到山貓的機會減少（被吃掉的概率降低），結果山兔的數量又逐漸增加，這樣山貓得到食物的機會也隨之增加，其數量又再一次增加，而山兔的數量又會再一次隨之減少，如此不斷迴圈。

2. 解析求解

設任意 $t$ 時刻山兔與山貓的數量分別是 $\phi$ 和 $\psi$ ，二者的變化服從下面動力學方程

\[\begin{aligned} \frac{\mathrm{d}\phi}{\mathrm{d}t} &= k_1 \phi - \mu\phi\psi \\ \frac{\mathrm{d}\psi}{\mathrm{d}t} &= \nu\phi\psi - k_2 \psi \end{aligned} \tag1 \]

其中，$k_1$，$k_2$，$\mu$ 和 $\nu$ 都是正常數。

在上述方程中有幾點需要注意：

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

方程組（1）因為含有乘積項，因此是非線性的。現採用線性化的特殊方法求解，即研究種群數量 $\phi$ 和 $\psi$ 在其穩定值附近的微小漲落。設方程組（1）的穩態解為 $\phi=\phi_0$，$\psi=\psi_0$，它們由下面條件決定

\[\begin{aligned} \left . \frac{\mathrm{d}\phi}{\mathrm{d}t} \right |_{\phi=\phi_0,\psi=\psi_0} &= 0 \\ \left . \frac{\mathrm{d}\psi}{\mathrm{d}t} \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 \]

代數方程（2）的解為

\[\begin{aligned} \phi_0 &= \frac{k_2}{\nu} \\ \psi_0 &=\frac{k_1}{\mu} \end{aligned} \]

現在，將方程組（1）的解寫為下面形式

\[\begin{aligned} \phi &= \phi_0+ \xi \\ \psi &= \psi_0 + \eta \end{aligned} \]

其中，$\xi$ 和 $\eta$ 與 $\phi_0$ 和 $\psi_0$ 相比都是小量。將上述解帶入方程組（1）中可以得到關於變數 $\xi$ 和 $\eta$ 的方程組

\[\begin{aligned} \frac{\mathrm{d}\xi}{\mathrm{d}t} &= k_1\xi-\mu\phi_0\eta-\mu\psi_0\xi-\mu\xi\eta\\ \frac{\mathrm{d}\eta}{\mathrm{d}t} &= \nu\phi_0\eta + \nu\psi_0\xi - k_2\eta+\nu\xi\eta \end{aligned} \tag3 \]

其中非線性項 $\mu\xi\eta$ 和 $\nu\xi\eta$ 為二階小量，可以忽略；再將穩態解代入可得線性化的耦合方程組

\[\begin{aligned} \frac{\mathrm{d}\xi}{\mathrm{d}t} &= -k_2\frac{\mu}{\nu}\eta\\ \frac{\mathrm{d}\eta}{\mathrm{d}t} &= k_1\frac{\nu}{\mu}\xi \end{aligned} \]

解耦後可得到

\[\begin{aligned} \frac{\mathrm{d}^2\xi}{\mathrm{d}t^2} +k_1k_2\xi&= 0\\ \frac{\mathrm{d}^2\eta}{\mathrm{d}t^2} +k_1k_2\eta&= 0 \end{aligned} \tag4 \]

可以知道，式（4）與 L-C 震盪電路及單擺問題同屬於相同的數學模型

\[\frac{\mathrm{d}^2y}{\mathrm{d}t^2} + k^2 y = 0 \]

其通解為

\[y(t) = E\sin(kt+\delta)\ \ \ \ 或\ \ \ \ y(t) = E\cos(kt+\delta) \]

其中，$E$ 和 $\delta$ 為振幅和初相位，與具體問題有關。

那麼我們也可以得到本問題的最終解的形式為

\[\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} \]

其中，每個公式中振幅與初相位取決於各自的初始條件。

3. 數值求解

從上一節可知，我們需要數值求解一個耦合的常微分方程組，可以用RungeKutta法$^{[2]}$。簡單推導過程如下：

\[\begin{aligned} \frac{\mathrm{d}\phi}{\mathrm{d}t} &= f_1\left( \phi,\psi \right) \\ \frac{\mathrm{d}\psi}{\mathrm{d}t} &= f_2\left( \phi,\psi \right) \\ \end{aligned} \]

其中，

\[\begin{aligned} f_1\left( \phi,\psi \right) &= k_1 \phi - \mu\phi\psi \\ f_2\left( \phi,\psi \right) &= \nu\phi\psi - k_2 \psi \\ \end{aligned} \]

四階Runge-Kutta方法可以表示為：

\[\begin{aligned} \phi^{k+1} &= \phi^{k} + \frac{\Delta t}{6} \left( f_{11} + 2f_{12} + 2f_{13} + f_{14} \right) \\ \psi^{k+1} &= \psi^{k} + \frac{\Delta t}{6} \left( f_{21} + 2f_{22} + 2f_{23} + f_{24} \right) \\ \end{aligned} \]

其中，

\[\begin{aligned} f_{i1} &= f_i \left( \phi_k, \psi_k \right) \\ f_{i2} &= f_i \left( \phi_k+\frac{\Delta t}{2}f_{11}, \psi_k+\frac{\Delta t}{2}f_{21} \right) \\ f_{i3} &= f_i \left( \phi_k+\frac{\Delta t}{2}f_{12}, \psi_k+\frac{\Delta t}{2}f_{22} \right) \\ f_{i4} &= f_i \left( \phi_k+{\Delta t}f_{11}, \psi_k+{\Delta t}f_{21} \right) \\ \end{aligned} \ \ \ \ i=1,2 \]

求解程式碼採用 Python 編寫，如下所示

#!/usr/bin/python3
# -*- coding:utf-8 -*-

import numpy as np

k1 = 0.7
k2 = 0.5
mu = 0.1
nu = 0.02

def f1(phi,psi):
    return k1*phi-mu*phi*psi

def f2(phi,psi):
    return nu*phi*psi-k2*psi

tStart = 0
tEnd   = 100.0
n      = 100000
deltaT = tEnd / n
halfDeltaT = deltaT / 2.0
Solution = np.ndarray([n+1,2])
Solution[0] = [30,20] 

for i in range(n):
    f11 = f1(Solution[i][0], Solution[i][1])
    f21 = f2(Solution[i][0], Solution[i][1])

    f12 = f1(Solution[i][0] + halfDeltaT * f11, Solution[i][1] + halfDeltaT * f21)
    f22 = f2(Solution[i][0] + halfDeltaT * f11, Solution[i][1] + halfDeltaT * f21)

    f13 = f1(Solution[i][0] + halfDeltaT * f12, Solution[i][1] + halfDeltaT * f22)
    f23 = f2(Solution[i][0] + halfDeltaT * f12, Solution[i][1] + halfDeltaT * f22)

    f14 = f1(Solution[i][0] + deltaT * f11, Solution[i][1] + deltaT * f21)
    f24 = f2(Solution[i][0] + deltaT * f11, Solution[i][1] + deltaT * f21)

    Solution[i+1][0] = Solution[i][0] + deltaT / 6.0 * (f11 + 2*f12 + 2*f13 + f14)
    Solution[i+1][1] = Solution[i][1] + deltaT / 6.0 * (f21 + 2*f22 + 2*f23 + f24)
    print((i+1)*deltaT,Solution[i+1][0],Solution[i+1][1])

4. OpenFOAM 求解

使用OpenFOAM 數值求解常微分方程（組）主要用到 ODESystem.H（構造微分方程系統）和 ODESolver.H（求解器）；此外，在 OpenFOAM 中需要對常微分方程（組）進行整理$^{[3]}$，進而方便編寫程式碼進行求解。

對於任意階常微分方程可以轉化為一系列一階常微分方程，這個過程稱為降階，一階常微分方程的個數與原方程的階數相等（對於耦合常微分方程組，其階數等於所有方程階數之和）。對於某個 $n$ 階常微分方程，可按下面形式降階

\[y^{(n)}(x) = f \left( x, y^{(0)}, y^{(1)},\ldots,y^{(n-1)} \right) \]

其中，$n$ 為階數，$y^{(0)}=y$ 。

進一步，引入符號 $\mathrm{D}$ 對各階導數重新定義，此過程稱為轉換

\[\mathrm{D}_j = y^{(j-1)}\ \ \ \ j=1,2,\ldots,n-1 \]

最終，使用新符號重新表達原系統，此過程稱為誘導

\[\begin{aligned} \mathrm{D}'_j &= \mathrm{D}_{j+1} \\ \mathrm{D}'_n = y^{(n)} &= f\left( x, \mathrm{D}_1, \mathrm{D}_2,\ldots,\mathrm{D}_n \right) \end{aligned} \]

在 OpenFOAM 中，存在另外一個過程，該過程僅與剛性系統求解器相關，這類求解器需要雅可比矩陣和對自變數的偏導數，即

\[J = \begin{bmatrix} \frac{\partial \mathrm{D}'_1}{\partial \mathrm{D}_1} & \frac{\partial \mathrm{D}'_1}{\partial \mathrm{D}_2} & \cdots & \frac{\partial \mathrm{D}'_1}{\partial \mathrm{D}_n}\\ \frac{\partial \mathrm{D}'_2}{\partial \mathrm{D}_1} & \frac{\partial \mathrm{D}'_2}{\partial \mathrm{D}_2} & \cdots & \frac{\partial \mathrm{D}'_2}{\partial \mathrm{D}_n}\\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial \mathrm{D}'_n}{\partial \mathrm{D}_1} & \frac{\partial \mathrm{D}'_n}{\partial \mathrm{D}_2} & \cdots & \frac{\partial \mathrm{D}'_n}{\partial \mathrm{D}_n}\\ \end{bmatrix} \ \ \ \ 和 \ \ \ \ \frac{\partial \mathrm{D}'_1}{\partial x},\frac{\partial \mathrm{D}'_2}{\partial x}, ,\ldots, \frac{\partial \mathrm{D}'_n}{\partial x} \]

接下來，我們看一下如何實現相關求解程式碼。首先看一下如何構造方程系統。系統程式碼需要繼承 Foam::ODESystem 抽象類，並且需要全部實現三個方法nEqns()、 derivatives() 和 jacobian()，其中 jacobian() 方法對於非剛性求解器可以將實現置空（空函數體）。

讓我們重新回顧一下公式（1），可知 nEqns() 應該返回 2；此外，定義 $Y=[\phi,\psi]^{\mathrm{T}}$ ，公式（1）可整理成如下向量形式

\[\frac{\mathrm{d}Y}{\mathrm{d}t} = \begin{bmatrix} k_1 & -\mu\phi \\ \nu\psi & -k_2 \\ \end{bmatrix} Y \]

因此，導數可按照公式（1）編寫即可，只不過需要注意是向量形式。最後，對應之前的描述的降階過程，可以知道

\[Y' = f\left( t, Y\right) \]

進而可以知道， $D_1 = Y, D'_1=Y'$，可得到雅可比矩陣和對自變數的偏導數分別為

\[\frac{\partial \mathrm{D}'_1}{\partial \mathrm{D}_1} = \frac{\partial Y'}{\partial Y} = \begin{bmatrix} k_1 & -\mu\phi \\ \nu\psi & -k_2 \\ \end{bmatrix},\ \ \ \ \frac{\partial \mathrm{D}'_1}{\partial t} = 0 \]

需要注意的是，雅可比矩陣只有一個元素 $\frac{\partial \mathrm{D}'_1}{\partial \mathrm{D}_1}$，只不過這個元素是一個塊的形式。

具體程式碼實現如下所示

#include "ODESystem.H"

class ODEs : public Foam::ODESystem
{
public:
    ODEs() {}
    ~ODEs() {}
    // 初始化引數
    ODEs(const Foam::scalar k1, const Foam::scalar mu, const Foam::scalar k2,
         const Foam::scalar nu)
    {
        k1_ = k1;
        mu_ = mu;
        k2_ = k2;
        nu_ = nu;
    }
    // 方程個數
    Foam::label nEqns() const override { return 2; }
    // 求導
    void derivatives(const Foam::scalar x, const Foam::scalarField& y,
                     Foam::scalarField& dydx) const override
    { // 兩個未知量存成向量，y[0] -> \phi, y[1] -> \psi
        dydx[0] = k1_ * y[0] - mu_ * y[0] * y[1];
        dydx[1] = nu_ * y[0] * y[1] - k2_ * y[1];
    }
    // 計算符號的雅可比矩陣和關於自變數的導數
    void jacobian(const Foam::scalar x, const Foam::scalarField& y, Foam::scalarField& dfdx,
                  Foam::scalarSquareMatrix& dfdy) const override
    {
        dfdx[0] = 0;
        dfdx[1] = 0;

        dfdy[0][0] = k1_;
        dfdy[0][1] = -mu_ * y[0];

        dfdy[1][0] = nu_ * y[1];
        dfdy[1][1] = -k2_;
    }

private:
    Foam::scalar k1_;
    Foam::scalar mu_;
    Foam::scalar k2_;
    Foam::scalar nu_;
};

對應的，我們實現下主函數

#include <iostream>
#include <memory>

#include "ODESystem.H"
#include "ODESolver.H"

class ODEs : public Foam::ODESystem
{
    // 這裡的程式碼在上邊已經介紹，此處省略
};

int main(int argc, char* argv[])
{
    const Foam::scalar startTime = 0.0;         // 開始時間
    const Foam::scalar endTime   = 100.0;       // 結束時間
    const Foam::scalar phi0      = 30;          // 山兔初始值
    const Foam::scalar psi0      = 20;          // 山貓初始值
    const Foam::label n          = 100000;      //
    const Foam::scalar deltaT    = endTime / n; // 步長
    // 係數，參考自文獻[4]
    const Foam::scalar k1 = 0.7;
    const Foam::scalar mu = 0.1;
    const Foam::scalar k2 = 0.5;
    const Foam::scalar nu = 0.02;
    // 構造物件
    ODEs odes(k1, mu, k2, nu);

    // 構造求解器，具體使用的演演算法通過引數傳遞
    Foam::dictionary dict;
    dict.add("solver", argv[1]);
    Foam::autoPtr<Foam::ODESolver> solver = Foam::ODESolver::New(odes, dict);

    // 初始化一些變數
    Foam::scalar tStart = startTime;
    Foam::scalarField PhiPsi(odes.nEqns()); // 因變數
    PhiPsi[0] = phi0;
    PhiPsi[1] = psi0;
    Foam::scalarField ddt(odes.nEqns()); // 儲存導數值

    // 計算過程
    for (Foam::label i = 0; i < n; ++i)
    {
        Foam::scalar dtEst = deltaT / 2;
        Foam::scalar tEnd  = tStart + deltaT;
        //
        odes.derivatives(tStart, PhiPsi, ddt);
        solver->solve(tStart, tEnd, PhiPsi, dtEst);
        //
        tStart = tEnd;
        //
        Foam::Info << tStart << "," << PhiPsi[0] << "," << PhiPsi[1] << Foam::endl;
    }

    return 0;
}

此外，CMakeLists.txt 檔案可參考筆者之前的隨筆，如 OpenFOAM程式設計 | Hello OpenFOAM 和 OpenFOAM 程式設計 | One-Dimensional Transient Heat Conduction，此處不再贅述。

5. 資料分析

筆者通過命令列引數分別採用RKCK45 演演算法和 seulex 演演算法（需要用到雅可比矩陣）對該問題進行求解，從下圖可見二者求解得到的結果是一致的。

同時執行筆者之前提到的 Python 程式碼後得到的數值結果與 OpenFOAM 計算結果繪製在同一張圖中，二者高度重合。

同時，解析解法（線性化的特殊解法）得到的結論是二者均按照 $\sqrt{k_1k_2}$ 圓頻率震盪，那麼對應的週期為 $T = 2\pi / \sqrt{k_1k_2} = 2 \pi / \sqrt{0.7*0.5} \approx 10.62 $，而數值解中得到的週期為 12.425，筆者認為在本文的條件假設下，其中的差距來自於線性解法中沒有考慮非線性，但這個解法仍然具有實際意義；讀者可以嘗試改用絕對值較小的係數來降低其非線性程度。

另外，感興趣的讀者可以嘗試使用 Matlab 或 GNU Octave 求解該問題。

參考文獻

[1] 顧樵. 數學物理方法[M]. 北京：科學出版社, 2012.
[2] Chenglin LI.數值計算（四十七）RungeKutta求解常微分方程組
[3] Hassan Kassem. How to solve ODE in OpenFOAM
[4] 捕食者與被捕食者模型——logistic-volterra

防止迷路，請關注筆者部落格部落格園@Fiatanium。
喜歡的朋友還請點贊、收藏、轉發，您的支援將是筆者創作的最大動力。