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 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}\]接下来，我们看一下如何实现相关求解代码 。首先看一下如何构造方程系统 。系统代码需要继承
经验总结扩展阅读
- 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特点是什么？