四阶非线性时间分数阶方程的混合虚拟元方法

The Mixed Virtual Element Method for Nonlinear Fourth-Order Time-Fractional Equations

摘要: 本文主要研究二维非线性四阶时间分数阶方程的混合虚拟元解法。通过引入辅助变量，将四阶时间分数阶方程转化为低阶的耦合方程。在时间方向上采用广义BDF2-θ得到半离散格式，进而在空间方向上应用混合虚拟元方法得到该方程的全离散格式。本文将给出具体的数值求解过程，并通过一个二维数值算例验证该算法的有效性和收敛性。

Abstract: This article primarily investigates the mixed virtual element method of two-dimensional nonlinear fourth-order time-fractional equations. By introducing an auxiliary variable, the fourth-order time-fractional equation is transformed into lower-order coupled equations. The semi-discrete scheme is obtained by employing the generalized BDF2-θ in the temporal direction. Subsequently, a fully discrete scheme for the equation is achieved by applying the mixed virtual element method in the spatial direction. This article will present the specific numerical process and validate the effectiveness and convergence of algorithm through a two-dimensional numerical example.

1. 引言

本文研究以下二维非线性四阶时间分数阶方程

$u_t-{}_0{}^{C}D{}_t^{\alpha }\text{Δ}u+{}_0{}^{C}D{}_t^{\beta }{\text{Δ}}^{2}u-\text{Δ}u+\epsilon \cdot \nabla f\left(u\right)=g\left(x,y,t\right),\left(x,y,t\right)\in \overline{\Omega }\times I,$ (1)

具有以下边界条件

$u\left(x,y,t\right)=\text{Δ}u\left(x,y,t\right)=0,\left(x,y,t\right)\in \partial \text{Ω}\times \overline{I},$ (2)

初始条件为

$u\left(x,y,0\right)=u_0\left(x,y\right),\left(x,y\right)\in \overline{\Omega }$(3)

其中$\overline{\text{Ω}}\subset {ℛ}^2$是空间区域，$I$ 是时间区域，$I=\left(0,T\right]$ ，$T\in \left(0,\infty \right)$ ，$\epsilon ={\left({\epsilon }_1,{\epsilon }_2\right)}^{\text{T}}$ ，${\epsilon }_1,{\epsilon }_2$ 均为正常数，$f\left(u\right)$ 是非线性项，$g\left(x,y,t\right)$ 是已知的源项。

与整数阶导数相比，分数阶导数在诸多科学领域展现出独特优势，能够更为精准地刻画系统的动态特性。这一特性使得分数阶模型在物理、生物、化学等诸多学科中发挥着关键作用，吸引了众多学者的目光[1]-[8]。以粘弹性材料为例，时间分数阶模型可有效描述其应力松弛与蠕变行为；在流体力学中，传统整数阶模型难以全面捕捉湍流的动态变化及多尺度流动现象，而分数阶模型则能很好地弥补这一不足，类似的应用实例不胜枚举。随着科学的不断进步，大量的数值计算方法被应用于求解四阶分数阶问题，例如有限差分方法[9] [10]，谱方法[11]，有限元方法[12] [13]，有限体积方法[14]，LDG方法[15]-[17]等。其中有限元方法是最为普遍的，但相较于传统有限元的三角形或四边形剖分以及单项式或样条函数作为基函数来说，虚拟元方法则是传统有限元方法的进一步推广。2013年L. beirao da Veiga等人首次提出虚拟元方法[18]，该方法弥补了有限元在剖分和基函数选取上的局限性，虚拟元方法可以采用凸或者非凸多边形作为单元剖分，而虚拟元空间的基函数也没有具体的形式，它的计算是利用自由度。这就使得虚拟元方法有很好的网格适应性和数值稳定性，降低了复杂几何区域上网格生成的难度[19]。迄今为止，虚拟元方法广泛应用于各种模型中，例如Navier-Stokes方程[20]，Sobolev方程[21]，Cahn-Hilliard方程[22]等。近年来虚拟元也被广泛应用于求解分数阶问题，Zhang和Feng在文献[23]中研究了具有非光滑数据的时间分数阶对流扩散反应方程，利用指数变换消去对流项，通过能量投影算子推导出最优误差估计。Li等人[24]利用非协调虚拟元方法求解了时间分数阶反应–亚扩散方程，详细证明了空间半离散和时空全离散下的最优误差估计。针对四阶问题，Zhang和Zhao在文献[25]中利用虚拟元研究了具有内罚的四阶奇异摄动问题，Zhang和Feng在[26]中利用混合虚拟元方法求解四阶线性时间分数阶亚扩散方程，基于L1格式证明变量在全离散情况下的无条件稳定。

本文在空间方向上基于混合虚拟元方法，不同于Zhang和Feng在[26]中的时间离散方式，本文在时间方向上采用广义BDF2-$\theta$ 方法[27]。由于时间分数阶导数中存在积分形式，直接进行数值计算较为困难，需要通过离散化将其转化为可数值求解的形式，广义BDF2-$\theta$ 方法可在其中发挥重要作用，形成多种计算格式。

文章内容具体安排如下：在第2节，利用Caputo与Riemann-Liouville导数之间的关系，即 ${}_0{}^{c}D{}_t^{\alpha }u={}_0{}^{R}D{}_t^{\alpha }\left(u-u_0\right)$ ，可将方程中的Caputo型导数转化为Riemann-Liouville导数进行计算。运用广义BDF2-$\theta$ 方法与混合虚拟元方法，构建模型的全离散格式。第3节，详细阐述具体的数值计算流程。第4节通过一个实际的数值算例，对所提算法的有效性予以验证。最后，在第5节对全文进行简要总结。

2. 数值格式

在本节中，我们将利用广义BDF2-$\theta$ 结合混合虚拟元方法对方程(1)进行时间和空间的离散从而得到全离散的格式。首先，用广义BDF2-$\theta$ 得到时间半离散格式，为此将时间区间$\left[0,T\right]$ 均匀剖分， $0=t_0<t_1<\cdots <t_N=T$ ，$N$ 为正整数，时间步长$\tau =T/N$ ，时间节点$t_n=n\tau$ ，$\left(n=1,2,\cdots ,N\right)$ 。

为了便于计算，令$\overline{u}:=u-u_0$ ，将其代入方程(1)~(3)可得

${\overline{u}}_t-{}_0^c{D}_t^{\alpha }\text{Δ}\overline{u}+{}_0^c{D}_t^{\beta }{\text{Δ}}^2\overline{u}-\text{Δ}\overline{u}+\epsilon \cdot \nabla f\left(\overline{u}\right)=\overline{g}\left(x,y,t\right),\left(x,y,t\right)\in \overline{\text{Ω}}\times \overline{I},$ (4)

边界条件为

$\overline{u}\left(x,y,t\right)=\text{Δ}\overline{u}\left(x,y,t\right)=0,\left(x,y,t\right)\in \partial \text{Ω}\times \overline{I},$(5)

初始条件为

$\overline{u}\left(x,y,0\right)=0,\left(x,y\right)\in \overline{\text{Ω}},$(6)

其中，$\overline{g}=g\left(x,y,t\right)+\text{Δ}{u}_0\left(x,y\right)$。Riemann-Liouville分数阶导数在$t_{n-\theta }$ 时刻的广义BDF2-$\theta$ 的逼近公式为[27]：

$\begin{array}{c}{}_0^R{D}_t^{\gamma }{v}^{n-\theta }={\tau }^{-\gamma }{\displaystyle \sum _{j=1}^n{\omega }_j^{\left(\gamma \right)}{v}^{n-j}}+{\tau }^{\gamma }{\displaystyle \sum _{j=1}^{\kappa }{\omega }_{n,j}^{\left(\gamma \right)}{v}^j}+O\left({\tau }^2\right)\\ \triangleq {\text{Ψ}}_{\tau }^{\gamma ,n}v+S_{\tau ,\kappa }^{\gamma ,n}v+O\left({\tau }^2\right),\end{array}$

其中${\text{Ψ}}_{\tau }^{\gamma ,n}v={\tau }^{-\gamma }{\displaystyle \sum _{j=1}^n{\omega }_j^{\left(\gamma \right)}{v}^{n-j}}$ 为离散卷积部分，它的卷积系数${\left\{{\omega }_j^{\left(\gamma \right)}\right\}}_{j=0}^{\infty }$ 由以下函数生成

${\omega }^{\left(\gamma \right)}\left(\xi \right)={\left(\frac{3\gamma -2\theta }{2\gamma }-\frac{2\gamma -2\theta }{\gamma }\xi +\frac{\gamma -2\theta }{2\gamma }{\xi }^2\right)}^{\gamma },0\le \theta \le \min \left\{\gamma ,\frac{1}{2}\right\}$

广义BDF2-$\theta$ 的卷积系数${\left\{{\omega }_j^{\left(\gamma \right)}\right\}}_{j=0}^{\infty }$ [27]为：

$\begin{array}{l}{\omega }_0^{\left(\gamma \right)}={\left(\frac{3\gamma -2\theta }{2\gamma }\right)}^{\gamma },{\omega }_1^{\left(\gamma \right)}=2\left(\theta -\gamma \right){\left(\frac{2\gamma }{3\gamma -2\theta }\right)}^{1-\gamma },\\ {\omega }_j^{\left(\gamma \right)}=\frac{2\gamma }{j\left(3\gamma -2\theta \right)}\left[2\left(\gamma -\theta \right)\left(\frac{j-1}{\gamma }-1\right){\omega }_{j-1}^{\left(\gamma \right)}+\left(\gamma -2\theta \right)\left(1-\frac{j-2}{2\gamma }\right){\omega }_{j-2}^{\left(\gamma \right)}\right],j\ge 2.\end{array}$ (7)

引入辅助变量$\sigma =\text{Δ}\overline{u}$，方程(4)可以写成如下两个低阶方程的耦合系统

$\{\begin{array}{l}{\overline{u}}_t-{}_0^c{D}_t^{\alpha }\text{Δ}\overline{u}+{}_0^c{D}_t^{\beta }\text{Δ}\sigma -\text{Δ}\overline{u}+\epsilon \cdot \nabla f\left(\overline{u}\right)=\overline{g}\left(x,y,t\right),\\ \sigma =\text{Δ}\overline{u},\end{array}$ (8)

在(8)的第一个方程两边同时乘以$v$ ，在第二个方程两边同时乘以$w$ ，并利用分部积分可得如下弱形式：

找到$\left(\overline{u},\sigma \right)\in H_0^1\times H_0^1$ ，使得对任意的$\left(v,w\right)\in H_0^1\times H_0^1$ ，有

$\{\begin{array}{l}\left({\overline{u}}_t,v\right)+\left({}_0^R{D}_t^{\alpha }\nabla \overline{u},\nabla v\right)-\left({}_0^R{D}_t^{\beta }\nabla \sigma ,\nabla v\right)+\left(\nabla \overline{u},\nabla v\right)+\left(\epsilon \cdot \nabla f\left(\overline{u}\right),v\right)=\left(\overline{g},v\right),\\ \left(\sigma ,w\right)+\left(\nabla \overline{u},\nabla w\right)=0,\end{array}$ (9)

令${\overline{u}}^n=\overline{u}\left(x,y,t_n\right),{\sigma }^n=\sigma \left(x,y,t_n\right)$，在$t_{n-\theta }$ 处，$u\left(t_{n-\theta }\right),g\left(t_{n-\theta }\right),f\left(u\left(t_{n-\theta }\right)\right)$ 分别有如下逼近公式[27] [28]：

$\begin{array}{c}u\left(t_{n-\theta }\right)=u^{n-\theta }+S_{\tau ,\kappa }^{0,n}u+O\left({\tau }^2\right)\\ \triangleq \left(1-\theta \right)u^n+\theta u^{n-1}+S_{\tau ,\kappa }^{0,n}u+O\left({\tau }^2\right),\end{array}$ (10)

$g\left(t_{n-\theta }\right)=g^{n-\theta }+O\left({\tau }^2\right),$ (11)

$f\left(u\left(t_{n-\theta }\right)\right)=f\left(u^{n-\theta }\right)+O\left({\tau }^2\right),$ (12)

其中$g^{n-\theta }\triangleq \left(1-\theta \right)g^n+\theta g^{n-1},f\left(u^{n-\theta }\right)\triangleq \left(2-\theta \right)f\left(u^n\right)-\left(1-\theta \right)f\left(u^{n-1}\right)$ 。针对非线性项$f\left(u\right)$ ，多数文献采用的处理方法是线性化和非线性牛顿迭代法。本文采用Liu等人在文献[29]中提出的在$n-\theta$ 层线性化处理的方法。用已知层来表示未知层，在数值算法上很大程度地提高了计算效率。使用牛顿迭代法虽然可以得到更高的收敛精度，但可能耗费的时间更久。在该模型中我们用线性化处理的方法已经可以达到所需的精度，因此无需再使用牛顿迭代法。

只考虑离散卷积部分，由广义BDF2-$\theta$ 逼近公式得到(9)在$t_{n-\theta }$ 处的时间半离散格式：

$\{\begin{array}{l}\left({\text{Ψ}}_{\tau }^{1,n}\overline{u},v\right)+\left({\text{Ψ}}_{\tau }^{\alpha ,n}\nabla \overline{u},\nabla v\right)-\left({\text{Ψ}}_{\tau }^{\beta ,n}\nabla \sigma ,\nabla v\right)+\left(\nabla {\overline{u}}^{n-\theta },\nabla v\right)=\left(\epsilon \cdot \nabla f\left({\overline{u}}^{n-\theta }\right),v\right)+\left({\overline{g}}^{n-\theta },v\right),\\ \left({\sigma }^{n-\theta },w\right)+\left(\nabla {\overline{u}}^{n-\theta },\nabla w\right)=0,\end{array}$ (13)

在此基础上，利用混合虚拟元方法对空间进行离散从而得到全离散格式。为此我们需要引入虚拟元

空间。${\mathcal{T}}_h$ 是区域$\Omega$ 的剖分，${ϵ}_h$ 是${\mathcal{T}}_h$ 的边$e$ 构成的集合。$h=\underset{E\in {\mathcal{T}}_h}{\max }{h}_E$ ，$h_E$ 是单元$E$ 的直径，对于任意的单元$E\in {\mathcal{T}}_h$ ，对于$k\in N$ ，令${\mathbb{P}}_k\left(E\right)$ 为单元$E$ 上次数不超过$k$ 的多项式，${\mathbb{P}}_{-1}=\left\{0\right\}$ 。如下给出两个假设[19]：

(A1) 单元$E$ 的形状为星形。

(A2) 单元$E$ 的任意两个顶点之间的距离$\ge \rho h_E$ 。

由文献[18]，考虑$k\ge 1$ 时的局部函数空间$V_h^E$ 定义为：

$V_h^E=\left\{v_h\in H^1\left(E\right):v_h\in C^0\left(\partial E\right),\text{Δ}{v}_h\in {\mathbb{P}}_k\left(E\right),{v_h|}_e\in {\mathbb{P}}_k\left(e\right),\forall e\subset \partial E\right\}.$

在函数空间$V_h^E$ 中选择如下自由度：

(d1) $v_h$ 在单元顶点处的取值。

(d2) 对于$k>1$ ，$v_h$ 在每条边$e$ 上的$k-1$ 个均匀间隔点处的取值。

(d3) 对于$k>1$ ，$v_h$ 在单元$E$ 上的$k-2$ 阶矩，即$\frac{1}{\left|E\right|}{\displaystyle {\int }_{E}m\left(x\right)v_h\left(x\right)\text{d}x}$ $\forall m\in {ℳ}_{k-2}\left(E\right)$ ，其中

${ℳ}_{k-2}=\left\{{\left(\frac{x-x_E}{h_E}\right)}^s,\left|s\right|\le k-2\right\}$

$x_E$ 为单元的重心坐标。定义投影算子${\text{Π}}_k^{\nabla ,E}:V_h^E\to {\mathbb{P}}_k\left(E\right)\subset V_h^E$ ：

$\{\begin{array}{l}{\left(\nabla \left(u-{\text{Π}}_k^{\nabla ,E}u\right),\nabla q_k\right)}_E=0,\forall q_k\in {\mathbb{P}}_k\left(E\right),\hfill \\ P_0\left({\text{Π}}_k^{\nabla ,E}{u}_h\right)=P_0{u}_h,u_h\in V_h^E,\hfill \end{array}$

其中$P_0$ 是从$L^2\left(E\right)$ 到${\mathbb{P}}_0\left(E\right)$ 的投影。定义一个标准的$L^2$ 投影算子${\text{Π}}_k^{0,E}:V_h^E\to {\mathbb{P}}_k\left(E\right)\subset V_h^E$ ：

${\left(u-{\text{Π}}_k^{0,E}u,q_k\right)}_E=0,\forall q_k\in {\mathbb{P}}_k\left(E\right).$

现在，我们来定义所需的局部虚拟元空间：

$W_h^E=\left\{w_h\in V_h^E:{\displaystyle {\int }_E\left({\text{Π}}_k^{\nabla ,E}{w}_h\right)q_k\text{d}x}={\displaystyle {\int }_E{w}_h{q}_k\text{d}x},\forall q_k\in {\mathbb{P}}_k\backslash {\mathbb{P}}_{k-2}\left(E\right)\right\},$

其中${\mathbb{P}}_k\backslash {\mathbb{P}}_{k-2}\left(E\right)$ 表示多项式等价类的商空间。将局部虚拟元空间拼接到一起构成全局虚拟元空间

$W_h=\left\{w\in H_0^1\left(\text{Ω}\right):{w|}_E\in W_h^E,\forall E\in {\mathcal{T}}_h\right\}.$

在虚拟元空间中定义全局双线性形式$a_h\left(\cdot ,\cdot \right):W_h\times W_h\to ℛ$ 和$m_h\left(\cdot ,\cdot \right):W_h\times W_h\to ℛ$

$a_h\left(u_h,v_h\right)={\displaystyle \sum _{E\in {\mathcal{T}}_h}{a}_h^E}\left(u_h,v_h\right),\forall u_h,v_h\in W_h^E,$

$m_h\left(u_h,v_h\right)={\displaystyle \sum _{E\in {\mathcal{T}}_h}{m}_h^E}\left(u_h,v_h\right),\forall u_h,v_h\in W_h^E,$

式中，$a_h^E\left(\cdot ,\cdot \right):W_h^E\times W_h^E\to ℛ$ 和$m_h^E\left(\cdot ,\cdot \right):W_h^E\times W_h^E\to ℛ$ 均为局部离散双线性形式：

$a_h^E\left(u_h,v_h\right)=a^E\left({\text{Π}}_k^{\nabla ,E}{u}_h,{\text{Π}}_k^{\nabla ,E}{v}_h\right)+s_a^E\left(\left(I-{\text{Π}}_k^{\nabla ,E}\right)u_h,\left(I-{\text{Π}}_k^{\nabla ,E}\right)v_h\right),\forall u_h,v_h\in W_h^E,$

$m_h^E\left(u_h,v_h\right)=m^E\left({\text{Π}}_k^{0,E}{u}_h,{\text{Π}}_k^{0,E}{v}_h\right)+s_m^E\left(\left(I-{\text{Π}}_k^{0,E}\right)u_h,\left(I-{\text{Π}}_k^{0,E}\right)v_h\right),\forall u_h,v_h\in W_h^E,$

$a_h^E\left(\cdot ,\cdot \right)$ 和$m_h^E\left(\cdot ,\cdot \right)$ 为了满足可计算以及分别逼近两个局部连续双线性形式$a^E\left(\cdot ,\cdot \right)$ 和$m^E\left(\cdot ,\cdot \right)$ ，需满足[30]：

(1) 多项式一致性：对所有的$h>0$ ，$E\in {\mathcal{T}}_h$

$a_h^E\left(p,\upsilon \right)=a^E\left(p,\upsilon \right)\forall p\in {\mathbb{P}}_k\left(E\right),\upsilon \in W_h^E$ ，

$m_h^E\left(p,\upsilon \right)=m^E\left(p,\upsilon \right)\forall p\in {\mathbb{P}}_k\left(E\right),\upsilon \in W_h^E$ ，

(2) 稳定性：存在两对不依赖于$h$ 但依赖网格正则性的常数$\left({\alpha }_{\ast },{\alpha }^{\ast }\right)$ 和$\left({\beta }_{\ast },{\beta }^{\ast }\right)$ ，$0<{\alpha }_{\ast }\le {\alpha }^{\ast }$ ， $0<{\beta }_{\ast }\le {\beta }^{\ast }$ ，有：

${\alpha }_{\ast }{a}^E\left(v,v\right)\le a_h^E\left(v,v\right)\le {\alpha }^{\ast }{a}^E\left(v,v\right)$ ，${\beta }_{\ast }{m}^E\left(v,v\right)\le m_h^E\left(v,v\right)\le {\beta }^{\ast }{m}^E\left(v,v\right)$ .

同时为处理右端项，定义$g_h\left(x,y,t\right)$ ：

$g_h\left(x,y,t\right)={\text{Π}}_k^{0}g\left(x,y,t\right)$ ，${\left({\text{Π}}_k^{0}g\right)|}_E={\text{Π}}_k^{0,E}\left({g|}_E\right)$ .

基于上述虚拟元空间$W_h$ ，得到的全离散格式：找$\left(u_h^{n-\theta },{\sigma }_h^{n-\theta }\right)\in W_h\times W_h$ ，使得对任意的$\left(v_h,w_h\right)\in W_h\times W_h$ 满足：

$\{\begin{array}{l}{m}_h\left({\text{Ψ}}_{\tau }^{1,n}{\overline{u}}_h,v_h\right)+a_h\left({\text{Ψ}}_{\tau }^{\alpha ,n}{\overline{u}}_h,v_h\right)-a_h\left({\text{Ψ}}_{\tau }^{\beta ,n}{\sigma }_h,v_h\right)+a_h\left({\overline{u}}_h^{n-\theta },v_h\right)=\left(\epsilon \cdot \nabla f\left({\overline{u}}_h^{n-\theta }\right),v_h\right)+\left({\overline{g}}_h^{n-\theta },v_h\right).\\ m_h\left({\sigma }_h^{n-\theta },w_h\right)+a_h\left({\overline{u}}_h^{n-\theta },w_h\right)=0.\end{array}$ (14)

3. 算法实现

在本节中，我们取$k=1$ 来具体阐述该算法的步骤。

当$k=1$ 时，此时虚拟元空间仅有顶点处的自由度，没有内部自由度。设${\phi }_i\in W_h$ ，$\left(i=1,2,\cdots ,m\right)$ ，$m=dim\left(W_h\right)$ ，是虚拟元空间的一组基。令$v_h={\phi }_j$ ，$w_h={\phi }_j$ ，$\left(j=1,2,\cdots ,m\right)$

$u_h^{n-\theta }={\displaystyle \sum _{i=1}^m{X}_i^{n-\theta }{\phi }_i},$ (15)

${\sigma }_h^{n-\theta }={\displaystyle \sum _{i=1}^m{Y}_i^{n-\theta }{\phi }_i}.$ (16)

将$v_h,w_h$ 和(15)、(16)代入(14)，并利用卷积系数(7)、逼近公式(12)、(11)、(10)对分数阶项、非线性项$f$ 、源项$\overline{g}$ 及${\overline{u}}_h^{n-\theta }$ 进行处理得到

$\begin{array}{l}{\tau }^{-1}{\omega }_0^{\left(1\right)}M{X}^n+{\tau }^{-\alpha }{\omega }_0^{\left(\alpha \right)}A{X}^n+\left(1-\theta \right)A{X}^n-{\tau }^{-\beta }{\omega }_0^{\left(\beta \right)}A{Y}^n\\ =F^{n-\theta }+G^{n-\theta }-{\tau }^{-1}{\displaystyle \sum _{k=1}^n{\omega }_k^{\left(1\right)}M{X}^{n-k}}-{\tau }^{-\alpha }{\displaystyle \sum _{k=1}^n{\omega }_k^{\left(\alpha \right)}A{X}^{n-k}}-\theta A{X}^{n-1}+{\tau }^{-\beta }{\displaystyle \sum _{k=1}^n{\omega }_k^{\left(\beta \right)}A{Y}^{n-k}}\end{array}$

$\left(1-\theta \right)A{X}^n+\left(1-\theta \right)M{Y}^n=-\theta A{X}^{n-1}-\theta M{Y}^{n-1}.$

其中，

$A={\left[a_h\left({\phi }_i,{\phi }_j\right)\right]}_{1\le i,j\le m}^{\text{T}},M={\left[m_h\left({\phi }_i,{\phi }_j\right)\right]}_{1\le i,j\le m}^{\text{T}},$

对于刚度矩阵$A$ 以及质量矩阵$M$ 中离散双线性形式的具体实现过程可见文章[31]。

$X^n={\left[X_1^n,X_2^n,\cdots ,X_m^n\right]}^{\text{T}}$ ，$Y^n={\left[Y_1^n,Y_2^n,\cdots ,Y_m^n\right]}^{\text{T}}$ ，

$F_1^{n-\theta }=\left[\begin{array}{c}{ϵ}_1\left(2-\theta \right)\left(f{\left({\overline{u}}_h^{n-1}\right)}_x,{\phi }_1\right)-{ϵ}_1\left(1-\theta \right)\left(f{\left({\overline{u}}_h^{n-2}\right)}_x,{\phi }_1\right)\\ {ϵ}_1\left(2-\theta \right)\left(f{\left({\overline{u}}_h^{n-1}\right)}_x,{\phi }_2\right)-{ϵ}_1\left(1-\theta \right)\left(f{\left({\overline{u}}_h^{n-2}\right)}_x,{\phi }_2\right)\\ ⋮\\ {ϵ}_1\left(2-\theta \right)\left(f{\left({\overline{u}}_h^{n-1}\right)}_x,{\phi }_m\right)-{ϵ}_1\left(1-\theta \right)\left(f{\left({\overline{u}}_h^{n-2}\right)}_x,{\phi }_m\right)\end{array}\right]$ ，

$F_2^{n-\theta }=\left[\begin{array}{c}{ϵ}_2\left(2-\theta \right)\left(f{\left({\overline{u}}_h^{n-1}\right)}_y,{\phi }_1\right)-{ϵ}_2\left(1-\theta \right)\left(f{\left({\overline{u}}_h^{n-2}\right)}_y,{\phi }_1\right)\\ {ϵ}_2\left(2-\theta \right)\left(f{\left({\overline{u}}_h^{n-1}\right)}_y,{\phi }_2\right)-{ϵ}_2\left(1-\theta \right)\left(f{\left({\overline{u}}_h^{n-2}\right)}_y,{\phi }_2\right)\\ ⋮\\ {ϵ}_2\left(2-\theta \right)\left(f{\left({\overline{u}}_h^{n-1}\right)}_y,{\phi }_m\right)-{ϵ}_2\left(1-\theta \right)\left(f{\left({\overline{u}}_h^{n-2}\right)}_y,{\phi }_m\right)\end{array}\right]$ ，

$G^{n-\theta }=\left[\begin{array}{c}\left(1-\theta \right)\left({\overline{g}}_h^n,{\phi }_1\right)+\theta \left({\overline{g}}_h^{n-1},{\phi }_1\right)\\ \left(1-\theta \right)\left({\overline{g}}_h^n,{\phi }_2\right)+\theta \left({\overline{g}}_h^{n-1},{\phi }_2\right)\\ ⋮\\ \left(1-\theta \right)\left({\overline{g}}_h^n,{\phi }_m\right)+\theta \left({\overline{g}}_h^{n-1},{\phi }_m\right)\end{array}\right]$ ，

$\begin{array}{l}{F}^{n-\theta }=F_1^{n-\theta }+F_2^{n-\theta },\\ K_1={\tau }^{-1}{\omega }_0^{\left(1\right)}M+{\tau }^{-\alpha }{\omega }_0^{\left(\alpha \right)}A+\left(1-\theta \right)A,\\ K_2=-{\tau }^{-\beta }{\omega }_0^{\left(\beta \right)}A,\\ K_3=\left(1-\theta \right)A,\\ K_4=\left(1-\theta \right)M,\\ b_1^1=F^{n-\theta }+G^{n-\theta }-\theta A{X}^{n-1},\end{array}$

$b_1^2={\tau }^{-1}{\displaystyle \sum _{k=1}^n{\omega }_k^{\left(1\right)}M{X}^{n-k}}-{\tau }^{-\alpha }{\displaystyle \sum _{k=1}^n{\omega }_k^{\left(\alpha \right)}A{X}^{n-k}}+{\tau }^{-\beta }{\displaystyle \sum _{k=1}^n{\omega }_k^{\left(\beta \right)}A{Y}^{n-k}}$ ，

$b_1=b_1^1-b_1^2$ ，

$b_2=-\theta A{X}^{n-1}-\theta M{Y}^{n-1}$ ，

可得

$\left[\begin{array}{cc}{K}_1& K_2\\ K_3& K_4\end{array}\right]\left[\begin{array}{c}{X}^n\\ Y^n\end{array}\right]=\left[\begin{array}{c}{b}_1\\ b_2\end{array}\right]$ ，

记作：

$KZ=b,Z={\left[X^n,Y^n\right]}^{\text{T}},b={\left[b_1,b_2\right]}^{\text{T}}$

处理边界条件和右端项后，得到全离散格式的矩阵形式：

$\widehat{K}Z=\widehat{b}$ ，

通过求解线性方程组即可得到数值解。

4. 数值算例

在本节中，我们将通过一个二维算例来验证算法的有效性并得到收敛阶。在文献[32]中，Yang与Liu利用混合有限元算法针对该模型的一维情况做出了稳定性分析并给出数值模拟结果，本文的模拟结果与文献[32]中的结果是一致的。

令$\overline{\text{Ω}}=\left[0,1\right]\times \left[0,1\right]$，时间区间$I=\left(0,1\right]$ ，精确解$u\left(x,y,t\right)=t^2\sin \left(2\pi x\right)\sin \left(2\pi y\right)$ ，非线性项$f\left(u\right)=u^2$ 和源项$g\left(x,y,t\right)$

$\begin{array}{l}g\left(x,y,t\right)=\sin \left(2\pi x\right)\sin \left(2\pi y\right)\left(2t+\frac{16{\pi }^2{t}^{2-\alpha }}{\text{Γ}\left(3-\alpha \right)}+\frac{128{\pi }^4{t}^{2-\beta }}{\text{Γ}\left(3-\beta \right)}+8{\pi }^2{t}^2\right)\\ +2\pi t^4\left({\epsilon }_1\sin \left(4\pi x\right){\left(\sin \left(2\pi y\right)\right)}^2+{\epsilon }_2\sin \left(4\pi y\right){\left(\sin \left(2\pi x\right)\right)}^2\right).\end{array}$

Table 1. The $L^2$ error and spatial convergence order with $t=1$ and $\alpha =0.2,\beta =0.8,\tau =1/100$

表1. $t=1$ 时刻，固定$\alpha =0.2,\beta =0.8,\tau =1/100$ 的$L^2$ 误差和空间收敛阶

$\theta$	$1/h$	$\Vert u_h-u\Vert$	收敛阶	$\Vert {\sigma }_h-\sigma \Vert$	收敛阶
0.2	10	5.5249E−02	-	3.9939E+00	-
	20	1.9186E−02	1.5259	1.0423E+00	1.9380
	40	5.1989E−03	1.8838	2.6361E−01	1.9833
	80	1.3186E−03	1.9792	6.5443E−02	2.0101
−0.5	10	5.5220E−02	-	3.9915E+00	-
	20	1.9152E−02	1.5277	1.0397E+00	1.9408
	40	5.1640E−03	1.8909	2.6094E−01	1.9944
	80	1.2835E−03	2.0084	6.2783E−02	2.0553
−1	10	5.5192E−02	-	3.9892E+00	-
	20	1.9120E−02	1.5294	1.0372E+00	1.9434
	40	5.1309E−03	1.8978	2.5842E−01	2.0049
	80	1.2503E−03	2.0370	6.0282E−02	2.0999

虚拟元相较于传统有限元的最大好处之一就是可以选择凸或者非凸多边形网格，于是通过此算例我们给出，固定时间步长$\tau =1/100$ ，在非凸多边形网格下的数值结果。空间步长分别为1/10，1/20，1/40，1/80。表1~3给出在$t=1$ 时，针对不同的分数阶参数$\alpha ,\beta$ 以及$\theta$ 得到的$L^2$ 误差和空间收敛阶，可以看到在非凸多边形网格下空间收敛阶为$O\left(h^2\right)$ ，该算例充分验证了算法的有效性和收敛性。

Table 2. The $L^2$ error and spatial convergence order with $t=1$ and $\alpha =0.5,\beta =0.5,\tau =1/100$

表2. $t=1$ 时刻，固定$\alpha =0.5,\beta =0.5,\tau =1/100$ 的$L^2$ 误差和空间收敛阶

$\theta$	$1/h$	$\Vert u_h-u\Vert$	收敛阶	$\Vert {\sigma }_h-\sigma \Vert$	收敛阶
0.5	10	5.5165E−02	-	3.9924E+00	-
	20	1.9131E−02	1.5279	1.0396E+00	1.9412
	40	5.1897E−03	1.8822	2.6331E−01	1.9812
	80	1.3239E−03	1.9709	6.5947E−02	1.9974
0.2	10	5.5161E−02	-	3.9921E+00	-
	20	1.9127E−02	1.5281	1.0393E+00	1.9415
	40	5.1856E−03	1.8830	2.6300E−01	1.9825
	80	1.3000E−03	1.9960	6.5600E−02	2.0033
−1	10	5.5077E−02	-	3.9853E+00	-
	20	1.9032E−02	1.5331	1.0319E+00	1.9493
	40	5.0871E−03	1.9035	2.5548E−01	2.0140
	80	1.2207E−03	2.0591	5.8179E−02	2.1347

Table 3. The $L^2$ error and spatial convergence order with $t=1$ and $\alpha =0.8,\beta =0.2,\tau =1/100$

表3. $t=1$ 时刻，固定$\alpha =0.8,\beta =0.2,\tau =1/100$ 的$L^2$ 误差和空间收敛阶

$\theta$	$1/h$	$\Vert u_h-u\Vert$	收敛阶	$\Vert {\sigma }_h-\sigma \Vert$	收敛阶
0.2	10	5.5049E−02	-	3.9894E+00	-
	20	1.9048E−02	1.5310	1.0351E+00	1.9464
	40	5.1645E−03	1.8830	2.6187E−01	1.9828
	80	1.3168E−03	1.9716	6.5533E−02	1.9986
0	10	5.5047E−02	-	3.9893E+00	-
	20	1.9047E−02	1.5311	1.0350E+00	1.9466
	40	5.1626E−03	1.8833	2.6173E−01	1.9834
	80	1.3149E−03	1.9731	6.5400E−02	2.0007
−1	10	5.4853E−02	-	3.9735E+00	-
	20	1.8827E−02	1.5428	1.0180E+00	1.9647
	40	4.9357E−03	1.9315	2.4451E−01	2.0577
	80	1.0885E−03	2.1809	4.8689E−02	2.3282

5. 总结

本文在时间方向上利用广义BDF2-$\theta$ 格式，空间方向上采用混合虚拟元方法求解二维非线性四阶时间分数阶方程，通过引入辅助变量将高阶方程转化为低阶耦合系统，很大程度地降低了计算的难度。并且相比于传统的有限元方法，虚拟元方法有很好的网格适应性和数值稳定性，降低了复杂几何区域上网格生成的难度。通过数值算例也验证了该算法的有效性。

NOTES

^*通讯作者。