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] -> \psidydx[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 算法（需要用到雅可比矩阵）对该问题进行求解，从下图可见二者求解得到的结果是一致的 。
经验总结扩展阅读
- Masked Label Prediction: Unified Message Passing Model for Semi-Supervised Classification
- DUCK 谣言检测《DUCK: Rumour Detection on Social Media by Modelling User and Comment Propagation Networks》
- 谣言检测《Data Fusion Oriented Graph Convolution Network Model for Rumor Detection》
- 【论文翻译】KLMo: Knowledge Graph Enhanced Pretrained Language Model with Fine-Grained Relationships
- 特斯拉modely和modelx有什么区别？
- 国产model3有氛围灯吗？
- modely和model3尺寸比较是什么？
- model3离地间隙是多少？
- 特斯拉modely底盘高度是多少？
- 特斯拉modely特点是什么？