含两个时间分数阶导数的二维反常扩散方程求解与微分阶数反演

王凤丹，孙春龙，李功胜，贾现正

(山东理工大学 理学院，山东 淄博 255049)

含两个时间分数阶导数的二维反常扩散方程求解与微分阶数反演

王凤丹，孙春龙，李功胜，贾现正

(山东理工大学 理学院，山东 淄博 255049)

对于一类含两个时间分数阶导数的二维反常扩散方程,基于对时间分数阶导数在Caputo意义下的离散,得到一个有限差分格式;利用分离变量法与Laplace变换得到该问题的解析解,并将两种方法得到的解进行数值比较.进一步,给定终值时刻数据,应用同伦正则化算法对扩散方程中的两个时间微分阶数进行数值反演,并给出反演算例.数值结果表明随着数据扰动水平的降低,解误差逐步变小,所用的反演算法对微分阶数反问题是有效的.

时间分数阶导数；二维分数阶扩散；解析解；微分阶数反问题；数值反演

近年来,反常扩散模型的应用越来越广泛.分数阶扩散方程在大气污染、水文地质学、金融学、物理学、力学、生物医学工程等诸多领域有了若干非常成功的应用[1],分数阶扩散模型相比经典的整数阶高斯扩散模型能更好的模拟扩散过程和重建试验数据.分数阶导数具有记忆性、遗传性和整体性,特别对于污染物长距离传输时表现出的非对称、非线性的整体性行为模式的研究,分数阶扩散模型可能是一个很有效的研究工具[2].对于含多个时间分数阶扩散方程正问题的研究,刘发旺和Meerschaert[3]等人,研究了含多个时间分数阶扩散-波动方程的数值计算方法,给出了隐式差分格式,并证明了隐式差分近似的无条件稳定性和收敛性. Luchko[4]应用分离变量法与傅里叶变换研究了一类含多个时间分数阶扩散方程的解析解,并证明了解的唯一性.Gejji[5]等利用分离变量法推导出了一维非齐次和二维齐次时间分数阶扩散-波动方程在齐次混合边值条件下的解析解.

对于分数阶扩散方程反问题的研究,越来越受到人们的关注. 郑光辉和魏婷[6-7]等应用傅里叶正则化,研究了带形区域上空间分数阶扩散方程的逆时问题,同时考虑了空间分数阶扩散方程中的源项反演问题,并且应用Tikhonov正则化对反问题数值求解.李功胜等[8]利用最佳摄动量算法对于一维时间分数阶扩散方程的微分阶数与扩散系数进行了数值联合反演. 刘继军等[9]应用准可逆性正则化方法探讨了一个时间分数阶扩散中由终值数据确定初始分布的逆时反问题.

目前,关于高维分数阶扩散方程相关反问题的研究还不多见.本文主要探讨含两个时间分数阶导数的二维时间分数阶齐次扩散方程的微分阶数反演问题.对于正问题的求解,通过对时间分数阶导数的数值离散,建立了一个隐式差分格式；同时利用分离变量法与Laplace变换得到该问题的解析解. 在正问题求解的基础上,给定终值时刻数据为附加数据，考虑一个确定扩散方程中两个时间微分阶数的反问题.应用同伦正则化算法进行数值反演,并讨论初始迭代、扰动水平等因素对反演算法的影响.

1 正问题及其数值解

考虑矩形域Ω=(0,l1)×(0,l2)上的二维时间分数阶齐次扩散方程的初边值问题：

(1)

其中(x,y)∈Ω,t>0,D>0是扩散系数,α,β∈(0,1)是时间分数阶导数的阶数,方程中的时间分数阶导数均用Caputo的定义，分别为：

下面给出数值求解二维问题(1)的一个隐式差分格式.

(2)

(3)

利用通常的差分离散方法,方程中的整数阶导数项可以进行以下近似：

则原方程可以离散为

(4)

(5)

如果令ck=dk-1-dk,k=1,2,…,n，则差分格式可以进一步整理为矩阵形式为

(6)

其中

(7)

i=1,2,…,M-1; I为(M-1)×(M-1)阶单位矩阵.

2 正问题的解析解

本节应用分离变量法和拉普拉斯变换,给出正问题(1)的基于Mittag-Leffler函数的解析解表达式,并进行数值模拟.

设u(x,y,t)=X(x)Y(y)T(t),代入方程可得

(8)

对式(8)整理可得

(9)

(10)

(11)

(12)

(13)

对t作Laplace变换可得

(14)

将式(14)代入式(13)中有：

(15)

(16)

令c=D(λn+μm)(常数).

(17)

由引理1，即有

于是，可得

(18)

3 正问题的数值模拟

取l1=l2=π,T=1,扩散系数D=0.5,初始函数u(x,y,0)=sinxsiny,微分阶数取为α=0.8,β=0.6,则由式(18),问题(1)的解析解可以表示为

应用差分格式(6)进行数值计算,得到的数值解记为u*(x,y,t),解析解与数值解的误差表示为Err=‖u(x,y,t)-u*(x,y,t)‖2/‖u(x,y,t)‖2.不同时间步长和空间步长的计算结果列于表1、表2.

表1 时间步长与解误差(T=1)

Δt1/501/1001/200Err3．26654×10-23．31394×10-23．34039×10-2

表2 空间步长与解误差(T=1)

hx=hyπ/10π/20π/50Err6．06776×10-23．31394×10-21．35194×10-3

表1、表2的计算结果表明，在当前算法下数值解能够较好地吻合精确解,时间步长变化对解误差影响不大,而随着空间步长的减小,解误差逐步减小.

进一步,取空间步长hx=hy=π/20,时间步长τ=1/100,数值解与解析解在t=T时的图像绘于图1.

图1 α=0.8,β=0.6,t=T时的数值解与解析解

4 确定微分阶数的反问题与数值反演

(19)

则由正问题(1)联合(19)式构成了一个确定微分阶数α和β的反问题.

(20)

其中λ∈(0,1)为同伦参数.当待确定的未知量为常数时,同伦正则化算法过程主要包括迭代逼近、线性化近似、梯度近似与同伦正则化逼近等,而影响算法实现的主要因素是初始迭代、梯度近似中的微分步长、同伦参数、迭代精度或迭代次数以及正问题的计算精度等,具体步骤参见文[11-12]等,且选取同伦参数为

(21)

其中n是迭代次数,n0是当前同伦参数取值下降至0.5时的预估迭代次数,而β0>0是校正参数.下面给出数值反演算例.

设微分阶数的真值为αtrue=[0.8,0.6].模型中其他参数取值同于第3节相应正问题数值算例中的取值, 且式(21)中的预估迭代次数与校正参数分别取为n0=5,β0=0.8.我们主要考察初始迭代与数据扰动对反演算法的影响.

首先看初始迭代对反演算法的影响,反演计算结果列于表3,其中a0表示初始迭代,ainv表示反演解,Err=‖ainv-atrue‖2/‖ainv‖2表示反演解与真解的误差.

表3 初始迭代对反演结果的影响

a0ainvErr(0，0)(0．80000000，0．59999999)3．134×10-13(0．2，0．1)(0．80000000，0．60000000)2．982×10-14(0．5，0．3)(0．80000000，0．60000000)8．385×10-14(1．0，1．0)(1．5，1．3)(2．0，2．0)(0．79999999，0．60000000)(0．60000000，0．80000000)发散2．698×10-132．828×10-1

从表3可以看出,初始迭代的选取对反演结果有一定的影响.当其逐渐远离真值时,反演解与真解的误差变大.

其次,考察数据扰动水平对算法的影响.设带扰动的附加数据表示为:

(22)

表4 扰动水平对反演结果的影响

δ-ainvErr5%(1．23458930，1．01325050)3．685×10-11%(0．81621652，0．61058542)6．693×10-20．5%(0．80227755，0．60374853)9．628×10-30．1%0(0．79872782，0．59549973)(0．80000000，0．59999999)6．697×10-33．134×10-13

从表4可以看出,随着数据扰动水平的减小,反演解与精确解的误差逐渐减小,在不加扰动时,反演解几乎接近真解,表明了反演算法的数值稳定性.

5 结束语

本文给出了含两个时间分数阶导数的时间分数阶二维齐次扩散方程正问题求解的差分格式与解析解,并应用同伦正则化算法进行了微分阶数的数值反演.反演结果表明,当数据扰动水平减小时,反演解与真解的误差变小;当没有扰动时,反演解与真解吻合.当扰动水平大于5%时,解误差变大,这表明所讨论的微分阶数反问题具有一定的病态性.在数据有较大扰动时,构建更有效的反演算法是今后的一项主要工作.

[1]曹建雄. 分数阶扩散方程的有限差方法及其应用[D]. 上海：上海大学, 2015.

[2]陈文, 孙洪广, 李西成,等.力学与工程问题的分数阶导数建模[M]. 北京：科学出版社, 2010.

[3]LIUF,MEERCHAERTMM,MCGOUGHRJ,etal.Numericalmethodsforsolvingthemulti-termtime-fractionalwave-diffusionequation[J].FractionalCalculus&AppliedAnalysis, 2013, 16(1):9-25.

[4]LUCHKOY.Initial-boundary-valueproblemsforthegeneralizedmulti-termtime-fractionaldiffusionequation[J].JournalofMathematicalAnalysis&Applications, 2011, 374(2):538-548.

[5]DAFTARDAR-GEJJIV.Positivesolutionsofasystemofnon-autonomousfractionaldifferentialequations[J].JournalofMathematicalAnalysis&Applications, 2005, 302(1):56-64.

[6]ZHENGGH,WEIT.SpectralregularizationmethodforaCauchyproblemofthetimefractionaladvection-dispersionequation[J].Mathematics&ComputersinSimulation, 2010, 233(10):2631-2640.

[7]ZHENGGH,WEIT.TworegularizationmethodsforsolvingaRiesz-Fellerspace-fractionalbackwarddiffusionproblem[J].InverseProblems, 2010, 26(11):115017.

[8]谷文娟, 李功胜, 殷凤兰,等. 一个时间分数阶扩散方程的参数反演问题[J]. 山东理工大学学报(自然科学版), 2010, 24(6):22-25.

[9]LIUJJ,YAMAMOTOM.Abackwardproblemforthetime-fractionaldiffusionequation[J].ApplicableAnalysis, 2010, 89(11):1 769-1 788.

[10]PODLUBNYI.FractionalDifferentialEquations[M].SanDiego:Academicpress,1999.

[11]孙春龙, 李功胜, 贾现正,等. 含三个时间分数阶导数的反常扩散方程求解与微分阶数反演[J]. 山东理工大学学报(自然科学版), 2015, 29(3):1-7.

[12]贾现正, 张大利, 李功胜,等. 空间-时间分数阶变系数对流扩散方程微分阶数的数值反演[J]. 计算数学, 2014, 36(2):113-132.

(编辑：姚佳良)

The solution to the 2D two-term time-fractional diffusion equation and numerical inversion for the fractional orders

WANG Feng-dan, SUN Chun-long, LI Gong-sheng, JIA Xian-zheng

(School of Science, Shandong University of Technology, Zibo 255049, China)

A finite difference scheme is introduced to solve the 2D two-term time-fractional diffusion equation based on Caputo′s discretization to the time fractional derivatives. Using the method of separation of variables and Laplace transform, analytical solution to the forward problem is obtained, and numerical test is presented to compare the finite difference solution with the analytical solution. Furthermore, the homotopy regularization algorithm is applied to solve the inverse problem of determining the fractional orders given additional data at the final time. Numerical inversions with noisy data are performed, and the inversion solutions error becomes small as the noise level goes to small demonstrating the effectiveness of the proposed algorithm.

time fractional derivative;2D fractional diffusion; analytical solution; inverse problem of fractional order; numerical inversion

2016-05-18

国家自然科学基金项目(11371231, 11071148)

王凤丹,女,949752533@qq.com; 通信作者：李功胜，男，lgs9901@163.com

1672-6197(2017)02-0001-07

O175

A