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

对应的,我们实现下主函数
#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 算法(需要用到雅可比矩阵)对该问题进行求解,从下图可见二者求解得到的结果是一致的 。

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

文章插图
同时运行笔者之前提到的 Python 代码后得到的数值结果与 OpenFOAM 计算结果绘制在同一张图中,二者高度重合 。
predator-prey model OpenFOAM 编程 | 求解捕食者与被捕食者模型问题(ODEs)

文章插图
同时,解析解法(线性化的特殊解法)得到的结论是二者均按照 \(\sqrt{k_1k_2}\) 圆频率震荡,那么对应的周期为 $T = 2\pi / \sqrt{k_1k_2} = 2 \pi / \sqrt{0.7*0.5} \approx 10.62 $,而数值解中得到的周期为 12.425,笔者认为在本文的条件假设下,其中的差距来自于线性解法中没有考虑非线性,但这个解法仍然具有实际意义 。
另外 , 感兴趣的读者可以尝试使用 MatlabGNU Octave 求解该问题 。
参考文献[1] 顾樵. 数学物理方法[M]. 北京:科学出版社, 2012.[2] Chenglin LI.数值计算(四十七)RungeKutta求解常微分方程组[3] Hassan Kassem. How to solve ODE in OpenFOAM[4] 捕食者与被捕食者模型——logistic-volterra
防止迷路,请关注笔者博客 博客园@Fiatanium 。喜欢的朋友还请点赞、收藏、转发,您的支持将是笔者创作的最大动力 。

推荐阅读