時間-空間雙邊分數(shù)階擴散方程微分階數(shù)的反演

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

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

時間-空間雙邊分數(shù)階擴散方程微分階數(shù)的反演

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

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

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

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

收稿日期：2015-01-29

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

通信作者：

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

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

中圖分類號：O175

文獻標志碼：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

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

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

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

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

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

(1)

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

(2)

空間雙邊分數(shù)階導數(shù)定義為：

(3)

(4)

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

1.1　正問題求解的差分格式

(5)

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

(6)

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

(7)

當k=0時,

(8)

當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的結構,利用引理1,可得如下命題:

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

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

(12)

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

(13)

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

(14)

1.3　數(shù)值試驗

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

u(x,t)=

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

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

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

2反問題與反演算法

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

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

(16)

則由式(1)聯(lián)合式(16)構成了確定微分階數(shù)α與γ的一個反問題.另外,如果給定區(qū)域內某一點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),設R=(γ,α)∈I,則上述問題可以轉化為如下極小泛函的求解:

(18)

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

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ù)值微分步長.上述泛函極小值問題相當于求解規(guī)范方程:

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

(23)

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

(24)

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

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

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

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

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ù)的反演)設f(x)=sinx,取微分階數(shù)真解Rtrue=[0.6,1.8],初始迭代K0=[0.2,1.4].首先考察微分階數(shù)對反演的影響,反演結果列于表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得出,時間分數(shù)階對反演算法的實施具有一定的影響,隨著時間微分階數(shù)的變大,迭代步數(shù)慢慢增加,而空間分數(shù)階對反演的影響不大.

進一步考察附加數(shù)據(jù)選取的多少對反演算法的影響,計算結果列于表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ù)的擾動對反演結果的影響.由于隨機擾動的影響,每次的反演計算解并不一樣.表4列出了連續(xù)10次反演計算的平均值.

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

表4附加數(shù)據(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ù)的反演)設f(x)=sin2x,微分階數(shù)真解與初始迭代與算例1一樣.本例僅考察附加數(shù)據(jù)的多少以及數(shù)據(jù)擾動對反演的影響,計算結果列于表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ù)的隨機擾動對反演算法的影響比算例1要小,反映出算例1的病態(tài)性較大.

4結束語

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

參考文獻

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

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

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

[4]郭柏靈，蒲學科，黃鳳輝.分數(shù)階偏微分方程及其數(shù)值解[M].北京：科學出版社, 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)正.非對稱分數(shù)階對流彌散的數(shù)值模擬及參數(shù)反演[J]. 高等學校計算數(shù)學學報, 2013, 35(4): 309-326.

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

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

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

[12]蘇麗娟,王文洽,雙邊空間分數(shù)階對流擴散方程的一種有限差分解法[J]. 山東大學學報：理學版，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ù)值解法及應用[M]. 西安:西北工業(yè)出版社, 1995.

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

[17]李功勝，姚德.擴散模型的源項反演及其應用[M].北京:科學出版社, 2014.

(編輯：姚佳良)