NLS-KDV方程的格子Boltzmann方法

Lattice Boltzmann Method for the NLS-KDV Equation

摘要: 建立一维格子Boltzmann模型的演化方程，运用Taylor展开和Chapman-Enskog多尺度分析技术，推导出能够恢复一类非线性耦合的NLS-KDV方程的平衡态分布函数和修正函数。最后，数值算例验证出该方法的计算结果与给出的精确解有很好的一致性。

Abstract: The evolution equations of the one-dimensional lattice Boltzmann model are established, and the equilibrium distribution function and the correction function that can recover a class of nonlinearly coupled NLS-KDV equations are derived by using Taylor expansion and Chapman-Enskog multiscale analysis techniques. Finally, numerical examples verify that the computational results of the method are in good agreement with the given exact solutions.

1. 引言

非线性Schrödinger方程是孤立光学方程的必要组成部分，有助于理解非线性光学、原子物理和量子场论等多种过程[1] [2]。KDV方程在流体动力学、等离子体物理学等领域发挥着重要作用，例如在等离子体物理学中产生了离子声学解，在地球物理流体动力学中表征了深海和潜水中的长波[3]。在实际物理系统中，不同波动模式会相互作用，因此将两个方程结合，得到耦合NLS-KDV方程。耦合NLS-KDV方程是非线性演化方程在等离子体物理学中的一个关键系统，控制着沿第三个非谐波介质的能量传递问题的长波极限[1]，用于在以离子声速度运动的坐标系统中模拟一维朗缪尔波和离子声波的非线性动力学[4]，描述了尘埃等离子体中的各种过程，如朗缪尔波、尘埃声波和电磁波的波动[5]。在数值计算方面，目前该方程的求解方法有有限元法[6] [7]、无网格法[8]、分解法[9]、同伦摄动法[10]和变分迭代法[11] [12]等。这些方法的迭代格式复杂，应用灵活性较差，不利于处理复杂的边界条件，也无法进行大规模并行计算。

格子Boltzmann方法(Lattice Boltzmann Method, LBM)已经广泛地被认为是描述流体运动与处理工程问题的有效手段[13] [14]。LBM已逐步应用于很多偏微分方程的求解，如长短波方程[15]，EFK方程[16]，STO方程[17]，Burgers方程[18]等。LBM具备易于并行计算、程序易于实施等优势。因此，笔者利用LBM求解一维NLS-KDV方程的初边值问题。

2. NLS-KDV方程的格子Boltzmann模型

在长波和短波之间的非线性相互作用中，基于弱非线性假设，两种典型的相互作用方程可以统一为具有如下标准形式的NLS-KDV方程：

$\{\begin{array}{l}i{u}_t+\varphi u_{xx}-\phi uv-\lambda {\left|u\right|}^{2}u=F\left(x,t\right),\\ v_t+\alpha v_{xxx}+{\left(\beta v^m+\gamma {\left|u\right|}^2\right)}_x=G\left(x,t\right).\end{array}$ (1)

其中：$\varphi ,\phi ,\lambda ,\alpha ,\beta$ 和$\gamma$ 为常数；$u\left(x,t\right)$ 是用来描述朗缪尔振荡电场的复值函数，$v\left(x,t\right)$ 是用来描述低频密度扰动的实值函数；$F\left(x,t\right)$ 和$G\left(x,t\right)$ 分别为已知的复函数和实函数。

对方程(1)采用D1Q5模型，构造格子Boltzmann演化方程为：

$\begin{array}{c}{f}_{si}\left(x+c_i\Delta t,t+\Delta t\right)-f_{si}\left(x,t\right)=-\frac{1}{{\tau }_s}\left[f_{si}\left(x,t\right)-f_{si}^{eq}\left(x,t\right)\right]\\ +\Delta t{h}_{si}\left(x,t\right)\text{+}\Delta t^2{\sigma }_{si}\left(x,t\right)s=1,2.\end{array}$ (2)

其中：$f_{si}\left(x,t\right)$ 为$\left(x,t\right)$ 处粒子的分布函数，$f_{si}^{eq}\left(x,t\right)$ 为$\left(x,t\right)$ 处粒子在平衡状态下的分布函数；$h_{si}\left(x,t\right)$ 和${\sigma }_{si}\left(x,t\right)$ 为修正函数；$\epsilon$ 为Kundsen (克努森)数，$c$ 为格子速度，$c_i=c\cdot e_i$ 表示不同方向的离散格子速度，${\tau }_s$ 为弛豫时间。

将(2)式左边利用Taylor公式展开，得到

$\begin{array}{l}\Delta t\left({\partial }_t+c_i\partial x\right)f_{si}+\frac{\Delta t^2}{2}{\left({\partial }_t+c_i\partial x\right)}^2{f}_{si}+\frac{\Delta t^3}{6}{\left({\partial }_t+c_i\partial x\right)}^3{f}_{si}+O\left(\Delta t^4\right)\\ =-\frac{1}{{\tau }_s}\left(f_{si}-f_{si}^{eq}\right)+\Delta t{h}_{si}\left(x,t\right)+\Delta t^2{\sigma }_{si}\left(x,t\right).\end{array}$ (3)

对(3)式进行Chapman-Enskog多尺度展开：

$\{\begin{array}{l}{\partial }_t={\partial }_{t_0}+\epsilon {\partial }_{t_1}+{\epsilon }^2{\partial }_{t_2}+O\left({\epsilon }^3\right),\\ f_{si}=f_{si}^{\left(0\right)}+\epsilon f_{si}^{\left(1\right)}+{\epsilon }^2{f}_{si}^{\left(2\right)}+{\epsilon }^3{f}_{si}^{\left(3\right)}+O\left({\epsilon }^4\right),\\ h_{si}={\epsilon }^2{h}_{si}^{\left(1\right)}+O\left({\epsilon }^3\right),\\ D={\partial }_{t_0}+c_i\partial x.\end{array}$ (4)

令$\Delta t=\epsilon$ ，将(4)式代入(3)式，并比较$\epsilon$ 同阶项的系数，可得

$O\left({\epsilon }^1\right):C_{1}D{f}_{si}^{\left(0\right)}=-\frac{1}{{\tau }_s}{f}_{si}^{\left(1\right)},$ (5)

$\begin{array}{l}O\left({\epsilon }^2\right):{\partial }_{t_1}{f}_{si}^{\left(0\right)}+D{f}_{si}^{\left(1\right)}+\frac{1}{2}{D}^2{f}_{si}^{\left(0\right)}=-\frac{1}{{\tau }_s}{f}_{si}^{\left(2\right)}+{\sigma }_{si}\left(x,t\right),\\ \end{array}$ (6)

$\begin{array}{l}O\left({\epsilon }^3\right):{\partial }_{t_2}{f}_{si}^{\left(0\right)}+{\partial }_{t_1}{f}_{si}^{\left(1\right)}+D{f}_{si}^{\left(2\right)}+\frac{1}{2}{D}^2{f}_{si}^{\left(1\right)}+D{\partial }_{t_1}{f}_{si}^{\left(0\right)}+\frac{1}{6}{D}^3{f}_{si}^{\left(0\right)}=-\frac{1}{{\tau }_s}{f}_{si}^{\left(3\right)}+h_{si}^{\left(1\right)},\\ \end{array}$ (7)

由(5)式和(6)式可以推导出

$O\left({\epsilon }^2\right):{\partial }_{t_1}{f}_{si}^{\left(0\right)}+C_2{D}^2{f}_{si}^{\left(0\right)}=-\frac{1}{{\tau }_s}{f}_{si}^{\left(2\right)}+{\sigma }_{si}.$ (8)

在(5)式和(8)式的基础上，用(7)式推导出

$\begin{array}{l}O\left({\epsilon }^3\right):{\partial }_{t_2}{f}_{si}^{\left(0\right)}+2C_{2}D{\partial }_{t_1}{f}_{si}^{\left(0\right)}+C_3{D}^3{f}_{si}^{\left(0\right)}+{\tau }_{s}D{\sigma }_{si}=-\frac{1}{{\tau }_s}{f}_{si}^{\left(3\right)}+h_{si}^{\left(1\right)}.\\ \end{array}$ (9)

对(5) (8) (9)式两边关于$i$ 求和，可得：

${\displaystyle \sum _{i}D{f}_{si}^{\left(0\right)}}=0.$ (10)

${\partial }_{t_1}{\displaystyle \sum _i{f}_{si}^{\left(0\right)}}+C_2{\displaystyle \sum _i{D}^2{f}_{si}^{\left(0\right)}}={\displaystyle \sum _i{\sigma }_{si}}.$ (11)

${\partial }_{t_2}{\displaystyle \sum _i{f}_{si}^{\left(0\right)}}+2C_2{\partial }_{t_1}{\displaystyle \sum _{i}D{f}_{si}^{\left(0\right)}}+C_3{\displaystyle \sum _i{D}^3{f}_{si}^{\left(0\right)}+}{\tau }_s{\displaystyle \sum _{i}D{\sigma }_{si}}={\displaystyle \sum h_{si}^{\left(1\right)}}.$ (12)

其中：$C_1=1;C_2=-{\tau }_s+\frac{1}{2};C_3={\tau }_s^2-{\tau }_s+\frac{1}{6}.$

2.1. LBM模型

在D1Q5模型的基础上，建立一维方程的格子玻尔兹曼模型，离散速度为

$e_i=\left(e_0,e_1,e_2,e_3,e_4\right)=\left(0,1,-1,2,-2\right).$

2.1.1. 恢复NLS方程

当$s=1$ 时，为了恢复NLS方程，选择满足

$\{\begin{array}{l}{\displaystyle \sum _i{f}_{1i}^{\left(0\right)}}=iu,\\ {\displaystyle \sum _i{c}_i{f}_{1i}^{\left(0\right)}}=0,\\ {\displaystyle \sum _i{c}_i^2{f}_{1i}^{\left(0\right)}}=\frac{\varphi u}{C_2\epsilon },\\ {\displaystyle \sum _i{c}_i^3{f}_{1i}^{\left(0\right)}}=0,\\ {\displaystyle \sum _i{c}_i^4{f}_{1i}^{\left(0\right)}}=0,\end{array}$ (13)

的平衡态分布函数$f_{1i}^{eq}\left(i=0,1,2,3,4\right)$ 和修正函数$h_{1i},{\sigma }_{1i},$

${\sigma }_{1i}=0,h_{1i}=\frac{1}{5}\left(F\left(x,t\right)+\phi uv+\lambda {\left|u\right|}^{2}u\right)i=0,1,2,3,4.$ (14)

将(13) (14)式代入(10)~(12)式，可得

$i{\partial }_{t_0}u=0,$ (15)

$i{\partial }_{t_1}u+\frac{\varphi }{\epsilon }{u}_{xx}=0,$ (16)

$i{\partial }_{t_2}u={\displaystyle \sum _i{h}_{1i}^{\left(1\right)}}.$ (17)

将$\left(15\right)式+\epsilon \times \left(16\right)式+{\epsilon }^2\times \left(17\right)式,$ 可以恢复出如下形式的NLS方程：

$i{u}_t+\varphi u_{xx}-\phi uv-\lambda {\left|u\right|}^{2}u=F\left(x,t\right).$

2.1.2. 恢复KDV方程

当$s=2$ 时，为了恢复KDV方程，选择满足

$\{\begin{array}{l}{\displaystyle \sum _i{f}_{2i}^{\left(0\right)}}=v,\\ {\displaystyle \sum _i{c}_i{f}_{2i}^{\left(0\right)}}=0,\\ {\displaystyle \sum _i{c}_i^2{f}_{2i}^{\left(0\right)}}=0,\\ {\displaystyle \sum _i{c}_i^3{f}_{2i}^{\left(0\right)}}=\frac{\alpha v}{{\epsilon }^2{C}_3},\\ {\displaystyle \sum _i{c}_i^4{f}_{2i}^{\left(0\right)}}=0,\end{array}$ (18)

$\{\begin{array}{l}{\displaystyle \sum _i{\sigma }_{2i}}=0,\\ {\displaystyle \sum _i{c}_i{\sigma }_{2i}}=\frac{\beta v^m+\gamma {\left|u\right|}^2}{{\tau }_2{\epsilon }^2},\\ {\displaystyle \sum _i{c}_i^2{\sigma }_{2i}}=0,\\ {\displaystyle \sum _i{c}_i^3{\sigma }_{2i}}=0,\\ {\displaystyle \sum _i{c}_i^4{\sigma }_{2i}}=0,\end{array}$ (19)

的平衡态分布函数$f_{2i}^{eq}\left(i=0,1,2,3,4\right)$ 和修正函数$h_{2i},{\sigma }_{2i},$

$h_{2i}=\frac{1}{5}G\left(x,t\right)i=0,1,2,3,4.$ (20)

将(18)~(20)式代入(10)~(12)式，可得

${\partial }_{t_0}v=0,$ (21)

${\partial }_{t_1}v=0,$ (22)

${\partial }_{t_2}v+\frac{\alpha }{{\epsilon }^2}{v}_{xxx}+\frac{1}{{\epsilon }^2}{\left(\beta v^m+\gamma {\left|u\right|}^2\right)}_x={\displaystyle \sum h_{2i}^{\left(1\right)}}.$ (23)

将$\left(21\right)式+\epsilon \times \left(22\right)式+{\epsilon }^2\times \left(23\right)式$ ，可以恢复出如下形式的KDV方程：

$v_t+\alpha v_{xxx}+{\left(\beta v^m+\gamma {\left|u\right|}^2\right)}_x=G\left(x,t\right).$

3. 数值算例

现给出一维NLS-KDV方程的两个数值算例，用来验证建立的格子Boltzmann模型对求解此方程的有效性。为了分析初边值问题LBM解和精确解之间的误差精度，引入全局相对误差(Global Relative Error，符号为${\delta }_{GRE}$ )和最大绝对误差(Maximum Absolute Error，符号为${\delta }_{MAE}$ )，定义如下：

$\begin{array}{l}{\delta }_{GRE}\left(u\right)=\frac{{\displaystyle \sum _k\left|u\left(x_k,t\right)-u^{*}\left(x_k,t\right)\right|}}{{\displaystyle \sum _k\left|u^{*}\left(x_k,t\right)\right|}},\\ {\delta }_{MAE}\left(u\right)=\underset{k}{\max }\left|u\left(x_k,t\right)-u^{*}\left(x_k,t\right)\right|,\end{array}$

其中：$u\left(x_k,t\right),u^{*}\left(x_k,t\right)$ 分别为$\left(x_k,t\right)$ 处的数值解和精确解。

3.1. 算例1

考虑方程(1)中参数$\varphi =0.1,\phi =1,\lambda =-2,m=1,\alpha =0.01,\beta =0.1,\gamma =1,$ 函数

$\{\begin{array}{l}F\left(x,t\right)=\left(0.4+t\sin 2x+2{\sin }^{2}x\right)\sin x{e}^{-\frac{it}{2}},\\ G\left(x,t\right)=-\sin 2x-0.12t\cos 2x+\sin 2x.\end{array}$

此时方程(1)有如下解析解：

$\{\begin{array}{l}u\left(x,t\right)=\sin x{e}^{-\frac{it}{2}},\\ v\left(x,t\right)=-t\sin 2x.\end{array}$

设置该例求解的区间范围为$\left[0,\pi \right]$ 。图1是在$N=120,\Delta t=0.00001,{\tau }_1=6.1,{\tau }_2=1,T=1$ 条件下，$\left|u\right|$ 和$v$ 的精确解和LBM解。由图1可知，LBM解与精确解吻合得较好，具有很好的一致性，说明模型具有较高的可靠性。图2分别给出了$\left|u\right|$ 和$v$ 在$T=1$ 时的LBM解的三维演化结果图以及精确解在时间和空间的分布图，可以看出数值解和精确解有很好的相似性。表1是$T=1$ 时不同网格数下$u$ 和$v$ 的误差，表2是$N=120$ 时不同时刻$u$ 和$v$ 的误差，可以看出全局相对误差和最大绝对误差都达到${10}^{-3}$ 精度，说明模型的求解精度较高。

(a) $\left|u\right|$ (b) $v$

Figure 1. Evolution of exact solutions and LBM numerical solutions

图1. 精确解和LBM数值解的演化

Table 1. Global relative error and maximum absolute error of u and v with different N when T = 1

表1. T = 1时不同网格数下u和v的全局相对误差和最大绝对误差

误差	N
	40	60	80	100	120
${\delta }_{GRE}\left(u\right)/{10}^{-3}$	67.946	27.716	15.773	9.6757	6.9809
${\delta }_{GRE}\left(v\right)/{10}^{-3}$	17.534	17.069	8.0779	6.3539	5.1506
${\delta }_{MAE}\left(u\right)/{10}^{-3}$	87.325	28.169	15.157	10.644	7.9288
${\delta }_{MAE}\left(v\right)/{10}^{-3}$	28.757	25.202	13.786	10.880	9.9770

Table 2. Global relative error and maximum absolute error of u and v with different T

表2. 不同时刻T下u和v的全局相对误差和最大绝对误差

误差	T
	0.2	0.4	0.6	0.8	1
${\delta }_{GRE}\left(u\right)/{10}^{-3}$	3.1279	5.2708	6.3663	7.2294	6.9809
${\delta }_{GRE}\left(v\right)/{10}^{-3}$	8.1168	7.0354	5.4325	4.6511	5.1506
${\delta }_{MAE}\left(u\right)/{10}^{-3}$	3.7756	5.2574	6.2732	6.7194	7.9288
${\delta }_{MAE}\left(v\right)/{10}^{-3}$	1.6345	3.9306	5.7468	7.8233	9.9770

3.2. 算例2

考虑方程(1)中参数$\varphi =-1,\phi =1,\lambda =-5,m=1,\alpha =0.1,\beta =1,\gamma =1,$ 函数

$\{\begin{array}{l}F\left(x,t\right)=\frac{\left(7+t\right)e^{-it}}{{\cosh }^{3}x},\\ G\left(x,t\right)=\frac{-{\cosh }^{3}x+\left(2.8t-2\right)\sinh x{\cosh }^{2}x-2.4t\sinh x}{{\cosh }^{5}x}.\end{array}$

此时方程(1)有如下解析解：

(a) (b)

(c) (d)

Figure 2. Evolution of exact solutions and LBM numerical solutions: (a) numerical solution of u; (b) exact solution of u; (c) numerical solution of v; (d) exact solution of v

图2. 精确解和LBM数值解的演化：(a) u的数值解；(b) u的精确解；(c) v的数值解；(d) v的精确解

(a) $\left|u\right|$ (b) $v$

Figure 3. Evolution of exact solutions and LBM numerical solutions

图3. 精确解和LBM数值解的演化

$\{\begin{array}{l}u\left(x,t\right)=\mathrm{sech}x{e}^{-it},\\ v\left(x,t\right)=-t{\mathrm{sech}}^{2}x.\end{array}$

设置该例求解的区间范围为$\left[-5,5\right]$ 。图3是在$N=180,\Delta t=0.00001,{\tau }_1=1.5,{\tau }_2=1,T=1$ 条件下，$\left|u\right|$ 和$v$ 的精确解和LBM解。由图3可知，LBM解与精确解吻合得较好，具有很好的一致性，说明模型具有较高的可靠性。表3是$T=1$ 时不同网格数下$u$ 和$v$ 的误差，表4是$N=180$ 时不同时刻$u$ 和$v$ 的误差，可以看出全局相对误差和最大绝对误差都达到${10}^{-3}$ 精度，说明模型的求解精度较高。

Table 3. Global relative error and maximum absolute error of u and v with different N when T = 1

表3. T = 1时不同网格数下u和v的全局相对误差和最大绝对误差

误差	N
	60	90	120	150	180
${\delta }_{GRE}\left(u\right)/{10}^{-3}$	50.049	23.905	14.502	9.9713	7.3721
${\delta }_{GRE}\left(v\right)/{10}^{-3}$	38.138	17.578	10.349	6.8545	5.2158
${\delta }_{MAE}\left(u\right)/{10}^{-3}$	26.647	12.817	8.6151	6.0742	4.5503
${\delta }_{MAE}\left(v\right)/{10}^{-3}$	45.782	21.559	11.564	7.4338	5.6430

Table 4. Global relative error and maximum absolute error of u and v with different T

表4. 不同时刻T下u和v的全局相对误差和最大绝对误差

误差	T
	0.2	0.4	0.6	0.8	1
${\delta }_{GRE}\left(u\right)/{10}^{-3}$	3.9149	5.3276	6.4918	7.4508	7.3721
${\delta }_{GRE}\left(v\right)/{10}^{-3}$	8.4048	6.2536	6.3615	5.4409	5.2158
${\delta }_{MAE}\left(u\right)/{10}^{-3}$	4.6479	3.2862	3.5933	4.0468	4.5503
${\delta }_{MAE}\left(v\right)/{10}^{-3}$	1.5953	2.7316	4.0137	4.6579	5.6430

4. 结语

笔者构建了双分布函数的一维格子Boltzmann模型，通过在模型中加入修正项函数来灵活地处理NLS-KDV方程中的非线性项和计算过程中产生的误差项，从而恢复出宏观方程。数值算例结果表明，数值解与精确解拟合较好，说明双分布函数的格子Boltzmann模型求解一维NLS-KDV方程的初边值问题是有效的。接下来，笔者将利用该格子Boltzmann模型求解高维情形的NLS-KDV方程，并将其推广到更高阶偏微分方程的数值求解中，以拓宽模型的应用范围。