二维非线性四阶分数阶波动方程的BDF2-WSGI有限元算法

BDF2-WSGI Finite Element Algorithm for a Two-Dimensional Nonlinear Fourth-Order Fractional Wave Equation

摘要: 本文主要研究了二维非线性四阶分数阶波动方程的有效数值算法。通过结合二阶BDF2-WSGI时间离散格式与有限元方法对二维非线性四阶分数阶方程进行求解。首先，引入辅助变量，将分数阶四阶波动问题转化为低阶耦合方程，然后利用Riemann-Liouville分数阶积分对所得方程进行积分，最后使用WSGI逼近公式逼近分数阶积分，形成二阶BDF2有限元格式。本文给出了详细的数值算法，并通过一个二维算例进行了数值试验，验证了算法的有效性和收敛性。

Abstract: This article mainly studies effective numerical algorithms for two-dimensional nonlinear fourth-order fractional wave equations. We combine the second-order BDF2-WSGI time discretization scheme with the finite element method to solve two-dimensional nonlinear fourth-order fractional equations. Firstly, introducing auxiliary variables transforms the fractional fourth-order wave problem into a low-order coupled equation. Then, the Riemann-Liouville fractional integration is used to integrate the resulting equation. Finally, the WSGI approximation formula approximates the fractional integration, forming a second-order BDF2 finite element scheme. This article provides a detailed numerical algorithm and conducts numerical experiments on a two-dimensional example to verify the effectiveness and convergence of the algorithm.

1. 引言

本文主要对以下二维非线性四阶分数阶波动方程进行了数值研究

$\begin{array}{l}{}_0{}^{C}D{}_t^{\gamma }u+{\Delta }^{2}u-\Delta f\left(u\right)=g\left(\mathrm{w},t\right),\left(\mathrm{w},t\right)\in \Omega \times (0,T],\\ u=\Delta u=0,\mathrm{w}\in \partial \Omega ,t\in \left[0,T\right],\\ u\left(\mathrm{w},0\right)={\bar{u}}_0\left(\mathrm{w}\right),u_t\left(\mathrm{w},0\right)={\hat{u}}_0\left(\mathrm{w}\right),\mathrm{w}\in \bar{\Omega }.\end{array}$ (1)

其中 $1\le \gamma \le 2$ ， $0\le T\le \infty$ ，二维空间 $\Omega \subset R^2$ 。 $f\left(u\right)$ 是非线性项， $g\left(\mathrm{w},t\right)$ 是源项， ${}_0{}^{C}D{}_t^{\gamma }u$ 是Caputo分数阶导数

${}_0{}^{C}D{}_t^{\gamma }u=\frac{1}{\Gamma \left(2-\gamma \right)}{\displaystyle {\int }_0^t\frac{{\partial }^{2}u\left(\mathrm{w},t\right)}{\partial s^2}}{\left(t-s\right)}^{1-\gamma }\text{d}s.$

分数阶波动模型 [1] - [7] 作为分数阶偏微分方程的一个重要模型，在化学、生物和物理领域扮演着重要的角色，受到了专家学者们的广泛关注，越来越多的人开始研究分数阶波动模型的解析解和数值解。在 [4] 中，Dahaghin和Hassani采用一种基于广义多项式的优化方法来求解分数阶波动模型。在 [7] 中，Ye等人研究了分数阶波动方程的紧致差分格式，并讨论了差分解的唯一可解性。近年来，具有空间四阶导数的分数阶波动方程 [8] [9] [10] [11] [12] 由于其复杂的结构，精确解难以求解，很多学者使用不同的数值方法进行有关这类方程的数值研究，主要方法有有限元法、谱方法、有限体积法、有限差分法等。在 [11] 中，Huang等人研究了具有空间四阶导数的时间分数阶非线性波动方程的两种线性化格式，并严格证明了这两种格式的收敛性和无条件稳定性。在 [12] 中，Li等人在空间上采用参数五次样条的方法来求解四阶导数的分数阶波动方程，并验证了此方法的可解性，稳定性和收敛性。在 [13] 中，Liu等人提出一种基于混合有限元方法的双网格算法，来求解具有Caputo型时间分数阶导数的非线性四阶反应扩散方程。在 [14] 中，Wang等人针对具有分数阶导数的四阶扩散波模型，提出了一种改进L1格式的混合单元算法，并分析了其误差估计。在 [15] 中，Wang等人在空间方向上使用Galerkin混合有限元方法，在时间方向上使用具有WSGD公式的二阶θ格式研究了二维非线性四阶时间分数阶波动模型。还有诸多研究，我们不能一一列举，但我们可以看到具有空间四阶导数的分数阶波动方程具有很好的研究空间。

求解分数阶模型的重要步骤之一是对分数阶导数项进行处理，许多学者使用WSGD [16] [17] 逼近公式对分数阶导数项进行处理，在 [18] 中Feng研究了一类二维非线性时间分数阶耦合亚扩散方程组模型，在时间上采用Crank-Nicolson格式和WSGD逼近，在空间上采用有限元方法。近年来，WSGI逼近公式作为一个求解分数阶偏微分方程的重要方法也受到了学者的广泛研究。在 [19] 中，Wang等人在时间上采用WSGI近似的Crank-Nicolson格式，在空间上采用有限元方法研究了二维非线性四阶时间分数阶波动模型的算法。在 [20] 中，Cao等人结合二阶Crank-Nicolson格式和WSGI逼近与有限元方法对多维时间分数阶波动方程进行求解。

在本文中，我们将结合BDF2-WSGI格式研究二维非线性四阶分数阶波动方程。为此，我们引入了一个非线性辅助变量，将分数阶四阶波动模型转换成一个低阶耦合系统，对其中一个方程进行积分，然后将时间上具有WSGI近似的BDF2格式与空间上的有限元方法相结合，得到完全离散格式。我们详细展示了数值算法以及如何执行数值计算。最后，通过数值算例验证了方法的可行性，并给出了收敛结果。

本文的其余部分安排如下：在第2节中，我们构造了二维非线性四阶分数阶波动方程完全离散格式。在第3节中，我们提供了一个详细的数值算法。在第4节中，给出了一个数值例子来验证我们算法的有效性。最后，在第5节中，我们将对本文做一个简要的总结。

2. 混合元全离散格式

在本节，我们构造式子(1)的全离散混合元格式。

通过引入辅助变量 $v=\Delta u-f\left(u\right)$ ，方程(1)可以改写为如下低阶耦合系统

$\begin{array}{l}{}_0{}^{C}D{}_t^{\gamma }u+\Delta v=g\left(\mathrm{w},t\right),\left(\mathrm{w},t\right)\in \Omega \times (0,T],\\ v-\Delta u+f\left(u\right)=0,\left(\mathrm{w},t\right)\in \Omega \times (0,T],\\ u=v=0,\mathrm{w}\in \partial \Omega ,t\in \left[0,T\right].\end{array}$ (2)

为了构造全离散格式，对时间区间 $\left[0,T\right]$ 进行均匀剖分使得 $0=t_0<t_1<\cdots <t_n=T$ ，时间节点 $t_n=n\tau$ ，这里 $\tau =\frac{T}{N}$ 为时间不长，N是一个正整数。为简便，我们记 $D_t{u}^1=\frac{u^1-u^0}{\tau }$ ， $u^n=u\left(t_n\right)$ 。

为了分析，我们需要引入Riemann-Liouville分数阶积分，定义如下

${}_0\mathcal{I}{}_t^{\alpha }f=\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle {\int }_0^t\frac{f\left(s\right)}{{\left(t-s\right)}^{1-\alpha }}\text{d}s},$ (3)

其中 $0<\alpha =\gamma -1<1$ ，然后通过此积分来积分(2)的第一个式子，以得到一下耦合系统

$\begin{array}{l}{u}_t\left(\mathrm{w},t\right)+{}_0\mathcal{I}{}_t^{\alpha }\Delta v={}_0\mathcal{I}{}_t^{\alpha }g\left(\mathrm{w},t\right)+{\hat{u}}_0\left(\mathrm{w}\right),\left(\mathrm{w},t\right)\in \Omega \times (0,T],\\ v-\Delta u+f\left(u\right)=0,\left(\mathrm{w},t\right)\in \Omega \times (0,T].\end{array}$ (4)

为了求解(4)，我们需要引入WSGI近似来近似算子 ${}_0\mathcal{I}{}_t^{\alpha }$ 。

${}_0\mathcal{I}{}_t^{\alpha }{u}^{n+1}={\tau }^{\alpha }{\displaystyle \underset{k=0}{\overset{n+1}{\sum }}{\lambda }_k^{\left(\alpha \right)}}{u}^{n+1-k}+\tilde{E}\triangleq {}_0\mathcal{I}{}_t^{\alpha }{u}^{n+1}+\tilde{E},$ (5)

其中 $\tilde{E}=O\left({\tau }^2\right)$ ，并且

${\lambda }_0^{\left(\alpha \right)}=\left(1-\frac{\alpha }{2}\right){\omega }_0^{\left(\alpha \right)},{\lambda }_k^{\left(\alpha \right)}=\left(1-\frac{\alpha }{2}\right){\omega }_k^{\left(\alpha \right)}+\frac{\alpha }{2}{\omega }_{k-1}^{\left(\alpha \right)},k\ge 1,$ (6)

这里的系数 ${\omega }_k^{\left(\alpha \right)}$ 定义如下

${\omega }_k^{\left(\alpha \right)}={\left(-1\right)}^k\left(\begin{array}{l}-\alpha \\ k\end{array}\right),{\omega }_0^{\left(\alpha \right)}=1,{\omega }_k^{\left(\alpha \right)}=\left(1+\frac{\alpha -1}{k}\right){\omega }_{k-1}^{\left(\alpha \right)},k\ge 1.$ (7)

我们使用WSGI近似和BDF2格式得到以下时间方向的半离散格式：

情况1： $n\ge 1$

$\{\begin{array}{l}\frac{3u^{n+1}-4u^n+u^{n+1}}{2\tau }+{}_0\mathcal{I}{}_t^{\alpha }\Delta v^{n+1}={}_0\mathcal{I}{}_t^{\alpha }{g}^{n+1}+{\hat{u}}_0+E^{n+1},\\ v^{n+1}-\Delta u^{n+1}+f\left(u^{n+1}\right)=E_1^{n+1},\end{array}$ (8)

情况2： $n=0$

$\{\begin{array}{l}{D}_t{u}^1+{}_0\mathcal{I}{}_t^{\alpha }\Delta v^1={}_0\mathcal{I}{}_t^{\alpha }{g}^1+{\hat{u}}_0+E^1,\\ v^1-\Delta u^1+f\left(u^1\right)=E_1^1,\end{array}$ (9)

其中 $f\left(u^{n+1}\right)=2f\left(u^n\right)-f\left(u^{n-1}\right)$ ， $E^{n+1}=\tilde{E}+E_c^{n+1}=O\left({\tau }^2\right)$ ， $E_c^{n+1}=D_t{u}^{n+1}-u_t\left(t_{n+1}\right)=O\left({\tau }^2\right)$ ， $E_1^{n+1}=f\left(u^{n+1}\right)-f\left(u\left(t_{n+1}\right)\right)$ 。

接下来我们给出(8)和(9)的弱格式，(8)和(9)的两个式子分别和 $\chi ,\psi \in H_0^1$ 做内积，然后对空间区域 $\Omega$ 积分, 我们得到：

情况1： $n\ge 1$

$\{\begin{array}{l}\left(\frac{3u^{n+1}-4u^n+u^{n+1}}{2\tau },\chi \right)-\left({}_0\mathcal{I}{}_t^{\alpha }\nabla v^{n+1},\nabla \chi \right)=\left({}_0\mathcal{I}{}_t^{\alpha }{g}^{n+1},\chi \right)+\left({\hat{u}}_0,\chi \right)+\left(E^{n+1},\chi \right),\\ \left(v^{n+1},\psi \right)+\left(\nabla u^{n+1},\nabla \psi \right)+\left(f\left(u^{n+1}\right),\psi \right)=\left(E_1^{n+1},\psi \right).\end{array}$ (10)

情况2： $n=0$

$\{\begin{array}{l}\left(D_t{u}^1,\chi \right)-\left({}_0\mathcal{I}{}_t^{\alpha }\nabla v^1,\nabla \chi \right)=\left({}_0\mathcal{I}{}_t^{\alpha }{g}^1,\chi \right)+\left({\hat{u}}_0,\chi \right)+\left(E^1,\chi \right),\\ \left(v^1,\psi \right)+\left(\nabla u^1,\nabla \psi \right)+\left(f\left(u^1\right),\psi \right)=\left(E_1^1,\psi \right).\end{array}$ (11)

取有限元空间 $V_h\subset H_0^1$ 并得到如下全离散格式：

情况1： $n\ge 1$

$\{\begin{array}{l}\left(\frac{3u_h^{n+1}-4u_h^n+u_h^{n-1}}{2\tau },{\chi }_h\right)-{\tau }^{\alpha }{\displaystyle \underset{k=0}{\overset{n+1}{\sum }}{\lambda }_k^{\left(\alpha \right)}}\left(\nabla v_h^{n+1},\nabla {\chi }_h\right)={\tau }^{\alpha }{\displaystyle \underset{k=0}{\overset{n+1}{\sum }}{\lambda }_k^{\left(\alpha \right)}}\left(g_h^{n+1},{\chi }_h\right)+\left({\hat{u}}_0,{\chi }_h\right),\forall {\chi }_h\in V_h,\\ \left(v_h^{n+1},{\psi }_h\right)+\left(\nabla u_h^{n+1},\nabla {\psi }_h\right)+\left(f\left(u_h^{n+1}\right),{\psi }_h\right)=0,\forall {\psi }_h\in V_h.\end{array}$ (12)

情况2： $n=0$

$\{\begin{array}{l}\left(D_t{u}_h^1,{\chi }_h\right)-{\tau }^{\alpha }{\displaystyle \underset{k=0}{\overset{1}{\sum }}{\lambda }_k^{\left(\alpha \right)}}\left(\nabla v_h^1,\nabla {\chi }_h\right)={\tau }^{\alpha }{\displaystyle \underset{k=0}{\overset{1}{\sum }}{\lambda }_k^{\left(\alpha \right)}}\left(g_h^1,{\chi }_h\right)+\left({\hat{u}}_0,{\chi }_h\right)+\left(E_h^1,{\chi }_h\right),\forall {\chi }_h\in V_h,\\ \left(v_h^1,{\psi }_h\right)+\left(\nabla u_h^1,\nabla {\psi }_h\right)+\left(f\left(u_h^1\right),{\psi }_h\right)=0,\forall {\psi }_h\in V_h.\end{array}$ (13)

3. 算法实现

通过取有限元空间 $V_h=\left\{v_h|v_h={\displaystyle \underset{i=0}{\overset{M}{\sum }}{v}_i{\varphi }_i}\right\}$ ，其中M是空间单元的数量，我们有 $u_h\left(t\right)={\displaystyle \underset{i=0}{\overset{M}{\sum }}{u}_i\left(t\right){\varphi }_i}$ ， $v_h\left(t\right)={\displaystyle \underset{j=0}{\overset{M}{\sum }}{v}_j\left(t\right){\varphi }_j}$ 。将这些代入到(12)中并且令 ${\chi }_h={\varphi }_m$ ， ${\psi }_h={\varphi }_l$ ，我们可以得到

$\{\begin{array}{l}{\displaystyle \underset{i=0}{\overset{M}{\sum }}\left(u_i{\left(t\right)}_t\left({\varphi }_i,{\varphi }_m\right)\right)-{\displaystyle \underset{j=0}{\overset{M}{\sum }}{}_0\mathcal{I}{}_t^{\alpha }{v}_j\left(t\right)\left(\nabla {\varphi }_j,\nabla {\varphi }_m\right)}}=\left({}_0\mathcal{I}{}_t^{\alpha }g\left(t\right),{\varphi }_m\right)+\left({\hat{u}}_0,{\varphi }_m\right),m=0,1,2,\cdots ,M,\\ {\displaystyle \underset{j=0}{\overset{M}{\sum }}{v}_j\left(t\right)\left({\varphi }_j,{\varphi }_l\right)-{\displaystyle \underset{i=0}{\overset{M}{\sum }}{u}_j\left(t\right)\left(\nabla {\varphi }_i,\nabla {\varphi }_l\right)}}+\left(f\left({\displaystyle \underset{i=0}{\overset{M}{\sum }}{u}_i\left(t\right){\varphi }_i}\right),{\varphi }_l\right)=0,l=0,1,2,\cdots ,M.\end{array}$ (14)

现在，考虑边界条件，我们得到(14)的半离散有限元格式的矩阵形式

$\{\begin{array}{l}A{u}_{ht}\left(t\right)-B{}_0\mathcal{I}{}_t^{\alpha }{v}_h\left(t\right)={}_0\mathcal{I}{}_t^{\alpha }G\left(t\right)+{\hat{u}}_0,\\ B{u}_h\left(t\right)+A{v}_h\left(t\right)+F\left(u_h\left(t\right)\right)=0,\end{array}$ (15)

其中

$A={\left(\left({\varphi }_i,{\varphi }_m\right)\right)}_{\left(M-1\right)\times \left(M-1\right)}$ ， $B={\left(\left(\nabla {\varphi }_i,\nabla {\varphi }_m\right)\right)}_{\left(M-1\right)\times \left(M-1\right)}$ ，

$F\left(u_h\left(t\right)\right)={\left(\left(f\left({\displaystyle \underset{i=0}{\overset{M}{\sum }}{u}_i\left(t\right){\varphi }_i}\right),{\varphi }_1\right),\cdots ,\left(f\left({\displaystyle \underset{i=0}{\overset{M}{\sum }}{u}_i\left(t\right){\varphi }_i}\right),{\varphi }_{M-1}\right)\right)}^{\text{T}}$ ，

$v_h\left(t\right)={\left(v_1,\cdots ,v_{M-1}\right)}^{\text{T}}$ ， $u_h\left(t\right)={\left(u_1,\cdots ,u_{M-1}\right)}^{\text{T}}$ ，

$G\left(t\right)={\left(\left(g\left(t\right),{\varphi }_1\right),\cdots ,\left(g\left(t\right),{\varphi }_{M-1}\right)\right)}^{\text{T}}$ ，

${\hat{u}}_0={\left(\left({\hat{u}}_0,{\varphi }_1\right),\cdots ,\left({\hat{u}}_0,{\varphi }_{M-1}\right)\right)}^{\text{T}}$ 。

接下来我们给出(15)的全离散格式

$\{\begin{array}{l}A\frac{3u_h^{n+1}-4u_h^n+u_h^{n-1}}{2\tau }-B{}_0\mathcal{I}{}_t^{\alpha }{v}_h^{n+1}={}_0\mathcal{I}{}_t^{\alpha }{G}^{n+1}+{\hat{u}}_0,\\ B{u}_h^{n+1}+A{v}_h^{n+1}+F^{n+1}=0,\end{array}$ (16)

其中

${}_0\mathcal{I}{}_t^{\alpha }{v}_h^{n+1}={\tau }^{\alpha }{\lambda }_0^{\left(\alpha \right)}{v}_h^{n+1}+{\tau }^{\alpha }{\displaystyle \underset{k=1}{\overset{n+1}{\sum }}{\lambda }_k^{\left(\alpha \right)}{v}_h^{n+1-k}}$ ，

$u_h^k={\left(u_1^k,\cdots ,u_{M-1}^k\right)}^{\text{T}}$ ， $v_h^k={\left(v_1^k,\cdots ,v_{M-1}^k\right)}^{\text{T}}$ ，

$F^k={\left(\left(f\left({\displaystyle \underset{i=0}{\overset{M}{\sum }}{u}_i^k{\varphi }_i}\right),{\varphi }_1\right),\cdots ,\left(f\left({\displaystyle \underset{i=0}{\overset{M}{\sum }}{u}_i^k{\varphi }_i}\right),{\varphi }_{M-1}\right)\right)}^{\text{T}}$ ， (17)

${\hat{u}}_0={\left(\left({\hat{u}}_0,{\varphi }_1\right),\cdots ,\left({\hat{u}}_0,{\varphi }_{M-1}\right)\right)}^{\text{T}}$ ，

$g^k={\left(\left(g^k,{\varphi }_1\right),\cdots ,\left(g^k,{\varphi }_{M-1}\right)\right)}^{\text{T}}$ ，

$F^{n+1}=2F^n-F^{n-1}$ ，

其中 $k=1,2,\cdots ,n$ 。

我们可以将(16)简化为

$\{\begin{array}{l}\frac{3A}{2\tau }{u}_h^{n+1}-{\tau }^{\alpha }{\lambda }_0^{\left(\alpha \right)}B{v}_h^{n+1}\\ =\frac{2A}{\tau }{u}_h^n-\frac{A}{2\tau }{u}_h^{n-1}+{\tau }^{\alpha }{\lambda }_0^{\left(\alpha \right)}{g}^{n+1}+{\tau }^{\alpha }{\displaystyle \underset{k=1}{\overset{n+1}{\sum }}{\lambda }_0^{\left(\alpha \right)}{v}^{n+1-k}}+{\tau }^{\alpha }{\displaystyle \underset{k=1}{\overset{n+1}{\sum }}{\lambda }_0^{\left(\alpha \right)}{g}^{n+1-k}+u_0}\\ =q_h^n,\\ B{u}_h^{n+1}+A{v}_h^{n+1}=-2F^n+F^{n-1}=w_h^n,\end{array}$ (18)

所以当 $n\ge 1$ 时，矩阵形式为

$\left(\begin{array}{cc}\frac{3A}{2\tau }& -{\tau }^{\alpha }{\lambda }_0^{\left(\alpha \right)}B\\ B& A\end{array}\right)\left(\begin{array}{c}{u}_h^{n+1}\\ v_h^{n+1}\end{array}\right)=\left(\begin{array}{c}{q}_h^n\\ w_h^n\end{array}\right)$ (19)

当 $n=0$ 时，我们可以通过求解下列线性化的方程组来求解 $u_h^1$ 和 $v_h^1$

$\{\begin{array}{l}A\frac{u_h^1-u_h^0}{\tau }-B{}_0\mathcal{I}{}_t^{\alpha }{v}_h^0={}_0\mathcal{I}{}_t^{\alpha }{G}^0+{\hat{u}}_0={\tau }^{\alpha }{\lambda }_0^{\left(\alpha \right)}{G}^0+{\hat{u}}_0,\\ B{u}_h^1+A{v}_h^1+F^0=0,\end{array}$ (20)

然后我们将(20)写作

$\{\begin{array}{l}\frac{A}{\tau }{u}_h^1={\tau }^{\alpha }{\lambda }_0^{\alpha }\left(B{v}^0+G^0\right)+\frac{A}{\tau }{u}_h^0+{\hat{u}}_0\triangleq q_h^0,\\ B{u}_h^1+A{v}_h^1=-F^0\triangleq w_h^0,\end{array}$ (21)

所以当 $n=0$ 时，我们可以得到

$\left(\begin{array}{cc}\frac{A}{\tau }& O\\ B& A\end{array}\right)\left(\begin{array}{c}{u}_h^1\\ v_h^1\end{array}\right)=\left(\begin{array}{c}{q}_h^0\\ w_h^0\end{array}\right)$ (22)

4. 数值实验

在本节中，我们将提供一个二维例子来验证算法的有效性和收敛阶。

$\begin{array}{l}{}_0{}^{C}D{}_t^{\gamma }u+{\Delta }^{2}u-\Delta f\left(u\right)=g\left(\mathrm{w},t\right),\left(\mathrm{w},t\right)\in \Omega \times (0,T],\\ u=\Delta u=0,\mathrm{w}\in \partial \Omega ,t\in \left[0,T\right],\\ u\left(\mathrm{w},0\right)={\bar{u}}_0\left(\mathrm{w}\right),u_t\left(\mathrm{w},0\right)={\hat{u}}_0\left(\mathrm{w}\right),\mathrm{w}\in \bar{\Omega }.\end{array}$

在本算例中，我们取 $\Omega =\left(0,1\right)\times \left(0,1\right)$ ， $T=\left(0,1\right)$ ，初值 ${\bar{u}}_0={\hat{u}}_0=0$ ，精确解 $u\left(x,y,t\right)=t^3\sin \left(2\pi x\right)\sin \left(2\pi y\right)$ ，非线性项 $f\left(u\right)=u^3-u$ 和源项

$\begin{array}{l}g\left(x,y,t\right)=6t^{3-\gamma }\sin \left(2\pi x\right)\sin \left(2\pi y\right)\Gamma \left(4-\gamma \right)+64{\pi }^4{t}^3\sin \left(2\pi x\right)\sin \left(2\pi y\right)\\ -24{\pi }^2{t}^9\sin \left(2\pi x\right)\sin \left(2\pi y\right)\left({\sin }^2\left(2\pi y\right){\cos }^2\left(2\pi x\right)+{\sin }^2\left(2\pi x\right)\cos \left(2\pi y\right)\right)\\ -8{\pi }^2{t}^3\sin \left(2\pi x\right)\sin \left(2\pi y\right).\end{array}$

表1和表2显示了算例的计算结果，在表1中，对每一个分数阶参数 $\gamma =1.1,1.3,1.5,1.7,1.9$ ，固定时间步长为 $\tau =1/200$ ，空间网格步长分别为 $h=\frac{1}{10},\frac{1}{20},\frac{1}{40}$ ，给出数值算法的误差和空间收敛阶。表一中的数值结果表明完全离散格式具有二阶的空间收敛阶。在表2中，我们还给出了参数变化时的时空收敛阶，其中 $\left(\tau ,h\right)=\left(\frac{1}{10},\frac{1}{10}\right),\left(\frac{1}{20},\frac{1}{20}\right),\left(\frac{1}{40},\frac{1}{40}\right)$ ，结果表明WSGI公式具有几乎不受分数阶参数 $\gamma$ 变化影响的二阶时空收敛阶。

Table 1. Fixed time step τ = 1 200 , spatial convergence result

表1. 固定时间步长 $\tau =\frac{1}{200}$ ，空间收敛结果

Table 2. When h = τ , the convergence result of spatiotemporal error

表2. 当 $h=\tau$ 时, 时空误差收敛结果

续表

5. 总结

本文中，我们将时间上具有WSGI近似的BDF2格式与空间上的有限元方法相结合研究二维非线性四阶分数阶波动方程。我们引入了一个非线性辅助变量，将分数阶四阶波动模型转换成一个低阶耦合系统，对含有分数阶的方程进行积分，使用二阶BDF2-WSGI时间离散格式与有限元方法得到完全离散格式。我们给出了数值算法以及如何执行数值计算，最后通过数值算例验证算法的有效性。

基金项目

本文受到国家自然科学基金项目(12061053, 12161063)、内蒙古自然科学基金项目(2021MS01018)、内蒙古自治区高校创新科研团队计划项目(NMGIRT2413, NMGIRT2207)资助。