\[\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 npk1 = 0.7k2 = 0.5mu = 0.1nu = 0.02def f1(phi,psi):return k1*phi-mu*phi*psidef f2(phi,psi):return nu*phi*psi-k2*psitStart = 0tEnd= 100.0n= 100000deltaT = tEnd / nhalfDeltaT = deltaT / 2.0Solution = 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] -> \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_;};