杜殿虎， 李功勝， 賈現(xiàn)正， 孫春龍

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

時間-空間雙邊分?jǐn)?shù)階擴(kuò)散方程微分階數(shù)的反演

杜殿虎， 李功勝， 賈現(xiàn)正， 孫春龍

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

摘要：研究了一維時間-空間雙邊分?jǐn)?shù)階擴(kuò)散方程的求解與微分階數(shù)的數(shù)值反演問題.基于Caputo意義下時間分?jǐn)?shù)階導(dǎo)數(shù)和Grünward-Letnikov意義下空間雙邊分?jǐn)?shù)階導(dǎo)數(shù)的離散,給出了一個有限差分求解格式,證明了其穩(wěn)定性和收斂性.分別基于終值數(shù)據(jù)及區(qū)域中點處的觀測值作為附加數(shù)據(jù),應(yīng)用同倫正則化算法對微分階數(shù)進(jìn)行數(shù)值反演.反演結(jié)果表明同倫正則化算法對于分?jǐn)?shù)階擴(kuò)散方程的微分階數(shù)反演是有效的.

關(guān)鍵詞：分?jǐn)?shù)階擴(kuò)散方程； 反問題； 差分格式； 穩(wěn)定性與收斂性； 同倫正則化； 數(shù)值反演

收稿日期：2015-01-29

基金項目：國家自然科學(xué)基金項目(11371231, 11071148)； 山東省自然科學(xué)基金項目(ZR2011AQ014)

通信作者：

作者簡介：杜殿虎,男,dudianhu@126.com；李功勝，男，lgs9901@163.com

文章編號：1672-6197(2016)01-0021-05

中圖分類號：O175

文獻(xiàn)標(biāo)志碼：A

Abstract:This paper deals with numerical solution to the time-space Caputo-Riesz fractional diffusion equation and inversion for the fractional orders. A finite difference scheme is introduced to solve the forward problem based on Caputo’s discretization to the time fractional derivative and Grünward-Letnikov’s discretization to the space fractional derivatives, and unconditional stability and convergence for the difference scheme are proved by fine estimation to the spectral radius of the coefficient matrix of the scheme. Furthermore, numerical inversions for determining the fractional orders are carried out with random noisy data at the final time and the inner point using the homotopy regularization algorithm. The numerical solutions give good approximations to the exact solution demonstrating the effectiveness of the proposed algorithm

Numerical determination of the fractional orders for

time-space two-sided fractional diffusion equation

DU Dian-hu, LI Gong-sheng, JIA Xian-zheng, SUN Chun-long

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

Key words: fractional diffusion equation; inverse problem; finite difference scheme; stability and convergence; homotopy regularization algorithm; numerical inversion

隨著近年來應(yīng)用學(xué)科的發(fā)展,對于分?jǐn)?shù)階反常擴(kuò)散模型的研究逐漸成為現(xiàn)代各門科學(xué)與數(shù)學(xué)相互驅(qū)動發(fā)展的一個主要方面.已有研究表明,與傳統(tǒng)的整數(shù)階微分方程比較,分?jǐn)?shù)階微分方程能夠更好地擬合某些自然物理過程和動態(tài)系統(tǒng)過程,因此其在工程、物理、地下水、金融和環(huán)境問題中得到了廣泛的應(yīng)用[1-4].不過從數(shù)學(xué)理論及其方法研究的角度,分?jǐn)?shù)階微分方程還有許多問題值得深入研究.

對于分?jǐn)?shù)階擴(kuò)散的相關(guān)反問題研究,程晉等[5]人首先研究了時間分?jǐn)?shù)階反常擴(kuò)散方程在齊次Neumann邊界及δ函數(shù)為初始條件下,基于邊界測量數(shù)據(jù)確定分?jǐn)?shù)階數(shù)和擴(kuò)散系數(shù)的聯(lián)合反演問題,并證明了反問題解的唯一性；隨后李功勝等[6]繼續(xù)研究了一般連續(xù)初值條件下對于時間微分階數(shù)與擴(kuò)散系數(shù)聯(lián)合反演的唯一性,并給出了數(shù)值反演模擬.最近,東京大學(xué)的山本昌弘教授等[7]考察了含多個時間階導(dǎo)數(shù)的擴(kuò)散方程中多個微分階數(shù)的反問題，也給出了反問題解的唯一性證明；李新潔等[8]研究了非對稱雙邊空間分?jǐn)?shù)階擴(kuò)散方程的全離散差分格式求解,并應(yīng)用最佳攝動量算法對微分階數(shù)與擴(kuò)散系數(shù)進(jìn)行了聯(lián)合反演模擬；賈現(xiàn)正等[9]探討了時間-空間左側(cè)分?jǐn)?shù)階變系數(shù)對流擴(kuò)散方程差分求解與微分階數(shù)的數(shù)值反演問題,通過對系數(shù)矩陣譜半徑的精細(xì)估計,證明了差分格式的無條件穩(wěn)定性和收斂性,同時應(yīng)用同倫正則化算法對時間-空間微分階數(shù)進(jìn)行了有效的數(shù)值反演.

本文將探討時間-空間雙邊分?jǐn)?shù)階擴(kuò)散方程微分階數(shù)的反演問題.對于這類同時含有時間分?jǐn)?shù)階和空間雙側(cè)分?jǐn)?shù)階導(dǎo)數(shù)的擴(kuò)散方程正問題的數(shù)值求解,已有不少研究[10-14].主要是差分方法,且其差分格式收斂性與穩(wěn)定性的證明是基于逐行估計的方法.本文繼續(xù)研究該正問題的差分求解,不過文中將采用一種新的證明方法,即通過對系數(shù)矩陣主對角元素的精細(xì)估計,獲得系數(shù)矩陣的譜半徑估計,進(jìn)而直接證明差分格式的無條件穩(wěn)定性和收斂性.這種方法相比已有方法顯得簡潔清晰.進(jìn)一步,本文還將探討利用終值數(shù)據(jù)和內(nèi)點觀測數(shù)據(jù)確定時間、空間微分階數(shù)的反問題.應(yīng)用同倫正則化算法,分別進(jìn)行精確數(shù)據(jù)和擾動數(shù)據(jù)條件下的數(shù)值反演模擬,反演解與真解吻合較好.結(jié)果表明對于所研究的時間-空間雙邊分?jǐn)?shù)階擴(kuò)散模型,同倫正則化算法可以有效實現(xiàn)對微分階數(shù)的數(shù)值確定.

1正問題及其數(shù)值求解

對于l>0,T>0,考慮時間-空間雙邊分?jǐn)?shù)階擴(kuò)散方程的初邊值問題：

(1)

其中0<γ<1為時間分?jǐn)?shù)階導(dǎo)數(shù),1<α<2為空間分?jǐn)?shù)階導(dǎo)數(shù),D>0為擴(kuò)散系數(shù),f(x)∈C([0,l])為初值函數(shù).這里的時間分?jǐn)?shù)階導(dǎo)數(shù)取Caputo的定義：

(2)

空間雙邊分?jǐn)?shù)階導(dǎo)數(shù)定義為：

(3)

(4)

下面先給出求解正問題(1)的一個隱式差分格式.

1.1　正問題求解的差分格式

(5)

其次,根據(jù)式(4)并利用移位的Grünward近似公式,空間雙邊分?jǐn)?shù)階導(dǎo)數(shù)離散為:

(6)

于是,略去高階項,原方程被離散為：

(7)

當(dāng)k=0時,

(8)

當(dāng)k>0時,

(i=1,2…,M-1)

(9)

用矩陣形式表示,即有：

(10)

f=(f(x1),f(x2),…,f(xM-1))′;

cj=2j1-γ-(j+1)1-γ-(j-1)1-γ,

j=1,2,…,k;bk=(k+1)1-γ-k1-γ,

k=1,2,…,N-1;

A=(ai,j),i,j=1,2…

M-1,aij定義為:

(11)

1.2　穩(wěn)定性和收斂性

本小節(jié)給出差分格式(10)的穩(wěn)定性和收斂性.

引理3對于矩陣A∈R(M-1)×(M-1)以及任意的ε>0,至少存在一種矩陣范數(shù)‖·‖*使得:

‖A‖*≤ρ(A)+ε.

注意到由式(10)給出的矩陣A的結(jié)構(gòu),利用引理1,可得如下命題:

由命題1可知,以A為系數(shù)矩陣的差分格式有唯一解.進(jìn)一步可得系數(shù)矩陣的譜半徑估計.

定理1對于式(10)定義的系數(shù)矩陣A,記‖a‖,

(12)

利用差分格式線性性,容易得到：

(13)

定理2隱式差分格式(10)無條件穩(wěn)定.

(14)

1.3　數(shù)值試驗

運(yùn)用分離變量法及Laplace變換,可得正問題(1)的解析解表達(dá)式:

u(x,t)=

分別利用上述解析解表達(dá)式與差分格式(10)進(jìn)行正問題的計算,數(shù)值結(jié)果繪于圖1.

圖1 　α=1.9,γ=0.9,t=0.5時解析解與數(shù)值解的比較,-:解析解,*:數(shù)值解

從圖1看出，除了在區(qū)域中點處解誤差稍大以外,數(shù)值解與解析解都吻合得較好,說明了差分格式的有效性.

2反問題與反演算法

當(dāng)微分階數(shù)α與γ未知時,給定時刻t=T的觀測值

u(x,T)=θ(x),0

(16)

則由式(1)聯(lián)合式(16)構(gòu)成了確定微分階數(shù)α與γ的一個反問題.另外,如果給定區(qū)域內(nèi)某一點x=x0處的時間觀測值作為附加數(shù)據(jù),即有附加條件:

u(x0,t)=θ(t),0≤t≤T

(17)

則由(1)聯(lián)合(17)也形成一個確定微分階數(shù)的反問題.

本文對γ,α的數(shù)值反演采用同倫正則化算法[15-17].記I=(0,1)×(1,2),設(shè)R=(γ,α)∈I,則上述問題可以轉(zhuǎn)化為如下極小泛函的求解:

(18)

其中λ>0是同倫(正則)參數(shù).根據(jù)最佳攝動量算法思想,上述極小泛函問題可以轉(zhuǎn)換成對于給定的初始迭代值R0求解最佳攝動量δRj,即設(shè)

Rj+1=Rj+δRj,j=0,1,2,…

(19)

δRj是下述泛函的極小解:

F(δRj)=(1-λ)‖u(x,T;Rj+δRj)-

(20)

將u(x,T;Rj+δRj)在Rj處作Taylor展開,略去高階項得到:

u(x,T;Rj+δRj)≈u(x,T;Rj)+

則上述誤差函數(shù)可表示為:

F(δRj)=(1-λ)‖Tu(x,T;Rj)·δRj+

(21)

(22)

其中:

B=(bks)K×S,bks=

ξ=(u(x1,T;Rj),u(x2,T;Rj)，…，u(xk,T;Rj)),η=(θ(x1),θ(x2),…,θ(xK))

τ是數(shù)值微分步長.上述泛函極小值問題相當(dāng)于求解規(guī)范方程:

((1-λ)BTB+λI)δRj=(1-λ)BT(η-ξ)

(23)

實際反演計算中，同倫參數(shù)取為擬神經(jīng)網(wǎng)絡(luò)函數(shù),即有：

(24)

其中N0,β為可調(diào)參數(shù),j為迭代步數(shù),λ(j)為第j個迭代步時同倫參數(shù)的取值.上述算法實施步驟為：

步驟1給定初始迭代值R0和數(shù)值微分步長τ,求向量η,ξ及矩陣B;

步驟2按照同倫算法進(jìn)行計算,由式(23)求出δRj;

步驟3對于給定的精度eps,判斷是否滿足‖δRj‖≤eps;若是則Rj即為所求,算法終止;否則由式(20)得到Rj+1,轉(zhuǎn)到步驟1繼續(xù)進(jìn)行.

3數(shù)值反演

以下算法實現(xiàn)中,若沒有特別說明,正問題計算中取l=π,T=1,D=0.5,M=N=100,反演參數(shù)分別取為τ=0.01,eps=le-5,N0=5，β=0.5.

算例1 (基于終值數(shù)據(jù)的反演)設(shè)f(x)=sinx,取微分階數(shù)真解Rtrue=[0.6,1.8],初始迭代K0=[0.2,1.4].首先考察微分階數(shù)對反演的影響,反演結(jié)果列于表1與2,其中γ,α分別表示時間與空間微分階數(shù),Rinv表示反演解,j為迭代次數(shù),Err表示反演解與真解的相對誤差.

表1時間微分階數(shù)對反演的影響

γRinvErrj0.1[0.100000,1.800000]5.7544×10-8180.2[0.200000,1.800000]8.3978×10-8180.3[0.300000,1.800000]1.3684×10-7180.4[0.400000,1.800000]2.6642×10-7180.5[0.500000,1.800000]2.1568×10-8190.6[0.600000,1.800000]5.3187×10-8190.7[0.700000,1.800000]6.6357×10-8210.8[0.800000,1.800000]1.5706×10-8250.9[0.900000,1.800000]2.4516×10-754

表2空間微分階數(shù)對反演的影響

αRinvErrj1.3[0.600000,1.300000]5.1302×10-7201.4[0.600000,1.400000]5.4895×10-8191.5[0.600000,1.500000]4.8791×10-8191.6[0.600000,1.600000]7.0265×10-8191.7[0.600000,1.700000]9.3387×10-8191.8[0.600000,1.800000]5.3187×10-8191.9[0.600000,1.900000]2.8764×10-719

由表1,2得出,時間分?jǐn)?shù)階對反演算法的實施具有一定的影響,隨著時間微分階數(shù)的變大,迭代步數(shù)慢慢增加,而空間分?jǐn)?shù)階對反演的影響不大.

進(jìn)一步考察附加數(shù)據(jù)選取的多少對反演算法的影響,計算結(jié)果列于表3.

表3附加數(shù)據(jù)選取對反演的影響

附加數(shù)據(jù)RinvErrj[0.3,0.4][0.600000,1.800000]1.0315×10-832[0.3,0.4,0.5][0.600000,1.800000]1.1956×10-730[0.3,0.4,0.5,0.6][0.600000,1.800000]9.7956×10-830[0.3,0.4,0.5,0.6,0.7][0.600000,1.800000]1.3468×10-729[0.3,0.4,0.5,0.6,0.7,0.8][0.600000,1.800000]3.0919×10-827

由表3看出,附加數(shù)據(jù)選取的多少對反演精度影響不大,只是附加數(shù)據(jù)選取較多時迭代步數(shù)稍微減少.最后,考察附加數(shù)據(jù)的擾動對反演結(jié)果的影響.由于隨機(jī)擾動的影響,每次的反演計算解并不一樣.表4列出了連續(xù)10次反演計算的平均值.

由表4得出,數(shù)據(jù)擾動對反演結(jié)果影響較大.不過，隨著數(shù)據(jù)擾動水平的減小,反演精度逐步提高,這反映了反演算法的數(shù)值穩(wěn)定性.

表4附加數(shù)據(jù)的隨機(jī)擾動對反演的影響

δRinvErrj　2%[0.698726,1.795743]5.2082×10-2211%[0.628209,1.794657]1.5132×10-2200.5%[0.613034,1.797258]7.0197×10-3190.1%[0.602487,1.799442]1.3434×10-3190.01%[0.600246,1.799944]1.3313×10-419

算例2(基于區(qū)域中點處觀測數(shù)據(jù)的反演)設(shè)f(x)=sin2x,微分階數(shù)真解與初始迭代與算例1一樣.本例僅考察附加數(shù)據(jù)的多少以及數(shù)據(jù)擾動對反演的影響,計算結(jié)果列于表5與6.

表5附加數(shù)據(jù)選取對反演的影響

附加數(shù)據(jù)RinvErrj[0.3,0.4][0.600000,1.800000]9.0495×10-930[0.3,0.4,0.5][0.600000,1.800000]1.0009×10-827[0.3,0.4,0.5,0.6][0.600000,1.800000]8.6784×10-925[0.3,0.4,0.5,0.6,0.7][0.600000,1.800000]4.1451×10-823[0.3,0.4,0.5,0.6,0.7,0.8][0.600000,1.800000]1.4238×10-822

表6附加數(shù)據(jù)的擾動對反演的影響

δRinvErrj　5%[0.606782,1.825859]1.4090×10-2151%[0.601370,1.805340]2.9055×10-3150.5%[0.600686,1.802681]1.4585×10-3150.1%[0.600137,1.800538]2.9259×10-4150.01%[0.600014,1.800054]2.9217×10-515

由表5、6可以看出,與算例1一樣,附加數(shù)據(jù)選取的多少對反演精度影響也不大,較多的附加數(shù)據(jù)所需的迭代步數(shù)較少;不過,本例中數(shù)據(jù)的隨機(jī)擾動對反演算法的影響比算例1要小,反映出算例1的病態(tài)性較大.

4結(jié)束語

本文探討了時間-空間雙邊分?jǐn)?shù)階擴(kuò)散方程微分階數(shù)的反演問題.文中建立了正問題求解的一種隱式差分格式,并給出了差分格式的穩(wěn)定性和收斂性;進(jìn)一步對于時間、空間兩個微分階數(shù)的反演問題,應(yīng)用一種同倫正則化算法給出了精確數(shù)據(jù)與隨機(jī)擾動數(shù)據(jù)情形下的數(shù)值反演.計算結(jié)果表明應(yīng)用適當(dāng)?shù)姆囱菟惴梢杂行崿F(xiàn)對于復(fù)雜的分?jǐn)?shù)階擴(kuò)散方程相關(guān)系數(shù)反問題的數(shù)值反演.

參考文獻(xiàn)

[1]PodlubnyI.Fractionaldifferentialequations[M].SanDiego:Academic, 1999.

[2]KilbasAA,SrivastavaHM,TrujilloJJ.Theoryandapplicationsoffractionaldifferentialequations[M].Amsterdam:Elsevier, 2006.

[3]陳文，孫洪廣，李西成,等.力學(xué)與工程問題的分?jǐn)?shù)階導(dǎo)數(shù)建模[M].科學(xué)出版社, 2010.

[4]郭柏靈，蒲學(xué)科，黃鳳輝.分?jǐn)?shù)階偏微分方程及其數(shù)值解[M].北京：科學(xué)出版社, 2011.

[5]ChengJ,NakagawaJ,YamamotoM，et al，Uniquenessinaninverseproblemforaone-dimensionalfractionaldiffusionequation[J].InverseProblems, 2009, 25: 115002.

[6]LiGS,ZhangDL,JiaXZ, et al.Simultaneousinversionforthespace-dependentdiffusioncoefficientandthefractionalorderinthetime-fractionaldiffusionequation[J].InverseProblems, 2013, 29: 065014.

[7]LiZY,YamamotoM.Uniquenessforinverseproblemsofdeterminingordersofmulti-termtime-fractionalderivativesofdiffusionequation[J].ApplicableAnalysis,2015,94(3):570-579.

[8]李新潔，李功勝，賈現(xiàn)正.非對稱分?jǐn)?shù)階對流彌散的數(shù)值模擬及參數(shù)反演[J]. 高等學(xué)校計算數(shù)學(xué)學(xué)報, 2013, 35(4): 309-326.

[9]賈現(xiàn)正，張大利，李功勝,等.空間-時間分?jǐn)?shù)階變系數(shù)對流擴(kuò)散方程微分階數(shù)的數(shù)值反演[J].計算數(shù)學(xué), 2014, 36(2): 113-132.

[10]LiuF,ZhuangP,AnhV,et al.Stabilityandconvergenceofthedifferencemethodsforthespace-timefractionaladvection-diffusionequation[J].AppliedMathematicsandComputation, 2007, 191: 12-20.

[11]覃平陽,張曉丹,空間-時間分?jǐn)?shù)階對流擴(kuò)散方程的數(shù)值解法[J].計算數(shù)學(xué), 2008, 30: 305-310.

[12]蘇麗娟,王文洽,雙邊空間分?jǐn)?shù)階對流擴(kuò)散方程的一種有限差分解法[J]. 山東大學(xué)學(xué)報：理學(xué)版，2009, 44(10): 26-29.

[13]ShenSJ,LiuFW,AnhV.Numericalapproximationsandsolutiontechniquesforthespace-timeRiesz-Caputofractionaladvection-diffusionequation[J].NumericalAlgorithms, 2011, 56: 383-403.

[14]ChenZQ,MeerschaertMM,NaneE.Space-timefractionaldiffusiononboundeddomains[J].JournalofMathematicalAnalysisandApplications, 2012, 393: 479-488.

[15]蘇超偉.偏微分方程逆問題的數(shù)值解法及應(yīng)用[M]. 西安:西北工業(yè)出版社, 1995.

[16]LiGS,ChengJ,YaoD, et al.One-dimensionalequilibriummodelandsourceparameterdeterminationforsoil-columnexperiment[J].AppliedMathematicsandComputation, 2007, 190: 1365-1374.

[17]李功勝，姚德.擴(kuò)散模型的源項反演及其應(yīng)用[M].北京:科學(xué)出版社, 2014.

(編輯：姚佳良)