带耗散项RLW方程的高阶有限差分算法

High-Order Finite Difference Algorithm for RLW Equation with Dissipation Term

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

Abstract: To conduct numerical simulation research on the initial-boundary value problem of the nonlinear RLW equation, a high-order conservative compact finite difference scheme is first constructed. The derivative in the time direction is discretized using the implicit mid-point scheme, and the derivatives of each order in the spatial direction are discretized using the inverse compact operator, achieving second-order accuracy in the time direction and sixth-order accuracy in the spatial direction. Further, the conservation property of the difference scheme is proved, a priori estimates are made, and the convergence analysis is carried out. Finally, the correctness of the theory and the effectiveness and reliability of the scheme are verified through numerical examples.

1. 引言

1966年，Peregrine [1]在研究水波问题时首次提出非线性正则长波(RLW)方程(也称为BBM方程)

$u_t+u_x+u{u}_x-u_{xxt}=0,$ (1)

该方程与非线性KdV方程

$u_t+c_0{u}_x+\beta u{u}_x+\gamma u_{xxx}=0$ (2)

所描述的运动相似，能够描述浅水域中的弱非线性孤立波的演化过程。目前求解非线性RLW型方程的数值方法主要包括物理信息神经网络(PINN)方法、有限元法及(紧致)有限差分法，成果丰富[2]-[10]。紧致有限差分格式计算量较小，仅用3~5点模板即可达到四到六阶精度，适合求解数值耗散误差敏感、求解区域规则的问题。

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

$\begin{array}{l}{u}_t+\left(c_0+\beta u+\sigma u^2\right)u_x+\epsilon u_{xxt}+\gamma u=0,x\in \left(a,b\right),t\in \left(0,T\right],\\ u\left(a,t\right)=u\left(b,t\right),u_x\left(a,t\right)=u_x\left(b,t\right)=0,t\in \left(0,T\right],\\ u\left(x,0\right)=u_0\left(x\right),x\in \left[a,b\right],\end{array}$ (3)

其中，$c_0$ 为相速度，$\beta ,\sigma$ 分别为低阶非线性项系数和高阶非线性项系数，$\epsilon$ 为涡动系数，$\gamma$ 为耗散项系数。此非线性初边值问题的解析解比较困难，用数值方法研究该问题的解显得尤为重要。

全文结构如下：第2节引入一些相关的算子符号，列出相关的定义和引理。第3节构造方程的有限差分格式，对时间方向一阶导数采用隐中点格式进行离散，分别对空间方向的各阶导数用六阶逆紧致算子进行离散验证差分格式的守恒律，对格式进行先验估计与收敛性分析。第4节给出一些数值实例。

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

2. 差分格式的构造

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

$D_t^{\frac{\gamma }{2}}{u}^n=\frac{e^{\frac{\gamma }{2}\tau }{u}^{n+1}-e^{-\frac{\gamma }{2}\tau }{u}^n}{\tau },A_t^{\frac{\gamma }{2}}{u}_j=\frac{e^{\frac{\gamma }{2}\tau }{u}^{n+1}+e^{-\frac{\gamma }{2}\tau }{u}^n}{2}$ (4)

当$\gamma =0$ 时，这些算子退化为标准算子。

定义算子$P_x$ 、$Q_x$ 使得

$P_x{u}_j^n=\frac{u_{j+1}^n+3u_j^n+u_{j-1}^n}{3}\text{, }{Q}_x{u}_j^n=\frac{2u_{j+1}^n+11u_j^n+2u_{j-1}^n}{11}\text{,}$ (5)

满足

$P_x{\left(\frac{\partial U}{\partial x}\right)}_j^n={\delta }_x^1{U}_j^n+O\left(h^6\right)\text{, }{Q}_x{\left(\frac{{\partial }^{2}U}{\partial x^2}\right)}_j^n={\delta }_x^2{U}_j^n+O\left(h^6\right).$ (6)

由式(5)可知$P_x$ 、$Q_x$ 是对角占优的，所以可逆，从而有

${\left(\frac{\partial U}{\partial x}\right)}_j^n=P_x^{-1}{\delta }_x^1{U}_j^n+O\left(h^6\right),{\left(\frac{{\partial }^{2}U}{\partial x^2}\right)}_j^n=Q_x^{-1}{\delta }_x^2{U}_j^n+O\left(h^6\right).$ (7)

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

${\left(\frac{\partial U}{\partial x}\right)}_j^n=D_1{U}^n+O\left(h^6\right),{\left(\frac{{\partial }^{2}U}{\partial x^2}\right)}_j^n=D_2{U}^n+O\left(h^6\right).$ (8)

其中

$U^n={\left(U_0^n,\cdots ,U_J^n\right)}^{\text{T}}\text{,}{D}_1=P_x^{-1}{\delta }_x^1,D_2=Q_x^{-1}{\delta }_x^2.$

基于Strange分裂，将问题(3)中的控制方程分裂为方程

$u_t=-\left(c_0+\beta u+\sigma u^2\right)u_x-\epsilon u_{xxt},x\in \left(a,b\right),t\in \left(0,T\right],$ (9)

和方程

$u_t=-\gamma u,x\in \left(a,b\right),t\in \left(0,T\right].$ (10)

对第一个方程，时间方向的一阶导数采用隐中点格式进行离散，空间方向的各阶导数用六阶紧致差分格式进行离散，记$u^{n+1}={\xi }_{\tau }\left(u^n\right)$ 。第二个方程容易求得解析解，记$u^{n+1}={\varphi }_{\tau }\left(u^n\right)=e^{-\gamma \tau }\left(u^n\right)$ ，两个方程解的组合${\varphi }_{\tau /2}\circ {\xi }_{\tau }\circ {\varphi }_{\tau /2}$ 即为原方程的解。

对由方程(9)和方程(10)构成的方程组在网格点$\left(x_j,t_{n\text{+}\frac{1}{2}}\right)$ 处进行离散，其离散方程组为

$\{\begin{array}{l}{u}_t\left(x_j,t_{n+\frac{1}{2}}\right)=-\left(c_0{u}_x\left(x_j,t_{n+\frac{1}{2}}\right)+\sigma \left(u^2{u}_x\right)\left(x_j,t_{n+\frac{1}{2}}\right)+\beta \left(u{u}_x\right)\left(x_j,t_{n+\frac{1}{2}}\right)+\epsilon u_{xxt}\left(x_j,t_{n+\frac{1}{2}}\right)\right)\hfill \\ u_t\left(x_j,t_{n+\frac{1}{2}}\right)=-\gamma u\left(x_j,t_{n+\frac{1}{2}}\right).\hfill \end{array}$ (11)

对方程组(11)式中的第一个方程在时间方向用隐中点格式离散，在空间方向用六阶紧致差分格式离散得差分方程

$r_j^n+D_t^{\frac{\gamma }{2}}{U}_j^{n+\frac{1}{2}}+c_0{\left(D_1{A}_t^{\frac{\gamma }{2}}{U}^{n+\frac{1}{2}}\right)}_j+\frac{1}{2}\beta {\left(D_1{\left(A_t^{\frac{\gamma }{2}}{U}^{n+\frac{1}{2}}\right)}^2\right)}_j+\frac{1}{3}\sigma {\left(D_1{\left(A_t^{\frac{\gamma }{2}}{U}^{n+\frac{1}{2}}\right)}^3\right)}_j+\epsilon {\left(D_2{D}_t^{\frac{\gamma }{2}}{U}^{n+\frac{1}{2}}\right)}_j=0,$ (12)

其中$r_j^n$ 为截断误差。再用数值解$u_j^n$ 近似代替精确解$U_j^n$ 并忽略高阶小项，即可得问题(3)的六阶紧致差分格式：

(13)

由对D₁，D₂，D₃及A_t和${\delta }_t$ 的构造易知$r_j^n$ 的误差阶为$O\left({\tau }^2+h^6\right)$ 。

下面给出相关的定义和引理。

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

$\left(f^n,g^n\right)=h{\displaystyle \sum _{j=1}^{J-1}{f}_j^n}{g}_j^n,\Vert f^n\Vert =\sqrt{\left(f^n,f^n\right)},{\Vert f^n\Vert }_{\infty }=\underset{1\le j\le J-1}{\max }\left|f_j^n\right|.$

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

引理1 对任意的两个网格函数$u^n$ ，${\rho }^n\in Z_h$ ，有

$\left(P_x^{-1}{\delta }_x^1{u}^n,{\rho }^n\right)=-\left(u^n,P_x^{-1}{\delta }_x^1{\rho }^n\right),$

特别地，

$\left(P_x^{-1}{\delta }_x^1{u}^n,u^n\right)=0.$

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

$\left(Q_x^{-1}{\delta }_x^2\left(u^n-{\rho }^n\right),u^n\text{+}{\rho }^n\right)=-\left({\Vert P_x^{-1}{\delta }_x^1{u}^n\Vert }^2-{\Vert P_x^{-1}{\delta }_x^1{\rho }^n\Vert }^2\right).$

3. 先验估计与收敛性分析

3.1. 离散格式的守恒律

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

定理1 定义能量和质量分别为

$E^n=e^{-2\gamma \tau }{\Vert u^n\Vert }^2-\epsilon e^{-2\gamma \tau }{\Vert D_1{u}^n\Vert }^2,Q^n=e^{-\gamma \tau }h{\displaystyle \sum _{j=0}^J{u}_j^n}.$ (14)

对于任意的$n\ge 0$ ，格式(14)满足离散的能量和质量衰减律，即

$E^{n\text{+}1}=e^{-2\gamma \tau }{E}^n,Q^{n\text{+}1}=e^{-\gamma \tau }{Q}^n.$ (15)

证明 用$A_t^{\frac{\gamma }{2}}{u}^{n+\frac{1}{2}}$ 与式(13)作内积可得

(16)

由定义1，引理1及引理2可得

$e^{\gamma \tau }{\Vert u^{n+1}\Vert }^2-\epsilon e^{\gamma \tau }{\Vert D_1{u}^{n+1}\Vert }^2=e^{-\gamma \tau }{\Vert u^n\Vert }^2-\epsilon e^{-\gamma \tau }{\Vert D_1{u}^n\Vert }^2,$

由$E^n$ 的定义可知，$E^{n\text{+}1}=e^{-2\gamma \tau }{E}^n$ 。

接下来证明$Q^{n\text{+}1}={\text{e}}^{-\gamma \tau }{Q}^n$ 。由式(13)可得

$\begin{array}{l}h{\displaystyle \sum _{j=0}^J\frac{e^{\frac{\gamma \tau }{2}}{u}_j^{n+1}-e^{-\frac{\gamma \tau }{2}}{u}_j^n}{\tau }}+c_{0}h{\displaystyle \sum _{j=0}^J{A}_t^{\frac{\gamma }{2}}{P}_x^{-1}{\delta }_x^1{u}_j^{n+\frac{1}{2}}}+\frac{1}{2}\beta h{\displaystyle \sum _{j=0}^J{A}_t^{\frac{\gamma }{2}}{P}_x^{-1}{\delta }_x^1{\left(u_j^{n+\frac{1}{2}}\right)}^2}\\ +\frac{1}{3}\sigma h{\displaystyle \sum _{j=0}^J{A}_t^{\frac{\gamma }{2}}{P}_x^{-1}{\delta }_x^1{\left(u_j^{n+\frac{1}{2}}\right)}^3}+\epsilon h{\displaystyle \sum _{j=0}^J{D}_t^{\frac{\gamma }{2}}{Q}_x^{-1}{\delta }_x^2}{u}_j^{n+\frac{1}{2}}=0,\end{array}$ (17)

再由边界条件整理可得

$h{\displaystyle \sum _{j=0}^J{u}_j^{n+1}}=h{\displaystyle \sum _{j=0}^J{e}^{-\gamma \tau }}{u}_j^n$ .

由$Q^n$ 的定义可知$Q^{n\text{+}1}=e^{-\gamma \tau }{Q}^n$ 。

3.2. 先验估计与收敛性分析

定理2 假设$U_0\in H_E^1\left[a,b\right]$ ，则问题(3)的精确解$U$ 满足

$\Vert U\Vert \le M,\Vert U_x\Vert \le M.$

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

定理3 假设$u_0\in H_E^1\left[a,b\right]$ ，则格式(14)的数值解$u^n$ 满足

$\Vert u^n\Vert \le M,\Vert D_t^{\frac{\gamma }{2}}{u}^n\Vert \le M,{\Vert u^n\Vert }_{\infty }\le M.$

证明 由定理1经过推导易得。

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

证明 因为$U^n$ 是初边值问题(3)的精确解，故令$A_t^{\frac{\gamma }{2}}{e}_j^n=A_t^{\frac{\gamma }{2}}{u}_j^n-A_t^{\frac{\gamma }{2}}{U}_j^n$ ，由式(13)减式(12)得

$\begin{array}{c}{r}_j^n=D_t^{\frac{\gamma }{2}}{e}_j^{n+\frac{1}{2}}+c_0{\left(D_1{A}_t^{\frac{\gamma }{2}}{e}^{n+\frac{1}{2}}\right)}_j+\frac{1}{2}\beta {\left(D_1{\left(A_t^{\frac{\gamma }{2}}{u}^{n+\frac{1}{2}}\right)}^2\right)}_j-\frac{1}{2}\beta {\left(D_1{\left(A_t^{\frac{\gamma }{2}}{U}^{n+\frac{1}{2}}\right)}^2\right)}_j\\ +\frac{1}{3}\sigma {\left(D_1{\left(A_t^{\frac{\gamma }{2}}{u}^{n+\frac{1}{2}}\right)}^3\right)}_j-\frac{1}{3}\sigma {\left(D_1{\left(A_t^{\frac{\gamma }{2}}{U}^{n+\frac{1}{2}}\right)}^3\right)}_j+\epsilon {\left(D_2{D}_t^{\frac{\gamma }{2}}{e}^{n+\frac{1}{2}}\right)}_j,\end{array}$ (18)

将式(18)与$A_t^{\frac{\gamma }{2}}{e}_j^{n+\frac{1}{2}}$ 作内积，可得

$\begin{array}{c}\left(r_j^n,A_t^{\frac{\gamma }{2}}{e}_j^{n+\frac{1}{2}}\right)=\left(D_t^{\frac{\gamma }{2}}{e}_j^{n+\frac{1}{2}},A_t^{\frac{\gamma }{2}}{e}_j^{n+\frac{1}{2}}\right)+c_0\left({\left(D_1{A}_t{e}^{n+\frac{1}{2}}\right)}_j,A_t^{\frac{\gamma }{2}}{e}_j^{n+\frac{1}{2}}\right)\\ +\frac{1}{2}\beta \left({\left(D_1{\left(A_t^{\frac{\gamma }{2}}{u}^{n+\frac{1}{2}}\right)}^2\right)}_j,A_t^{\frac{\gamma }{2}}{e}_j^{n+\frac{1}{2}}\right)-\frac{1}{2}\beta \left({\left(D_1{\left(A_t^{\frac{\gamma }{2}}{U}^{n+\frac{1}{2}}\right)}^2\right)}_j,A_t^{\frac{\gamma }{2}}{e}_j^{n+\frac{1}{2}}\right)\\ +\frac{1}{3}\sigma \left({\left(D_1{\left(A_t^{\frac{\gamma }{2}}{u}^{n+\frac{1}{2}}\right)}^3\right)}_j,A_t^{\frac{\gamma }{2}}{e}_j^{n+\frac{1}{2}}\right)-\frac{1}{3}\sigma \left({\left(D_1{\left(A_t^{\frac{\gamma }{2}}{U}^{n+\frac{1}{2}}\right)}^3\right)}_j,A_t^{\frac{\gamma }{2}}{e}_j^{n+\frac{1}{2}}\right)\\ +\epsilon \left({\left(D_2{D}_t^{\frac{\gamma }{2}}{e}^{n+\frac{1}{2}}\right)}_j,A_t^{\frac{\gamma }{2}}{e}_j^{n+\frac{1}{2}}\right).\end{array}$ (19)

考虑到$D_1,D_3$ 是反对称矩阵，令矩阵$D_2$ 的最小特征值为$M$ ，则有

$\left(r_j^n,A_t^{\frac{\gamma }{2}}{e}_j^{n+\frac{1}{2}}\right)\le \Vert r^n\Vert \frac{e^{\gamma \tau }\Vert e^{n+1}\Vert +e^{-\gamma \tau }\Vert e^{n+1}\Vert }{2},$ (20)

而

$\left(r_j^n,A_t^{\frac{\gamma }{2}}{e}_j^{n+\frac{1}{2}}\right)=\frac{1}{\tau }\left(e^{\gamma \tau }{\Vert e^{n+1}\Vert }^2-e^{-\gamma \tau }{\Vert e^n\Vert }^2\right)+\left({\left(D_2{D}_t^{\frac{\gamma }{2}}{e}^{n+\frac{1}{2}}\right)}_j,2A_t^{\frac{\gamma }{2}}{e}_j^{n+\frac{1}{2}}\right).$ (21)

整理得

$\left(1+\epsilon M\right)\frac{1}{\tau }\left(e^{\gamma \tau }\Vert e^{n+1}\Vert -e^{-\gamma \tau }\Vert e^n\Vert \right)\le \frac{1}{2}\Vert r^n\Vert ,$ (22)

对式(22)中的n从0到N-1求和，则有

$e^{\gamma \tau }\Vert e^N\Vert \le e^{\text{-}\gamma \tau }\Vert e^0\Vert +{\displaystyle \sum _{n=0}^{N-1}\frac{\tau }{2\left(1+\epsilon M\right)}}\Vert r^n\Vert .$ (23)

假设$\Vert e^0\Vert =0$ ，根据

$\tau {\displaystyle \sum _{n=0}^{N-1}\Vert r^n\Vert }\le N\tau \underset{0\le n\le N-1}{\max }\Vert r^n\Vert =T\cdot O\left(h^6+{\tau }^2\right),$ (24)

可得

$\Vert e^N\Vert \le \frac{T}{2\left(1+\epsilon M\right)}\cdot O\left(h^6+{\tau }^2\right)$ . (25)

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

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

4. 数值实例

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

$E^n=e^{-2\gamma \tau }h{\displaystyle \sum _{j=0}^{N\text{-1}}\left({\left(u^n\right)}_j^2-\epsilon {\left(D_1{u}^n\right)}_j^2\right)}$ , $Q^n=e^{-\gamma \tau }h{\displaystyle \sum _{j=0}^{N\text{-1}}{u}_j^n}$ .

时间方向上与空间方向上的收敛阶如下

$Orde{r}_t={\log }_{\frac{{\tau }_1}{{\tau }_2}}\frac{err\left({\tau }_1\right)}{err\left({\tau }_2\right)},Orde{r}_h={\log }_{\frac{h_1}{h_2}}\frac{err\left(h_1\right)}{err\left(h_2\right)},$

其中

$err\left(\tau ,h\right)=\underset{0\le n\le N}{\max }\Vert u^n-U^n\Vert$ .

算例4.1 考虑用RLW方程

$\begin{array}{l}{u}_t+c_0{u}_x+\beta u{u}_x+\sigma u^2{u}_x+\epsilon u_{xxt}=0,x\in \left(a,b\right),t\in \left(0,T\right],\\ u\left(a,t\right)=u\left(b,t\right),u_x\left(a,t\right)=u_x\left(b,t\right)=0,t\in \left(0,T\right],\\ u\left(x,0\right)=u_0\left(x\right),x\in \left[a,b\right]\end{array}$ (26)

的初边值问题描述单波。

取$c_0=1$ ，$\beta =1$ ，$\epsilon =-1$ ，$\sigma =0$ ，方程退化为RLW方程(2)。RLW方程有解析解

$u\left(x,t\right)={\mathrm{sech}}^2\left(\left[\frac{1}{4}x-\frac{1}{3}t\right]\right).$

取$x\in \left[-50,50\right]$ ，$T=1$ 。时间方向上的收敛阶如表1所示，空间方向上的的收敛阶如表2所示。表3为不同格式下的误差比较。

Table 1. Numerical solution error and time convergence order (h = 0.5)

表1. 数值解误差与时间收敛阶(h = 0.5)

$\tau$	${\Vert \cdot \Vert }_{\infty }$ 误差	误差阶
0.5	0.0027	——
0.25	6.917e−04	1.9674
0.125	1.74164−04	1.9897
0.0625	4.3906e−05	1.9879

Table 2. Numerical solution error and spatial convergence order (τ = 0.001)

表2. 数值解误差与空间收敛阶(τ = 0.001)

h	${\Vert \cdot \Vert }_{\infty }$ 误差	误差阶
0.8	8.3389e−06	——
0.625	1.9453e−06	5.9028
0.5	4.8963e−07	6.1758
0.4	1.3000e−7	5.5429

Table 3. The comparison results of numerical errors of this paper’s format with those of references [11]-[13] at T = 1 moment

表3. 在T = 1时刻本文格式与文献[11]-[13]的数值误差比较结果

文献	h	$\tau$	${\Vert \cdot \Vert }_{\infty }$ 误差
格式(13)	0.4	0.01	1.0179e−06
文献[11]	0.1	0.00625	1.2671e−06
文献[12]	0.1	0.00625	1.4205e−06
文献[13]	0.1	0.00625	1.1220e−04

由表1和表2可见，时间方向具有二阶精度，空间方向上具有六阶精度，即数值结果与理论分析相一致，验证了算法的有效性。由表3可见，不需要太小步长即可达到更高精度，说明格式(13)在计算效率方面有明显的优势。

算例4.2 考虑带耗散项的RLW方程

$\begin{array}{l}{u}_t+c_0{u}_x+\beta u{u}_x+\sigma u^2{u}_x+\epsilon u_{xxt}+\gamma u=0,x\in \left(a,b\right),t\in \left(0,T\right],\\ u\left(a,t\right)=u\left(b,t\right),u_x\left(a,t\right)=u_x\left(b,t\right)=0,t\in \left(0,T\right],\\ u\left(x,0\right)=u_0\left(x\right),x\in \left[a,b\right]\end{array}$ (27)

的初边值问题模拟海洋内孤立波。

假设海水分为上下两层，上层厚度$h_1=100\text{m}$ ，下层厚度$h_2=300\text{m}$ ，上层密度${\rho }_1=1025\text{kg}/{\text{m}}^3$ ，下层密度${\rho }_2=1025\text{kg}/{\text{m}}^3$ 。参照文献[14]，取$\epsilon =-2\times {10}^{-6}{\text{km}}^2$ ，各系数计算如下：

$c_0=\sqrt{\frac{g{h}_1{h}_2\left({\rho }_2-{\rho }_1\right)}{{\rho }_1{h}_2+{\rho }_2{h}_1}}=1.621\times {10}^{-3}\text{km}/\text{s}$ ，$\beta =-\frac{3c_0}{2}\frac{{\rho }_1{h}_2^2-{\rho }_2{h}_1^2}{{\rho }_1{h}_1{h}_2^2+{\rho }_2{h}_2{h}_1^2}=-0.022{\text{s}}^{-1}$ ，

$\sigma =-3\frac{h_1^2-h_1{h}_2+h_2^2}{h_1{h}_2}{c}_0=-4.462\times {10}^{-4}{\text{km}}^{-1}\cdot {\text{s}}^{-1}$ 。

取初始波

$u\left(x,t\right)=a_0{\mathrm{sech}}^2\left(\frac{1}{4}x\right)$ ，

初始振幅$a_0$ 为20 m，$a=-50\text{km}$ ，$b=50\text{km}$ ，$T=100$ ，空间步长$h=0.5$ ，时间步长$\tau =0.25$ 。分别取$\gamma =0.01$ 和$\gamma =0.005$ ，海洋内孤立波的演化过程及质量与能量随时间的变化如图1~8所示。

Figure 1. The propagation and evolution diagram of a single solitary wave (γ = 0.01)

图1. 单孤立波传播演变图(γ = 0.01)

Figure 2. The propagation and evolution diagram of a single solitary wave (γ = 0.005)

图2. 单孤立波传播演变图(γ = 0.005)

Figure 3. The numerical quality from t = 0 to t = 100 (γ = 0.01)

图3. t从0到100时的质量(γ = 0.01)

Figure 4. The numerical quality from t = 0 to t = 100 (γ = 0.005)

图4. t从0到100时的质量(γ = 0.005)

Figure 5. The numerical energy from t = 0 to t = 100 (γ = 0.01)

图5. t从0到100时的能量(γ = 0.01)

Figure 6. The numerical energy from t = 0 to t = 100 (γ = 0.005)

图6. t从0到100时的能量(γ = 0.005)

如图1至图2所示，随着耗散系数$\gamma$ 的减小，海洋内孤立波的振幅在不断地增大。图3所示质量规律近似表示为$y=0.16e^{-0.01t}$ ，图4所示质量规律近似表示为$y=0.16e^{-0.005t}$ ，由图3和图4可见全局质量随时间呈指数衰减。图5所示能量规律近似表示为$y=0.021e^{-0.02t}$ ，图6所示能量规律近似表示为$y=0.021e^{-0.01t}$ ，表明全局能量随时间呈指数衰减。

算例4.3 考虑用RLW方程

$\begin{array}{l}{u}_t+c_0{u}_x+\beta u{u}_x+\sigma u^2{u}_x+\epsilon u_{xxt}+\gamma u=0,x\in \left(a,b\right),t\in \left(0,T\right],\\ u\left(a,t\right)=u\left(b,t\right),u_x\left(a,t\right)=u_x\left(b,t\right)=0,t\in \left(0,T\right],\\ u\left(x,0\right)=u_0\left(x\right),x\in \left[a,b\right]\end{array}$

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

$u\left(x,0\right)={\displaystyle \sum _{j=1}^{2}3c_j{\mathrm{sech}}^2\left(k_j\left(x-x_j\right)\right)}$ ，

其中$C_j=4k_j^2/\left(1-4k_j^2\right)$ ，$k_1=0.4$ ，$k_2=0.3$ ，$x_1=15$ ，$x_2=35$ 。

(a) γ = 0.001 (b) γ = 0.01

Figure 7. Three-dimensional graph of double waves varying with the dissipation term

图7. 双波随耗散项变化的三维图

(a) γ = 0.001 (b) γ = 0.01

Figure 8. Two-dimensional graph of double waves varying with the dissipation term

图8. 双波随耗散项变化的二维图

(a) (b)

(c)

Figure 9. Numerical simulation of double isolated waves. (a) Evolution of double solitary waves; (b) The numerical quality from t = 0 to t = 40; (c) The numerical energy from t = 0 to t = 40

图9. 双孤立波的数值模拟。(a) 双孤立波演化；(b) t从0到100时的质量；(c) t从0到100时的能量

取$a=0$ ，$b=150$ ，$T=40$ ，固定空间步长$h=0.5$ ，时间步长$\tau =0.125$ ，$\alpha =1$ ，$\beta =1$ ，$\sigma =0.01$ ，$\epsilon =-1$ 。下面分别讨论耗散项对双波结果的影响。

观察图7和图8可见，当耗散项$\gamma$ 增加时，双波振幅减小，传播速度减慢。为了便于直观观察，取$\sigma =0.005$ ，耗散项系数为$\gamma =1\times {10}^{-7}$ ，两波相撞前后波的演化过程、能量及质量随时间变化如图9所示。

由图9(a)可见，大振幅波传播速度快，小振幅波速度慢，大振幅波追上小振幅波后发生碰撞。碰撞瞬间，大振幅波和小振幅波的振幅叠加。碰撞后，两者均恢复至碰撞前的波形，并继续维持各自的传播方向与特性。观察图9(b)，由于耗散项的存在，质量随时间的变化规律呈现$y=37.916533e^{-{10}^{-7}t}$ ，表明质量呈现指数衰减。观察图9(c)，由于耗散项的存在，能量也呈现衰减趋势，在衰减的过程中，由于两波相撞出现了能量干扰，之后能量继续呈现指数衰减。

5. 小结

主要利用紧致差分格式研究带耗散项的高阶非线性RLW方程，数值解在时间上达到二阶、在空间上达到六阶收敛精度。目前，尚未发现带耗散项高阶RLW方程的初边值问题的精确解，所以采用Strange分裂把原方程分裂为两个方程，将两个解进行组合即可得到该方程的数值解，在保证原有精度的基础上，有效提升计算稳定性与可实施性。

通过带耗散项的RLW方程模拟浅水域海洋内孤立波，结果显示：当耗散项逐渐增大，经过长时间传播演化，振幅逐渐减小，波的传播速度减慢，符合实际的物理规律[15]，验证了算法的有效性与可靠性。同时，模拟结果显示质量与能量随时间变化呈现指数耗散的情况，这为刻画实际海洋环境中内孤立波的能量衰减、波形演变提供了依据。

事实上，没有耗散且波形保持不变的理想孤立波在实际问题中是比较少见的，该研究具有一定理论价值和实际意义。

基金项目

国家自然基金(62161045)；教育部重点实验室开放课题项目(2023KFZDO1)；内蒙古自治区自然基金(2023LHMS01007)；内蒙古自治区“揭榜挂帅”项目(2024JBGS0010-1)；内蒙古师范大学数学一流拔尖培育学科建设项目(2024YLKY19)。

NOTES

^*通讯作者。