非线性时间分数阶四阶混合次扩散和扩散波动方程的混合有限元算法

Mixed Finite Element Algorithm for a Nonlinear Time-Fractional Fourth-Order Mixed Sub-Diffusion and Diffusion Wave Equation

摘要: 本文数值求解了一个二维非线性时间分数阶四阶混合次扩散和扩散波动方程，在时间方向上采用L1-CN格式，在空间上通过混合有限元方法进行离散，并且在此基础上，给出了它的全离散格式。最后针对该数值格式提供了算法过程和数值算例，以及详细的收敛结果。

Abstract: In this article, the two-dimensional nonlinear time-fractional fourth-order mixed sub-diffusion and diffusion wave equation is numerically solved, where L1-CN scheme is used in the time direction, and the mixed finite element method is applied in space. Further, the fully discrete scheme is provided. Finally, the algorithm process, numerical example, and detailed convergence results are provided for this numerical scheme.

1. 引言

分数阶扩散方程作为分数阶偏微分方程的一个重要组成部分，常被用于解释不同物质之间的扩散现象，随着研究深入，含扩散项的方程不断演变，其中带有单个时间分数阶的扩散类方程最为常见，时间分数阶为Caputo型 [1] [2] [3] [4] 或Riemann-Liouville型 [5] [6] 居多。

随着研究更加深入，人们发现带有两项时间分数阶的扩散或扩散波动方程更能准确地把握物质的属性以及现象发生的过程。在 [7] 中，Feng等利用有限元方法首次研究了二维时间分数阶混合次扩散和扩散波动方程。近些年混合次扩散和扩散波动方程的研究理论不断完善，其中，Chen等 [8] 利用空间两网格有限元方法结合L1-CN格式解决了一个二维非线性时间分数阶两项混合次扩散和扩散波动方程。在 [3] 中，Sun等将有限差分法与L1格式相结合数值离散了一个分数阶混合次扩散和扩散波动方程。在 [9] 中，Shen等将单项时间分数阶扩散波动方程的H2N2方法扩展到带有非齐次原项的时间分数阶混合次扩散和扩散波动方程上。由此可知，L1-CN格式及其改进格式能有效求解混合次扩散和扩散方程。除此之外，不同形式的双分数阶方程多年来同样被大量学者重视。在 [10] 中，Yin等人在卷积求积框架下，利用Caputo型分数阶导数与Riemann-Liouville型分数阶导数的关系化解了系数为(1, 2)的分数阶导数卷积求解难的问题。纵然双分数阶问题被大量研究，但多年来该问题的四阶情况从未被考虑，因此选取合适的数值方法来处理含有四阶项的双分数阶问题具有非常重要的意义。据了解目前处理高阶问题常用的方法有混合有限元方法 [11] [12] ，有限差分法 [13] ，间断有限元方法 [14] 等。

在本文中，我们在空间上选用混合有限元方法，在时间上利用L1-CN格式考虑以下二维非线性四阶时间分数阶混合次扩散和扩散波动方程：

${}_0{}^{C}D{}_t^{\alpha }u+{}_0{}^{C}D{}_t^{\beta }u+{\Delta }^{2}u-\Delta u_t-\Delta u+f\left(u\right)=g\left(x,y,t\right),\left(x,y,t\right)\in \Omega \times J,$ (1)

满足边界条件

$u\left(x,y,t\right)=\Delta u\left(x,y,t\right)=0,\left(x,y\right)\in \partial \bar{\Omega },t\in \bar{J},$ (2)

和初始条件

$u\left(x,y,0\right)=u_0,u_t\left(x,y,0\right)=u_1,\left(x,y\right)\in \bar{\Omega }.$ (3)

其中 $\left(x,y\right)\in \Omega$ ， $\Omega$ 是一个凸多边形区域， $\partial \bar{\Omega }$ 是 $\bar{\Omega }$ 的边界， $J=\left(0,T\right]$ 是时间区间，并且 $0<T<\infty$ 。

${}_0{}^{C}D{}_t^{\alpha }u$ 和 ${}_0{}^{C}D{}_t^{\beta }u$ 均为Caputo型分数阶导数，有如下定义：

${}_0{}^{C}D{}_t^{\alpha }u=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^t\frac{\partial u\left(x,y,s\right)}{\partial s}}{\left(t-s\right)}^{\alpha }\text{d}s,\alpha \in \left(0,1\right),$

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

本文的重点如下：

1) 二维非线性时间分数阶四阶混合次扩散和扩散波动方程首次被研究。

2) 本文以混合有限元和L1-CN格式相结合的方法处理四阶混合次扩散和扩散波动方程。同时还呈现了该方程详细的混合有限元算法过程，以及算例和数据结果。

本文结构如下：

在第二节中，通过引入辅助变量 $q=\Delta u$ ，将四阶问题转化为二阶耦合系统，然后时间方向上结合L1-CN近似公式给出相对应的全离散格式。在第三节中，基于耦合系统的全离散格式，给出方程(1)~(3)的详细算法过程，并呈现了算例的误差和收敛结果。在第四节中，得出了一些结论。

2. 数值格式

在这一节中，我们在时间方向上利用L1-CN近似公式的引理结合空间混合有限元方法得出方程(1)~(3)的弱形式和全离散格式。

为了得到全离散格式，我们将时间区间 $\left[0,T\right]$ 均匀剖分为 $0=t_0<t_1<t_2<\cdots <t_N=T$ ，N为正整数，其中 $t_n=n\tau \left(n=1,2,\cdots ,N\right)$ ， $\tau =T/N$ 为时间步长。定义 $u^n=u\left(\cdot ,t_n\right)$ ，我们给出 $t_{n-\frac{1}{2}}$ 处的几个离散形式：

$u^{n-\frac{1}{2}}=\frac{u^n+u^{n-1}}{2}+O\left(\Delta t^2\right),u_t^{n-\frac{1}{2}}=\frac{u^n-u^{n-1}}{\tau }+O\left(\Delta t^2\right),$

${\partial }_t{u}^{n-\frac{1}{2}}=\frac{u^n-u^{n-1}}{\tau }+O\left(\Delta t^2\right).$

引理1

Caputo型导数 ${}_0{}^{C}D{}_t^{\alpha }u$ 和 ${}_0{}^{C}D{}_t^{\beta }u$ 的L1-CN近似公式 [4] [8] 如下所示：

$D_t^{\alpha }{u}^{n-\frac{1}{2}}=\frac{{\tau }^{1-\alpha }}{2\Gamma \left(2-\alpha \right)}\left({\displaystyle \underset{k=1}{\overset{n}{\sum }}{a}_{n-k}^{\left(\alpha \right)}}{\partial }_t{u}^{k-\frac{1}{2}}+{\displaystyle \underset{k=1}{\overset{n-1}{\sum }}{a}_{n-1-k}^{\left(\alpha \right)}}{\partial }_t{u}^{k-\frac{1}{2}}\right),$

$D_t^{\beta }{u}^{n-\frac{1}{2}}=\frac{{\tau }^{1-\beta }}{\Gamma \left(3-\beta \right)}\left(b_0^{\left(\beta \right)}{\partial }_t{u}^{n-\frac{1}{2}}-{\displaystyle \underset{k=1}{\overset{n-1}{\sum }}\left(b_{n-1-k}^{\left(\beta \right)}-b_{n-k}^{\left(\beta \right)}\right){\partial }_t{u}^{k-\frac{1}{2}}}-b_{n-1}^{\left(\beta \right)}{u}_1\right),$

其中 $a_k^{\left(\alpha \right)}={\left(k+1\right)}^{1-\alpha }-k^{1-\alpha }$ ， $b_k^{\left(\beta \right)}={\left(k+1\right)}^{2-\beta }-k^{2-\beta }$ ， $k=0,1,\cdots ,n-1$ 。

接下来引入一个辅助变量 $q=\Delta u$ ，方程(1)可以改写为：

${}_0{}^{C}D{}_t^{\alpha }u+{}_0{}^{C}D{}_t^{\beta }u+\Delta q-\Delta u_t-\Delta u+f\left(u\right)=g\left(x,y,t\right),$ (4)

$q=\Delta u.$ (5)

取 $\left(u,q\right)\in H_0^1\times H_0^1$ ， $\left(v,w\right)\in H_0^1\times H_0^1$ ，(4)~(5)在 $t_{n-\frac{1}{2}}$ 处相应的弱形式如下：

$\begin{array}{l}\left({}_0{}^{C}D{}_t^{\alpha }{u}^{n-\frac{1}{2}},v\right)+\left({}_0{}^{C}D{}_t^{\beta }{u}^{n-\frac{1}{2}},v\right)-\left(\nabla q^{n-\frac{1}{2}},\nabla v\right)+\left(\nabla u_t^{n-\frac{1}{2}},\nabla v\right)\\ +\left(\nabla u^{n-\frac{1}{2}},\nabla v\right)+\left(f\left(u^{n-\frac{1}{2}}\right),v\right)=\left(g,v\right),\end{array}$

$\left(q^{n-\frac{1}{2}},w\right)+\left(\nabla u^{n-\frac{1}{2}},\nabla w\right)=0.$

定义有限元空间 $V_h\in H_0^1$ 和 $W_h\in H_0^1$

$V_h=\left\{v_h|v_h\in H_0^1,v_h\in P_1\left(K\right),\forall K\in {\mathcal{K}}_h\right\},$

$W_h=\left\{w_h|w_h\in H_0^1,w_h\in P_1\left(I\right),\forall I\in {\mathcal{I}}_h\right\}.$

利用引理1，取 $u_h^n\in V_h$ ， $q_h^n\in W_h$ ，其中 $n=1,\cdots ,N$ ，可以得到以下全离散格式：

$\begin{array}{l}\left(D_t^{\alpha }{u}_h^{n-\frac{1}{2}},v_h\right)+\left(D_t^{\beta }{u}_h^{n-\frac{1}{2}},v_h\right)-\left(\nabla q_h^{n-\frac{1}{2}},\nabla v_h\right)+\left(\nabla {\partial }_t{u}_h^{n-\frac{1}{2}},\nabla v_h\right)\\ +\left(\nabla u_h^{n-\frac{1}{2}},\nabla v_h\right)+\left(f\left(u_h^{n-\frac{1}{2}}\right),v_h\right)=\left(g,v_h\right),\forall v_h\in V_h,\end{array}$ (6)

$\left(q_h^{n-\frac{1}{2}},w_h\right)+\left(\nabla u_h^{n-\frac{1}{2}},\nabla w_h\right)=0,\forall w_h\in W_h.$ (7)

3. 数值模拟

3.1. 数值算法

对 $\forall u_h^n\in V_h$ ， $q_h^n\in W_h$ ，且 $v_h\in V_h$ ， $w_h\in W_h$ ，我们记，

$V_h=span{\left\{{\varphi }_i\right\}}_{i=1}^{N_h},W_h=span{\left\{{\psi }_j\right\}}_{j=1}^{N_h},$

然后有

$u_h^n={\displaystyle \underset{i=1}{\overset{N_h}{\sum }}{u}_i^n}{\varphi }_i,q_h^n={\displaystyle \underset{i=1}{\overset{N_h}{\sum }}{q}_i^n}{\psi }_i,$

故令

$M={\left[\left({\varphi }_i,{\varphi }_j\right)\right]}_{1\le i,j\le N_h}^{\text{T}},A={\left[\left(\nabla {\varphi }_i,\nabla {\varphi }_j\right)\right]}_{1\le i,j\le N_h}^{\text{T}},$

$B={\left[\left(\nabla {\psi }_i,\nabla {\varphi }_j\right)\right]}_{1\le i,j\le N_h}^{\text{T}},C={\left[\left(\nabla {\varphi }_i,\nabla {\psi }_j\right)\right]}_{1\le i,j\le N_h}^{\text{T}},D={\left[\left({\psi }_i,{\psi }_j\right)\right]}_{1\le i,j\le N_h}^{\text{T}}.$

进而对于(6)~(7)有以下算法方程式成立：

当 $n=1$ 时，

$\{\begin{array}{l}\left(\mu M+\frac{A}{\tau }+\frac{A}{2}\right)U^1-\frac{B}{2}{Q}^1=H_1,\\ C{U}^1+D{Q}^1=-C{U}^0-D{Q}^0.\end{array}$

其中

$U^n={\left[u_1^n,u_2^n,\cdots ,u_{N_h}^n\right]}^{\text{T}},Q^n={\left[q_1^n,q_2^n,\cdots ,q_{N_h}^n\right]}^{\text{T}},$

$G^n={\left[\left(g^n,{\varphi }_1\right),\left(g^n,{\varphi }_2\right),\cdots ,\left(g^n,{\varphi }_{N_h}\right)\right]}^{\text{T}},F^1=\left(f\left(u_h^0\right),v_h\right),$

$\mu =\frac{{\tau }^{-\alpha }}{2\Gamma \left(2-\alpha \right)}{a}_0^{\left(\alpha \right)}+\frac{{\tau }^{-\beta }}{\Gamma \left(3-\beta \right)}{b}_0^{\left(\beta \right)},$

$H_1=\left(\mu M+\frac{A}{\tau }-\frac{A}{2}\right)U^0+\frac{B}{2}{Q}^0+\frac{M{\tau }^{1-\beta }}{\Gamma \left(3-\beta \right)}{b}_0^{\left(\beta \right)}{u}_1-F^1+G^1,$

当 $n>1$ 时，

$\{\begin{array}{l}\left(\mu M+\frac{A}{\tau }+\frac{A}{2}\right)U^n-\frac{B}{2}{Q}^n=H_n,\\ C{U}^n+D{Q}^n=-C{U}^{n-1}-D{Q}^{n-1}.\end{array}$

其中

${\sigma }_{n-k}^{\left(\alpha \right)}=\frac{{\tau }^{1-\alpha }}{2\Gamma \left(2-\alpha \right)}\left(a_{n-1-k}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)}\right),{\sigma }_{n-k}^{\left(\beta \right)}=\frac{{\tau }^{1-\beta }}{\Gamma \left(3-\beta \right)}\left(b_{n-1-k}^{\left(\beta \right)}-b_{n-k}^{\left(\beta \right)}\right),$

$\begin{array}{c}{H}_n=\left(\mu M+\frac{A}{\tau }-\frac{A}{2}\right)U^{n-1}+\frac{B}{2}{Q}^{n-1}+{\displaystyle \underset{k=1}{\overset{n-1}{\sum }}\left({\sigma }_{n-k}^{\left(\beta \right)}-{\sigma }_{n-k}^{\left(\alpha \right)}\right)}M{\partial }_t{U}^{k-\frac{1}{2}}\\ +\frac{M{\tau }^{1-\beta }}{\Gamma \left(3-\beta \right)}{b}_{n-1}^{\left(\beta \right)}{u}_1-F^n+G^n,\end{array}$

$F^n=\left(\frac{3f\left(u_h^{n-1}\right)-f\left(u_h^{n-2}\right)}{2},v_h\right).$

3.2. 数值算例

在本节中，通过数值算例验证数值格式的有效性。

以下方程的精确解表示为 $u\left(x,y,t\right)=t^{3+\alpha +\beta }\sin \left(\pi x\right)\sin \left(\pi y\right)$ ，

$\{\begin{array}{l}{}_0{}^{C}D{}_t^{\alpha }u+{}_0{}^{C}D{}_t^{\beta }u+{\Delta }^{2}u-\Delta u_t-\Delta u+u^2=g_1\left(x,y,t\right),\left(x,y,t\right)\in \Omega \times J,\\ u\left(x,y,t\right)=\Delta u\left(x,y,t\right)=0,\left(x,y\right)\in \partial \Omega ,t\in \bar{J},\\ u\left(x,y,0\right)=0,u_t\left(x,y,0\right)=0,\left(x,y\right)\in \bar{\Omega },\end{array}$

其中 $\bar{\Omega }=\left[0,1\right]\times \left[0,1\right]$ ， $J=\left(0,1\right]$ ， $g_1$ 是已知源项

$\begin{array}{c}{g}_1\left(x,y,t\right)=\sin \left(\pi x\right)\sin \left(\pi y\right)\left(\frac{t^{3+\beta }\Gamma \left(4+\alpha +\beta \right)}{\Gamma \left(4+\beta \right)}+\frac{t^{3+\alpha }\Gamma \left(4+\alpha +\beta \right)}{\Gamma \left(4+\alpha \right)}+4{\pi }^4{t}^{3+\alpha +\beta }\right)\\ +\sin \left(\pi x\right)\sin \left(\pi y\right)\left(2{\pi }^2\left(3+\alpha +\beta \right)t^{2+\alpha +\beta }+2{\pi }^2{t}^{3+\alpha +\beta }\right)\\ +t^{6+2\alpha +2\beta }{\sin }^2\left(\pi x\right){\sin }^2\left(\pi y\right).\end{array}$

在表1中，固定时间步长 $\tau =1/400$ ，选择空间网格步长 $H=h_x\times h_y=1/100,1/400,1/900,1/1600$ ，得到了基于不同参数 $\alpha =0.25$ ， $\beta =1.75$ ， $\alpha =0.5$ ， $\beta =1.5$ 以及 $\alpha =0.75$ ， $\beta =1.25$ 下的空间收敛结果。在表2中，固定空间网格步长 $H=1/180\times 1/180$ ，选择时间步长 $\tau =1/4,1/8,1/12$ ，得到了基于不同参数 $\alpha =0.25$ ， $\beta =1.75$ ， $\alpha =0.5$ ， $\beta =1.5$ 以及 $\alpha =0.75$ ， $\beta =1.25$ 下的时间收敛结果。表1，表2的收敛结果表明数值格式(6)~(7)对于求解方程(1)~(3)是非常有效的。

Table 1. Spatial convergence results with τ = 1 / 400

表1. 当 $\tau =1/400$ 时，空间收敛结果

Table 2. Time convergence results with H = 1 / 32400

表2. 当 $H=1/32400$ 时，空间收敛结果

为了更加直观清楚地观察数值模拟的效果，给出了一些数值解以及误差的图像。图1~3展示了精确解u和数值解 $u_h$ 基于不同参数 $\alpha$ 和 $\beta$ 的图像。在图4，图5中，给出了不同参数 $\alpha$ 和 $\beta$ 下的 $\Vert u-u_h\Vert$ 和 $\Vert q-q_h\Vert$ 的等高线图。最后在图6，图7中，我们结合表1，表2的数据分别画出了 $\Vert u-u_h\Vert$ 和 $\Vert q-q_h\Vert$ 的空间收敛阶图和时间收敛阶图。

Figure 1. The images of numerical solution $u_h$ and exact solution u with $\alpha =0.25$ , $\beta =1.75$

图1. 当 $\alpha =0.25$ ， $\beta =1.75$ 数值解 $u_h$ 与精确解u的图像

Figure 2. The images of numerical solution $u_h$ and exact solution u with $\alpha =0.5$ , $\beta =1.5$

图2. 当 $\alpha =0.5$ ， $\beta =1.5$ 数值解 $u_h$ 与精确解u的图像

Figure 3. The images of numerical solution $u_h$ and exact solution u with $\alpha =0.75$ , $\beta =1.25$

图3. 当 $\alpha =0.75$ ， $\beta =1.25$ 数值解 $u_h$ 与精确解u的图像

Figure 4. Error contour of $\Vert u-u_h\Vert$ with different parameters $\alpha$ and $\beta$

图4. 不同参数 $\alpha$ 和 $\beta$ 下， $\Vert u-u_h\Vert$ 的误差等高线图

Figure 5. Error contour of $\Vert q-q_h\Vert$ with different parameters and $\beta$

图5. 不同参数 $\alpha$ 和 $\beta$ 下， $\Vert q-q_h\Vert$ 的误差等高线图

Figure 6. Spatial convergence order of $\Vert u-u_h\Vert$ and $\Vert q-q_h\Vert$ with different parameters $\alpha$ and $\beta$

图6. 不同参数 $\alpha$ 和 $\beta$ 下， $\Vert u-u_h\Vert$ 和 $\Vert q-q_h\Vert$ 的空间收敛阶图

Figure 7. Spatial convergence order of $\Vert u-u_h\Vert$ and $\Vert q-q_h\Vert$ with different parameters $\alpha$ and $\beta$

图7. 不同参数 $\alpha$ 和 $\beta$ 下， $\Vert u-u_h\Vert$ 和 $\Vert q-q_h\Vert$ 的时间收敛阶图

4. 结论

本文提出了一种数值求解二维非线性时间分数阶四阶混合次扩散和扩散波动方程的方法，具体是在空间方向上采用混合有限元方法，在时间方向上利用L1-CN格式。进一步给出了它的全离散格式以及混合有限元算法过程。最后通过算例及其结果表明我们的方法能够有效地求解含有四阶导数项的混合次扩散和扩散波动方程。