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

杜殿虎， 李功勝， 賈現正， 孫春龍

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

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

杜殿虎， 李功勝， 賈現正， 孫春龍

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

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

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

收稿日期：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

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

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

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

1正問題及其數值求解

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

(1)

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

(2)

空間雙邊分數階導數定義為：

(3)

(4)

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

1.1　正問題求解的差分格式

(5)

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

(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　穩定性和收斂性

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

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

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

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

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

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

(12)

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

(13)

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

(14)

1.3　數值試驗

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

u(x,t)=

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

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

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

2反問題與反演算法

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

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

(16)

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

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

(17)

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

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

(18)

其中λ>0是同倫(正則)參數.根據最佳攝動量算法思想,上述極小泛函問題可以轉換成對于給定的初始迭代值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)+

則上述誤差函數可表示為:

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))

τ是數值微分步長.上述泛函極小值問題相當于求解規范方程:

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

(23)

實際反演計算中，同倫參數取為擬神經網絡函數,即有：

(24)

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

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

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

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

3數值反演

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

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

表1時間微分階數對反演的影響

γ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空間微分階數對反演的影響

α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得出,時間分數階對反演算法的實施具有一定的影響,隨著時間微分階數的變大,迭代步數慢慢增加,而空間分數階對反演的影響不大.

進一步考察附加數據選取的多少對反演算法的影響,計算結果列于表3.

表3附加數據選取對反演的影響

附加數據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看出,附加數據選取的多少對反演精度影響不大,只是附加數據選取較多時迭代步數稍微減少.最后,考察附加數據的擾動對反演結果的影響.由于隨機擾動的影響,每次的反演計算解并不一樣.表4列出了連續10次反演計算的平均值.

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

表4附加數據的隨機擾動對反演的影響

δ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(基于區域中點處觀測數據的反演)設f(x)=sin2x,微分階數真解與初始迭代與算例1一樣.本例僅考察附加數據的多少以及數據擾動對反演的影響,計算結果列于表5與6.

表5附加數據選取對反演的影響

附加數據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附加數據的擾動對反演的影響

δ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一樣,附加數據選取的多少對反演精度影響也不大,較多的附加數據所需的迭代步數較少;不過,本例中數據的隨機擾動對反演算法的影響比算例1要小,反映出算例1的病態性較大.

4結束語

本文探討了時間-空間雙邊分數階擴散方程微分階數的反演問題.文中建立了正問題求解的一種隱式差分格式,并給出了差分格式的穩定性和收斂性;進一步對于時間、空間兩個微分階數的反演問題,應用一種同倫正則化算法給出了精確數據與隨機擾動數據情形下的數值反演.計算結果表明應用適當的反演算法可以有效實現對于復雜的分數階擴散方程相關系數反問題的數值反演.

參考文獻

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

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

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

[4]郭柏靈，蒲學科，黃鳳輝.分數階偏微分方程及其數值解[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]李新潔，李功勝，賈現正.非對稱分數階對流彌散的數值模擬及參數反演[J]. 高等學校計算數學學報, 2013, 35(4): 309-326.

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

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

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

[12]蘇麗娟,王文洽,雙邊空間分數階對流擴散方程的一種有限差分解法[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]蘇超偉.偏微分方程逆問題的數值解法及應用[M]. 西安:西北工業出版社, 1995.

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

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

(編輯：姚佳良)