非结构网格下2D Riesz分数阶方程的Galerkin有限元方法卜玮平【摘要】讨论了2D Riesz分数阶扩散方程的Galerkin有限元方法. 基于非结构网格,采用Lagrange线性分片多项式作为基函数,详细描述了分数阶扩散方程的有限元实现. 与现有方法相比, 该方法有效地降低了计算成本, 提高了刚度矩阵的精度.最后,数值算例验证了所提方法的有效性.【期刊名称】《纯粹数学与应用数学》【年(卷),期】2019(035)002【总页数】13页(P169-181)【关键词】Riesz分数阶扩散方程;有限元方法;非结构网格【作者】卜玮平【作者单位】湘潭大学数学与计算科学学院,湖南湘潭 411105【正文语种】中文【中图分类】O1781 引言最近,分数阶微分方程数值方法的研究越来越受关注.有限差分方法(FDM)是求解分数阶微分方程最常见的数值方法之一.文献[1]是研究FDM求解分数阶微分方程数值方法的先驱.随后,利用有限差分方法人们对各种分数阶微分方程进行了求解[2-5].谱方法(SM)和谱元方法(SEM)是求解分数阶微分方程的重要方法,它们通常具有很高的收敛精度.目前,SM和SEM求解分数阶微分方程的工作有文献[6-8]等.相比于FDM方法和SM(或SEM)方法,有限元方法(FEM)的主要特征是它能够较易处理复杂区域且对解的光滑性要求较低.文献[9]首次尝试用有限元方法求解分数阶微分方程,并给出了有限元逼近的理论框架.此后,文献[10]建立了有限元求解一维和二维分数阶对流弥散方程的理论框架.随后,在有限元求解分数阶微分方程方面涌现了越来越多的工作,其中包括文献[11-15]等.近年来,关于二维分数阶微分方程的有限元方法也有一些研究成果.文献[10]考虑了基于分数阶方向导数的二维分数阶对流-弥散方程的有限元方法.通过使用矩阵转换技术,文献[16]讨论了基于分数阶拉普拉斯算子的二维分数阶扩散方程有限元方法.文献[17]发展了基于分数拉普拉斯算子的分数阶扩散方程自适应有限元方法.为了求解二维Risez/Riemann-Liouville分数阶扩散方程(2DRFDE/2DRLFDE),基于一致三角网格剖分,文献[18-19]考虑了Galerkin有限元方法.然而,对于不规则区域上2DRFDE/2DRLFDE的有限元方法,通常需要采用非结构网格.基于非结构网格,文献[20]利用三角形上定义的Lagrange多项式,建立了2DRLFDE的间断Galerkin方法.文献[21]考虑了非线性2DRFDE的有限元方法,并描述了有限元方法的实现.文献[22]讨论了二维时空分数阶波方程在非规则凸域上的有限元方法.尽管在文献[20-22]中已有非结构网格下有限元求解分数阶微分方程的工作,但是这些工作使用的是相同的有限元实现技巧.值得注意,上述文中提到的有限元实现方法有不足之处.因此,在这篇文章中将提出一些新的技巧,以改进现有的有限元实现方法.本文的主要贡献如下:首先,对于刚度矩阵的计算降低了计算花销,原来的计算花销为O(Ne3),现在的花销为O(Ne2),Ne为总剖分单元数;其次,提高了刚度矩阵的元素的精度,对于三角形单元的高斯积分,不同于现有的计算方法,高斯积分区域为被积函数的非零的区域,这比由以往方法得到的刚度矩阵元素更精确;第三,将Riemann-Liouville导数转化为Caputo导数,简化了内积的计算.本文的结构安排如下:在第2节,给出了分数阶导数的定义、模型问题及有限元全离散格式的推导;在第3节,首先介绍了现有的刚度矩阵计算方法,然后详细描述了有限元的实现方法,并与现有方法进行了比较;在第4节,给出了数值实验来证明方法的有效性;最后,对本文进行了总结.2 准备工作令Ω⊂R2,则x,y方向的左右Riemann-Liouville分数阶导数定义如下:其中γ,n−1<γ≤n,n∈N,a(y),b(y),c(x),d(x)定义如图1.图1 关于Ω上a(y),b(y),c(x),d(x)的定义进一步,定义x,y方向的γ(γ1,3,···)阶 Riesz分数阶导数如下类似地,x,y方向γ阶右Caputo分数阶导数定义如下在文献[23]中已经讨论了Riemann-Liouville分数阶导数和Caputo分数阶导数的等价关系,如果u(a(y),y)=0,u(x,c(x))=0,0<γ<1,则本文考虑如下分数阶扩散方程其中,β<1,P,Q,A为常数.对上述方程,首先考虑其变分形式和有限元全离散格式.对于任意γ≥0,记γ(Ω)为Hγ(Ω)的子集且其元素在Ω 外的零扩张属于Hγ(R2).令V=α(Ω)∩ β(Ω).于是,由文献[18]中引理5可得如下变分问题:寻找u∈V使得其中(·,·):=(·,·)L2(Ω),F(v):=(f,v),且为了得到上述变分问题的全离散格式,先将Ω进行剖分.令{Th}是Ω的一个正则的三角剖分,h为所有三角形单元的最大直径.定义如下有限元空间这里PK(x,y)为K次多项式,接下来,定义(10)式的全离散格式:寻找uh∈XhK使得3 有限元方法的实现为了计算全离散格式(14),需要计算刚度矩阵和荷载向量.由于分数阶导数为非局部算子,与整数阶微分方程相比,其主要区别在于分数阶微分方程刚度矩阵的计算非常复杂.因此,方法实现的重点将放在刚度矩阵的计算上.令 ,则uh可写成,其中ϕi(x,y)是属于剖分节点i的Lagrange线性基函数,N为剖分节点总数目.如果三角形单元e包含顶点i且它的三个顶点逆时针排列依次为i,j,k,则有这里∆e为三角形单元e的面积,(xi,yi),(xj,yj),(xk,yk)是对应顶点i,j,k的坐标.接下来,考虑B(uh,vh)的计算.令这里S={Sij}N×N为刚度矩阵,其元素为3.1 已有的实现方法目前,文献[20-22]已经考虑了Sij的计算.然而,应该注意这些文章关于Sij的计算方法是类似的.这里简单对文献[20-22]中Sij的计算方法进行描述.由于(17)中四个内积的计算具有相似性,因此仅仅以为例进行讨论.考虑分数阶算子的非局部性并利用高斯积分可得这里Ge为三角形单元e上的高斯点,ωl是高斯点(xl,yl)对应的权重.上面(18)式中的方法可以用来计算Sij,然而该方法有三点不足之处.令Ne为剖分三角形单元的总数.首先,对于Sij的计算需要考虑如下积分在每个三角形单元的高斯积分因此,为了得到刚度矩阵的一个元素Sij,上述积分需要进行Ne次,而为了获得整个刚度矩阵,需要计算次,这将使得刚度矩阵的计算花销随Ne的增涨而快速增涨.其次,对于下列积分考虑被积函数在三角形单元上的非零区域.由于分数阶算子具有非局部性,显然存在满足下列情形的三角形单元:被积函数在三角形的某些部分为零,在其他部分非零.因此,被积函数在满足上述条件的三角形单元上为间断函数.如果对这些单元,在整个三角形上运用高斯积分势必达不到数值积分的相应精度.图2给了描述上述情形的例子(图形(a)描述了的非零区域,其中ϕi为节点i的基函数,ϕk为节点k的基函数;图形(b)描述了被积函数的非零区域,它为支集的交集).图2 的非零区域和的非零区域最后,对于被积函数 ,将说明它在某些区域光滑性差的特点.假设ϕi(x,y)和点(xp,yp)由图3给出在(xp,yp)的计算涉及的三角形的边界点p0,p1,p2,它们为积分路径与三角形单元e1,e2,e3边界的交点,它们在x轴方向的坐标分别设为x0,x1,x2).图3 在(xp,yp)的计算结果计算ϕi(x,y)的左Riemann-Liouville分数阶导数在(xp,yp)点的值从公式(21)可以看出, 在支集ϕi(x,y)的某些点处光滑性较差.类似地,也可以证明存在光滑性较差的区域.因此,为了保持(19)式的精度,计算时需要取大量的高斯点. 3.2 实现方法针对已有方法的不足,设计一种新的求解二维Riesz分数阶扩散方程的有限元实现方法.因为,显然因此,为了计算(16)中刚度矩阵S,根据对称性仅需要计算Sx和Sy.下面以为例来描述计算方法.假设ϕi(x,y),ϕj(x,y)的支集定义如图4(Lagrange线性基函数ϕi(x,y),ϕj(x,y)的支集如图4,其中为ϕi(x,y)支集上的三角形单元).图4 Lagrange线性基函数ϕi(x,y),ϕj(x,y)的支集这里是ϕi(x,y)关于变量x的系数.比较 (26)式与 (18)式,显然的计算花销将要下降.为了得到刚度矩阵,此时仅仅需要计算次积分.为了计算积分Ii,已有的方法是在整个三角形单元ei上使用高斯积分.然而,上面已经提到这样做将降低高斯积分的效果,因为分数阶算子具有非局部性, 的非零区域在某些单元可能只占有一部分(即在该单元为间断函数),如图5.图5 的支集的支集包含三角形单元因此,I1,I2能够在上选择高斯点来计算.然而, 仅仅含有支集的一部分.因此,可得其中为在上的支集.由于是一个三角形区域,因此在其上可以选用高斯点来计算I3.由于是四边形区域,因此要计算I4,I5,一种方法是在四边形上选用高斯点,另一种方法是将四边形分解成两个三角形然后分别作高斯积分.考虑图6所描述的情形,即是一个多边形区域.因此对I5的计算,一种方法是将该区域剖分成三角形与四边形然后进行高斯积分,另一种方法可以将积分区域全部剖分成三角形在每个三角形上进行高斯积分.图6 的支集注意到,对于ϕi(x,y), 由 (9)式可得于是利用高斯积分有这里e为的非零区域,ai,∆e定义在(15)式中.假设ϕj(x,y)及其积分路径被定义在图3中,对于高斯点(xp,yp),有显然(30)式比(21)式更容易计算.4 数值实验本节给出一个数值例子来验证方法的有效性.例 4.1在模型方程(10)中,取P=Q=2,A=4.考虑如下两种情形:(a)假设考虑问题区域为[0,1]×[0,1],其精确解为u=100x2(1−x)2y2(1−y)2;(b)假设考虑问题区域为,其精确解为,相应的右端函数分别定义如下当选择不同的α和β时,表1-表4分别列出了情形(a)和情形(b)所得的数值结果. 可以看出,所得误差的收敛率是最优收敛率.表1 基于Lagrange线性多项式与α=0.6,β=0.6计算情形(a)的数值误差与收敛率h‖uh −u (x,y) ‖ 0 收敛率1 4 2.193 4e-2 –1 8 4.730 6e-3 2.213 1 1 16 1.106 0e-3 2.096 7 1 32 2.489 8e-4 2.151 2表2 基于Lagrange线性多项式与α=0.7,β=0.8计算情形(a)的数值误差与收敛率h‖uh −u (x,y) ‖ 0 收敛率1 4 2.259 1 e− 2 –1 8 4.998 3 e− 3 2.176 2 1 16 1.173 6 e− 3 2.090 5 1 32 2.906 8 e− 4 2.013 4表3 基于Lagrange线性多项式与α=0.7,β=0.6计算情形(b)的数值误差与收敛率h‖uh −u (x,y) ‖ 0 收敛率1 4 4.567 0 e− 3 –1 8 1.106 5 e− 3 2.045 2 1 16 2.374 1 e− 4 2.220 6 1 32 5.610 1 e− 5 2.081 3表4 基于Lagrange线性多项式与α=0.6,β=0.8计算情形(b)的数值误差与收敛率h‖uh −u (x,y) ‖ 0 收敛率1 4 4.712 8 e− 3 –1 8 1.133 7 e− 3 2.055 5 1 16 2.467 7 e− 4 2.199 8 1 32 5.767 4 e− 5 2.097 25 总结本文研究了非结构网格下利用Lagrange线性基函数求解2D Riesz分数阶扩散方程的有限元方法实现.首先,描述了现有有限元全离散格式的实现方法,并指出了现有方法的不足之处.随后,针对这些缺点设计了一种新的实现方法,提高了有限元方法的计算效率和刚度矩阵的精度.最后,给出了数值算例,数值结果验证了本文所提方法的有效性.参考文献【相关文献】[1]Lubich C.Discretized fractional calculus[J].SIAM J.Math.Ana.,1986,17(3):704-719.[2]Liu F,Zhuang P,Anh V,et al.Stability and convergence of the difference methods for the space-time fractional advection-diffusion equation[J]pu.,2007,191(1):12-20.[3]Sun Z,Wu X.A fully discrete difference scheme for a diffusion-wave system[J].Applied Numerical Mathematics,2006,56(2):193-209.[4]Ding Z,Xiao A,Li M.Weighted finite difference methods for a class of space fractional partial differential equations with variablecoefficients[J].Appl.Math.,2010,233(8):1905-1914.[5]Tian W Y,Zhou H,Deng W.A class of second order difference approximations for solving space fractional diffusion equations[J].Mathematics of Computation,2015,84(294):1703-1727.[6]Zayernouri M,Karniadakis G E.Fractional Sturm-Liouville eigen-problems:theory and numerical approximation[J].Journal of Computational Physics,2013,252:495-517.[7]Li X,Xu C.Existence and uniqueness of the weak solution of the space-time fractional diffusion equation and a spectral method approximation[J].Phy.,2010,8(5):1016-1021.[8]Zeng F,Liu F,Li C,et al.A Crank–Nicolson ADI spectral method for a two-dimensional Riesz space fractional nonlinear reaction-diffusion equation[J].SIAMJ.Num.Ana.,2014,52(6):2599-2622.[9]Fix G J,Roof J P.Least squares finite-element solution of a fractional order two-point boundary value problem[J].Computers&Mathematics with Applications,2004,48(7-8):1017-1033.[10]Roop J P.Variational solution of the fractional advection dispersion equation[D].South Carolina:Clemson University,2004.[11]Wang H,Yang D.Wellposedness of variable-coefficient conservative fractional elliptic differential equations[J].SIAM Journal on Numerical Analysis,2013,51(2):1088-1107. [12]Jin B,Lazarov R,Pasciak J,et al.Error analysis of a finite element method for the space-fractional parabolic equation[J].SIAM Journal on Numerical Analysis,2014,52(5):2272-2294.[13]Deng W.Finite element method for the space and time fractional Fokker-Planck equation[J].SIAM journal on numerical analysis,2008,47(1):204-226.[14]Li C,Zhao Z,Chen Y Q.Numerical approximation of nonlinear fractional differential equations with subdiffusion and superdiffusion[J].Compu.Math.Appl.,2011,62(3):855-875.[15]Ma J,Liu J,Zhou Z.Convergence analysis of moving finite element methods for space fractional differential equations[J].Journal of Computational and Applied Mathematics,2014,255:661-670.[16]Yang Q,Turner I,Liu F,et al.Novel numerical methods for solving the time-space fractional diffusion equation in two dimensions[J].SIAM .,2011,33(3):1159-1180.[17]Chen L,Nochetto R H,Otrola E,et al.A PDE approach to fractional diffusion:a posteriori error analysis[J].Journal of Computational Physics,2015,293:339-358.[18]Bu W,Tang Y,Yang J.Galerkin finite element method for two-dimensional Riesz space fractional diffusion equations[J].Journal of Computational Physics,2014,276:26-38. [19]Zhao Y,Bu W,Huang J,et al.Finite element method for two-dimensional space-fractional advection-dispersion equations[J].Applied Mathematics andComputation,2015,257:553-565.[20]Qiu L,Deng W,Hesthaven J S.Nodal discontinuous Galerkin methods for fractional diffusion equations on 2D domain with triangular meshes[J].Phy.,2015,298:678-694.[21]Yang Z,Yuan Z,Nie Y,et al.Finite element method for nonlinear Riesz space fractionaldiffusion equations on irregular domains[J].Journal of ComputationalPhysics,2017,330:863-883.[22]Fan W,Liu F,Jiang X,et al.A novel unstructured mesh finite element method for solving the time-space fractional wave equation on a two-dimensional irregular convex domain[J].Fractional Calculus and Applied Analysis,2017,20(2):352-383.[23]Bu W,Tang Y,Wu Y,et al.Finite difference/ finite element method for two-dimensional space and time fractional Bloch-Torrey equations[J].Journal of Computational Physics,2015,293:264-279.。
![[硕士]题解无网格法及其与有限元法的比较研究初步_pdf](https://img.taocdn.com/s1/m/8b58c6db51e79b8968022675.png)
