非线性RLW-EKdV方程的高精度紧致有限差分算法

High-precision Compact Finite Difference Algorithm for Nonlinear RLW-EKdV Equation

摘要

为了利用非线性RLW-EKdV方程的初边值问题进行数值模拟研究，首先构造高精度守恒型紧致有限差分格式，对时间方向的导数采用隐中点格式进行离散，对空间方向的各阶导数用逆紧致算子进行离散，使其在时间方向上的精度达到二阶，在空间方向上的精度达到六阶。进一步，对差分格式进行守恒性证明、先验估计及收敛性分析。最后，通过数值算例验证理论的正确性和格式的有效性与可靠性。

Abstract

To numerically simulate the initial-boundary value problems of the nonlinear RLW-EKdV equation， a high-precision conservative compact finite difference scheme was first constructed. The time derivative was discretized using the implicit midpoint scheme， and the spatial derivatives of various orders were discretized using the inverse compact operator， achieving second-order accuracy in time and sixth-order accuracy in space. Furthermore， the conservation proof， a priori estimates， and convergence analysis of the difference scheme were conducted. Finally， numerical examples were used to verify the correctness of the theory and the effectiveness and reliability of the scheme.

Graphical abstract

关键词

RLW-EKdV方程 / 紧致差分格式 / 守恒性 / 收敛性 / 误差估计

Key words

RLW-EKdV equation / compact difference scheme / conservation / convergence / error estimation

水波是自由面或层结流体在重力与浮力共同作用下产生的波动，是一种常见的物理现象，小至微风轻拂水面而起的细微涟漪，大至海洋内部振幅达数十米、绵延数百公里的内孤立波，均属该范畴，水波也因此成为流体力学及工程海洋学等领域的核心研究对象。1966年，Peregrine^［1］在研究水波问题时提出非线性正则长波（RLW）方程（也称为BBM方程）

u_t+u_x+u{u}_x-u_{xxt}=\mathrm{0}，（1）

该方程能够描述浅水域中弱非线性孤立波的演化过程。在此基础上发展起来的非线性广义BBM-KdV（g-BBM-KdV）方程

u_t+\alpha u_x+\sigma u^p{u}_x+\gamma u_{xxx}+\epsilon u_{xxt}=\mathrm{0}，（2）

可作为描述包含高阶非线性作用的小振幅长波模型^［2］，通常情况下，取p=1或p=2。

目前求解非线性KdV型方程及RLW型方程的数值方法主要包括物理信息神经网络（PINN）方法、有限元法及（紧致）有限差分法，研究成果丰富^［3⁃11］。紧致有限差分格式计算量较小，仅用三到五点模板即可达到四到六阶精度，适合求解数值耗散误差敏感、求解区域规则的问题。Hu等^［12］对RLW方程提出来一个结合Richardson外推法的三层线性化隐式有限差分方案，其时间精度为二阶，空间精度为四阶。Yang等^［13］分别采用四阶紧致差分格式和padé四阶紧致差分格式离散RLW方程，并采用托马斯（Thomas）算法对离散后的线性方程组进行求解。邓雅清等^［14］对RLW-KdV方程提出一种在空间上四阶收敛且在时间上二阶收敛的紧致有限差分格式，用离散能量法证明了差分格式的守恒性。Polwang等^［15］对KdV-BBM方程提出了3种高效高阶线性数值方案，其时间收敛阶均为二阶，空间收敛阶均为四阶，同时保持质量和能量守恒。

本文研究RLW-EKdV方程的初边值问题

\begin{array}{l}{u}_t+\alpha u_x+\beta u{u}_x+\sigma u^{\mathrm{2}}{u}_x+\gamma u_{xxx}+\epsilon u_{xxt}=\mathrm{0},x\in (a,b),t\in (\mathrm{0},T],\\ u(a,t)=u(b,t),u_x(a,t)=u_x(b,t),u_{xx}(a,t)=u_{xx}(b,t),t\in (\mathrm{0},T],\\ u(x,\mathrm{0})=u_{\mathrm{0}}(x),x\in [a,b]\mathrm{。}\end{array}（3）

此非线性初边值问题的解析解比较困难，用数值方法研究该问题的解显得尤为重要。

如果没有特殊说明，本文对求解区域x\in [a,b]，t\in [\mathrm{0},T]进行网格剖分，取时间步长\tau =T/N，空间步长h=(b-a)/J，其中J和N为正整数。记网格点x_j=a+jh，(\mathrm{0}\le j\le J)，t_n=n\tau，(\mathrm{0}\le n\le N)。用u_j^n表示u(x,t)在网格点(x_j,t_n)处的近似值，I_h=\left\{x_{\mathrm{0}},x_{\mathrm{1}},\dots ,x_J\right\}表示区间I上的节点集合，Z_h表示在离散网格点上定义的函数集合，U_j^n表示网格点(x_j,t_n)处的精确值，H_E^{\mathrm{1}}表示定义在I_h上的一组网格函数u满足周期边界条件。约定M、M₁、M₂、C₁及C₂为一般正常数，即在不同处有不同的取值。

1 记号和引理

为了方便后续内容的叙述，引入记号

\begin{array}{l}{\delta }_t{u}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}=\frac{u_j^{n+\mathrm{1}}-u_j^n}{\tau },{\delta }_x{u}_j^n=\frac{u_{j+\mathrm{1}}^n-u_j^n}{\mathrm{\hslash }},A_t{u}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}=\frac{u_j^{n+\mathrm{1}}+u_j^n}{\mathrm{2}},\\ {\delta }_x^{\mathrm{1}}{u}_j^n=\frac{u_{j+\mathrm{2}}^n+\mathrm{28}{u}_{j+\mathrm{1}}^n-\mathrm{28}{u}_{j-\mathrm{1}}^n-u_{j-\mathrm{2}}^n}{\mathrm{36}h},{\delta }_x^{\mathrm{2}}{u}_j^n=\frac{\mathrm{3}{u}_{j+\mathrm{2}}^n+\mathrm{48}{u}_{j+\mathrm{1}}^n-\mathrm{102}{u}_j^n+\mathrm{48}{u}_{j-\mathrm{1}}^n+\mathrm{3}{u}_{j-\mathrm{2}}^n}{\mathrm{44}{h}^{\mathrm{2}}},\\ {\delta }_x^{\mathrm{3}}{u}_j^n=\frac{-u_{j+\mathrm{3}}^n+\mathrm{64}{u}_{j+\mathrm{2}}^n-\mathrm{125}{u}_{j+\mathrm{1}}^n+\mathrm{125}{u}_{j-\mathrm{1}}^n-\mathrm{64}{u}_{j-\mathrm{2}}^n+u_{j-\mathrm{3}}^n}{\mathrm{64}{h}^{\mathrm{3}}}\mathrm{。}\end{array}

定义算子P_x、Q_x和S_x

P_x{u}_j^n=\frac{u_{j+\mathrm{1}}^n+\mathrm{3}{u}_j^n+u_{j-\mathrm{1}}^n}{\mathrm{3}},Q_x{u}_j^n=\frac{\mathrm{2}{u}_{j+\mathrm{1}}^n+\mathrm{11}{u}_j^n+\mathrm{2}{u}_{j-\mathrm{1}}^n}{\mathrm{11}},S_x{u}_j^n=\frac{\mathrm{7}{u}_{j+\mathrm{1}}^n+\mathrm{16}{u}_j^n+\mathrm{7}{u}_{j-\mathrm{1}}^n}{\mathrm{16}},（4）

满足^［16］

{\mathit{P}}_x{\left(\frac{\partial U}{\partial x}\right)}_j^n={\delta }_x^{\mathrm{1}}{U}_j^n+O(h^{\mathrm{6}}),{\mathit{Q}}_x{\left(\frac{{\partial }^{\mathrm{2}}U}{\partial x^{\mathrm{2}}}\right)}_j^n={\delta }_x^{\mathrm{2}}{U}_j^n+O(h^{\mathrm{6}}),{\mathit{S}}_x{\left(\frac{{\partial }^{\mathrm{3}}U}{\partial x^{\mathrm{3}}}\right)}_j^n={\delta }_x^{\mathrm{3}}{U}_j^n+O(h^{\mathrm{6}})\mathrm{。}

由式（4）可知{\mathit{P}}_x、{\mathit{Q}}_x及{\mathit{S}}_x是对角占优的，所以可逆，从而有

{\left(\frac{\partial U}{\partial x}\right)}_j^n={\mathit{P}}_x^{-\mathrm{1}}{\delta }_x^{\mathrm{1}}{U}_j^n+O(h^{\mathrm{6}}),{\left(\frac{{\partial }^{\mathrm{2}}U}{\partial x^{\mathrm{2}}}\right)}_j^n={\mathit{Q}}_x^{-\mathrm{1}}{\delta }_x^{\mathrm{2}}{U}_j^n+O(h^{\mathrm{6}}),{\left(\frac{{\partial }^{\mathrm{3}}U}{\partial x^{\mathrm{3}}}\right)}_j^n={\mathit{S}}_x^{-\mathrm{1}}{\delta }_x^{\mathrm{3}}{U}_j^n+O(h^{\mathrm{6}})\mathrm{。}

U关于x的各阶导数可表示为如下的矩阵形式

{\left(\frac{\partial U}{\partial x}\right)}_j^n={\mathit{D}}_{\mathrm{1}}{U}^n+O(h^{\mathrm{6}}),{\left(\frac{{\partial }^{\mathrm{2}}U}{\partial x^{\mathrm{2}}}\right)}_j^n={\mathit{D}}_{\mathrm{2}}{U}^n+O(h^{\mathrm{6}}),{\left(\frac{{\partial }^{\mathrm{3}}U}{\partial x^{\mathrm{3}}}\right)}_j^n={\mathit{D}}_{\mathrm{3}}{U}^n+O(h^{\mathrm{6}})\mathrm{。}

其中

U^n={\left(U_{\mathrm{0}}^n,U_{\mathrm{1}}^n,\dots ,U_J^n\right)}^{\mathrm{T}},{\mathit{D}}_{\mathrm{1}}=P_x^{-\mathrm{1}}{\delta }_x^{\mathrm{1}},{\mathit{D}}_{\mathrm{2}}=Q_x^{-\mathrm{1}}{\delta }_x^{\mathrm{2}},{\mathit{D}}_{\mathrm{3}}=S_x^{-\mathrm{1}}{\delta }_x^{\mathrm{3}}\mathrm{。}

下面给出几个定义和引理。

定义1 对于任意的两个网格函数f^n,g^n\in Z_h，定义如下的离散内积和范数：

(f^n,g^n)=h{\displaystyle \sum _{j=\mathrm{1}}^{J-\mathrm{1}}}{f}_j^n{g}_j^n,||f^n||=\sqrt[]{(f^n,f^n)},||f^n||{}_{\mathrm{\infty }}=\underset{\mathrm{1}\le j\le J-\mathrm{1}}{\mathrm{m}\mathrm{a}\mathrm{x}}|f_j^n|\mathrm{。}

经过简单推导很容易得到下列结论。

引理1^［17］ 对任意的两个网格函数u^n，{\rho }^n\in Z_h，有\left({\mathit{P}}_x^{-\mathrm{1}}{\delta }_x^{\mathrm{1}}{u}^n,{\rho }^n\right)=-\left(u^n,{\mathit{P}}_x^{-\mathrm{1}}{\delta }_x^{\mathrm{1}}{\rho }^n\right), 特别地，有

\left({\mathit{P}}_x^{-\mathrm{1}}{\delta }_x^{\mathrm{1}}{u}^n,u^n\right)=\mathrm{0}\mathrm{。}

引理2^［17］ 对任意的两个网格函数u^n，{\rho }^n\in Z_h，有

\left({\mathit{Q}}_x^{-\mathrm{1}}{\delta }_x^{\mathrm{2}}\left(u^n-{\rho }^n\right),u^n+{\rho }^n\right)=-\left({‖{\mathit{P}}_x^{-\mathrm{1}}{\delta }_x^{\mathrm{1}}{u}^n‖}^{\mathrm{2}}-{‖{\mathit{P}}_x^{-\mathrm{1}}{\delta }_x^{\mathrm{1}}{\rho }^n‖}^{\mathrm{2}}\right)\mathrm{。}

引理3^［17］ 对任意的网格函数u^n\in Z_h，存在两个正常数C_{\mathrm{1}}和C_{\mathrm{2}}，使得

{‖u^n‖}_{\mathrm{\infty }}\le C_{\mathrm{1}}‖u^n‖+C_{\mathrm{2}}‖{\delta }_x{u}^n‖。

引理4^［17］ 对任意网格函数u^n\in Z_h，有

{\displaystyle \sum _{j=\mathrm{1}}^J}{\mathit{P}}_x{u}_j={\displaystyle \sum _{j=\mathrm{1}}^J}{u}_j,{\displaystyle \sum _{j=\mathrm{1}}^J}{\mathit{S}}_x{u}_j={\displaystyle \sum _{j=\mathrm{1}}^J}{u}_j,{\displaystyle \sum _{j=\mathrm{1}}^J}{\mathit{Q}}_x{u}_j={\displaystyle \sum _{j=\mathrm{1}}^J}{u}_j,

相似地，{{\mathit{P}}_x}^{-\mathrm{1}},{{\mathit{S}}_x}^{-\mathrm{1}},{{\mathit{Q}}_x}^{-\mathrm{1}}也有以下性质：

{\displaystyle \sum _{j=\mathrm{1}}^J}{\mathit{P}}_x^{-\mathrm{1}}{u}_j={\displaystyle \sum _{j=\mathrm{1}}^J}{u}_j,{\displaystyle \sum _{j=\mathrm{1}}^J}{\mathit{S}}_x^{-\mathrm{1}}{u}_j={\displaystyle \sum _{j=\mathrm{1}}^J}{u}_j,{\displaystyle \sum _{j=\mathrm{1}}^J}{\mathit{Q}}_x^{-\mathrm{1}}{u}_j={\displaystyle \sum _{j=\mathrm{1}}^J}{u}_j\mathrm{。}

引理5^［17］ 对于矩阵{\mathit{P}}_x^{-\mathrm{1}}有{\mathit{P}}_x^{-\mathrm{1}}={\mathit{R}}_{\mathrm{1}}{\mathit{R}}_{\mathrm{1}}^{\mathrm{T}}，且C_{\mathrm{1}}‖u^n‖\le ‖R_{\mathrm{1}}{u}^n‖\le C_{\mathrm{2}}‖u^n‖，其中R_{\mathrm{1}}=\mathrm{C}\mathrm{h}\mathrm{o}\mathrm{l}(P_x^{-\mathrm{1}})，C_{\mathrm{1}}，C_{\mathrm{2}}是正常数，\mathrm{C}\mathrm{h}\mathrm{o}\mathrm{l}(\cdot )表示\mathrm{C}\mathrm{h}\mathrm{o}\mathrm{l}\mathrm{e}\mathrm{s}\mathrm{k}\mathrm{y}分解。

2 差分格式的构造

2.1 格式构造

对问题（3）在网格点\left(x_j,t_{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)处进行离散，在时间方向用隐中点格式离散，在空间方向用六阶紧致差分格式离散，得差分格式：

（5）

2.2 守恒律

下面给出方程满足的能量守恒律和质量守恒律。

定理1 设U(x,t)为问题（3）的精确解，记能量

E(t)={‖U‖}^{\mathrm{2}}-\epsilon {‖U_x‖}^{\mathrm{2}}，（6）

则有E(t)=E(\mathrm{0})，

证明 用u(x,\mathrm{t})与问题（3）的控制方程两边作内积，得

{\int }_a^{b}u\left(x,t\right)u_t\left(x,t\right)\mathrm{d}x+\alpha {\int }_a^{b}u\left(x,t\right)u_x\left(x,t\right)\mathrm{d}x+\beta {\int }_a^b{u}^{\mathrm{2}}\left(x,t\right)u_x\left(x,t\right)\mathrm{d}x+\sigma {\int }_a^b{u}^{\mathrm{3}}\left(x,t\right)u_x\left(x,t\right)\mathrm{d}x+\gamma {\int }_a^{b}u\left(x,t\right)u_{xxx}\left(x,t\right)\mathrm{d}x+\epsilon {\int }_a^{b}u\left(x,t\right)u_{xxt}\left(x,t\right)\mathrm{d}x=\mathrm{0},（7）

经简单计算得

\frac{\mathrm{d}}{\mathrm{d}t}\left[{\int }_a^b{u}^{\mathrm{2}}\left(x,t\right)\mathrm{d}x-\epsilon {{\int }_a^b\left(u_x\right)}^{\mathrm{2}}\mathrm{d}x\right]=\mathrm{0}。（8）

即

E(t)=E(\mathrm{0})，t\in [\mathrm{0},T]。（9）

该定理通常称为方程的能量守恒律，下面证明质量守恒律。

定理2 设u(x,t)为问题（3）的精确解，记质量

Q(t)={\int }_a^{b}u\left(x,t\right)\mathrm{d}x,

则有Q(t)=Q(\mathrm{0})，。

证明 对问题（3）的控制方程的每一项进行积分，即

\begin{array}{l}{\int }_a^b{u}_t\left(x,t\right)\mathrm{d}x+\alpha {\int }_a^b{u}_x\left(x,t\right)\mathrm{d}x+\beta {\int }_a^{b}u\left(x,t\right)u_x\left(x,t\right)\mathrm{d}x+\\ \sigma {\int }_a^b{u}^{\mathrm{2}}\left(x,t\right)u_x\left(x,t\right)\mathrm{d}x+\gamma {\int }_a^b{u}_{xxx}\left(x,t\right)\mathrm{d}x+\epsilon {\int }_a^b{u}_{xxt}\left(x,t\right)\mathrm{d}x=\mathrm{0},\end{array}

经过简单计算可知\frac{\mathrm{d}}{\mathrm{d}t}{\int }_a^{b}u\left(x,t\right)\mathrm{d}x=\mathrm{0}。即Q(t)=Q(\mathrm{0})，t\in [\mathrm{0},T]。

下面给出离散格式满足的能量守恒律和质量守恒律。

定理3 格式（5）满足离散的能量和质量守恒律，定义能量和质量分别为

E^n={‖u^n‖}^{\mathrm{2}}-\epsilon {‖{\mathit{D}}_{\mathrm{1}}{u}^n‖}^{\mathrm{2}},Q^n=h{\displaystyle \sum _{j=\mathrm{0}}^J}{u}_j^n。（10）

对于任意的n\ge \mathrm{0}，有E^n=E^{n+\mathrm{1}},Q^n=Q^{n+\mathrm{1}}\mathrm{。}

证明 首先证明E^n=E^{n+\mathrm{1}}。用A_t{u}^{n+\frac{\mathrm{1}}{\mathrm{2}}}与式（5）的控制方程作内积，得

\frac{\mathrm{1}}{\tau }\left({‖u^{n+\mathrm{1}}‖}^{\mathrm{2}}-{‖u^n‖}^{\mathrm{2}}\right)+\alpha \left({\left({\mathit{D}}_{\mathrm{1}}{A}_t{u}^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)}_j,A_t{u}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)+\frac{\mathrm{1}}{\mathrm{3}}\sigma \left({\left({\mathit{D}}_{\mathrm{1}}{\left(A_t{u}^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)}^{\mathrm{3}}\right)}_j,A_t{u}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)+\frac{\mathrm{1}}{\mathrm{2}}\beta \left({\left({\mathit{D}}_{\mathrm{1}}{\left(A_t{u}^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)}^{\mathrm{2}}\right)}_j,A_t{u}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)+\left(\gamma {\left({\mathit{D}}_{\mathrm{3}}\left(A_t{u}^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)\right)}_j,A_t{u}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)+\epsilon \left({\delta }_t{\left({\mathit{D}}_{\mathrm{2}}{u}^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)}_j,A_t{u}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)=\mathrm{0},（11）

整理上式得{‖u^{n+\mathrm{1}}‖}^{\mathrm{2}}-\epsilon {‖{\mathit{D}}_{\mathrm{1}}{u}^{n+\mathrm{1}}‖}^{\mathrm{2}}={‖u^n‖}^{\mathrm{2}}-\epsilon {‖{\mathit{D}}_{\mathrm{1}}{u}^n‖}^{\mathrm{2}}。由E^n的定义可知，E^{n+\mathrm{1}}=E^n。

接下来证明Q^n=Q^{n+\mathrm{1}}\mathrm{。}由式（5）的控制方程可得

h{\displaystyle \sum _{j=\mathrm{0}}^J}\frac{u_j^{n+\mathrm{1}}-u_j^n}{\tau }+\alpha h{\displaystyle \sum _{j=\mathrm{0}}^J}{A}_t{P}_x^{-}{\delta }_x^{\mathrm{1}}{u}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}+\frac{\mathrm{1}}{\mathrm{3}}\sigma h{\displaystyle \sum _{j=\mathrm{0}}^J}{A}_t{P}_x^{-}{\delta }_x^{\mathrm{1}}{\left(u_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)}^{\mathrm{3}}+\frac{\mathrm{1}}{\mathrm{2}}\beta h{\displaystyle \sum _{j=\mathrm{0}}^J}{A}_t{P}_x^{-}{\delta }_x^{\mathrm{1}}{\left(u_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)}^{\mathrm{2}}+\gamma h{\displaystyle \sum _{j=\mathrm{0}}^J}{A}_t{S}_x^{-}{\delta }_x^{\mathrm{3}}{\left(u_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)}^{\mathrm{2}}+\epsilon h{\displaystyle \sum _{j=\mathrm{0}}^J}{\delta }_t{Q}_x^{-}{\delta }_x^{\mathrm{2}}{u}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}=\mathrm{0}。（12）

再由边界条件及引理4有h{\displaystyle \sum _{j=\mathrm{0}}^J}{u}_j^n=h{\displaystyle \sum _{j=\mathrm{0}}^J}{u}_j^{n+\mathrm{1}}。由Q^n的定义可知Q^n=Q^{n+\mathrm{1}}\mathrm{。}

3 先验估计与收敛分析

利用先验估计进行收敛性分析。

定理4 假设U_{\mathrm{0}}\in H_E^{\mathrm{1}}\left[a,b\right]，则问题 （3） 的精确解U满足‖U‖\le M,‖U_x‖\le M\mathrm{。}

证明 结合定理1经简单推导可得。

定理5 假设u_{\mathrm{0}}\in H_E^{\mathrm{1}}\left[a,b\right]，则格式（5）的数值解u^n满足‖u^n‖\le M,‖{\delta }_x^{\mathrm{1}}{u}^n‖\le M,{‖u^n‖}_{\mathrm{\infty }}\le M\mathrm{。}

证明 由定理3和引理5经过简单推导可得。

定理6 设初始条件u^{\mathrm{0}}\in H_E^{\mathrm{1}}\left[a,b\right]，则差分格式（5）的数值解u^n以||\cdot ||收敛到初边值问题（3）的精确解U^n，且收敛阶为O\left({\tau }^{\mathrm{2}}+h^{\mathrm{6}}\right)\mathrm{。}

证明 因为U^n是初边值问题（3）的精确解，故令e_j^n=u_j^n-U_j^n，误差表示为

r_j^n={\delta }_t{e}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}+\alpha {\left({\mathit{D}}_{\mathrm{1}}{A}_t{e}^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)}_j+\frac{\mathrm{1}}{\mathrm{2}}\beta {\left({\mathit{D}}_{\mathrm{1}}{\left(u^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)}^{\mathrm{2}}\right)}_j-\frac{\mathrm{1}}{\mathrm{2}}\beta {\left({\mathit{D}}_{\mathrm{1}}{\left(U^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)}^{\mathrm{2}}\right)}_j+\frac{\mathrm{1}}{\mathrm{3}}\sigma {\left({\mathit{D}}_{\mathrm{1}}{\left(A_t{u}^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)}^{\mathrm{3}}\right)}_j-\frac{\mathrm{1}}{\mathrm{3}}\sigma \left({\mathit{D}}_{\mathrm{1}}{\left(A_t{U}^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)}^{\mathrm{3}}\right)+\gamma {\left({\mathit{D}}_{\mathrm{3}}\left(A_t{e}^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)\right)}_j+\epsilon {\left({\mathit{D}}_{\mathrm{2}}{\delta }_t{e}^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)}_j\mathrm{。}（13）

为了得到误差函数e^N的范数估计，将式（13）与\mathrm{2}{e}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}作内积，可得

\begin{array}{l}\left(r_j^n,\mathrm{2}{e}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)=\left({\delta }_t{e}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}},\mathrm{2}{e}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)+\alpha \left({\left({\mathit{D}}_{\mathrm{1}}{A}_t{e}^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)}_j,\mathrm{2}{e}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)+\frac{\mathrm{1}}{\mathrm{2}}\beta \left({\left({\mathit{D}}_{\mathrm{1}}{\left(A_t{u}^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)}^{\mathrm{2}}\right)}_j,\mathrm{2}{e}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)-\frac{\mathrm{1}}{\mathrm{2}}\beta \left({\left({\mathit{D}}_{\mathrm{1}}{\left(A_t{U}^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)}^{\mathrm{2}}\right)}_j,\mathrm{2}{e}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)+\frac{\mathrm{1}}{\mathrm{3}}\sigma \left({\left({\mathit{D}}_{\mathrm{1}}{\left(A_t{u}^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)}^{\mathrm{3}}\right)}_j,\mathrm{2}{e}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)-\\ \frac{\mathrm{1}}{\mathrm{3}}\sigma \left({\left({\mathit{D}}_{\mathrm{1}}{\left(A_t{U}^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)}^{\mathrm{3}}\right)}_j,\mathrm{2}{e}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)+\gamma \left({\mathit{D}}_{\mathrm{3}}\left(A_t{e}^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right),\mathrm{2}{e}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)+\epsilon \left({\left({\mathit{D}}_{\mathrm{2}}{d}_t{e}^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)}_j,\mathrm{2}{e}_j^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)\mathrm{。}\end{array}（14）

由引理1及Cauchy-Schwarz不等式可知，

\left(\mathrm{1}+\epsilon M_{\mathrm{2}}\right)\frac{\mathrm{1}}{\tau }\left(‖e^{n+\mathrm{1}}‖-‖e^n‖\right)\le \frac{\mathrm{1}}{\mathrm{2}}‖r^n‖,（15）

对式（15）中的n从0到N-1求和，则有

‖e^N‖\le ‖e^{\mathrm{0}}‖+{\displaystyle \sum _{n=\mathrm{0}}^{N-\mathrm{1}}}\frac{\tau }{\mathrm{2}(\mathrm{1}+\epsilon M_{\mathrm{2}})}‖r^n‖，（16）

考虑‖e^{\mathrm{0}}‖=\mathrm{0}，且由第2节对D₁，D₂，D₃及A_t 和{\delta }_t的构造易知误差r_j^n的阶为O\left({\tau }^{\mathrm{2}}+h^{\mathrm{6}}\right)，有

\tau {\displaystyle \sum _{\mathrm{n}=\mathrm{0}}^{N-\mathrm{1}}}‖{\mathrm{r}}^n‖\le N\tau \underset{\mathrm{0}\le n\le N-\mathrm{1}}{\mathrm{m}\mathrm{a}\mathrm{x}}‖r^n‖=T\cdot O\left(h^{\mathrm{6}}+{\tau }^{\mathrm{2}}\right),（17）

则可得

‖e^N‖\le \frac{T}{\mathrm{2}(\mathrm{1}+\epsilon M_{\mathrm{2}})}\cdot O\left(h^{\mathrm{6}}+{\tau }^{\mathrm{2}}\right)。（18）

即，差分格式 （5） 的数值解u^n以||\cdot ||收敛到初边值问题 （3） 的精确解U^n，且收敛阶为O\left({\tau }^{\mathrm{2}}+h^{\mathrm{6}}\right)。

由||\cdot ||和||\cdot |{|}_{\mathrm{\infty }}的关系可知，数值解u^n以||\cdot |{|}_{\mathrm{\infty }}收敛到精确解U^n，且收敛阶为O\left({\tau }^{\mathrm{2}}+h^{\mathrm{6}}\right)。

4 数值实例

定义2 记能量和质量分别为

E^n=h{\displaystyle \sum _{j=\mathrm{0}}^{N-\mathrm{1}}}\left({\left(u^n\right)}_j^{\mathrm{2}}-\epsilon {\left({\mathit{D}}_{\mathrm{1}}{u}^n\right)}_j^{\mathrm{2}}\right)，Q^n=h{\displaystyle \sum _{j=\mathrm{0}}^{N-\mathrm{1}}}{u}_j^n，

在t=n\tau处的误差分别记为

R_{\epsilon }(t)=\left|E^n-E^{\mathrm{0}}\right|，Q_{\epsilon }(t)=\left|Q^n-Q^{\mathrm{0}}\right|\mathrm{。}

时间方向上与空间方向上的收敛阶为

\mathrm{O}\mathrm{r}\mathrm{d}\mathrm{e}{\mathrm{r}}_t=\mathrm{l}\mathrm{o}{\mathrm{g}}_{\frac{{\tau }_{\mathrm{1}}}{{\tau }_{\mathrm{2}}}}\frac{\mathrm{e}\mathrm{r}\mathrm{r}\left({\tau }_{\mathrm{1}}\right)}{\mathrm{e}\mathrm{r}\mathrm{r}\left({\tau }_{\mathrm{2}}\right)},\mathrm{O}\mathrm{r}\mathrm{d}\mathrm{e}{\mathrm{r}}_h=\mathrm{l}\mathrm{o}{\mathrm{g}}_{\frac{h_{\mathrm{1}}}{h_{\mathrm{2}}}}\frac{\mathrm{e}\mathrm{r}\mathrm{r}\left(h_{\mathrm{1}}\right)}{\mathrm{e}\mathrm{r}\mathrm{r}\left(h_{\mathrm{2}}\right)},

其中

\mathrm{e}\mathrm{r}\mathrm{r}\left(h,\tau \right)=\underset{\mathrm{0}\le n\le N}{\mathrm{m}\mathrm{a}\mathrm{x}}{‖u^n-U^n‖}_{\mathrm{\infty }}。

算例1 考虑用RLW⁃EKdV方程

\begin{array}{l}{u}_t+\alpha u_x+\beta u{u}_x+\sigma u^{\mathrm{2}}{u}_x+\gamma u_{xxx}+\epsilon u_{xxt}=\mathrm{0},x\in (a,b),t\in (\mathrm{0},T],\\ u(a,t)=u(b,t),u_x(a,t)=u_x(b,t),u_{xx}(a,t)=u_{xx}(b,t),t\in (\mathrm{0},T],\\ u(x,\mathrm{0})=u_{\mathrm{0}}(x),x\in [a,b]\end{array}（19）

的初边值问题描述单波。

取\alpha =\mathrm{1},\beta =\mathrm{1},\gamma =\mathrm{1},\epsilon =-\mathrm{1},\sigma =\mathrm{0}，则RLW-EKdV方程退化为RLW-KdV方程。取初始条件

u\left(x,\mathrm{0}\right)=\mathrm{0.9}\mathrm{s}\mathrm{e}\mathrm{c}{\mathrm{h}}^{\mathrm{2}}(\sqrt[]{\frac{\mathrm{3}}{\mathrm{92}}}x)\mathrm{。}

取x\in [-\mathrm{60,60}]，T=\mathrm{1}时间方向上的收敛阶如表1所示，空间方向上的收敛阶如表2所示。

由表1和表2可见，数值格式的收敛阶与理论结果一致，时间方向具有二阶收敛，空间方向上具有六阶精度，即数值结果与理论分析相一致，验证了算法的有效性。

取x\in [-\mathrm{100,100}]，T=\mathrm{40}，\sigma =\mathrm{0}，RLW⁃KdV方程的能量与质量误差见图1。由图1可见，RLW-KdV方程的能量和质量的误差至少达到O（10^-10）。

RLW-EKdV模型更具有通用性，当\sigma =\mathrm{0}，\gamma =\mathrm{0}时，模型（19）退化为RLW方程；当\sigma =\mathrm{0}，\beta =\mathrm{0}，时，退化为KdV方程；当\epsilon =\mathrm{0}时，退化为EKdV方程；当\sigma =\mathrm{0}时，退化为RLW-KdV方程。KdV型方程数值解模拟图见图2，图2（a）为KdV，EKdV，RLW，RLW-KdV方程数值解的比较。由图2（b）与图2（c）可见，RLW-EKdV方程的数值模拟结果与水槽实验观测到的波形特征更接近。

算例2 考虑用RLW-EKdV方程

\begin{array}{l}{u}_t+\alpha u_x+\beta u{u}_x+\sigma u^{\mathrm{2}}{u}_x+\gamma u_{xxx}+\epsilon u_{xxt}=\mathrm{0},x\in (a,b),t\in (\mathrm{0},T],\\ u(a,t)=u(b,t),u_x(a,t)=u_x(b,t),u_{xx}(a,t)=u_{xx}(b,t),t\in (\mathrm{0},T],\\ u(x,\mathrm{0})=u_{\mathrm{0}}(x),x\in [a,b]\end{array}

的初边值问题描述双孤立波。选取初值条件为

u(x,\mathrm{0})={\displaystyle \sum _{j=\mathrm{1}}^{\mathrm{2}}}\mathrm{3}{c}_j\mathrm{s}\mathrm{e}\mathrm{c}{\mathrm{h}}^{\mathrm{2}}\left(k_j\left(x-x_j\right)\right)，

其中C_j=\mathrm{4}{k}_j^{\mathrm{2}}/\left(\mathrm{1}-\mathrm{4}{k}_j^{\mathrm{2}}\right)，k_{\mathrm{1}}=\mathrm{0.4}，k_{\mathrm{2}}=\mathrm{0.3}，x_{\mathrm{1}}=\mathrm{15}，x_{\mathrm{2}}=\mathrm{35}。

取a=\mathrm{0}，b=\mathrm{140}，T=40，固定空间步长h=0.5，时间步长\tau =\mathrm{0.125}\mathrm{。}\mathrm{取}\alpha =\mathrm{1}，\beta =\mathrm{1}，\epsilon =-\mathrm{1}，\gamma =\mathrm{0.001}，\sigma =\mathrm{0.01}，为了便于直观观察两波相撞前后波的演化过程、能量及质量随时间变化，数值模拟如图3所示。

由图3可见，大振幅波传播速度快，小振幅波传播速度慢，大振幅波追上小振幅波后发生碰撞。在碰撞瞬间，大振幅波和小振幅波的振幅叠加。碰撞后，两者均恢复至碰撞前的波形，并继续维持各自的传播方向与特性，两波碰撞结束后与碰撞前能量及质量均守恒。

5 结论

本文主要利用紧致差分格式研究RLW-EKdV方程，数值格式保持守恒律，数值解在时间上达到二阶、在空间上达到六阶收敛精度，模拟结果可以描述内波在水槽及海洋中传播的某些特性。

（1）RLW-EKdV方程更具有通用性，当系数发生变化时，它可以退化为KdV方程、EKdV方程、RLW方程及RLW-KdV方程等。

（2）目前，尚未发现RLW-EKdV方程初边值问题的精确解，故数值实验中利用RLW-KdV初边值问题的精确解获取初始值进行数值模拟，当结果稳定后由图可见有尾波出现，这与内波在水槽中的某些传播特性更接近。且算法的质量误差阶可达到O（10^-12），能量误差阶至少可达到O（10^-10）。

模拟结果显示紧致差分方法在保持质量和动量守恒方面表现良好，这有利于对相应物理问题和物理现象的准确理解。事实上，没有尾波且波形保持不变的理想孤立波在实际问题中是比较少见的，RLW-EKdV方程可以描述内孤立波在水槽中和海洋的传播问题，该研究具有一定价值和实际意义。

参考文献

[1] PEREGRINE D H.Calculations of the development of an unduiar bore［J］.Journal of Fluid Mechanics，1966，25：321-330．

[2] 吕秀敏，葛倩，李金.重心插值配点法求解小振幅长波广义BBM-KdV方程［J］.山东大学学报（理学版），2024，59（8）：67-76.

[3] 邱天威，魏光美，宋禹欣，等.基于PINN方法的KdV类方程新孤子解的研究［J］.应用数学和力学，2025，46（1）：105-113.

[4] 王一辰，王桂霞，李骞.基于物理信息神经网络算法的海洋内孤立波研究［J］.内蒙古师范大学学报（自然科学版），2025，54（2）：198-206.

[5] 栗雪娟，刘瑜欣.基于PINN及其改进算法求解KdV-mKdV方程［J］.浙江大学学报（理学版），2024，51（6）：702-711.

[6] 李烨.非线性BBM-KdV方程的局部间断有限元方法［D］.哈尔滨：哈尔滨工业大学，2022.

[7] ASOKAN R， VINODH D.Soliton and exact solutions for the KdV-BBM type equations by tanh-coth and transformed rational function methods［J］.International Journal of Applied and Computational Mathematics，2018，4（4）： 100-107.

[8] 文宗川，郭彦，梁静国，等.KdV方程的特征线混合间断有限元方法（英文）［J］.应用数学，2007，20（4）：791-800.

[9] LI S G.Numerical study of a conservative weighted compact difference scheme for the symmetric regularized long wave equations［J］.Numerical Methods for Partial Differential Equations，2019，35（1）：60-83.

[10] HE Y Y， WANG X F， CHENG H， et al.Numerical analysis of a high-order accurate compact finite difference scheme for the SRLW equation［J］.Applied Mathematics and Computation，2022，418：126837.

[11] HE Y Y， WANG X F， ZHONG R H. A new linearized fourth-order conservative compact difference scheme for the SRLW equations［J］. Advances in Computational Mathematics， 2022， 48（3）： 27.

[12] HU J S， LI J J， WANG X.New high-order conservative difference scheme for regularized long wave equation with Richardson extrapolation［J］.Thermal Science，2019，23（3）：737-745.

[13] YANG X J， ZHANG L， GE Y B.High-order compact finite difference schemes for solving the regularized long-wave equation［J］.Applied Numerical Mathematics，2023，185：165-187.

[14] 邓雅清，王晓峰，何育宇，等.RLW-KdV方程的高阶紧致有限差分格式［J］.数学进展，2022，51（2）：360-374.

[15] POLWANG A， POOCHINAPAN K， WONGSAIJAI B.Numerical simulation of wave flow： Integrating the BBM-KdV equation using compact difference schemes［J］.Mathematics and Computers in Simulation，2025，236： 70-89.

[16] 王晨逦，王桂霞，萨和雅.EKdV方程的高阶多辛保结构算法及孤立波解的数值模拟［J］.内蒙古师范大学学报（自然科学汉文版），2023，52（2）：119-126.

[17] 高晶英，何斯日古楞，青梅，等.对称正则长波方程的高效紧致差分格式［J］.应用数学和力学，2025， 46（3）： 412-424.

[18] 郅长红，徐双东，韩盼盼，等.高阶单向传播内孤立波理论模型适用性［J］.物理学报，2022，71（17）： 254-263.

基金资助

国家自然科学基金资助项目“全海深内孤立波统一表征及其在深⁃浅过渡区遥感方法反演应用研究”(62161045)

内蒙古师范大学数学一流拔尖培育学科建设资助项目“哈密顿保结构算法在海洋内波问题中的应用研究”(2024YLKY19)

内蒙古师范大学学校专项基金资助项目“非线性海洋内波方程的数值方法研究”(KJJYZX2026098)