¹呼和浩特民族学院，数学与大数据学院，内蒙古 呼和浩特

²呼和浩特民族学院，计算机科学与信息工程学院，内蒙古 呼和浩特

收稿日期：2022年7月11日；录用日期：2022年8月4日；发布日期：2022年8月12日

摘要

本文基于应力佳点对偶网格剖分以及分片三次Lagrange插值试探函数空间和分片常数检验函数空间的三次有限体积元法和Caputo导数的L1-逼近公式构造数值格式求解一维时间分数阶扩散方程，并证明了格式的L2范数在时间和空间方向分别 $2-\alpha$ 阶和四阶收敛误差估计。通过数值实验验证了理论分析结果以及所提格式的有效性。

关键词

时间分数阶扩散方程，三次有限体积元法，L1-公式

A Cubic Finite Volume Element Method for the Time Fractional Diffusion Equation

Siriguleng He¹, Chelimuge², Jingying Gao¹

¹School of Mathematics and Big Data, Hohhot Minzu College, Hohhot Inner Mongolia

²School of Computer Science and Information Engineering, Hohhot Minzu College, Hohhot Inner Mongolia

Received: Jul. 11^th, 2022; accepted: Aug. 4^th, 2022; published: Aug. 12^th, 2022

ABSTRACT

In this article, the one-dimensional time fractional diffusion equation is solved by a numerical scheme which is constructed by the L1-formula of approximating the Caputo fractional derivative and a cubic finite volume element method. The cubic finite volume element method is based on the optimal stress points dual partition, and the trial function space of piecewise cubic Lagrange interpolation and the test function space of piecewise constant. The L2-norm error estimate of fourth order convergence in space and $2-\alpha$ order convergence in time is proved. Numerical experiments are given to verify the effectiveness of the theoretical analysis results and the proposed scheme.

Keywords:Time Fractional Diffusion Equation, Cubic Finite Volume Element Method, L1-Formula

Copyright © 2022 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

有限体积元法具有局部积分守恒的性质，能保持局部的质量、动量或能量的守恒性，与有限元方法有同阶精度但要比有限元方法便于计算，因此有限体积元法被认为是求解偏微分方程最有效的计算方法之一 [1]。李荣华教授于1982年提出了基于变分形式定义的有限体积元法并奠定了理论基础 [2]。此后，线性元等低次有限体积元法得到了很好的发展，形成了比较完善的理论框架。高次有限体积元法对离散空间的光滑度要求高，其格式构造和理论分析都比较复杂 [3] [4] [5] [6] [7]。基于插值导数超收敛点(在力学上称为应力佳点)的方法，王同科 [8] [9] [10] 等和李永海 [11] [12] 等各自研究了两点边值问题、椭圆和抛物问题的基于应力佳点的高阶有限体积元法，并给出了数值格式的最优L2-范数等误差估计。文献 [13] [14] 结合文献 [9] 的思想针对对流扩散方程分别构造了Crank-Nicolson三次有限体积元格式和时间间断时空有限体积元格式，并分析了所提格式的L2模最优收敛误差估计。这些方法在不增加计算量的情况下，可达到更高的计算精度且具有天然的梯度超收敛性。

本文基于文献 [9] 的体积元构造思想，针对一维时间分数扩散方程 [15] 构造了一类三次有限体积元格式，并分析了格式的L2模最优收敛误差估计。文献 [15] 利用有限元法研究了时间分数阶扩散方程，而本文研究了该方程的三次有限体积元方法。本文方法具有局部积分守恒的性质，与三次有限元方法有同阶精度但要比三次有限元方法便于计算。

2. 三次有限体积元格式

考虑下面的时间分数阶扩散方程

$\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}-\frac{{\partial }^{2}u}{\partial x^2}=f\left(x,t\right),\left(x,t\right)\in \left(a,b\right)\times \left(0,T\right]$， (1)

以及初边值条件

$\{\begin{array}{l}u\left(a,t\right)=0,u\left(b,t\right)=0,t\in \left(0,T\right],\hfill \\ u\left(x,0\right)=\phi \left(x\right),x\in \left(a,b\right),\hfill \end{array}$ (2)

其中 $0<\alpha <1$，$f,\phi$ 是给定的光滑函数，Caputo分数阶导数定义为

$\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^t\frac{\partial u\left(x,s\right)}{\partial s}}\frac{\text{d}s}{{\left(t-s\right)}^{\alpha }}$。 (3)

对正整数N，令 $\tau =T/N$ 是时间步长， $u^m=u\left(x,t_m\right)\left(t_m=m\tau ,0\le m\le N\right)$ 和 ${\bar{u}}^m=\left(u^m+u^{m-1}\right)/2$。分数阶导数(3)在点 $t_m$ 处有 [15]

$\frac{{\partial }^{\alpha }u\left(x,t_m\right)}{\partial t^{\alpha }}=\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left[u^m-{\displaystyle \underset{k=1}{\overset{m-1}{\sum }}{\beta }_k{u}^{m-k}-{\beta }_{m-1}{u}^0}\right]+R^m$，

其中 ${\beta }_k=\left(d_{k-1}-d_k\right),{\beta }_{m-1}=d_{m-1},d_k={\left(k+1\right)}^{1-\alpha }-k^{1-\alpha }$ 且 $R^m=o\left({\tau }^{2-\alpha }\right)$。于是方程(1)的时间半离散形式为

$u^m-{\displaystyle \underset{k=1}{\overset{m-1}{\sum }}{\beta }_k{u}^{m-k}}-{\beta }_{m-1}{u}^0-\Gamma \left(2-\alpha \right){\tau }^{\alpha }{u}_{xx}^m=\Gamma \left(2-\alpha \right){\tau }^{\alpha }{f}^m+\Gamma \left(2-\alpha \right){\tau }^{\alpha }{R}^m$。

对区间 $I=\left[a,b\right]$ 做剖分 $I_h$，节点为 $a=x_0<x_1<\cdots <x_n<x_{n+1}<\cdots <x_{2n}<x_{2n+1}<\cdots <x_{3n}=b$，称 $I_i=\left[x_{3i-3},x_{3i}\right]$ 为单元，单元内的三等分点为 $x_{3i-2},x_{3i-1}$，令 $h_i$ 为每等分长度， $h={\max }_{1\le i\le n}{h}_i$。在区间 $I_i$ 上，经过 $x_{3i-3},x_{3i-2},x_{3i-1},x_{3i}$ 四点的三次插值应力佳点分别为

$x_{3i,l}=x_{3i}-\frac{3+\sqrt{5}}{2}{h}_i,x_{3i,r}=x_{3i}-\frac{3-\sqrt{5}}{2}{h}_i,x_{3i,c}=x_{3i}-\frac{3}{2}{h}_i$。

在这三点处，函数的三次插值导数具有四阶精度，而其余点处仅有三阶精度 [9]。记 $I_i^{*}=\left[x_{3i,l},x_{3i,c}\right]$，$I_i^{**}=\left[x_{3i,c},x_{3i,r}\right]$，$I_i^{***}=\left[x_{3i,r},x_{3\left(i+1\right),l}\right]$，$\left(i=1,2,\cdots ,n-1\right)$，$I_0^{***}=\left[x_0,x_{3,l}\right]$，$I_n^{***}=\left[x_{3n,r},x_{3n}\right]$。所有的 $I_i^{*},I_i^{**},I_i^{***}$ 构成剖分 $I_h$ 的对偶剖分 $I_h^{*}$，称 $I_i^{*},I_i^{**},I_i^{***}$ 分别为对应节点 $x_{3i-2},x_{3i-1},x_{3i}$ 的控制体积。

这里除了用标准Sobolev空间及L2范数 ${\Vert \cdot \Vert }_0$ 和H1半范数 ${\left|\cdot \right|}_1$，需要如下离散空间

$U_h=\left\{u_h\in H_0^1\left(I\right):{u_h|}_{I_i}\in P_3\left(I_i\right),\forall I_i\in I_h,u_h\left(a,t\right)=u_h\left(b,t\right)=0\right\}$，

$V_h=\left\{v_h\in L_2\left(I\right):{v_h|}_{I_i^{*}\left(或I_i^{**}或I_i^{***}\right)}=常数,\forall I_i^{*}\left(或I_i^{**}或I_i^{***}\right)\in I_h^{*},{v_h|}_{I_0^{***}}=0,{v_h|}_{I_n^{***}}=0\right\}$，

和如下离散范数

${\Vert u_h\Vert }_{0,h}^2=\frac{3}{8}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{h}_i\left(u_{3i-3}^2+3u_{3i-2}^2+3u_{3i-1}^2+u_{3i}^2\right)},\forall u_h\in U_h$，

${\left|u_h\right|}_{1,h}^2={\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{1}{h_i}\left[{\left(u_{3i-2}-u_{3i-3}\right)}^2+{\left(u_{3i-1}-u_{3i-2}\right)}^2+{\left(u_{3i}-u_{3i-1}\right)}^2\right]},\forall u_h\in U_h$，

其中令 $u_i=u_h\left(x_i\right)$。定义插值算子 ${\pi }_h^{*}:U_h\to V_h$ 使得

${\pi }_h^{*}{u}_h={\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(u_h\left(x_{3i-2}\right){\psi }_{3i-2}+u_h\left(x_{3i-1}\right){\psi }_{3i-1}+u_h\left(x_{3i}\right){\psi }_{3i}\right)}$，

其中 ${\psi }_{3i-2},{\psi }_{3i-1},{\psi }_{3i}$ 分别是控制体积 $I_i^{*},I_i^{**},I_i^{***}$ 上的特征函数。

问题(1)~(2)的有限体积元弱形式为：对 $\forall v\in V_h$，求 $u\in H_0^1\left(I\right)$ 使得

$\begin{array}{l}\left(u^m,v\right)+\Gamma \left(2-\alpha \right){\tau }^{\alpha }a\left(u^m,v\right)\\ ={\displaystyle \underset{k=1}{\overset{m-1}{\sum }}{\beta }_k}\left(u^{m-k},v\right)+{\beta }_{m-1}\left(u^0,v\right)+\Gamma \left(2-\alpha \right){\tau }^{\alpha }\left(f^m,v\right)+\Gamma \left(2-\alpha \right){\tau }^{\alpha }\left(R^m,v\right)\end{array}$。

由于 ${\pi }_h^{*}$ 是 $U_h$ 到 $V_h$ 的迁移算子，上述弱形式等价于：对 $\forall v\in U_h$，求 $u\in H_0^1\left(I\right)$ 使得

$\begin{array}{l}\left(u^m,{\pi }_h^{*}v\right)+\Gamma \left(2-\alpha \right){\tau }^{\alpha }a\left(u^m,{\pi }_h^{*}v\right)\\ ={\displaystyle \underset{k=1}{\overset{m-1}{\sum }}{\beta }_k}\left(u^{m-k},{\pi }_h^{*}v\right)+{\beta }_{m-1}\left(u^0,{\pi }_h^{*}v\right)+\Gamma \left(2-\alpha \right){\tau }^{\alpha }\left(f^m,{\pi }_h^{*}v\right)+\Gamma \left(2-\alpha \right){\tau }^{\alpha }\left(R^m,{\pi }_h^{*}v\right)\end{array}$。 (4)

相应的有限体积元格式为：对 $\forall v\in U_h$，求 $u_h\in U_h$ 使得

$\begin{array}{l}\left(u_h^m,{\pi }_h^{*}v\right)+\Gamma \left(2-\alpha \right){\tau }^{\alpha }a\left(u_h^m,{\pi }_h^{*}v\right)\\ ={\displaystyle \underset{k=1}{\overset{m-1}{\sum }}{\beta }_k}\left(u_h^{m-k},{\pi }_h^{*}v\right)+{\beta }_{m-1}\left(u_h^0,{\pi }_h^{*}v\right)+\Gamma \left(2-\alpha \right){\tau }^{\alpha }\left(f^m,{\pi }_h^{*}v\right)\end{array}$， (5)

其中

$a\left(u_h,{\pi }_h^{*}v\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}\left[v_{3i-2}a\left(u_h,{\psi }_{3i-2}\right)+v_{3i-1}a\left(u_h,{\psi }_{3i-1}\right)\right]}+{\displaystyle \underset{i=1}{\overset{n-1}{\sum }}{v}_{3i}a\left(u_h,{\psi }_{3i}\right)}$，

$\left(u_h,{\pi }_h^{*}v\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}\left[v_{3i-2}\left(u_h,{\psi }_{3i-2}\right)+v_{3i-1}\left(u_h,{\psi }_{3i-1}\right)\right]}+{\displaystyle \underset{i=1}{\overset{n-1}{\sum }}{v}_{3i}\left(u_h,{\psi }_{3i}\right)}$，

且

$a\left(u_h,{\psi }_{3i-2}\right)={u^{\prime }}_h\left(x_{3i,l}\right)-{u^{\prime }}_h\left(x_{3i,c}\right),a\left(u_h,{\psi }_{3i-1}\right)={u^{\prime }}_h\left(x_{3i,c}\right)-{u^{\prime }}_h\left(x_{3i,r}\right),a\left(u_h,{\psi }_{3i}\right)={u^{\prime }}_h\left(x_{3i,r}\right)-{u^{\prime }}_h\left(x_{3\left(i+1\right),l}\right)$

$\left(u_h,{\psi }_{3i-2}\right)={\displaystyle {\int }_{I_i^{*}}{u}_h\text{d}x},\left(u_h,{\psi }_{3i-1}\right)={\displaystyle {\int }_{I_i^{**}}{u}_h\text{d}x},\left(u_h,{\psi }_{3i}\right)={\displaystyle {\int }_{I_i^{***}}{u}_h\text{d}x}$。

注：格式(5)是便于使用有限元框架进行误差分析。在实际计算中，我们可使用与格式(5)等价的下述有限体积格式：求 $u_h\in U_h$ 使得

${\displaystyle {\int }_{I_i^{*}}{u}^m\text{d}x}-{\displaystyle \underset{k=1}{\overset{m-1}{\sum }}{\beta }_k{\displaystyle {\int }_{I_i^{*}}{u}^{m-k}\text{d}x}}-{\beta }_{m-1}{\displaystyle {\int }_{I_i^{*}}{u}^0\text{d}x}-\Gamma \left(2-\alpha \right){\tau }^{\alpha }{u_x^m|}_{x_{3i,l}}^{x_{3i,c}}=\Gamma \left(2-\alpha \right){\tau }^{\alpha }{\displaystyle {\int }_{I_i^{*}}{f}^m\text{d}x},i=1,\cdots ,n$，

${\displaystyle {\int }_{I_i^{**}}{u}^m\text{d}x}-{\displaystyle \underset{k=1}{\overset{m-1}{\sum }}{\beta }_k{\displaystyle {\int }_{I_i^{**}}{u}^{m-k}\text{d}x}}-{\beta }_{m-1}{\displaystyle {\int }_{I_i^{**}}{u}^0\text{d}x}-\Gamma \left(2-\alpha \right){\tau }^{\alpha }{u_x^m|}_{x_{3i,c}}^{x_{3i,r}}=\Gamma \left(2-\alpha \right){\tau }^{\alpha }{\displaystyle {\int }_{I_i^{*\text{*}}}{f}^m\text{d}x},i=1,\cdots ,n$，

${\displaystyle {\int }_{I_i^{***}}{u}^m\text{d}x}-{\displaystyle \underset{k=1}{\overset{m-1}{\sum }}{\beta }_k{\displaystyle {\int }_{I_i^{***}}{u}^{m-k}\text{d}x}}-{\beta }_{m-1}{\displaystyle {\int }_{I_i^{***}}{u}^0\text{d}x}-\Gamma \left(2-\alpha \right){\tau }^{\alpha }{u_x^m|}_{x_{3i,r}}^{x_{3\left(i+1\right),l}}=\Gamma \left(2-\alpha \right){\tau }^{\alpha }{\displaystyle {\int }_{I_i^{*\text{**}}}{f}^m\text{d}x},i=1,\cdots ,n-1$，

其中各个积分项的数值计算可参考文献 [9]。该格式具有与三次有限元方法有同阶精度但要比三次有限元方法更加便于实现。

3. 误差分析

首先给出误差分析所需一些基本引理和定义。

引理1 [9] 对任意的 $u_h\in U_h$，半范数 ${\left|\cdot \right|}_{\text{1,}h}$ 与 ${\left|\cdot \right|}_1$ 等价以及范数 ${\Vert \cdot \Vert }_{\text{0,}h}$ 与 ${\Vert \cdot \Vert }_0$ 等价，即

${\left|u_h\right|}_{\text{1,}h}\le {\left|u_h\right|}_1\le \frac{9\sqrt{10}}{20}{\left|u_h\right|}_{1,h},0.59{\Vert u_h\Vert }_{0,h}\le {\Vert u_h\Vert }_0\le 1.16{\Vert u_h\Vert }_{0,h}$。

引理2 [9] 对充分小的h， $a\left(u_h,{\pi }_h^{*}{u}_h\right)$ 是正定的，即存在正常数 $\sigma$，使得

$a\left(u_h,{\pi }_h^{*}{u}_h\right)\ge \sigma {\left|u_h\right|}_{\text{1}}^2,\forall u_h\in U_h$。

定义1 [9] 椭圆投影算子 $P_h:H^2\left(I\right)\cap H_0^1\left(I\right)\to U_h$ 使得

$a\left(P_{h}u,{\pi }_h^{*}v\right)=a\left(u,{\pi }_h^{*}v\right),\forall v\in U_h$。

引理3 [9] 设 $P_h$ 是由式(9)所定义的椭圆投影算子，则对 $\forall u\in H^5\left(I\right)\cap H_0^1\left(I\right)$，有

${\left|u-P_{h}u\right|}_{\text{1}}\le C{h}^3{\left|u\right|}_4,{\Vert u-P_{h}u\Vert }_0\le C{h}^4{\Vert u\Vert }_5$。

引理4 [9] 对 $\forall u_h,v_h\in U_h$，迁移算子 ${\pi }_h^{*}$ 满足下列不等式

${\Vert {\pi }_h^{\text{*}}{u}_h\Vert }_{\text{0}}\le 2.9214{\Vert u_h\Vert }_0,\left|\left(u_h,{\pi }_h^{*}{v}_h\right)-\left(u_h,v_h\right)\right|\le \frac{0.0133068}{\sigma }{h}^2{\Vert u\Vert }_{0,h}^2+\sigma {\left|v\right|}_{1,h}^2$。

引理5 [15] 设 ${\epsilon }^k\ge 0,k=0,1,\cdots ,N$ 且 $\gamma >0$，则当

${\epsilon }^m\le {\displaystyle \underset{k=1}{\overset{m-1}{\sum }}{\beta }_k{\epsilon }^{m-k}+\gamma }$，

时有 ${\epsilon }^m\le C{\tau }^{-\alpha }\gamma ,\left(m=1,2,\cdots ,N\right)$，其中C只依赖于 $T,\alpha$ 和u。

其次，利用上述引理和弱形式(4)和格式(5)得到如下误差分析结果。

定理1 设 $u\left(t_m\right)$ 和 $u_h^m$ 分别是问题(1)和格式(5)的解，并 ${\tau }^{\alpha }=O\left(h^2\right)$，$u_h^0=P_h\phi \left(x\right)$，则存在与剖分步长h和 $\tau$ 无关的正常数C，使得

${\Vert u\left(t_m\right)-u_h^m\Vert }_0\le C\left(h^4+{\tau }^{2-\alpha }\right)$。

证 令 $u^m-u_h^m=u^m-P_h{u}^m+P_h{u}^m-u_h^m={\rho }^m+{\theta }^m$，则由引理3知 ${\Vert {\rho }^m\Vert }_0\le C{h}^4{\Vert u^m\Vert }_5$。将方程(4)和(5)相减，得误差方程

$\begin{array}{l}\left({\theta }^m,{\pi }_h^{*}v\right)+\Gamma \left(2-\alpha \right){\tau }^{\alpha }a\left({\theta }^m,{\pi }_h^{*}v\right)\\ =\left({\displaystyle \underset{k=1}{\overset{m-1}{\sum }}{\beta }_k{\theta }^{m-k}}+{\beta }_{m-1}{\theta }^0,{\pi }_h^{*}v\right)+\Gamma \left(2-\alpha \right){\tau }^{\alpha }\left(R^m,{\pi }_h^{*}v\right)-\left({\rho }^m-{\displaystyle \underset{k=1}{\overset{m-1}{\sum }}{\beta }_k{\rho }^{m-k}}-{\beta }_{m-1}{\rho }^0,{\pi }_h^{*}v\right)\end{array}$ (6)

上式中取 $v={\theta }^m$，并令

$\left({\theta }^m,{\pi }_h^{*}{\theta }^m\right)=\left({\theta }^m,{\theta }^m\right)+\left({\theta }^m,{\pi }_h^{*}{\theta }^m\right)-\left({\theta }^m,{\theta }^m\right)=T_1+T_2$。 (7)

利用Cauchy-Schwarz不等式、引理1和引理4，有

$\left|T_2\right|=\left|\left({\theta }^m,{\pi }_h^{*}{\theta }^m\right)-\left({\theta }^m,{\theta }^m\right)\right|\le \frac{C{h}^2}{\Gamma \left(2-\alpha \right){\tau }^{\alpha }}{\Vert {\theta }^m\Vert }_0^2+\sigma \Gamma \left(2-\alpha \right){\tau }^{\alpha }{\left|{\theta }^m\right|}_1^2$， (8)

$\left|\left({\theta }^{m-k},{\pi }_h^{*}{\theta }^m\right)\right|\le C{\Vert {\theta }^{m-k}\Vert }_0{\Vert {\theta }^m\Vert }_0,a\left({\theta }^m,{\pi }_h^{*}{\theta }^m\right)\ge \sigma {\left|{\theta }^m\right|}_1^2$。 (9)

令 $S_{\rho }=\frac{1}{\Gamma \left(2-\alpha \right){\tau }^{\alpha }}\left({\rho }^m-{\displaystyle \underset{k=1}{\overset{m-1}{\sum }}{\beta }_k{\rho }^{m-k}}-{\beta }_{m-1}{\rho }^0\right)$，则经过展开计算和重新合并后有

$S_{\rho }=\frac{{\tau }^{1-\alpha }}{\Gamma \left(2-\alpha \right)}{\displaystyle \underset{k=0}{\overset{m-1}{\sum }}{d}_k\left(\frac{{\rho }^{m-k}-{\rho }^{m-k-1}}{\tau }\right)}$，

又由于

${\Vert \frac{{\rho }^{m-k}-{\rho }^{m-k-1}}{\tau }\Vert }_0={\Vert \frac{1}{\tau }{\displaystyle {\int }_{t_{m-k-1}}^{t_{m-k}}{\rho }_t\left(x,s\right)\text{d}s}\Vert }_0\le C\underset{t\in \left[0,T\right]}{\max }{\Vert u_t\Vert }_5$。

因此

${\Vert S_{\rho }\Vert }_{\text{0}}\le C{T}^{1-\alpha }{h}^4\underset{t\in \left[0,T\right]}{\max }{\Vert u_t\Vert }_5$。 (10)

于是由不等式(6)~(10)，并取 $u_h^0=P_h\phi \left(x\right)$，利用引理5得 ${\Vert {\theta }^m\Vert }_{\text{0}}\le C\left(h^4+{\tau }^{2-\alpha }\right)$。再与 ${\Vert {\rho }^m\Vert }_0\le C{h}^4{\Vert u^m\Vert }_5$ 结合可得结论，证毕。

4. 数值算例

本节用所提格式在区间 $\left(x,t\right)\in \left[0,1\right]\times \left[0,1\right]$ 上数值求解具有解析解 $u=t^2\sin \left(2\pi x\right)$ 的问题(1)~(2)。此时，问题(1)~(2)的初始值 $\phi =0$。空间和时间均匀网格步长分别取为 $h=1/3n,\tau =1/N$，其中 $n,N$ 分别是空间和时间剖分个数。时间和空间收敛阶计算公式分别为

$Rat{e}_t=\left|\frac{\ln \left(\Vert 细网格误差\Vert \right)/\ln \left(\Vert 粗网格误差\Vert \right)}{\ln \left(细网格剖分N\right)/\ln \left(粗网格剖分N\right)}\right|$，$Rat{e}_s=\left|\frac{\ln \left(\Vert 细网格误差\Vert \right)/\ln \left(\Vert 粗网格误差\Vert \right)}{\ln \left(细网格剖分n\right)/\ln \left(粗网格剖分n\right)}\right|$。

当分数阶导数参数 $\alpha =0.2$ 时，表1左半部分给出了取定充分小的时间步长 $\tau =0.5\text{e}-4$ 而空间步长 $h=1/20,1/25,1/30,1/35$ 时的空间方向收敛阶 $Rat{e}_s$，右半部分给出了取定充分小的空间步长 $h=0.5\text{e}-4$ 而时间步长 $h=1/75,1/80,1/90,1/200$ 时的时间方向收敛阶 $Rat{e}_t$。表2和表3中分别给出了分数阶导数参数 $\alpha =0.5,0.8$ 的数值计算结果，表格的结构与表1类似。表1~3中数据表明时间方向收敛阶接近 $2-\alpha$ 、空间方向收敛阶接近4阶与定理1的理论分析结果吻合。

表1. 当取 $\alpha =0.2$ 时，用格式计算算例的数值结果

表2. 当取 $\alpha =0.5$ 时，用格式计算算例的数值结果

表3. 当取 $\alpha =0.8$ 时，用格式计算算例的数值结果

5. 结论

针对一维时间分数扩散方程初边值问题，本文基于Lagrange插值多项式的应力佳点对偶剖分和L1-逼近公式构建了一种三次有限体积元格式，并详细推导了格式的L2-模 $O\left({\tau }^{2-\alpha },h^4\right)$ 阶误差估计，其中 $\tau$ 和h分别表示时间剖分步长和空间剖分步长。数值实验验证了该方法的有效性和理论分析结果，并且数值结果表明本文方法拥有较高的计算精度。

基金项目

2022年度自治区直属高校基本科研业务费项目，呼和浩特民族学院校级科学研究项目(HM-ZD-202101)资助。

文章引用

何斯日古楞,澈力木格,高晶英. 时间分数阶扩散方程的一类三次有限体积元方法
A Cubic Finite Volume Element Method for the Time Fractional Diffusion Equation[J]. 应用数学进展, 2022, 11(08): 5529-5535. https://doi.org/10.12677/AAM.2022.118582

参考文献