时间分数阶扩散方程的一类三次有限体积元方法
A Cubic Finite Volume Element Method for the Time Fractional Diffusion Equation
DOI: 10.12677/AAM.2022.118582, PDF, HTML, XML,    科研立项经费支持
作者: 何斯日古楞, 高晶英：呼和浩特民族学院，数学与大数据学院，内蒙古 呼和浩特；澈力木格：呼和浩特民族学院，计算机科学与信息工程学院，内蒙古 呼和浩特
关键词: 时间分数阶扩散方程；三次有限体积元法；L1-公式；Time Fractional Diffusion Equation； Cubic Finite Volume Element Method； L1-Formula
摘要: 本文基于应力佳点对偶网格剖分以及分片三次Lagrange插值试探函数空间和分片常数检验函数空间的三次有限体积元法和Caputo导数的L1-逼近公式构造数值格式求解一维时间分数阶扩散方程，并证明了格式的L2范数在时间和空间方向分别 阶和四阶收敛误差估计。通过数值实验验证了理论分析结果以及所提格式的有效性。
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 inter-polation and the test function space of piecewise constant. The L2-norm error estimate of fourth or-der convergence in space and order convergence in time is proved. Numerical experiments are given to verify the effectiveness of the theoretical analysis results and the proposed scheme.
文章引用：何斯日古楞, 澈力木格, 高晶英. 时间分数阶扩散方程的一类三次有限体积元方法[J]. 应用数学进展, 2022, 11(8): 5529-5535. https://doi.org/10.12677/AAM.2022.118582

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. 三次有限体积元格式

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

∂ α u ( x , t ) ∂ t α − ∂ 2 u ∂ x 2 = f ( x , t ) ,   ( x , t ) ∈ ( a , b ) × ( 0 , T ] ， (1)

以及初边值条件

{ u ( a , t ) = 0 ,   u ( b , t ) = 0 ,   t ∈ ( 0 , T ] , u ( x , 0 ) = φ ( x ) ,   x ∈ ( a , b ) , (2)

其中 0 < α < 1 ， f , φ 是给定的光滑函数，Caputo分数阶导数定义为

∂ α u ( x , t ) ∂ t α = 1 Γ ( 1 − α ) ∫ 0 t ∂ u ( x , s ) ∂ s d s ( t − s ) α 。 (3)

对正整数N，令 τ = T / N 是时间步长， u m = u ( x , t m )   ( t m = m τ , 0 ≤ m ≤ N ) 和 u ¯ m = ( u m + u m − 1 ) / 2 。分数阶导数(3)在点 t m 处有 [15]

∂ α u ( x , t m ) ∂ t α = τ − α Γ ( 2 − α ) [ u m − ∑ k = 1 m − 1 β k u m − k − β m − 1 u 0 ] + R m ，

其中 β k = ( d k − 1 − d k ) ,   β m − 1 = d m − 1 ,   d k = ( k + 1 ) 1 − α − k 1 − α 且 R m = o ( τ 2 − α ) 。于是方程(1)的时间半离散形式为

u m − ∑ k = 1 m − 1 β k u m − k − β m − 1 u 0 − Γ ( 2 − α ) τ α u x x m = Γ ( 2 − α ) τ α f m + Γ ( 2 − α ) τ α R m 。

对区间 I = [ a , b ] 做剖分 I h ，节点为 a = x 0 < x 1 < ⋯ < x n < x n + 1 < ⋯ < x 2 n < x 2 n + 1 < ⋯ < x 3 n = b ，称 I i = [ x 3 i − 3 , x 3 i ] 为单元，单元内的三等分点为 x 3 i − 2 , x 3 i − 1 ，令 h i 为每等分长度， h = max 1 ≤ i ≤ n h i 。在区间 I i 上，经过 x 3 i − 3 , x 3 i − 2 , x 3 i − 1 , x 3 i 四点的三次插值应力佳点分别为

x 3 i , l = x 3 i − 3 + 5 2 h i , x 3 i , r = x 3 i − 3 − 5 2 h i , x 3 i , c = x 3 i − 3 2 h i 。

在这三点处，函数的三次插值导数具有四阶精度，而其余点处仅有三阶精度 [9]。记 I i * = [ x 3 i , l , x 3 i , c ] ， I i * * = [ x 3 i , c , x 3 i , r ] ， I i * * * = [ x 3 i , r , x 3 ( i + 1 ) , l ] ， ( i = 1 , 2 , ⋯ , n − 1 ) ， I 0 * * * = [ x 0 , x 3 , l ] ， I n * * * = [ x 3 n , r , x 3 n ] 。所有的 I i * , I i * * , I i * * * 构成剖分 I h 的对偶剖分 I h * ，称 I i * , I i * * , I i * * * 分别为对应节点 x 3 i − 2 , x 3 i − 1 , x 3 i 的控制体积。

这里除了用标准Sobolev空间及L2范数 ‖   ⋅   ‖ 0 和H1半范数 |   ⋅   | 1 ，需要如下离散空间

U h = { u h ∈ H 0 1 ( I ) : u h | I i ∈ P 3 ( I i ) ,   ∀ I i ∈ I h ,   u h ( a , t ) = u h ( b , t ) = 0 } ，

V h = { v h ∈ L 2 ( I ) :   v h | I i * ( 或 I i * * 或 I i * * * ) = 常 数 ,   ∀ I i * ( 或 I i * * 或 I i * * * ) ∈ I h * ,   v h | I 0 * * * = 0 , v h | I n * * * = 0 } ，

和如下离散范数

‖ u h ‖ 0 , h 2 = 3 8 ∑ i = 1 n h i ( u 3 i − 3 2 + 3 u 3 i − 2 2 + 3 u 3 i − 1 2 + u 3 i 2 )   ,   ∀ u h ∈ U h ，

| u h | 1 , h 2 = ∑ i = 1 n 1 h i [ ( u 3 i − 2 − u 3 i − 3 ) 2 + ( u 3 i − 1 − u 3 i − 2 ) 2 + ( u 3 i − u 3 i − 1 ) 2 ]   ,   ∀ u h ∈ U h ，

其中令 u i = u h ( x i ) 。定义插值算子 π h * : U h → V h 使得

π h * u h = ∑ i = 1 n ( u h ( x 3 i − 2 ) ψ 3 i − 2 + u h ( x 3 i − 1 ) ψ 3 i − 1 + u h ( x 3 i ) ψ 3 i ) ，

其中 ψ 3 i − 2 , ψ 3 i − 1 , ψ 3 i 分别是控制体积 I i * , I i * * , I i * * * 上的特征函数。

问题(1)~(2)的有限体积元弱形式为：对 ∀ v ∈ V h ，求 u ∈ H 0 1 ( I ) 使得

( u m , v ) + Γ ( 2 − α ) τ α a ( u m , v ) = ∑ k = 1 m − 1 β k ( u m − k , v ) + β m − 1 ( u 0 , v ) + Γ ( 2 − α ) τ α ( f m , v ) + Γ ( 2 − α ) τ α ( R m , v ) 。

由于 π h * 是 U h 到 V h 的迁移算子，上述弱形式等价于：对 ∀ v ∈ U h ，求 u ∈ H 0 1 ( I ) 使得

( u m , π h * v ) + Γ ( 2 − α ) τ α a ( u m , π h * v ) = ∑ k = 1 m − 1 β k ( u m − k , π h * v ) + β m − 1 ( u 0 , π h * v ) + Γ ( 2 − α ) τ α ( f m , π h * v ) + Γ ( 2 − α ) τ α ( R m , π h * v ) 。 (4)

相应的有限体积元格式为：对 ∀ v ∈ U h ，求 u h ∈ U h 使得

( u h m , π h * v ) + Γ ( 2 − α ) τ α a ( u h m , π h * v ) = ∑ k = 1 m − 1 β k ( u h m − k , π h * v ) + β m − 1 ( u h 0 , π h * v ) + Γ ( 2 − α ) τ α ( f m , π h * v ) ， (5)

其中

a ( u h , π h * v ) = ∑ i = 1 n [ v 3 i − 2 a ( u h , ψ 3 i − 2 ) + v 3 i − 1 a ( u h , ψ 3 i − 1 ) ] + ∑ i = 1 n − 1 v 3 i a ( u h , ψ 3 i ) ，

( u h , π h * v ) = ∑ i = 1 n [ v 3 i − 2 ( u h , ψ 3 i − 2 ) + v 3 i − 1 ( u h , ψ 3 i − 1 ) ] + ∑ i = 1 n − 1 v 3 i ( u h , ψ 3 i ) ，

且

a ( u h , ψ 3 i − 2 ) = u ′ h ( x 3 i , l ) − u ′ h ( x 3 i , c ) ,   a ( u h , ψ 3 i − 1 ) = u ′ h ( x 3 i , c ) − u ′ h ( x 3 i , r ) ,   a ( u h , ψ 3 i ) = u ′ h ( x 3 i , r ) − u ′ h ( x 3 ( i + 1 ) , l )

( u h , ψ 3 i − 2 ) = ∫ I i * u h   d x ,   ( u h , ψ 3 i − 1 ) = ∫ I i * * u h   d x ,   ( u h , ψ 3 i ) = ∫ I i * * * u h   d x 。

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

∫ I i * u m d x − ∑ k = 1 m − 1 β k ∫ I i * u m − k d x − β m − 1 ∫ I i * u 0 d x − Γ ( 2 − α ) τ α u x m | x 3 i , l x 3 i , c = Γ ( 2 − α ) τ α ∫ I i * f m d x ,   i = 1 , ⋯ , n ，

∫ I i * * u m d x − ∑ k = 1 m − 1 β k ∫ I i * * u m − k d x − β m − 1 ∫ I i * * u 0 d x − Γ ( 2 − α ) τ α u x m | x 3 i , c x 3 i , r = Γ ( 2 − α ) τ α ∫ I i * * f m d x ,   i = 1 , ⋯ , n ，

∫ I i * * * u m d x − ∑ k = 1 m − 1 β k ∫ I i * * * u m − k d x − β m − 1 ∫ I i * * * u 0 d x − Γ ( 2 − α ) τ α u x m | x 3 i , r x 3 ( i + 1 ) , l = Γ ( 2 − α ) τ α ∫ I i * ** f m d x ,   i = 1 , ⋯ , n − 1 ，

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

3. 误差分析

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

引理1 [9] 对任意的 u h ∈ U h ，半范数 |   ⋅   | 1, h 与 |   ⋅   | 1 等价以及范数 ‖   ⋅   ‖ 0, h 与 ‖   ⋅   ‖ 0 等价，即

| u h | 1, h ≤ | u h | 1 ≤ 9 10 20 | u h | 1 , h ,   0.59 ‖ u h ‖ 0 , h ≤ ‖ u h ‖ 0 ≤ 1.16 ‖ u h ‖ 0 , h 。

引理2 [9] 对充分小的h， a ( u h , π h * u h ) 是正定的，即存在正常数 σ ，使得

a ( u h , π h * u h ) ≥ σ | u h | 1 2 ,   ∀ u h ∈ U h 。

定义1 [9] 椭圆投影算子 P h : H 2 ( I ) ∩ H 0 1 ( I ) → U h 使得

a ( P h u , π h * v ) = a ( u , π h * v ) ,   ∀ v ∈ U h 。

引理3 [9] 设 P h 是由式(9)所定义的椭圆投影算子，则对 ∀ u ∈ H 5 ( I ) ∩ H 0 1 ( I ) ，有

| u − P h u | 1 ≤ C h 3 | u | 4 ,   ‖ u − P h u ‖ 0 ≤ C h 4 ‖ u ‖ 5 。

引理4 [9] 对 ∀ u h ,   v h ∈ U h ，迁移算子 π h * 满足下列不等式

‖ π h * u h ‖ 0 ≤ 2.9214 ‖ u h ‖ 0 ,   | ( u h , π h * v h ) − ( u h , v h ) | ≤ 0.0133068 σ h 2 ‖ u ‖ 0 , h 2 + σ | v | 1 , h 2 。

引理5 [15] 设 ε k ≥ 0 ,   k = 0 , 1 , ⋯ , N 且 γ > 0 ，则当

ε m ≤ ∑ k = 1 m − 1 β k ε m − k + γ ，

时有 ε m ≤ C τ − α γ , ( m = 1 , 2 , ⋯ , N ) ，其中C只依赖于 T , α 和u。

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

定理1 设 u ( t m ) 和 u h m 分别是问题(1)和格式(5)的解，并 τ α = O ( h 2 ) ， u h 0 = P h φ ( x ) ，则存在与剖分步长h和 τ 无关的正常数C，使得

‖ u ( t m ) − u h m ‖ 0 ≤ C ( h 4 + τ 2 − α ) 。

证 令 u m − u h m = u m − P h u m + P h u m − u h m = ρ m + θ m ，则由引理3知 ‖ ρ m ‖ 0 ≤ C h 4 ‖ u m ‖ 5 。将方程(4)和(5)相减，得误差方程

( θ m , π h * v ) + Γ ( 2 − α ) τ α a ( θ m , π h * v ) = ( ∑ k = 1 m − 1 β k θ m − k + β m − 1 θ 0 , π h * v ) + Γ ( 2 − α ) τ α ( R m , π h * v ) − ( ρ m − ∑ k = 1 m − 1 β k ρ m − k − β m − 1 ρ 0 , π h * v ) (6)

上式中取 v = θ m ，并令

( θ m , π h * θ m ) = ( θ m , θ m ) + ( θ m , π h * θ m ) − ( θ m , θ m ) = T 1 + T 2 。 (7)

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

| T 2 | = | ( θ m , π h * θ m ) − ( θ m , θ m ) | ≤ C h 2 Γ ( 2 − α ) τ α ‖ θ m ‖ 0 2 + σ   Γ ( 2 − α ) τ α | θ m | 1 2 ， (8)

| ( θ m − k , π h * θ m ) | ≤ C ‖ θ m − k ‖ 0 ‖ θ m ‖ 0 ,   a ( θ m , π h * θ m ) ≥ σ | θ m | 1 2 。 (9)

令 S ρ = 1 Γ ( 2 − α ) τ α ( ρ m − ∑ k = 1 m − 1 β k ρ m − k − β m − 1 ρ 0 ) ，则经过展开计算和重新合并后有

S ρ = τ 1 − α Γ ( 2 − α ) ∑ k = 0 m − 1 d k ( ρ m − k − ρ m − k − 1 τ ) ，

又由于

‖ ρ m − k − ρ m − k − 1 τ ‖ 0 = ‖ 1 τ ∫ t m − k − 1 t m − k ρ t ( x , s ) d s ‖ 0 ≤ C max t ∈ [ 0 , T ] ‖ u t ‖ 5 。

因此

‖ S ρ ‖ 0 ≤ C T 1 − α h 4 max t ∈ [ 0 , T ] ‖ u t ‖ 5 。 (10)

于是由不等式(6)~(10)，并取 u h 0 = P h φ ( x ) ，利用引理5得 ‖ θ m ‖ 0 ≤ C ( h 4 + τ 2 − α ) 。再与 ‖ ρ m ‖ 0 ≤ C h 4 ‖ u m ‖ 5 结合可得结论，证毕。

4. 数值算例

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

R a t e t = | ln ( ‖ 细 网 格 误 差 ‖ ) / ln ( ‖ 粗 网 格 误 差 ‖ ) ln ( 细 网 格 剖 分 N ) / ln ( 粗 网 格 剖 分 N ) | ， R a t e s = | ln ( ‖ 细 网 格 误 差 ‖ ) / ln ( ‖ 粗 网 格 误 差 ‖ ) ln ( 细 网 格 剖 分 n ) / ln ( 粗 网 格 剖 分 n ) | 。

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

Table 1. Numerical results with α = 0.2 for the example that calculated by using the scheme

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

Table 2. Numerical results with α = 0.5 for the example that calculated by using the scheme

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

Table 3. Numerical results with α = 0.8 for the example that calculated by using the scheme

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

5. 结论

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

基金项目

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

参考文献

[1] 李荣华, 陈仲英. 微分方程的广义差分法[M]. 长春: 吉林大学出版社, 1994.
[2] 李荣华, 祝丕奇. 二阶椭圆偏微分方程的广义差分方法(I)-三角网情形[J]. 高等学校计算数学学报, 1982, 4(2): 140-152.
[3] Tian, M. and Chen, Z. (1991) Quadratical Element Generalized Differential Methods for Elliptic Equations. Numerical Mathematics: A Journal of Chinese Universities, 13, 99-113.
[4] Liebau, F. (1996) The Finite Volume Element Method with Quadratic Basis Functions. Computing, 57, 281-299.
https://doi.org/10.1007/BF02252250
[5] Chen, Z., Wu, J. and Xu, Y. (2012) Higher-Order Finite Volume Meth-ods for Elliptic Boundary Value Problems. Advances in Computational Mathematics, 37, 191-253.
https://doi.org/10.1007/s10444-011-9201-8
[6] Wang X. and Li, Y.H. (2016) L2 Error Estimates for High Order Finite Volume Methods on Triangular Meshes. SIAM Journal on Numerical Analysis, 54, 2729-2749.
https://doi.org/10.1137/140988486
[7] Zou, Q. (2017) An Unconditionally Stable Quadritic Finite Volume Scheme Over Triangular Meshes for Elliptic Equations. Journal of Scientific Computing, 70, 112-124.
https://doi.org/10.1007/s10915-016-0244-3
[8] Wang, T.K. and Gu, Y. (2010) Superconvergent Biquadratic Fi-nite Volume Element Method for Two-Dimensional Poisson's Equations. Journal of Computational and Applied Math-ematics, 234, 447-460.
https://doi.org/10.1016/j.cam.2009.12.036
[9] Gao, G.H. and Wang, T.K. (2010) Cubic Superconvergent Finite Volume Element Method for One-dimensional Elliptic and Parabolic Equations. Journal of Computational and Applied Mathematics, 233, 2285-2301.
https://doi.org/10.1016/j.cam.2009.10.013
[10] 王星, 高广花, 王同科. 半线性抛物问题的一类三次有限体积元方法[J]. 辽宁工程技术大学学报(自然科学版), 2015, 34(2): 281-284.
[11] 于长华, 王晓玲, 李永海. 解两点边值问题的一类修改的三次有限体积元方法[J]. 计算数学, 2010, 32(4): 385-398.
[12] Yu, C.H. and Li, Y.H. (2011) Biquadratic Finite Volume Element Methods Based on Optimal Stress Points for Parabolic Problem. Journal of Compu-tational and Applied Mathematics, 236, 1055-1068.
https://doi.org/10.1016/j.cam.2011.07.030
[13] 杨凯丽, 何斯日古楞, 肖宇宇. 对流扩散方程的三次有限体积元法[J]. 内蒙古大学学报(自然科学版), 2021, 52(3): 250-256.
[14] 肖宇宇, 何斯日古楞, 杨凯丽. 对流扩散方程的时间间断时空有限体积元法[J]. 高校应用数学学报, 2021, 36(2): 179-192.
[15] Jiang, Y. J. and Ma, J.T. (2011) High-Order Finite Element Methods for Time-Fractional Partial Differential Equations. Journal of Computational Applied Mathematics, 235, 3285-3290.
https://doi.org/10.1016/j.cam.2011.01.011

为你推荐