分数阶守恒Swift-Hohenberg方程的三种显隐Runge-Kutta方法

Three Explicit and Implicit Runge-Kutta Methods for Fractional Conservation Swift-Hohenberg Equation

摘要: 分数阶守恒Swift-Hohenberg (SH)方程是材料学中模拟凝固微观组织的基本模型。由于分数阶导算子及非局部守恒项的影响，许多求解经典整数阶SH方程行之有效的数值方法在解决此类问题时存在严重困难。本文针对分数阶守恒SH方程，研究其高效数值逼近算法。首先，在时间方向采用显隐Runge-Kutta方法，空间方向采用傅里叶谱方法，构造分数阶守恒SH方程的数值格式；其次，进一步给出所建立格式质量守恒的理论分析；最后，通过数值实验验证了格式的收敛阶和能量递减性，同时对长时间动力行为进行模拟，验证了算法的有效性。

Abstract: Fractional conservation Swift-Hohenberg (SH) equation is a basic model for simulating solidification microstructure in material science. Due to the influence of fractional derivative operators and nonlocal conservation terms, many effective numerical methods for solving classical integer order SH equations have severe difficulties in solving such problems. In this paper, an efficient numerical approximation algorithm for fractional conservative SH equation is studied. Firstly, the explicit and implicit Runge-Kutta method is used in the time direction and Fourier spectrum method is used in the space direction to construct the numerical schemes of fractional order conserved SH equation. Secondly, the theoretical analysis of mass conservation of these schemes is given. Finally, the convergence order and energy decrease of these schemes are verified by numerical experiments, and the validity of these algorithms are verified by long-term dynamic behavior simulation.

1. 引言

相场晶体模型，是一种能够有效模拟凝固微观组织的新方法 [1]，可用于模拟裂纹萌生、相分离、外延生产机制、晶界迁移和位错运动等，因此相场晶体模型成为科研学者进行微观组织研究的重要工具 [2]，而Swift-Hohenberg (SH)方程就是一种重要的相场晶体模型。SH方程是由Swift和Hohenberg [3] 在1976年研究卷波的Rayleigh-Bénard不稳定性的简单模型时提出的，随后有大量学者对SH方程进行了数值求解 [4] [5] [6]。

经典的SH模型是四阶整数阶的偏微分方程，整数阶微分方程很难刻画具有记忆效应的长时间动力行为，无法模拟许多现实生活中的自然现象。而采用分数阶微分方程求解的数学模型能获得更精确、更符合实际的结果，因此国内外众多科研人员将分数阶微分思想应用到SH方程的求解上。例如：Li和Pang [7] 研究了一种迭代方法，将其应用于具有初始值的时间分数阶SH方程，得到了具有数值图形的初值问题的近似解析解；Zahra等 [8] 提出了一种基于指数拟合技术的非线性时间分数阶SH方程求解新方案；Hussam等 [9] 使用Laplace Adomian分解方法给出了分数阶SH方程的近似解析解；Prakasha等 [10] 借助于残差幂级数方法研究了分数阶SH方程的近似解析解；Weng等 [11] 给出了具有非局部非线性的SH方程的一种高阶精确的快速显式算子分裂方案。

上述方法都是在经典的SH方程上研究的，经典的SH方程具有梯度流动结构，有能量递减特性，但是存在质量不守恒的问题。为了改善这个缺陷，2013年，Elsey和Wirth [12] 引入了非局部Lagrange乘子，提出了质量守恒的SH方程，且证明了该乘子不会影响原始的能量定律，因此该方程具有重要的科研价值。在此之后，守恒型SH方程相关数值算法的实现得到广泛的发展与应用。Lee [13] 研究了质量守恒型SH方程的高效解法及其数值实现；Wang和Zhai [14] 基于算子分裂法和谱方法，提出了一种快速高效的SH方程数值算法；Zhang和Yang [15] 通过将IEQ方法与稳定技术相结合，推导出了质量守恒型SH方程。

由于SH方程求解的复杂性与应用的广泛性，国内外已有许多学者对其数值算法进行研究。例如Kenichi等 [16] 利用Painleve分析，Hirota多线性方法和直接分析技巧得到了(1 + 1)维复三次和五次SH方程的一些精确孤子解；王志林和纪峰波 [17] 利用齐次平衡原则和F-展开式方法，得到一维复系数的SH方程的精确孤立波解；林国广等 [18] 在非线性正核 $G\left(r\right)$ 的有界性及光滑性条件下，获得了非局部二维SH方程的惯性流形存在性的结果；兰天柱等 [19] 运用多尺度方法研究了一类推广的SH方程的规范型；Lee基于傅里叶谱方法，提出了该方程的半解析方法 [20] 和具有二次三次非线性项的能量稳定方法 [21]。

基于前人的工作，本文将在时间方向上采用Runge-Kutta方法，在空间方向上采用傅里叶谱方法求解SH方程，以验证分数阶守恒SH方程数值格式能量递减和质量守恒的特性，最后通过数值实验验证格式的高精度和有效性。

2. 数值格式的构造

SH方程可以视为 Lyapunov 能量泛函的 $L^2$ 梯度流，设基本能量泛函为 $E\left(\varphi \right)$，即

$E\left(\varphi \right)={\displaystyle {\int }_{\Omega }\left(F\left(\varphi \right)-{\left|\Delta \varphi \right|}^2+\frac{1}{2}{\left(\Delta \varphi \right)}^2\right)\text{d}x}$ (1)

其中， $F\left(\varphi \right)=\frac{1}{4}{\varphi }^4-\frac{g}{3}{\varphi }^3-\frac{\epsilon }{2}{\varphi }^2$，Ω是 ${ℝ}^d$ 中的一个有界域 $\left(d=1,2,3\right)$，且能量泛函 $E\left(\varphi \right)$ 满足耗散性质。

本文研究具有周期边界条件的分数阶守恒SH方程，如下所示：

$\{\begin{array}{l}\frac{\partial \varphi }{\partial t}=-\left(F^{\prime }\left(\varphi \right)+2{\Delta }^{\frac{\alpha }{2}}\varphi +{\Delta }^{\alpha }\varphi \right)+{\left(F\left(\varphi \right)\right)}^r\beta \left(t\right),\left(x,t\right)\in {\left[0,L\right]}^2\times \left[0,T\right],\\ {\varphi }_0=\varphi \left(x,0\right).\end{array}$ (2)

其中， $\varphi$ 是密度场，参数 $g\ge 0$ 和 $\epsilon \ge 0$ 为具有物理意义的正常数， $\alpha \in \left(1,2\right]$， $\beta \left(t\right)=\frac{{\displaystyle {\int }_{\Omega }{F}^{\prime }\left(\varphi \left(x,t\right)\right)\text{d}x}}{{\displaystyle {\int }_{\Omega }{\left(F\left(\varphi \left(x,t\right)\right)\right)}^r\text{d}x}}$。

SH方程的主要优点是添加了非局部Lagrange乘子，使得方程具有质量守恒性。本文在空间方向上，采用傅里叶谱方法：

令 ${\varphi }_m^n$ 为 $\varphi \left(x_m,t^n\right)$ 的近似值，其中 $x_i=mh,\left(m=0,1,\cdots ,M-1\right)$， $t^n=n\Delta t$，离散傅里叶变换为

${\hat{\varphi }}_k={\displaystyle \underset{m=0}{\overset{M-1}{\sum }}{\varphi }_m{\text{e}}^{-i{x}_m{\xi }_k}}$ (3)

其逆变换为

${\varphi }_m=\frac{1}{M}{\displaystyle \underset{k=0}{\overset{M-1}{\sum }}{\hat{\varphi }}_m{\text{e}}^{-i{x}_m{\xi }_k}}$ (4)

其中， ${\xi }_k=\frac{2\pi k}{L}$。

进一步结合一阶向前差分、二阶和三阶Runge-Kutta方法，我们就能得到方程(2)的如下三种求解格式：

S1: ${\varphi }_m^{n+1}={\varphi }_m^n+\Delta t\left(p\left({\varphi }_m^{n+1}\right)+q\left({\varphi }_m^n\right)\right)$

S2: $\{\begin{array}{l}{\varphi }_m^{\left(1\right)}={\varphi }_m^n+\Delta t\left(\gamma p\left({\varphi }_m^{\left(1\right)}\right)+\gamma q\left({\varphi }_m^n\right)\right),\\ {\varphi }_m^{n+1}={\varphi }_m^n+\Delta t\left(\gamma p\left({\varphi }_m^{n+1}\right)+\left(1-\gamma \right)p\left({\varphi }_m^{\left(1\right)}\right)+\left(1-\delta \right)q\left({\varphi }_m^{\left(1\right)}\right)+\delta q\left({\varphi }_m^n\right)\right)\end{array}$

S3: $\{\begin{array}{l}{\varphi }_m^{\left(1\right)}={\varphi }_m^n+\Delta t\left(\frac{1}{2}p\left({\varphi }_m^{\left(1\right)}\right)+\frac{1}{2}q\left({\varphi }_m^n\right)\right),\\ {\varphi }_m^{\left(2\right)}={\varphi }_m^n+\Delta t\left(\frac{1}{2}p\left({\varphi }_m^{\left(1\right)}\right)+\frac{1}{6}p\left({\varphi }_m^{\left(1\right)}\right)+\frac{1}{18}q\left({\varphi }_m^{\left(1\right)}\right)+\frac{11}{18}q\left({\varphi }_m^n\right)\right),\\ {\varphi }_m^{\left(3\right)}={\varphi }_m^n+\Delta t\left(\frac{1}{2}p\left({\varphi }_m^{\left(3\right)}\right)+\frac{1}{2}p\left({\varphi }_m^{\left(2\right)}\right)-\frac{1}{2}p\left({\varphi }_m^{\left(1\right)}\right)+\frac{1}{2}q\left({\varphi }_m^{\left(2\right)}\right)-\frac{5}{6}q\left({\varphi }_m^{\left(1\right)}\right)+\frac{5}{6}q\left({\varphi }_m^n\right)\right),\\ {\varphi }_m^{n+1}={\varphi }_m^n+\Delta t(\frac{1}{2}p\left({\varphi }_m^{n+1}\right)+\frac{1}{2}p\left({\varphi }_m^{\left(3\right)}\right)-\frac{3}{2}p\left({\varphi }_m^{\left(2\right)}\right)+\frac{3}{2}p\left({\varphi }_m^{\left(1\right)}\right)-\frac{7}{4}q\left({\varphi }_m^{\left(3\right)}\right)+\frac{3}{4}q\left({\varphi }_m^{\left(2\right)}\right)\\ +\frac{7}{4}q\left({\varphi }_m^{\left(1\right)}\right)+\frac{1}{4}q\left({\varphi }_m^n\right)),\end{array}$

其中 $\{\begin{array}{l}p\left({\varphi }_m^n\right)=-\left(s{\varphi }_m^n+2\frac{{\partial }^{\alpha }{\varphi }_m^n}{\partial x^{\alpha }}+\frac{{\partial }^{2\alpha }{\varphi }_m^n}{\partial x^{2\alpha }}\right)\\ q\left({\varphi }_m^n\right)=-\left(F^{\prime }\left({\varphi }_m^n\right)-s{\varphi }_m^n\right)+{\left(F\left({\varphi }_m^n\right)\right)}^r{\beta }^n\end{array}$， ${\beta }^n=\frac{{\displaystyle \underset{m=0}{\overset{M-1}{\sum }}{F}^{\prime }\left({\varphi }_m^n\right)}}{{\displaystyle \underset{m=0}{\overset{M-1}{\sum }}{\left(F\left({\varphi }_m^n\right)\right)}^r}}$， $\gamma =\frac{2-\sqrt{2}}{2}$， $\delta =1-\frac{1}{2\gamma }$。

3. 守恒性

定理1：格式S1满足质量守恒，即 ${\displaystyle \underset{m=0}{\overset{M-1}{\sum }}{\varphi }_m^{n+1}}={\displaystyle \underset{m=0}{\overset{M-1}{\sum }}{\varphi }_m^n}$。

证明：对S1采用谱变换后可得

${\hat{\varphi }}_k^{n+1}=\frac{{\hat{\varphi }}_k^n+\Delta t{\hat{q}}_k^n}{1+\Delta t\left(s-2{\xi }_k^{\alpha }+{\xi }_k^{2\alpha }\right)}$ (5)

结合公式(3)可得

${\displaystyle \underset{m=0}{\overset{M-1}{\sum }}{\varphi }_m^{n+1}}=\frac{{\hat{\varphi }}_0^n+\Delta t{\hat{q}}_0^n}{1+s\Delta t}$ (6)

由于

$\begin{array}{c}{\hat{q}}_0^{\left(n\right)}=-{\displaystyle \underset{m=0}{\overset{M-1}{\sum }}\left(F^{\prime }\left({\varphi }_m^{\left(n\right)}\right)-s{\varphi }_m^{\left(n\right)}\right)}+{\displaystyle \underset{m=0}{\overset{M-1}{\sum }}{\left(F\left({\varphi }_m^{\left(n\right)}\right)\right)}^r}\frac{{\displaystyle \underset{m=0}{\overset{M-1}{\sum }}\left(F^{\prime }\left({\varphi }_m^{\left(n\right)}\right)\right)}}{{\displaystyle \underset{m=0}{\overset{M-1}{\sum }}{\left(F\left({\varphi }_m^{\left(n\right)}\right)\right)}^r}}\\ =s{\displaystyle \underset{m=0}{\overset{M-1}{\sum }}{\varphi }_m^{\left(n\right)}}=s{\hat{\varphi }}^{\left(n\right)}\end{array}$ (7)

因此

${\displaystyle \underset{m=0}{\overset{M-1}{\sum }}{\varphi }_m^{n+1}}={\hat{\varphi }}_0^{n+1}=\frac{{\hat{\varphi }}_0^n+\Delta t{\hat{q}}_0^n}{1+s\Delta t}=\frac{{\hat{\varphi }}_0^n+s\Delta t{\hat{\varphi }}_0^n}{1+s\Delta t}={\hat{\varphi }}_0^n={\displaystyle \underset{m=0}{\overset{M-1}{\sum }}{\varphi }_m^n}$ (8)

得证方法S1质量守恒。

定理2：格式S2满足质量守恒，即 ${\displaystyle \underset{m=0}{\overset{M-1}{\sum }}{\varphi }_m^{n+1}}={\displaystyle \underset{m=0}{\overset{M-1}{\sum }}{\varphi }_m^n}$。

证明：由S1可知

${\hat{\varphi }}_0^{\left(1\right)}=\frac{{\hat{\varphi }}_0^n+rs\Delta t{\hat{\varphi }}_0^n}{1+rs\Delta t}={\hat{\varphi }}_0^n$ (9)

对S2的第二个公式采用谱变换，并整理可得

${\hat{\varphi }}_0^{n+1}=\frac{{\hat{\varphi }}_0^n+\Delta t\left[-\left(1-r\right)s{\hat{\varphi }}_0^{\left(1\right)}+\left(1-\delta \right){\hat{\varphi }}_0^{\left(1\right)}+\delta s{\hat{\varphi }}_0^n\right]}{1+rs\Delta t}=\frac{{\hat{\varphi }}_0^n+rs\Delta t{\hat{\varphi }}_0^n}{1+rs\Delta t}={\hat{\varphi }}_0^n$ (10)

因此

${\displaystyle \underset{m=0}{\overset{M-1}{\sum }}{\varphi }_m^{n+1}=}{\hat{\varphi }}_0^{n+1}={\hat{\varphi }}_0^n={\displaystyle \underset{m=0}{\overset{M-1}{\sum }}{\varphi }_m^n}$ (11)

得证方法S2质量守恒。

定理3：格式S3满足质量守恒，即 ${\displaystyle \underset{m=0}{\overset{M-1}{\sum }}{\varphi }_m^{n+1}}={\displaystyle \underset{m=0}{\overset{M-1}{\sum }}{\varphi }_m^n}$。

证明：由S1与S2推导的结果得

${\hat{\varphi }}_m^{\left(2\right)}={\hat{\varphi }}_m^{\left(1\right)}={\hat{\varphi }}_0^n$ (12)

对S2的第三个公式采用谱变换，并整理可得：

$\begin{array}{c}{\hat{\varphi }}_0^{\left(3\right)}=\frac{{\hat{\varphi }}_0^n+\Delta t\left[-\frac{1}{2}s{\hat{\varphi }}_0^{\left(2\right)}+\frac{1}{2}s{\hat{\varphi }}_0^{\left(1\right)}+\frac{1}{2}s{\hat{\varphi }}_0^{\left(2\right)}-\frac{5}{6}s{\hat{\varphi }}_0^{\left(1\right)}+\frac{5}{6}s{\hat{\varphi }}_0^n\right]}{1+\frac{1}{2}s\Delta t}\\ =\frac{{\hat{\varphi }}_0^n+\frac{1}{2}s\Delta t{\hat{\varphi }}_0^n}{1+\frac{1}{2}s\Delta t}={\hat{\varphi }}_0^n\end{array}$ (13)

因此

$\begin{array}{c}{\displaystyle \underset{m=0}{\overset{M-1}{\sum }}{\varphi }_m^{n+1}}={\hat{\varphi }}_0^{n+1}\\ =\frac{{\hat{\varphi }}_0^n+\Delta t\left(-\frac{s}{2}{\hat{\varphi }}_0^{\left(3\right)}+\frac{3s}{2}{\hat{\varphi }}_0^{\left(2\right)}-\frac{3s}{2}{\hat{\varphi }}_0^{\left(1\right)}-\frac{7s}{4}{\hat{\varphi }}_0^{\left(3\right)}+\frac{3s}{4}{\hat{\varphi }}_0^{\left(2\right)}+\frac{s}{4}{\hat{\varphi }}_0^{\left(3\right)}\right)}{1+\frac{1}{2}s\Delta t}\\ =\frac{{\hat{\varphi }}_0^n+\frac{1}{2}s\Delta t{\hat{\varphi }}_0^n}{1+\frac{1}{2}s\Delta t}\\ ={\displaystyle \underset{m=0}{\overset{M-1}{\sum }}{\varphi }_m^n}\end{array}$ (14)

得证方法S3质量守恒。

4. 数值实验

4.1. 收敛性测试

考虑区域 $\Omega ={\left[0,32\right]}^2$ 上的分数阶守恒SH方程，初值如下：

$\varphi \left(x,0\right)={\left[0.07-0.02\cos \left(\frac{2\pi \left(x-12\right)}{32}\right)+0.02{\cos }^2\left(\frac{\pi \left(x+10\right)}{32}\right)-0.01{\sin }^2\left(\frac{4\pi x}{32}\right)\right]}^2$

参数选取如下：

$T=5,\epsilon =0.1,g=0,r=0,s=2,h=\frac{1}{3}$

为方便分析，定义 $e_i^k={\varphi }_i^k-{\hat{\varphi }}_i^k$，其中， $e_i^k$ 表示时间剖分为N时，对应位置上函数 $\varphi$ 的值； ${\hat{\varphi }}_i^k$ 表示时间剖分为2N时对应位置上函数 $\varphi$ 的值，定义最大误差 $Er{r}_{\infty }$ 和 $L_2$ 误差 $Er{r}_2$ 分别为： $Er{r}_{\infty }=\underset{1\le i\le N}{\max }\left|e_i^k\right|$ 和 $Er{r}_2=h\sqrt{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left(e_i^k\right)}^2}}$。

下面我们计算 $T=5$ 时不同时间步长( $\Delta t$ 分别取 $\frac{1}{200}$， $\frac{1}{400}$， $\frac{1}{800}$)、不同分数阶( $\alpha$ 分别取1.2、1.5、1.8、2.0)时的最大误差 $Er{r}_{\infty }$ 及相应的收敛阶，结果如下：

Table 1. Maximum Error and Convergence order of S1, S2 and S3 at T = 5

表1. 格式S1、S2、S3在T = 5的最大误差及收敛阶

Table 2. L2 error and Convergence order of S1, S2 and S3 at T = 5

表2. 格式S1、S2、S3在T = 5的L₂误差及收敛阶

如表1、表2所示：随着网格剖分精度的增加，S1、S2、S3的 $Er{r}_{\infty }$ 和 $Er{r}_2$ 均减小，且收敛精度逐渐接近理论结果；在达到相同的误差精度时，S3对应的网格剖分精度比S2更低，S2对应的网格剖分精度比S1更低。

4.2. 二维相图仿真

定义能量函数 $E\left(\varphi \right)$ 为

$E\left(\varphi \right)={\displaystyle {\int }_{\Omega }\left(\frac{\varphi }{2}{\left(\Delta +1\right)}^{\alpha }\varphi +F\left(\varphi \right)\right)\text{d}x}$

取定义域 $\Omega ={\left[0,32\right]}^2$ 上的初值为

$\varphi \left(x,0\right)=\bar{\varphi }+\text{rand}\left(x,0\right),$

其中， $\bar{\varphi }=0.02$，rand是网格点处介于−0.1和0.1之间的随机数，其他相应的参数为 $\epsilon =0.1$， $g=0$， $r=0$，

$s=2$， $h=\frac{1}{3}$。

对于三阶显隐向前差分谱格式S3，为了节约计算时间，我们选择不同的时间步长：

$\Delta t=\{\begin{array}{l}\frac{1}{4},T\in \left(1,128\right],\\ 4,T\in \left(128,2176\right].\end{array}$

在上述条件下， $\varphi$ 随时间变化图和能量递减图，如下图：

Figure 1. The phase diagram of the equation $\varphi \left(x,t\right)$ is solved by S3 at $\left(\bar{\varphi },\epsilon \right)=\left(0.02,0.1\right)$ and $\alpha =1.2,1.5,1.8,2.0$

图1. $\left(\bar{\varphi },\epsilon \right)=\left(0.02,0.1\right)$ 且 $\alpha =1.2,1.5,1.8,2.0$，运用S3格式求解出方程 $\varphi \left(x,t\right)$ 的相图

图1给出了密度场 $\varphi$ 在不同 $\alpha$ 下，不同时刻的数值解图像。由图1可知：随着 $\alpha$ 的增大，相图到达稳态的时间变短。图2展示了不同 $\alpha$ 取值下的数值解能量图，进一步验证了算法满足能量递减规律。

Figure 2. Energy consumption at $\alpha =1.2,1.5,1.8,2.0$

图2. $\alpha =1.2,1.5,1.8,2.0$ 的能量耗损

4.3. 二维过冷液体中的晶体生长

在下面的例子中，在定义域 $\left[0,80\right]\times \left[0,80\right]$ 上考虑多晶在2维过冷液体中的生长。我们通过在三个小正方形块中设置三个具有不同取向的完美微晶来表示一种恒定构型，使用以下表达式来定义初始微晶:

${\varphi }_0\left(x_l,y_l\right)=\bar{\varphi }+C_1\left(\cos \left(C_2{x}_l\right)\cos \left(\frac{C_2}{\sqrt{3}}{y}_l\right)-0.5\cos \left(\frac{2C_2}{\sqrt{3}}{y}_l\right)\right),l=1,2,3,$

其中 $\bar{\varphi }=0.285$， $x_l$ 和 $y_l$ 定义了一个局部的笛卡尔坐标系，其常数参数分别为 $C_1=0.446$， $C_2=0.66$。三个长度为4的小正方形块的中心分别位于(32, 22)，(62, 40)，(32, 62)处。为了生成具有不同取向的微晶，

我们采用以下角度变换来产生两个不同角度 $\theta =\frac{\pi }{4},-\frac{\pi }{3}$ 的旋转：

$x_l\left(x,y\right)=x\sin \left(\theta \right)+y\cos \left(\theta \right),y_l\left(x,y\right)=-x\cos \left(\theta \right)+y\sin (\; \theta \; )$

其他相应参数为 $\epsilon =0.1$， $g=0$， $r=0$， $s=2$， $h=\frac{80}{256}$， $\Delta t=\frac{1}{4}$。

在上述条件下， $\varphi$ 随时间变化图和能量递减图，如下图：

Figure 3. Numerical solution images of different moments is solved by S3 at $\alpha =1.5,1.8$

图3. $\alpha =1.5,1.8$ 时，运用S3求解出的不同时刻的数值解图像

Figure 4. Energy consumption at $\alpha =1.5,1.8$

图4. $\alpha =1.5,1.8$ 的能量耗损

图3显示了密度场 $\varphi$ 在 $\alpha =1.5$ 和 $\alpha =1.8$ 时，不同时刻的数值解图像。由图可知，随着时间的推移，微晶相互碰撞并形成晶界，最后充满整个区域并趋于稳态。图4展示了不同 $\alpha$ 取值下的数值解能量图，验证了算法满足能量递减规律。

5. 结论

本文利用显隐Runge-Kutta方法和傅里叶谱方法求解分数阶守恒Swift-Hohenberg方程，建立了求解该方程的三种高效数值格式，并给出了质量守恒的理论分析。最后，通过数值实验验证了所提出格式的收敛阶和能量递减性。

致谢

由衷感谢华侨大学数学科学学院翟术英老师对本论文提供的支持与帮助！

基金项目

华侨大学2021年大学生创新创业训练项目(S202010385038)；华侨大学实验教学与管理改革课题 (SY2021J12)。

NOTES

^*通讯作者。