王鳳丹，孫春龍，李功勝，賈現(xiàn)正

(山東理工大學(xué) 理學(xué)院，山東 淄博 255049)

含兩個(gè)時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的二維反常擴(kuò)散方程求解與微分階數(shù)反演

王鳳丹，孫春龍，李功勝，賈現(xiàn)正

(山東理工大學(xué) 理學(xué)院，山東 淄博 255049)

對(duì)于一類(lèi)含兩個(gè)時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的二維反常擴(kuò)散方程,基于對(duì)時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)在Caputo意義下的離散,得到一個(gè)有限差分格式;利用分離變量法與Laplace變換得到該問(wèn)題的解析解,并將兩種方法得到的解進(jìn)行數(shù)值比較.進(jìn)一步,給定終值時(shí)刻數(shù)據(jù),應(yīng)用同倫正則化算法對(duì)擴(kuò)散方程中的兩個(gè)時(shí)間微分階數(shù)進(jìn)行數(shù)值反演,并給出反演算例.數(shù)值結(jié)果表明隨著數(shù)據(jù)擾動(dòng)水平的降低,解誤差逐步變小,所用的反演算法對(duì)微分階數(shù)反問(wèn)題是有效的.

時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)；二維分?jǐn)?shù)階擴(kuò)散；解析解；微分階數(shù)反問(wèn)題；數(shù)值反演

近年來(lái),反常擴(kuò)散模型的應(yīng)用越來(lái)越廣泛.分?jǐn)?shù)階擴(kuò)散方程在大氣污染、水文地質(zhì)學(xué)、金融學(xué)、物理學(xué)、力學(xué)、生物醫(yī)學(xué)工程等諸多領(lǐng)域有了若干非常成功的應(yīng)用[1],分?jǐn)?shù)階擴(kuò)散模型相比經(jīng)典的整數(shù)階高斯擴(kuò)散模型能更好的模擬擴(kuò)散過(guò)程和重建試驗(yàn)數(shù)據(jù).分?jǐn)?shù)階導(dǎo)數(shù)具有記憶性、遺傳性和整體性,特別對(duì)于污染物長(zhǎng)距離傳輸時(shí)表現(xiàn)出的非對(duì)稱(chēng)、非線性的整體性行為模式的研究,分?jǐn)?shù)階擴(kuò)散模型可能是一個(gè)很有效的研究工具[2].對(duì)于含多個(gè)時(shí)間分?jǐn)?shù)階擴(kuò)散方程正問(wèn)題的研究,劉發(fā)旺和Meerschaert[3]等人,研究了含多個(gè)時(shí)間分?jǐn)?shù)階擴(kuò)散-波動(dòng)方程的數(shù)值計(jì)算方法,給出了隱式差分格式,并證明了隱式差分近似的無(wú)條件穩(wěn)定性和收斂性. Luchko[4]應(yīng)用分離變量法與傅里葉變換研究了一類(lèi)含多個(gè)時(shí)間分?jǐn)?shù)階擴(kuò)散方程的解析解,并證明了解的唯一性.Gejji[5]等利用分離變量法推導(dǎo)出了一維非齊次和二維齊次時(shí)間分?jǐn)?shù)階擴(kuò)散-波動(dòng)方程在齊次混合邊值條件下的解析解.

對(duì)于分?jǐn)?shù)階擴(kuò)散方程反問(wèn)題的研究,越來(lái)越受到人們的關(guān)注. 鄭光輝和魏婷[6-7]等應(yīng)用傅里葉正則化,研究了帶形區(qū)域上空間分?jǐn)?shù)階擴(kuò)散方程的逆時(shí)問(wèn)題,同時(shí)考慮了空間分?jǐn)?shù)階擴(kuò)散方程中的源項(xiàng)反演問(wèn)題,并且應(yīng)用Tikhonov正則化對(duì)反問(wèn)題數(shù)值求解.李功勝等[8]利用最佳攝動(dòng)量算法對(duì)于一維時(shí)間分?jǐn)?shù)階擴(kuò)散方程的微分階數(shù)與擴(kuò)散系數(shù)進(jìn)行了數(shù)值聯(lián)合反演. 劉繼軍等[9]應(yīng)用準(zhǔn)可逆性正則化方法探討了一個(gè)時(shí)間分?jǐn)?shù)階擴(kuò)散中由終值數(shù)據(jù)確定初始分布的逆時(shí)反問(wèn)題.

目前,關(guān)于高維分?jǐn)?shù)階擴(kuò)散方程相關(guān)反問(wèn)題的研究還不多見(jiàn).本文主要探討含兩個(gè)時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的二維時(shí)間分?jǐn)?shù)階齊次擴(kuò)散方程的微分階數(shù)反演問(wèn)題.對(duì)于正問(wèn)題的求解,通過(guò)對(duì)時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的數(shù)值離散,建立了一個(gè)隱式差分格式；同時(shí)利用分離變量法與Laplace變換得到該問(wèn)題的解析解. 在正問(wèn)題求解的基礎(chǔ)上,給定終值時(shí)刻數(shù)據(jù)為附加數(shù)據(jù)，考慮一個(gè)確定擴(kuò)散方程中兩個(gè)時(shí)間微分階數(shù)的反問(wèn)題.應(yīng)用同倫正則化算法進(jìn)行數(shù)值反演,并討論初始迭代、擾動(dòng)水平等因素對(duì)反演算法的影響.

1 正問(wèn)題及其數(shù)值解

考慮矩形域Ω=(0,l1)×(0,l2)上的二維時(shí)間分?jǐn)?shù)階齊次擴(kuò)散方程的初邊值問(wèn)題：

(1)

其中(x,y)∈Ω,t>0,D>0是擴(kuò)散系數(shù),α,β∈(0,1)是時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的階數(shù),方程中的時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)均用Caputo的定義，分別為：

下面給出數(shù)值求解二維問(wèn)題(1)的一個(gè)隱式差分格式.

(2)

(3)

利用通常的差分離散方法,方程中的整數(shù)階導(dǎo)數(shù)項(xiàng)可以進(jìn)行以下近似：

則原方程可以離散為

(4)

(5)

如果令ck=dk-1-dk,k=1,2,…,n，則差分格式可以進(jìn)一步整理為矩陣形式為

(6)

其中

(7)

i=1,2,…,M-1; I為(M-1)×(M-1)階單位矩陣.

2 正問(wèn)題的解析解

本節(jié)應(yīng)用分離變量法和拉普拉斯變換,給出正問(wèn)題(1)的基于Mittag-Leffler函數(shù)的解析解表達(dá)式,并進(jìn)行數(shù)值模擬.

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

(8)

對(duì)式(8)整理可得

(9)

(10)

(11)

(12)

(13)

對(duì)t作Laplace變換可得

(14)

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

(15)

(16)

令c=D(λn+μm)(常數(shù)).

(17)

由引理1，即有

于是，可得

(18)

3 正問(wèn)題的數(shù)值模擬

取l1=l2=π,T=1,擴(kuò)散系數(shù)D=0.5,初始函數(shù)u(x,y,0)=sinxsiny,微分階數(shù)取為α=0.8,β=0.6,則由式(18),問(wèn)題(1)的解析解可以表示為

應(yīng)用差分格式(6)進(jìn)行數(shù)值計(jì)算,得到的數(shù)值解記為u*(x,y,t),解析解與數(shù)值解的誤差表示為Err=‖u(x,y,t)-u*(x,y,t)‖2/‖u(x,y,t)‖2.不同時(shí)間步長(zhǎng)和空間步長(zhǎng)的計(jì)算結(jié)果列于表1、表2.

表1 時(shí)間步長(zhǎng)與解誤差(T=1)

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

表2 空間步長(zhǎng)與解誤差(T=1)

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

表1、表2的計(jì)算結(jié)果表明，在當(dāng)前算法下數(shù)值解能夠較好地吻合精確解,時(shí)間步長(zhǎng)變化對(duì)解誤差影響不大,而隨著空間步長(zhǎng)的減小,解誤差逐步減小.

進(jìn)一步,取空間步長(zhǎng)hx=hy=π/20,時(shí)間步長(zhǎng)τ=1/100,數(shù)值解與解析解在t=T時(shí)的圖像繪于圖1.

圖1 α=0.8,β=0.6,t=T時(shí)的數(shù)值解與解析解

4 確定微分階數(shù)的反問(wèn)題與數(shù)值反演

(19)

則由正問(wèn)題(1)聯(lián)合(19)式構(gòu)成了一個(gè)確定微分階數(shù)α和β的反問(wèn)題.

(20)

其中λ∈(0,1)為同倫參數(shù).當(dāng)待確定的未知量為常數(shù)時(shí),同倫正則化算法過(guò)程主要包括迭代逼近、線性化近似、梯度近似與同倫正則化逼近等,而影響算法實(shí)現(xiàn)的主要因素是初始迭代、梯度近似中的微分步長(zhǎng)、同倫參數(shù)、迭代精度或迭代次數(shù)以及正問(wèn)題的計(jì)算精度等,具體步驟參見(jiàn)文[11-12]等,且選取同倫參數(shù)為

(21)

其中n是迭代次數(shù),n0是當(dāng)前同倫參數(shù)取值下降至0.5時(shí)的預(yù)估迭代次數(shù),而β0>0是校正參數(shù).下面給出數(shù)值反演算例.

設(shè)微分階數(shù)的真值為αtrue=[0.8,0.6].模型中其他參數(shù)取值同于第3節(jié)相應(yīng)正問(wèn)題數(shù)值算例中的取值, 且式(21)中的預(yù)估迭代次數(shù)與校正參數(shù)分別取為n0=5,β0=0.8.我們主要考察初始迭代與數(shù)據(jù)擾動(dòng)對(duì)反演算法的影響.

首先看初始迭代對(duì)反演算法的影響,反演計(jì)算結(jié)果列于表3,其中a0表示初始迭代,ainv表示反演解,Err=‖ainv-atrue‖2/‖ainv‖2表示反演解與真解的誤差.

表3 初始迭代對(duì)反演結(jié)果的影響

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)發(fā)散2．698×10-132．828×10-1

從表3可以看出,初始迭代的選取對(duì)反演結(jié)果有一定的影響.當(dāng)其逐漸遠(yuǎn)離真值時(shí),反演解與真解的誤差變大.

其次,考察數(shù)據(jù)擾動(dòng)水平對(duì)算法的影響.設(shè)帶擾動(dòng)的附加數(shù)據(jù)表示為:

(22)

表4 擾動(dòng)水平對(duì)反演結(jié)果的影響

δ-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可以看出,隨著數(shù)據(jù)擾動(dòng)水平的減小,反演解與精確解的誤差逐漸減小,在不加擾動(dòng)時(shí),反演解幾乎接近真解,表明了反演算法的數(shù)值穩(wěn)定性.

5 結(jié)束語(yǔ)

本文給出了含兩個(gè)時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的時(shí)間分?jǐn)?shù)階二維齊次擴(kuò)散方程正問(wèn)題求解的差分格式與解析解,并應(yīng)用同倫正則化算法進(jìn)行了微分階數(shù)的數(shù)值反演.反演結(jié)果表明,當(dāng)數(shù)據(jù)擾動(dòng)水平減小時(shí),反演解與真解的誤差變小;當(dāng)沒(méi)有擾動(dòng)時(shí),反演解與真解吻合.當(dāng)擾動(dòng)水平大于5%時(shí),解誤差變大,這表明所討論的微分階數(shù)反問(wèn)題具有一定的病態(tài)性.在數(shù)據(jù)有較大擾動(dòng)時(shí),構(gòu)建更有效的反演算法是今后的一項(xiàng)主要工作.

[1]曹建雄. 分?jǐn)?shù)階擴(kuò)散方程的有限差方法及其應(yīng)用[D]. 上海：上海大學(xué), 2015.

[2]陳文, 孫洪廣, 李西成,等.力學(xué)與工程問(wèn)題的分?jǐn)?shù)階導(dǎo)數(shù)建模[M]. 北京：科學(xué)出版社, 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]谷文娟, 李功勝, 殷鳳蘭,等. 一個(gè)時(shí)間分?jǐn)?shù)階擴(kuò)散方程的參數(shù)反演問(wèn)題[J]. 山東理工大學(xué)學(xué)報(bào)(自然科學(xué)版), 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]孫春龍, 李功勝, 賈現(xiàn)正,等. 含三個(gè)時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的反常擴(kuò)散方程求解與微分階數(shù)反演[J]. 山東理工大學(xué)學(xué)報(bào)(自然科學(xué)版), 2015, 29(3):1-7.

[12]賈現(xiàn)正, 張大利, 李功勝,等. 空間-時(shí)間分?jǐn)?shù)階變系數(shù)對(duì)流擴(kuò)散方程微分階數(shù)的數(shù)值反演[J]. 計(jì)算數(shù)學(xué), 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

國(guó)家自然科學(xué)基金項(xiàng)目(11371231, 11071148)

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

1672-6197(2017)02-0001-07

O175

A