非局部粘性水波模型的直接间断Galerkin算法

Direct Discontinuous Galerkin Algorithm for a Nonlocal Viscous Water Wave Model

摘要: 本文基于直接间断有限元(DDG)方法数值求解非局部粘性水波模型。该算法结合L1近似公式与BDF2方法，系统构建了非线性时间分数阶偏微分方程的高效数值算法。首先，运用分部积分法对模型的弱形式进行降阶处理。其次，通过引入边界项和选用合适的数值通量，确保离散格式的稳定性。最后，针对时间导数项应用L1近似公式与BDF2时间差分的离散方法，建立全离散DDG格式。文中详细给出数值格式的构造过程并严格证明该算法的稳定性。数值实验部分选取无已知解析解的水波模型，验证该算法在时空离散上的高精度特性。

Abstract: This paper numerically solves the nonlocal viscous water wave model based on the Direct Discontinuous Galerkin (DDG) method. This algorithm combines the L1 approximation formula with the BDF2 method to systematically construct an efficient numerical algorithm for nonlinear time-fractional partial differential equations. Firstly, the integration by parts method is used to reduce the order of the weak formula. Secondly, the stability of the discrete scheme is ensured by introducing the boundary term and constructing a stable numerical flux. Finally, for the time derivative term, the discrete methods of the L1 approximation formula and the BDF2 time difference are applied to construct the fully discrete DDG scheme. This paper gives a detailed description of the construction process of the numerical scheme and strictly proves the stability of this algorithm. In the numerical experiment part, a water wave model without a known analytical solution is selected to verify the high-precision characteristics of this algorithm in space-time discretization.

1. 引言

本文研究以下具有周期边界条件的非局部粘性水波模型

$\begin{array}{l}{u}_t+u_x+\beta u_{xxx}+\frac{\sqrt{\nu }}{\sqrt{\pi }}{{\displaystyle \int }}_0^t\frac{\partial f\left(s\right)}{\partial s}\frac{\text{d}s}{\sqrt{t-s}}+\alpha u{u}_x=\gamma u_{xx},\left(x,t\right)\in \overline{\Omega }\times \left(0,T\right],\\ u\left(x,0\right)=u_0\left(x\right),x\in \overline{\Omega }.\end{array}$ (1)

其中$\overline{\Omega }=\left[a,b\right]$ 代表空间区间，$\beta ,T,\alpha ,\gamma$ 是给定的正常数。

分数阶水波模型是一类用于刻画复杂水波现象的数学模型，与整数阶水波模型相比，分数阶导数可以更准确地刻画水波的传播和衰减。该类模型在流体力学、海洋工程和环境科学等领域具有重要的研究价值。在[1]中，Dutykh与Dias通过对流体边界层推导，得出了非局部粘性势自由表面流动方程。2010年，Chen等人在[2]中从有限水深条件下的自由表面问题出发，成功获得了带有非局部粘性色散项的水波模型。鉴于分数阶导数与非线性项的存在，此类模型难以获得精确解。因此，对这类模型的数值方法展开研究具有极为重要的意义。Lin和Xu在文献[3]中研究了时间分数阶水波模型的谱逼近法。在[4]中，Liu等人提出一种快速时间双网格有限元算法求解分数阶水波模型。在[5]中，Li和Zhao运用差分格式求解具有空间分数阶导数的水波模型。Wang等学者在[6]中利用LDG (局部间断有限元)方法对非局部粘性水波模型展开研究。本文受该篇文章启发，计划在空间上运用间断有限元方法进行研究。

间断有限元(DG)方法于1973年由Reed和Hill在[7]中首次提出，随后学者们不断探索其改进形式以克服原始方法在热传导方程求解时的不足，相继发展了LDG [8]、Baumann-Oden DG [9]、dGRP (扩散广义黎曼问题) DG [10]、DDG [11]、UWDG (超弱间断有限元) [12]等方法。其中，DDG方法在处理高阶偏微分方程时展现出显著优势：(1) 无需引入辅助变量，计算效率高；(2) 仅需单次分部积分，格式构造简洁；(3) 保留了传统DG方法处理复杂边界条件和高精度求解的特性。Yan和Liu在文献[11]中首次将DDG方法应用于热传导方程求解，Yi等学者在[13]中进一步发展了能量守恒的DDG格式用于求解广义KdV方程。近年来，用DDG方法数值求解含有分数阶导数的偏微分方程成为研究热点：Huang等学者在[14]中发展了空间DDG方法，时间GMMP格式的混合算法求解含扩散项的时间分数阶方程；Liao等学者在[15]中结合L1近似公式建立了空间DDG方法误差估计。值得注意的是，同时包含扩散项与色散项的时间分数阶水波模型尚未有系统研究。

本文计划利用DDG方法与L1近似公式构造非局部粘性水波模型的数值算法。首先，对该模型弱形式中空间导数项进行一次分部积分后，在边界处引入稳定项，并选取合适的数值通量以保证格式的稳定性，特别地，对于含有空间导数的非线性项采用非线性迭代策略进行高效处理；时间导数离散运用BDF2方法与L1近似公式，从而构建全离散数值格式。其次，详细阐述算法执行流程及计算实现细节，并进行稳定性分析。最后，通过数值算例验证方法的可行性与收敛性，结果表明该算法在保持高精度的同时具备良好的计算效率。本文结构安排如下：第2节系统构造非局部水波模型的全离散格式；第3节给出详细的数值算法实现步骤；第4节进行稳定性分析；第5节通过典型算例验证算法有效性；第6节总结研究成果。

2. 全离散格式

在本节中，我们将采用L1近似公式以及BDF2方法对时间变量进行离散，同时在空间维度上运用DDG方法，构建全离散的数值格式。

下面对本文将要使用的一些符号进行说明。给定区间$\overline{\Omega }=\left[a,b\right]$ ，将区间分成N份，区间端点分别记为$a=x_{1/2}<x_{3/2}<\cdots <x_{N+1/2}=b$ ，区间中点记为$x_j=\frac{1}{2}\left(x_{j-1/2}+x_{j+1/2}\right)$ 。定义$h=\underset{1\le j\le N}{\max }{h}_j$ ，$u_{j+1/2}^{+}=\underset{x\to x_{j+1/2}^{+}}{\lim }u\left(x\right)$ ，$u_{j+1/2}^{-}=\underset{x\to x_{j+1/2}^{-}}{\lim }u\left(x\right)$ ，其中$1\le j\le N$ 。$\left[u\right]=u^{+}-u^{-}$ 代表u的跃度，$\overline{u}=\frac{1}{2}\left(u^{+}+u^{-}\right)$ 代表u的均值。

为了在时间方向上进行离散，将时间区间等分成M份，其中$0=t_0<t_1<\cdots <t_M=T$ 。时间步长$\tau =T/M$ 。为了简便用$u^{n+1}$ 代表$u$ 在$\left(x,t_n\right)$ 处的值。给出Caputo型时间分数阶导数分数阶为0.5阶的L1近似公式

${}_C{D}_{0,t}^{\frac{1}{2}}u=\frac{1}{\text{Γ}\left(1/2\right)}{\displaystyle {\int }_0^t\frac{\partial u\left(s\right)}{\partial s}}\frac{\text{d}s}{\sqrt{t-s}}=\frac{1}{{\tau }^{1/2}\text{Γ}\left(3/2\right)}\left(u^{n+1}-{\displaystyle \sum }_{i=1}^n\left(b_{i-1}-b_i\right)-b_n{u}^0\right)+R_0^{n+1}.$ (2)

其中$R_0^{n+1}=O\left({\tau }^{3/2}\right)$ ，$b_i={\left(i+1\right)}^{1/2}-i^{1/2}$ 。为了简便引入如下的符号：

${\delta }_t{u}^{n+1}=\{\begin{array}{l}\frac{u^{n+1}-u^n}{\tau },n=0,\\ \frac{3u^{n+1}-4u^n+u^{n-1}}{2\tau },n\ge 1.\end{array}$ (3)

则(1)在$t_{n+1}$ 处的弱形式为

当$n=0$ 时

$\begin{array}{l}\left({\delta }_t{u}^1,v\right)+g\left(u^1,v\right)+{\displaystyle \sum }_{j=1}^N\left({\left(u^1{v}^{-}\right)}_{j+\frac{1}{2}}-{\left(u^1{v}^{+}\right)}_{j-\frac{1}{2}}\right)-\left(u^1,v_x\right)\\ +\alpha \left({\displaystyle \sum }_{j=1}^N\left({\left(f\left(u^1\right)v^{-}\right)}_{j+\frac{1}{2}}-{\left(f\left(u^1\right)v^{+}\right)}_{j-\frac{1}{2}}\right)-\left(f\left(u^1\right),v_x\right)\right)\\ +\beta \left({\displaystyle \sum }_{j=1}^N\left({\left(u_{xx}^1{v}^{-}\right)}_{j+\frac{1}{2}}-{\left(u_{xx}^1{v}^{+}\right)}_{j-\frac{1}{2}}\right)-\left(u_{xx}^1,v_x\right)\right)\\ -\gamma \left({\displaystyle \sum }_{j=1}^N\left({\left(u_x^1{v}^{-}\right)}_{j+\frac{1}{2}}-{\left(u_x^1{v}^{+}\right)}_{j-\frac{1}{2}}\right)+\left(u_x^1,v_x\right)\right)\\ =g\left(u^0,v\right)+\left(R_0^1,v\right)+\left(R_1^1,v\right)\text{.}\end{array}$ (4)

当$n\ge 1$ 时

$\left({\delta }_t{u}^{n+1},v\right)+g\left(u^{n+1},v\right)+{\displaystyle \sum }_{j=1}^N\left({\left(u^{n+1}{v}^{-}\right)}_{j+\frac{1}{2}}-{\left(u^{n+1}{v}^{+}\right)}_{j-\frac{1}{2}}\right)-\left(u^{n+1},v_x\right)$

$\begin{array}{l}+\alpha \left({\displaystyle \sum }_{j=1}^N\left({\left(f\left(u^{n+1}\right)v^{-}\right)}_{j+\frac{1}{2}}-{\left(f\left(u^{n+1}\right)v^{+}\right)}_{j-\frac{1}{2}}\right)-\left(f\left(u^{n+1}\right),v_x\right)\right)\\ +\beta \left({\displaystyle \sum }_{j=1}^N\left({\left(u_{xx}^{n+1}{v}^{-}\right)}_{j+\frac{1}{2}}-{\left(u_{xx}^{n+1}{v}^{+}\right)}_{j-\frac{1}{2}}\right)-\left(u_{xx}^{n+1},v_x\right)\right)\\ -\gamma \left({\displaystyle \sum }_{j=1}^N\left({\left(u_x^{n+1}{v}^{-}\right)}_{j+\frac{1}{2}}-{\left(u_x^{n+1}{v}^{+}\right)}_{j-\frac{1}{2}}\right)+\left(u_x^{n+1},v_x\right)\right)\\ =g{b}_n\left(u^0,v\right)+g\left({\displaystyle \sum }_{i=1}^n\left(b_{i-1}-b_i\right)\left(u^{n+1-i},v\right)\right)+\left(R_0^{n+1},v\right)+\left(R_{n+1}^2,v\right).\end{array}$ (5)

其中$g=\frac{\sqrt{\upsilon }}{{\tau }^{1/2}\text{Γ}\left(3/2\right)}$ ，并且

$\begin{array}{l}{R}_0^{n+1}={}_C{D}_{0,t}^{1/2}{u}^{n+1}-\frac{1}{{\tau }^{1/2}\text{Γ}\left(3/2\right)}\left(u^{n+1}-{\displaystyle \sum }_{i=1}^n\left(b_{i-1}-b_i\right)-b_n{u}^0\right)=O\left({\tau }^{3/2}\right),\\ R_1^1={\delta }_t{u}^1-u_t^1=O\left(\tau \right),R_2^{n+1}={\delta }_t{u}^{n+1}-u_t^{n+1}=O\left({\tau }^2\right)\text{.}\end{array}$ (6)

为了进一步得到全离散数值格式，定义有限元空间$V_h^k=\left\{v\in L^2\left(\overline{\Omega }\right):{V|}_{I_j}\in p^k\left(I_j\right),j=1,\cdots ,N\right\}$ 。设$u_h\in V_h^k$ 是$u$ 的近似解。为了保证稳定性，需要对色散项增加额外的边界项，故全离散格式为

当$n=0$ 时

$\begin{array}{l}\left({\delta }_t{u}_h^1,v_h\right)+g\left(u_h^1,v_h\right)+{\displaystyle \sum }_{j=1}^N\left[{\left({\stackrel{⌣}{u}}_h^1{v}_h^{-}\right)}_{j+1/2}-{\left({\stackrel{⌣}{u}}_h^1{v}_h^{+}\right)}_{j-1/2}\right]-\left(u_h^1,v_{hx}\right)\\ +\alpha \left[{\displaystyle \sum }_{j=1}^N\left({\left(\widehat{f}\left(u_h^1\right)v_h^{-}\right)}_{j+1/2}-{\left(\widehat{f}\left(u_h^1\right)v_h^{+}\right)}_{j-1/2}\right)-\left(f\left(u_h^1\right),v_{hx}\right)\right]\\ +\beta [{\displaystyle \sum }_{j=1}^N\left({\left({\left(u_h^1\right)}_{xx}{v}_h^{-}\right)}_{j+1/2}-{\left({\left(u_h^1\right)}_{xx}{v}_h^{+}\right)}_{j-1/2}\right)-\left({\left(u_h^1\right)}_{xx},v_{hx}\right)\\ +{\displaystyle \sum }_{j=1}^N\left({\left(\left({\left(u_h^1\right)}_x-{\left({\widehat{u}}_h^1\right)}_x\right){\left(v_h\right)}_x^{-}\right)}_{j+1/2}-{\left(\left({\left(u_h^1\right)}_x-{\left({\widehat{u}}_h^1\right)}_x\right){\left(v_h\right)}_x^{-}\right)}_{j-1/2}\right)\\ +{\displaystyle \sum }_{j=1}^N\left({\left(\left(\left({\widehat{u}}_h^1\right)-\left(u_h^1\right)\right){\left(v_h\right)}_{xx}^{-}\right)}_{j+\frac{1}{2}}-{\left(\left(\left({\widehat{u}}_h^1\right)-\left(u_h^1\right)\right){\left(v_h\right)}_{xx}^{-}\right)}_{j-\frac{1}{2}}\right)]\\ -\gamma \left[{\displaystyle \sum }_{j=1}^N{\left(\widetilde{u_{hx}^1}{\left(u_h^1\right)}^{-}\right)}_{j+\frac{1}{2}}-{\left(\widetilde{u_{hx}^1}{\left(u_h^1\right)}^{+}\right)}_{j-\frac{1}{2}}+\left(u_{hx}^1,u_{hx}^1\right)-{\left(\left[u_h^1\right]\widetilde{v_x}\right)}_{j+\frac{1}{2}}-{\left(\left[u_h^1\right]\widetilde{v_x}\right)}_{j-\frac{1}{2}}\right].\end{array}$ (7)

当$n\ge 1$ 时

$\begin{array}{l}\left({\delta }_t{u}_h^{n+1},v_h\right)+g\left(u_h^{n+1},v_h\right)+{\displaystyle \sum }_{j=1}^N\left[{\left({\stackrel{⌣}{u}}_h^{n+1}{v}_h^{-}\right)}_{j+\frac{1}{2}}-{\left({\stackrel{⌣}{u}}_h^{n+1}{v}_h^{+}\right)}_{j-\frac{1}{2}}\right]-\left(u_h^{n+1},v_{hx}\right)\\ +\alpha \left[{\displaystyle \sum }_{j=1}^N\left({\left(\widehat{f}\left(u_h^{n+1}\right)v_h^{-}\right)}_{j+\frac{1}{2}}-{\left(\widehat{f}\left(u_h^{n+1}\right)v_h^{+}\right)}_{j-\frac{1}{2}}\right)-\left(f\left(u_h^{n+1}\right),v_{hx}\right)\right]\\ +\beta [{\displaystyle \sum }_{j=1}^N\left({\left({\left(u_h^{n+1}\right)}_{xx}{v}_h^{-}\right)}_{j+\frac{1}{2}}-{\left({\left(u_h^{n+1}\right)}_{xx}{v}_h^{+}\right)}_{j-\frac{1}{2}}\right)-\left({\left(u_h^{n+1}\right)}_{xx},v_{hx}\right)\\ +{\displaystyle \sum }_{j=1}^N\left({\left(\left({\left(u_h^{n+1}\right)}_x-{\left({\widehat{u}}_h^{n+1}\right)}_x\right){\left(v_h\right)}_x^{-}\right)}_{j+\frac{1}{2}}-{\left(\left({\left(u_h^{n+1}\right)}_x-{\left({\widehat{u}}_h^{n+1}\right)}_x\right){\left(v_h\right)}_x^{-}\right)}_{j-\frac{1}{2}}\right)\end{array}$

$\begin{array}{l}+{\displaystyle \sum }_{j=1}^N\left({\left(\left({\widehat{u}}_h^{n+1}-u_h^{n+1}\right){\left(v_h\right)}_{xx}^{-}\right)}_{j+\frac{1}{2}}-{\left(\left({\widehat{u}}_h^{n+1}-u_h^{n+1}\right){\left(v_h\right)}_{xx}^{-}\right)}_{j-\frac{1}{2}}\right)]\\ -\gamma \left[{\displaystyle \sum }_{j=1}^N\left({\left({\left(\widetilde{u_h^{n+1}}\right)}_x{v}_h^{-}\right)}_{j+\frac{1}{2}}-{\left({\left(\widetilde{u_h^{n+1}}\right)}_x{v}_h^{+}\right)}_{j-\frac{1}{2}}+{\left(\left[u_h^{n+1}\right]\widetilde{v_x}\right)}_{j+\frac{1}{2}}+{\left(\left[u_h^{n+1}\right]\widetilde{v_x}\right)}_{j-\frac{1}{2}}\right)+\left({\left(u_h^{n+1}\right)}_x,v_{hx}\right)\right]\\ =g{b}_n\left(u_h^0,v_h\right)+g\left[{\displaystyle \sum }_{i=1}^n\left(b_{i-1}-b_i\right)\left(u_h^{n+1-i},v_h\right)\right]\text{.}\end{array}$ (8)

(7)~(8)中数值通量选取如下[11] [13] [16]

$\begin{array}{l}{\stackrel{⌣}{u}}_h^{n+1}={\left(u_h^{n+1}\right)}^{-},{\widehat{u}}_h^{n+1}=\theta {\left(u_h^{n+1}\right)}^{-}+\left(1-\theta \right){\left(u_h^{n+1}\right)}^{+};\\ {\left({\widehat{u}}_h^{n+1}\right)}_x=\left(1-\eta \right){\left(u_h^{n+1}\right)}_x^{-}+\eta {\left(u_h^{n+1}\right)}_x^{+},\eta =\frac{1}{2};{\left({\widehat{u}}_h^{n+1}\right)}_{xx}=\left(1-\theta \right){\left(u_h^{n+1}\right)}_x^{-}+\theta {\left(u_h^{n+1}\right)}_x^{+};\\ \widehat{f}\left(u_h^{n+1}\right)=\widehat{f}\left({\left(u_h^{n+1}\right)}^{-},{\left(u_h^{n+1}\right)}^{+}\right),{\left(\widetilde{u_h^{n+1}}\right)}_x={\beta }_0\frac{\left[u_h^{n+1}\right]}{h}+\overline{{\left(u_h^{n+1}\right)}_x}+{\beta }_{1}h\left[{\left(u_h^{n+1}\right)}_{xx}\right],\end{array}$ (9)

其中$\theta \in \left[0,1\right]$ ，$\widehat{f}\left({\left(u_h^{n+1}\right)}^{-},{\left(u_h^{n+1}\right)}^{+}\right)$ 是一个与$f\left(u_h^{n+1}\right)$ 一致的单调数值通量，并且该通量关于${\left(u_h^{n+1}\right)}^{-},{\left(u_h^{n+1}\right)}^{+}$ 是Lipschitz连续的。例如对于特殊的非线性项$f\left(\omega \right)={\omega }^{p+1}$ 使用平方最优通量

$\widehat{f}\left({\omega }^{-},{\omega }^{+}\right)=\frac{1}{p+2}{\displaystyle \sum }_{j=0}^{p+1}{\left({\omega }^{-}\right)}^j{\left({\omega }^{+}\right)}^{p+1-j}\text{.}$ (10)

注1：上述数值通量中，非线性项还可选用Lax-Friedrich通量。经数值实验验证，该通量的选用不会对方法的精度造成影响。

3. 算法实现

选取有限元空间为$V_h^k$ ，在第j个区间上有$u_{h,j}\left(x,t_n\right)={\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_n\right){\varphi }_{i,j}$ ，将其代入到(7)~(8)中并且令$v_h={\varphi }_{m,j}$ ，可得以下全离散格式

当$n=0$ 时

$\begin{array}{l}\left({\delta }_t{\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_1\right){\varphi }_{i,j},{\varphi }_{m,j}\right)+g\left({\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_1\right){\varphi }_{i,j},{\varphi }_{m,j}\right)\\ +\left[{\left({\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_1\right){\varphi }_{i,j}{\varphi }_{m,j}^{-}\right)}_{j+1/2}-{\left({\displaystyle \sum }_{i=0}^k{u}_{i,j-1}\left(t_1\right){\varphi }_{i,j-1}{\varphi }_{m,j}^{+}\right)}_{j-1/2}\right]-\left({\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_1\right){\varphi }_{i,j},{\varphi }_{m,j,x}\right)\\ +\alpha {\text{H}}_j^1\left(u_{h,j},{\varphi }_{m,j}\right)+\beta {\text{B}}_j^1\left(u_{h,j},{\varphi }_{m,j}\right)-\gamma {\text{G}}_j^1\left(u_{h,j},{\varphi }_{m,j}\right)=g{b}_n\left({\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_0\right){\varphi }_{i,j},{\varphi }_{m,j}\right),\end{array}$ (11)

当$n\ge 1$ 时

$\begin{array}{l}\left({\delta }_t{\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_{n+1}\right){\varphi }_{i,j},{\varphi }_{i,j}\right)+g\left({\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_{n+1}\right){\varphi }_{i,j},{\varphi }_{i,j}\right)\\ +\left[{\left({\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_{n+1}\right){\varphi }_{i,j}{\varphi }_{i,j}^{-}\right)}_{j+1/2}-{\left({\displaystyle \sum }_{i=0}^k{u}_{i,j-1}\left(t_{n+1}\right){\varphi }_{i,j-1}{\varphi }_{i,j}^{+}\right)}_{j-1/2}\right]-\left({\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_{n+1}\right){\varphi }_{i,j},{\varphi }_{i,j,x}\right)\\ +\alpha {\text{H}}_j^{n+1}\left(u_{h,j},{\varphi }_{m,j}\right)+\beta {\text{B}}_j^{n+1}\left(u_{h,j},{\varphi }_{m,j}\right)-\gamma {\text{G}}_j^{n+1}\left(u_{h,j},{\varphi }_{m,j}\right)\\ =g{b}_n\left({\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_0\right){\varphi }_{i,j},{\varphi }_{i,j}\right)+g\left[{\displaystyle \sum }_{i=1}^n\left(b_{i-1}-b_i\right)\left({\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_{n+1-i}\right){\varphi }_{i,j,x},{\varphi }_{i,j}\right)\right],\end{array}$ (12)

其中，

$\begin{array}{c}{\text{B}}_j^n\left(u_{h,j},{\varphi }_{m,j}\right)=\left({\left({\left({\displaystyle \sum }_{i=0}^k\widehat{u_{i,j}\left(t_n\right)}{\varphi }_{i,j}\right)}_{xx}{\varphi }_{m,j}^{-}\right)}_{j+1/2}-{\left({\left({\displaystyle \sum }_{i=0}^k\widehat{u_{i,j}\left(t_n\right)}{\varphi }_{i,j}\right)}_{xx}{\varphi }_{m,j}^{+}\right)}_{j-1/2}\right)\\ -\left({\left({\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_n\right){\varphi }_{i,j}\right)}_{xx},{\varphi }_{m,j,x}\right)+{\left(\left({\left({\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_n\right){\varphi }_{i,j}\right)}_x-{\left({\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_n\right){\varphi }_{i,j}\right)}_x\right){\left({\varphi }_{m,j}\right)}_x^{-}\right)}_{j+1/2}\\ -{\left(\left({\left({\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_n\right){\varphi }_{i,j}\right)}_x-{\left({\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_n\right){\varphi }_{i,j}\right)}_x\right){\left({\varphi }_{m,j}\right)}_x^{-}\right)}_{j-1/2}\\ +{\left[\left(\left({\displaystyle \sum }_{i=0}^k\widehat{u_{i,j}\left(t_n\right)}{\varphi }_{i,j}\right)-\left({\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_n\right){\varphi }_{i,j}\right)\right){\left({\varphi }_{m,j}\right)}_{xx}^{-}\right]}_{j+1/2}\\ -{\left[\left(\left({\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_n\right){\varphi }_{i,j}\right)-\left({\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_n\right){\varphi }_{i,j}\right)\right){\left({\varphi }_{m,j}\right)}_{xx}^{-}\right]}_{j-1/2},\end{array}$

$\begin{array}{c}{\text{H}}_j^n\left(u_{h,j},{\varphi }_{m,j}\right)=\left({\left(\widehat{f}\left({\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_n\right){\varphi }_{i,j}\right){\varphi }_{m,j}^{-}\right)}_{j+1/2}-{\left(\widehat{f}\left({\displaystyle \sum }_{i=0}^k{u}_{i,j}\left(t_n\right){\varphi }_{i,j}\right){\varphi }_{m,j}^{+}\right)}_{j-1/2}\right)\\ -\left(f\left({\displaystyle \sum }_{i=0}^k\left(u_{i,j}\left(t_n\right){\varphi }_{i,j}\right)\right),{\varphi }_{m,j,x}\right),\end{array}$

$\begin{array}{c}{\text{G}}_j^n\left(u_{h,j},{\varphi }_{m,j}\right)=\left({\left({\left({\displaystyle \sum _{i=0}^k\widetilde{u_{i,j}\left(t_n\right)}}{\varphi }_{i,j}\right)}_x{\varphi }_{m,j}^{-}\right)}_{j+1/2}-{\left({\left({\displaystyle \sum _{i=0}^k\widetilde{u_{i,j}\left(t_n\right)}}{\varphi }_{i,j}\right)}_x{\varphi }_{m,j}^{+}\right)}_{j-1/2}\right)\\ +\left({\left({\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_i\right)}_x,{\varphi }_{m,j,x}\right).\end{array}$

为了构建整个区间上全离散格式的矩阵形式需要在每个区间构建如下的单元矩阵

$A_j={\left({\varphi }_{i,j},{\varphi }_{m,j}\right)}_{\left(k+1\right)\times \left(k+1\right)},{\left(U_{xx}{V}_x\right)}_j={\left({\varphi }_{i,j,xx},{\varphi }_{i,j,x}\right)}_{\left(k+1\right)\times \left(k+1\right)},$

${\left(U{V}_x\right)}_j={\left({\varphi }_{i,j},{\varphi }_{i,j,x}\right)}_{\left(k+1\right)\times \left(k+1\right)},{\left(U_x{V}_x\right)}_j={\left({\varphi }_{i,j,x},{\varphi }_{i,j,x}\right)}_{\left(k+1\right)\times \left(k+1\right)},$

$TU{V}_{xx,1,j}=\left(\begin{array}{c}-\left(1-\theta \right){\displaystyle \sum _{i=0}^k{u}_{i,j-1}}\left(t_n\right){\varphi }_{i,j-1}\left(x_{j-1/2}\right){\varphi }_{0,j,xx}\left(x_{j-1/2}\right)\\ ⋮\\ -\left(1-\theta \right){\displaystyle \sum _{i=0}^k{u}_{i,j-1}}\left(t_n\right){\varphi }_{i,j-1}\left(x_{j-1/2}\right){\varphi }_{k,j,xx}\left(x_{j-1/2}\right)\end{array}\right),$

$TU{V}_{xx,2,j}=\left(\begin{array}{c}\left(1-\theta \right){\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j-1}\left(x_{j-1/2}\right){\varphi }_{0,j}\left(x_{j+1/2}\right)-\theta {\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j-1}\left(x_{j-1/2}\right){\varphi }_{0,j,xx}\left(x_{j-1/2}\right)\\ ⋮\\ \left(1-\theta \right){\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j-1}\left(x_{j-1/2}\right){\varphi }_{k,j}\left(x_{j+1/2}\right)-\theta {\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j-1}\left(x_{j-1/2}\right){\varphi }_{k,j,xx}\left(x_{j-1/2}\right)\end{array}\right),$

$TU{V}_{xx,3,j}=\left(\begin{array}{c}\theta {\displaystyle \sum _{i=0}^k{u}_{i,j+1}}\left(t_n\right){\varphi }_{i,j+1}\left(x_{j+1/2}\right){\varphi }_{0,j,xx}\left(x_{j+1/2}\right)\\ ⋮\\ \theta {\displaystyle \sum _{i=0}^k{u}_{i,j+1}}\left(t_n\right){\varphi }_{i,j+1}\left(x_{j+1/2}\right){\varphi }_{k,j,xx}\left(x_{j+1/2}\right)\end{array}\right),TU{V}_{xx,j}=\left[TU{V}_{xx,1,j}TU{V}_{xx,2,j}TU{V}_{xx,3,j}\right],$

$T{U}_x{V}_{x,1,j}=\left(\begin{array}{c}\left(1-\eta \right){\displaystyle \sum _{i=0}^k{u}_{i,j-1}}\left(t_n\right){\varphi }_{i,j-1,x}\left(x_{j-1/2}\right){\varphi }_{0,j,x}\left(x_{j-1/2}\right)\\ ⋮\\ \left(1-\eta \right){\displaystyle \sum _{i=0}^k{u}_{i,j-1}}\left(t_n\right){\varphi }_{i,j-1,x}\left(x_{j-1/2}\right){\varphi }_{k,j,x}\left(x_{j-1/2}\right)\end{array}\right),$

$T{U}_x{V}_{x,2,j}=\left(\begin{array}{c}-\left(1-\eta \right){\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j-1,x}\left(x_{j-1/2}\right){\varphi }_{0,j,x}\left(x_{j+1/2}\right)+\eta {\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j-1,x}\left(x_{j-1/2}\right){\varphi }_{0,j,x}\left(x_{j-1/2}\right)\\ ⋮\\ -\left(1-\eta \right){\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j-1,x}\left(x_{j-1/2}\right){\varphi }_{k,j,x}\left(x_{j+1/2}\right)+\eta {\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j-1,x}\left(x_{j-1/2}\right){\varphi }_{k,j,x}\left(x_{j-1/2}\right)\end{array}\right),$

$T{U}_x{V}_{x,3,j}=\left(\begin{array}{c}-\eta {\displaystyle \sum _{i=0}^k{u}_{i,j+1}}\left(t_n\right){\varphi }_{i,j+1,x}\left(x_{j+1/2}\right){\varphi }_{0,j,x}\left(x_{j+1/2}\right)\\ ⋮\\ +\eta {\displaystyle \sum _{i=0}^k{u}_{i,j+1}}\left(t_n\right){\varphi }_{i,j+1,x}\left(x_{j+1/2}\right){\varphi }_{k,j,x}\left(x_{j+1/2}\right)\end{array}\right),T{U}_x{V}_{x,j}=\left[T{U}_x{V}_{x,1,j}T{U}_x{V}_{x,2,j}T{U}_x{V}_{x,3,j}\right],$

$T{U}_{xx}{V}_{1,j}=\left(\begin{array}{c}-\theta {\displaystyle \sum _{i=0}^k{u}_{i,j-1}}\left(t_n\right){\varphi }_{i,j-1,xx}\left(x_{j-1/2}\right){\varphi }_{0,j}\left(x_{j-1/2}\right)\\ ⋮\\ -\theta {\displaystyle \sum _{i=0}^k{u}_{i,j-1}}\left(t_n\right){\varphi }_{i,j-1,xx}\left(x_{j-1/2}\right){\varphi }_{k,j}\left(x_{j-1/2}\right)\end{array}\right),$

$T{U}_{xx}{V}_{2,j}=\left(\begin{array}{c}\theta {\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j,xx}\left(x_{j-1/2}\right){\varphi }_{0,j}\left(x_{j+1/2}\right)-\left(1-\eta \right){\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j,xx}\left(x_{j-1/2}\right){\varphi }_{0,j}\left(x_{j-1/2}\right)\\ ⋮\\ \theta {\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j,xx}\left(x_{j-1/2}\right){\varphi }_{k,j}\left(x_{j+1/2}\right)-\left(1-\eta \right){\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j,xx}\left(x_{j-1/2}\right){\varphi }_{k,j}\left(x_{j-1/2}\right)\end{array}\right),$

$T{U}_{xx}{V}_{3,j}=\left(\begin{array}{c}\left(1-\eta \right){\displaystyle \sum _{i=0}^k{u}_{i,j+1}}\left(t_n\right){\varphi }_{i,j+1,xx}\left(x_{j+1/2}\right){\varphi }_{0,j}\left(x_{j+1/2}\right)\\ ⋮\\ \left(1-\eta \right){\displaystyle \sum _{i=0}^k{u}_{i,j+1}}\left(t_n\right){\varphi }_{i,j+1,xx}\left(x_{j+1/2}\right){\varphi }_{k,j}\left(x_{j+1/2}\right)\end{array}\right),$

$T{U}_{xx}{V}_{x,j}=\left[T{U}_{xx}{V}_{1,j}T{U}_{xx}{V}_{2,j}T{U}_{xx}{V}_{3,j}\right]$

$TU{V}_{1,j}=\left(\begin{array}{c}-{\displaystyle \sum _{i=0}^k{u}_{i,j-1}}\left(t_n\right){\varphi }_{i,j-1}\left(x_{j-1/2}\right){\varphi }_{0,j}\left(x_{j-1/2}\right)\\ ⋮\\ -{\displaystyle \sum _{i=0}^k{u}_{i,j-1}}\left(t_n\right){\varphi }_{i,j-1}\left(x_{j-1/2}\right){\varphi }_{k,j}\left(x_{j-1/2}\right)\end{array}\right),TU{V}_{2,j}=\left(\begin{array}{c}{\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j}\left(x_{j-1/2}\right){\varphi }_{0,j}\left(x_{j+1/2}\right)\\ ⋮\\ {\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j}\left(x_{j-1/2}\right){\varphi }_{k,j}\left(x_{j+1/2}\right)\end{array}\right),$

$TU{V}_j=\left[\begin{array}{cc}TU{V}_{1,j}& TU{V}_{2,j}\end{array}\right],$

$T_{1,j}=\left(\begin{array}{c}\left(\frac{{\beta }_0}{h}{\displaystyle \sum _{i=0}^k{u}_{i,j-1}}\left(t_n\right){\varphi }_{i,j-1}\left(x_{j-1/2}\right)-\frac{1}{2}{\displaystyle \sum _{i=0}^k{u}_{i,j-1}}\left(t_n\right){\varphi }_{i,j-1,x}\left(x_{j-1/2}\right)+{\beta }_1{\displaystyle \sum _{i=0}^k{u}_{i,j-1}}\left(t_n\right){\varphi }_{i,j-1,xx}\left(x_{j-1/2}\right)\right){\varphi }_{0,j}\left(x_{j-1/2}\right)\\ ⋮\\ \left(\frac{{\beta }_0}{h}{\displaystyle \sum _{i=0}^k{u}_{i,j-1}}\left(t_n\right){\varphi }_{i,j-1}\left(x_{j-1/2}\right)-\frac{1}{2}{\displaystyle \sum _{i=0}^k{u}_{i,j-1}}\left(t_n\right){\varphi }_{i,j-1,x}\left(x_{j-1/2}\right)+{\beta }_1{\displaystyle \sum _{i=0}^k{u}_{i,j-1}}\left(t_n\right){\varphi }_{i,j-1,xx}\left(x_{j-1/2}\right)\right){\varphi }_{k,j}\left(x_{j-1/2}\right)\end{array}\right),$

$T_{2,j}=\left(\begin{array}{c}\left(-\frac{{\beta }_0}{h}{\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j}\left(x_{j-1/2}\right)+\frac{1}{2}{\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j,x}\left(x_{j-1/2}\right)-{\beta }_1{\displaystyle \sum _{i=0}^k{u}_{i,j-1}}\left(t_n\right){\varphi }_{i,j,xx}\left(x_{j-1/2}\right)\right){\varphi }_{0,j}\left(x_{j-1/2}\right)\\ -\left(\frac{{\beta }_0}{h}{\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j}\left(x_{j-1/2}\right)+\frac{1}{2}{\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j,x}\left(x_{j-1/2}\right)+{\beta }_1{\displaystyle \sum _{i=0}^k{u}_{i,j-1}}\left(t_n\right){\varphi }_{i,j,xx}\left(x_{j-1/2}\right)\right){\varphi }_{0,j}\left(x_{j+1/2}\right)\\ ⋮\\ \left(-\frac{{\beta }_0}{h}{\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j}\left(x_{j-1/2}\right)+\frac{1}{2}{\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j,x}\left(x_{j-1/2}\right)-{\beta }_1{\displaystyle \sum _{i=0}^k{u}_{i,j-1}}\left(t_n\right){\varphi }_{i,j-1,xx}\left(x_{j-1/2}\right)\right){\varphi }_{k,j}\left(x_{j-1/2}\right)\\ -\left(\frac{{\beta }_0}{h}{\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j}\left(x_{j-1/2}\right)+\frac{1}{2}{\displaystyle \sum _{i=0}^k{u}_{i,j}}\left(t_n\right){\varphi }_{i,j,x}\left(x_{j-1/2}\right)+{\beta }_1{\displaystyle \sum _{i=0}^k{u}_{i,j-1}}\left(t_n\right){\varphi }_{i,j,xx}\left(x_{j-1/2}\right)\right){\varphi }_{k,j}\left(x_{j+1/2}\right)\end{array}\right),$

$T_{3,j}=\left(\begin{array}{c}\left(\frac{{\beta }_0}{h}{\displaystyle \sum _{i=0}^k{u}_{i,j+1}}\left(t_n\right){\varphi }_{i,j+1}\left(x_{j+1/2}\right)+\frac{1}{2}{\displaystyle \sum _{i=0}^k{u}_{i,j+1}}\left(t_n\right){\varphi }_{i,j+1,x}\left(x_{j+1/2}\right)+{\beta }_1{\displaystyle \sum _{i=0}^k{u}_{i,j-1}}\left(t_n\right){\varphi }_{i,j+1,xx}\left(x_{j+1/2}\right)\right){\varphi }_{0,j}\left(x_{j+1/2}\right)\\ ⋮\\ \left(\frac{{\beta }_0}{h}{\displaystyle \sum _{i=0}^k{u}_{i,j+1}}\left(t_n\right){\varphi }_{i,j-1}\left(x_{j+1/2}\right)+\frac{1}{2}{\displaystyle \sum _{i=0}^k{u}_{i,j-1}}\left(t_n\right){\varphi }_{i,j+1,x}\left(x_{j+1/2}\right)+{\beta }_1{\displaystyle \sum _{i=0}^k{u}_{i,j+1}}\left(t_n\right){\varphi }_{i,j+1,xx}\left(x_{j+1/2}\right)\right){\varphi }_{k,j}\left(x_{j+1/2}\right)\end{array}\right),$

$T_j=\left[\begin{array}{ccc}{T}_{1,j}& T_{2,j}& T_{3,j}\end{array}\right],$

$F\left(u_{h,j}\left(t\right)\right)=\left(\begin{array}{c}\left(f\left({\displaystyle \sum _{i=0}^k{M}_i}{u}_{i,j}\left(t\right){\varphi }_i\right),{\varphi }_{0,x}\right)\\ ⋮\\ \left(f\left({\displaystyle \sum _{i=0}^k{M}_i}{u}_{i,j}\left(t\right){\varphi }_i\right),{\varphi }_{k,x}\right)\end{array}\right).$

根据周期边界条件构建整个区间上的矩阵，以下矩阵的规格均为$\left(k+1\right)N\times \left(k+1\right)N$

$A=\left[\begin{array}{ccc}{A}_1& & \\ & \ddots & \\ & & A_N\end{array}\right],U_{xx}{V}_x=\left[\begin{array}{ccc}{\left(U_{xx}{V}_x\right)}_1& & \\ & \ddots & \\ & & {\left(U_x{V}_x\right)}_N\end{array}\right],U_x{V}_x=\left[\begin{array}{ccc}{\left(U_x{V}_x\right)}_1& & \\ & \ddots & \\ & & {\left(U_x{V}_x\right)}_N\end{array}\right],$

$TU{V}_{xx}=\left[\begin{array}{cccc}TU{V}_{xx,2,1}& TU{V}_{xx,3,1}& \dots & TU{V}_{xx,1,1}\\ TU{V}_{xx,1,1}& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & TU{V}_{xx,3,N}\\ TU{V}_{xx,3,N}& \dots & TU{V}_{xx,1,N}& TU{V}_{xx,2,N}\end{array}\right],T{U}_x{V}_x=\left[\begin{array}{cccc}T{U}_{xx}{V}_{2,1}& T{U}_{xx}{V}_{3,1}& \dots & T{U}_{xx}{V}_{1,1}\\ T{U}_{xx}{V}_{1,2}& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & T{U}_{xx}{V}_{3,N}\\ T{U}_{xx}{V}_{3,N}& \dots & T{U}_{xx}{V}_{1,N}& T{U}_{xx}{V}_{2,N}\end{array}\right],$

$T{U}_{xx}V=\left[\begin{array}{cccc}T{U}_{xx}{V}_{2,1}& T{U}_{xx}{V}_{3,1}& \dots & T{U}_{xx}{V}_{1,1}\\ T{U}_{xx}{V}_{1,2}& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & T{U}_{xx}{V}_{3,N}\\ T{U}_{xx}{V}_{3,N}& \dots & T{U}_{xx}{V}_{1,N}& T{U}_{xx}{V}_{2,N}\end{array}\right],U{V}_x=\left[\begin{array}{ccc}{\left(U{V}_x\right)}_1& & \\ & \ddots & \\ & & {\left(U{V}_x\right)}_N\end{array}\right],\widehat{F^n}=\left(\begin{array}{c}\widehat{F}\left(u_{h,1}\left(t\right)\right)\\ ⋮\\ \widehat{F}\left(u_{h,N}\left(t\right)\right)\end{array}\right),$

$TUV=\left[\begin{array}{cccc}TU{V}_{2,1}& 0& \dots & T{U}_{xx}{V}_{1,1}\\ TU{V}_{1,2}& TU{V}_{2,2}& \ddots & ⋮\\ ⋮& \ddots & \ddots & 0\\ 0& \dots & T{U}_{xx}{V}_{1,N}& T{U}_{xx}{V}_{2,N}\end{array}\right],T=\left[\begin{array}{cccc}{T}_{2,1}& T_{3,1}& \dots & T{V}_{1,1}\\ T_{1,2}& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & T{V}_{3,N}\\ T_{3,N}& \dots & T_{1,N}& T_{2,N}\end{array}\right]，F^n=\left(\begin{array}{c}F\left(u_{h,1}\left(t_n\right)\right)\\ ⋮\\ F\left(u_{h,N}\left(t_n\right)\right)\end{array}\right).$

记$T=\beta \left(T{U}_{xx}V+T{U}_x{V}_x+TU{V}_{xx}-U_{xx}{V}_x\right)+TUV-U{V}_x-\alpha \left(T-U_x{V}_x\right)$ ，我们可以通过解以下非线性方程得到该模型的数值解

$\{\begin{array}{l}A\frac{u_h^1-u_h^0}{\tau }+A{u}_h^1+T{u}_h^1+{\widehat{F}}^1-F^1=g{b}_{n}A{u}_h^0,\hfill \\ A\frac{3u_h^{n+1}-4u_h^n+u_h^{n-1}}{2\tau }+A{u}_h^{n+1}+T{u}_h^{n+1}+{\widehat{F}}^{n+1}-F^{n+1}=g{b}_{n}A{u}_h^0+gA{\displaystyle \sum _{i=1}^n\left(b_{i-1}-b_i\right)u_h^{n+1-i}}.\hfill \end{array}$ (13)

4. 稳定性

在证明稳定性之前，引入一些算子和引理。

定义函数$\Delta {\displaystyle {\int }_{u^{-}}^{u^{+}}f\left(u\right)\text{d}u}=:\frac{{\displaystyle {\int }_{u^{-}}^{u^{+}}f}\left(u\right)\text{d}u}{\left[u\right]}$ 。引入变量$K\left({\alpha }_1\right)=\underset{v\in P^{k-1}\left[-1,1\right]}{\sup }\frac{{\left(v\left(1\right)-2{\alpha }_1{\partial }_{\xi }v\left(1\right)\right)}^2}{\Delta {\displaystyle {\int }_{-1}^1{v}^2\left(\xi \right)\text{d}\xi }}$ ，其中 $\xi =2\left(x-x_j\right)/h_j$ ，$K$ 是关于${\alpha }_1$ 和多项式次数$k$ 的函数。

定义1. [17]对于扩散项定义如下形式的算子

$A\left(\phi ,\psi \right)={\displaystyle \sum _{j=1}^N{\displaystyle {\int }_{I_j}{\phi }_x}}{\psi }_x\text{d}x+{\displaystyle \sum _{j=1}^N{\left(\widehat{{\phi }_x}\left[{\psi }_x\right]\right)}_{j+1/2}}+{\displaystyle \sum _{j=1}^N{\left(\widehat{{\psi }_x}\left[{\phi }_x\right]\right)}_{j+1/2}}.$

引理1. [17]本文数值通量的选择满足${\beta }_0>K\left({\beta }_1\right)$ ，故$\exists \zeta \in \left(0,1\right)$ 使得

$A\left(v,v\right)\ge \zeta {\Vert v\Vert }_E^2,\forall v\in V_h^k,$ (14)

其中${\Vert v_h\Vert }_E^2:={\displaystyle \sum _{j=1}^N{\displaystyle \underset{I_j}{\int }\left|{\left(v_h\right)}_x\right|\text{d}x}}+{\displaystyle \sum _{j=1}^N\frac{{\beta }_0}{h}}{\left[v_h\right]}^2{\mid }_{I_j}$ 。

引理2. 考虑特殊的非线性项，分部积分一次后为[16]

$\left({\left(u^{p+1}\right)}_x,v\right)=-\left(u^{p+1},v_x\right)-{\displaystyle \sum }_{j=1}^N{\left[u^{p+1}v\right]}_{j-1/2}\text{.}$

定义2. [16]定义非线性算子$\mathcal{N}:H^1\left(I_j\right)\to V_h^q$

$\left(\mathcal{N}\left(u\right),v\right)=-\left(f\left(u\right),v_x\right)-{\displaystyle \sum _{j=1}^N\widehat{f}}\left(u_{j-1/2}^{+},u_{j-1/2}^{-}\right){\left[v\right]}_{j-1/2}.$ (15)

引理3. [16]在(15)中定义的非线性算子，对于$\forall u\in C^1\left(\left[0,1\right]\right)$ 且边界条件为周期边界条件时，满足

(16)

定义3. [13]定义算子$D_{\theta }:H^3\left(I_j\right)\to V_h^q$

$\begin{array}{l}\left(D_{\theta }\left(u\right),v\right)=-\left(u_{xx},v_x\right)-{\displaystyle \sum }_{j=1}^N{\left[u_x\right]}_{j-\frac{1}{2}}{\left\{v_x\right\}}_{j-\frac{1}{2}}\\ -{\displaystyle \sum }_{j=1}^N\left({\left(\theta u_{xx}^{+}+\left(1-\theta \right)u_{xx}^{-}\right)}_{j-\frac{1}{2}}{\left[v\right]}_{j-\frac{1}{2}}-{\left[u\right]}_{j-\frac{1}{2}}{\left(\theta v_{xx}^{+}+\left(1-\theta \right)v_{xx}^{-}\right)}_{j-\frac{1}{2}}\right).\end{array}$ (17)

引理4. [16] [18]由(17)定义的算子与色散项是一致的且是斜伴随的，即

$\left({\mathcal{D}}_{\theta }\left(u\right),v\right)=\left(u_{xxx},v\right),\forall v\in H^3\left(I_j\right),$

$\left({\mathcal{D}}_{\theta }\left(u\right),u\right)=0,\forall v\in H^3\left(I_j\right).$ (18)

定理1. 具有周期边界条件的全离散系统(7)和(8)，其中通量满足(9)~(10)，则：

$\Vert u_h^{n+1}\Vert \le C\Vert u_h^0\Vert$

证明：(Ⅰ) 对于$n\ge 1$ 的情况，令$v_h=u_h^{n+1},v_{h,x}={\left(u_h^{n+1}\right)}_x,v_{h,xx}={\left(u_h^{n+1}\right)}_{xx}$ 代入(8)得

$\begin{array}{l}\left({\delta }_t{u}_h^{n+1},u_h^{n+1}\right)+g\left(u_h^{n+1},u_h^{n+1}\right)+\frac{1}{2}{\displaystyle \sum _{j=1}^N{〚u_h^{n+1}〛}_{j-1/2}^2}+\alpha N\left(u_h^{n+1},u_h^{n+1}\right)\\ +\beta \left({\mathcal{D}}_{\theta }\left(u_h^{n+1}\right),u_h^{n+1}\right)+\gamma A\left(u_h^{n+1},u_h^{n+1}\right)=g{b}_n\left(u_h^0,u_h^{n+1}\right)+g\left[{\displaystyle \sum _{i=1}^n\left(b_{i-1}-b_i\right)}\left(u_h^{n+1-i},u_h^{n+1}\right)\right].\end{array}$ (19)

其中，

${\displaystyle \sum _{j=1}^N\left[{\left({\stackrel{⌢}{u}}_h^{n+1}{\left(u_h^{n+1}\right)}^{-}\right)}_{j+1/2}-{\left({\stackrel{⌢}{u}}_h^{n+1}{\left(u_h^{n+1}\right)}^{+}\right)}_{j-1/2}\right]}-\left(u_h^{n+1},{\left(u_h^{n+1}\right)}_x\right)=\frac{1}{2}{\displaystyle \sum _{j=1}^N〚{u_h^{n+1}〛}_{j-1/2}^2}.$ (20)

设$\mathbb{B}\left(u_h^{n+1},u_h^n\right)\dot{=}{\Vert u_h^{n+1}\Vert }^2+{\Vert 2u_h^{n+1}-u_h^n\Vert }^2$ ，

$\left({\delta }_t{u}_h^{n+1},u_h^{n+1}\right)=\frac{1}{4\tau }\left[\mathbb{B}\left(u_h^{n+1},u_h^n\right)-\mathbb{B}\left(u_h^n,u_h^{n-1}\right)\right]+\frac{1}{4}\tau {\Vert u_h^{n+1}-2u_h^n+u_h^{n-1}\Vert }^2.$

将(14) (16) (18) (20)代入(19)利用Cauchy-Schwarz不等式和Young不等式可得

$\begin{array}{c}\frac{1}{4\tau }\mathbb{B}\left(u_h^{n+1},u_h^n\right)\le \frac{1}{4\tau }\mathbb{B}\left(u_h^n,u_h^{n-1}\right)-g\left[{\Vert u_h^{n+1}\Vert }^2-{\displaystyle \sum _{j=1}^n\left(b_{i-1}-b_i\right)}\left(u_h^{n+1-i},u_h^{n+1}\right)-b_n\left(u_h^0,u_h^{n+1}\right)\right]\\ \le \frac{1}{4\tau }\mathbb{B}\left(u_h^n,u_h^{n-1}\right)-\frac{g}{2}{\displaystyle \sum _{j=0}^n{b}_i}{\Vert u_h^{n+1-i}\Vert }^2+\frac{g}{2}{\displaystyle \sum _{j=0}^{n-1}{b}_i}{\Vert u_h^{n-i}\Vert }^2+\frac{g}{2}{b}_n{\Vert u_h^0\Vert }^2.\end{array}$ (21)

为了简便，引入以下符号

$\Xi \left(u_{h}n+1\right)\dot{=}\mathbb{B}\left(u_h^{n+1},u_h^n\right)+2g\left(\nu \right){\tau }^{\frac{1}{2}}{\displaystyle \sum _{j=0}^n{b}_j}{\Vert u_h^{n+1-j}\Vert }^2.$ (22)

其中$g\left(\nu \right)=\frac{\sqrt{\nu }}{\Gamma \left(3/2\right)}$ ，将(22)代入(21)并且将(21)乘$4\tau$ ，得

$\begin{array}{c}\Xi \left(u_h^{n+1}\right)\le \Xi \left(u_h^n\right)+2g\left(\nu \right){\tau }^{1/2}{b}_n{\Vert u_h^0\Vert }^2\\ \le \Xi \left(u_h^{n-1}\right)+2g\left(\nu \right){\tau }^{1/2}{b}_{n-1}{\Vert u_h^0\Vert }^2+2g\left(\nu \right){\tau }^{1/2}{b}_n{\Vert u_h^0\Vert }^2\\ \le \cdots \le \Xi \left(u_h^1\right)+2g\left(\nu \right){\tau }^{1/2}{\displaystyle \sum _{j=1}^n{b}_j}{\Vert u_h^0\Vert }^2\le \Xi \left(u_h^1\right)+2g\left(\nu \right)T^{1/2}{\Vert u_h^0\Vert }^2.\end{array}$ $$ (23)

下面估计$\Xi \left(u_h^1\right)$ 。

(Ⅱ) 对于$n=0$ 的情形，在(7)中取$v_h=u_h^1$ ，$v_{hx}=u_{hx}^1$ ，$v_{hxx}=u_{hxx}^1$ ，

(24)

故

$\left({\delta }_t{u}_h^1,u_h^1\right)+g\left(u_h^1,u_h^1\right)+\frac{1}{2}{\displaystyle \sum _{j=1}^N{〚u_h^1〛}_{j-\frac{1}{2}}^2}+\alpha \mathcal{N}\left(u_h^1,u_h^1\right)+\beta \mathcal{D}\left(u_h^1,u_h^1\right)+\gamma A\left(u_h^1,u_h^1\right)=g{b}_n\left(u_h^0,u_h^1\right).$

利用Cauchy-Schwarz不等式和Young不等式，得

$\begin{array}{c}\left(\frac{1}{\tau }+g\right){\Vert u_h^1\Vert }^2+\frac{1}{2}{\displaystyle \sum _{j=1}^N{〚u_h^1〛}_{j-\frac{1}{2}}^2}\le \frac{1}{\tau }\Vert u_h^0\Vert \Vert u_h^1\Vert +g{b}_0\Vert u_h^0\Vert \Vert u_h^1\Vert \\ \le \left(\frac{1}{2\tau }+\frac{g{b}_0}{2}\right)\left({\Vert u_h^0\Vert }^2+{\Vert u_h^1\Vert }^2\right).\end{array}$ (25)

用$\tau$ 乘(25)，可得

$\left(\frac{1}{2}+\frac{g\tau }{2}\right){\Vert u_h^1\Vert }^2\le \left(\frac{1}{2}+\frac{g\tau }{2}\right){\Vert u_h^0\Vert }^2.$ (26)

$n=0$ 的情况得证。

(Ⅲ) 注意到$\Xi \left(u_h^1\right)\le C{\Vert u_h^1\Vert }^2\le C{\Vert u_h^0\Vert }^2$ 。结合上面的证明，稳定性得证。

5. 数值实验

例：数值求解以下问题

$\begin{array}{l}\sqrt{\nu }{}_C{D}_{0,t}^{\xi }+u_t+u_x+\beta u_{xxx}+\alpha f{\left(u\right)}_x=\gamma u_{xx},\left(x,t\right)\in \overline{\Omega }\times \left(0,T\right],\\ u_0\left(x\right)=0.32{\mathrm{sech}}^2\left(0.4\left(x-x_0\right)\right),t=0.\end{array}$

边界条件为周期边界条件，其中$\overline{\Omega }=[a,b]$ ，$x_0=\left(b-a\right)/2$ ，$T=10$ 。

在表1中，令$\overline{\Omega }=\left[0,40\right]$ ，$\xi =0.5$ ，时间步长固定为$\tau =0.01$ ，分别取$N=20,40,80$ ，改变算例中空间导数项前的系数计算空间收敛阶。数值结果表明，对于空间收敛阶，扩散项前的系数对其影响较大。无论方程中各项系数如何变化，该方法均具有$O\left(h^k\right)$ 的空间收敛阶。在表2中，令$\overline{\Omega }=\left[0,40\right]$ ，$\xi =0.5$ 空间步长固定为$h=0.04$ ，分别取$M=20,40,80$ ，改变算例中空间导数项前的系数，给出数值算法的误差以及时间收敛阶。表2中的数值结果表明完全离散格式具有$O\left({\tau }^{2-\xi }\right)$ 的时间收敛阶。观察表中数据可以看出时间分数阶导数项前面的系数对时间收敛阶具有较大的影响。表3改变初始条件，令$\overline{\Omega }=\left[0,400\right]$ ，进行数值实验，可以看出初值的变化对时间收敛阶没有显著影响。

Table 1. Error and spatial convergence order of $\overline{\Omega }=\left[0,40\right],T=10,M=1000,k=2$

表1. $\overline{\Omega }=\left[0,40\right],T=10,M=1000,k=2$ 的误差和空间收敛阶

	$N$	$\Vert u-u_h\Vert$	收敛阶	${\Vert u-u_h\Vert }_{\infty }$	收敛阶
$\nu =\alpha =\beta =\gamma =1$	20	3.1451E−03	-	6.8540E−03	-
	40	6.8278E−04	2.2036	1.4804E−03	2.2109
	80	1.6506E−04	2.0484	3.5454E−04	2.0620
$\nu =\alpha =\beta =1,\gamma =0.01$	20	1.7447E−03		6.8540E−03
	40	2.0532E−04	3.0871	1.4804E−03	2.2109
	80	2.5384E−05	3.0158	3.5454E−04	2.0620
$\gamma =\alpha =\beta =1,\nu =10$	20	2.3484E−03	-	5.1577E−03	-
	40	4.4933E−04	2.3858	9.9242E−04	2.3777
	80	1.0372E−04	2.1151	2.2807E−04	2.1215
$\gamma =\alpha =\nu =1,\beta =10$	20	2.2270E−03	-	4.4947E−03	-
	40	3.8236E−04	2.5421	7.7061E−04	2.5441
	80	8.3689E−05	2.1918	1.8146E−04	2.0863
$\gamma =\beta =\nu =1,\alpha =0.1$	20	3.1419E−03	-	6.9893E−03	-
	40	6.7899E−04	2.2102	1.5253E−03	2.1961
	80	1.6389E−04	2.0506	3.6925E−04	2.0464

Table 2. Error and temporal convergence order of $\overline{\Omega }=\left[0,40\right],T=10,M=1000,k=2$

表2. $\overline{\Omega }=\left[0,40\right],T=10,M=1000,k=2$ 的误差和时间收敛阶

	$M$	$\Vert u-u_h\Vert$	收敛阶	${\Vert u-u_h\Vert }_{\infty }$	收敛阶
	20	2.2954E−02	-	2.3214E−05	-
$\nu =\alpha =\beta =\gamma =1$	40	5.6166E−03	2.0310	9.3187E−06	2.0817
	80	1.4384E−03	1.9653	3.6484E−06	1.9888
$\nu =\alpha =\beta =1,\gamma =0.01$	20	1.7447E−03	-	4.3869E−05	-
	40	2.0532E−04	1.4165	1.7213E−05	1.3497
	80	2.5384E−05	1.3894	6.7702E−06	1.3462
$\gamma =\alpha =\beta =1,\nu =10$	20	4.1222E−03	-	2.3705E−04	-
	40	1.8289E−03	1.1724	1.0385E−04	1.1907
	80	7.9512E−04	1.2018	4.4755E−05	1.2144
$\gamma =\beta =\nu =1,\alpha =0.1$	20	3.8475E−03	-	3.8475E−03	-
	40	1.5567E−03	1.3055	1.5567E−03	1.3055
	80	6.1394E−04	1.3423	6.1394E−04	1.3423
$\gamma =\alpha =\nu =1,\beta =10$	20	1.5871E−03	-	1.7547E−06	-
	40	6.2888E−04	1.3355	6.9266E−06	1.3410
	80	2.4643E−04	1.3516	2.7665E−06	1.3241

Table 3. Error and temporal convergence order of $\overline{\Omega }=\left[0,400\right],T=10,M=1000,k=2$

表3. $\overline{\Omega }=\left[0,400\right],T=10,M=1000,k=2$ 的误差和时间收敛阶

	$M$	$\Vert u-u_h\Vert$	收敛阶	${\Vert u-u_h\Vert }_{\infty }$	收敛阶
$\nu =\alpha =\beta =\gamma =1$	20	1.5000E−03	−	2.3214E−04	−
	40	6.0224E−04	1.3166	9.3187E−05	1.3168
	80	2.3653E−04	1.3483	3.6484E−05	1.3529

图1为$T=100,200,300,400$ ，$\overline{\Omega }=\left[0,200\right]$ 时数值解$u_h$ 的图像，其中$\nu =\alpha =\beta =\gamma =1$ 。通过观察图像中不同时间点的波形，可以直观看出波的衰减程度。

Figure 1. Numerical solution of u

图1. u的数值解图像

6. 总结

本文数值求解了非局部粘性水波模型。首先，在时间方向上使用BDF2方法和L1近似公式分别离散分数阶导数和整数阶导数，在空间方向上采用DDG方法，得到全离散格式。其次，给出了数值算法以及具体的矩阵计算方法并进行稳定性分析。最后，通过数值算例验证该算法的有效性。在对时间分数阶导数进行离散时，我们选用L1近似公式，此方法尽管在精度方面相对有限，但其对解的正则性要求较低，也不会因初始条件的变化产生较大波动。在以后的研究中，我们将考虑运用更高精度的方法进行离散。在空间方向上，与有限差分和有限元方法相比，DDG方法可以随着基函数个数的选取达到任意收敛阶，并且处理复杂的边界条件时更具灵活性。

NOTES

^*通讯作者。