池光勝,李功勝,2

(1.內(nèi)蒙古工業(yè)大學 理學院,呼和浩特 010051; 2.山東理工大學 理學院,淄博 255049)

分數(shù)階擴散方程中多參數(shù)聯(lián)合數(shù)值反演

池光勝1,李功勝1,2

(1.內(nèi)蒙古工業(yè)大學 理學院,呼和浩特 010051; 2.山東理工大學 理學院,淄博 255049)

考慮一維空間分數(shù)階擴散方程中同時確定空間微分階數(shù)、擴散系數(shù)和源項的多參數(shù)聯(lián)合反演問題.基于空間分數(shù)階導數(shù)離散系數(shù)的性質(zhì)和快速傅里葉變換,給出求解正問題的一種新的差分格式.利用同倫正則化方法對所提的空間微分階數(shù)、擴散系數(shù)和源項的多參數(shù)聯(lián)合反演問題進行精確數(shù)據(jù)與擾動數(shù)據(jù)情形下的數(shù)值反演.結(jié)果表明： 同倫正則化算法對空間分數(shù)階擴散中的多參數(shù)聯(lián)合反演是有效的.

分數(shù)階擴散方程； 差分格式； 多參數(shù)聯(lián)合反問題； 數(shù)值反演

1 研究背景

作為整數(shù)階的推廣,分數(shù)階具有非局部性,可以更好地描述生活中整數(shù)階無法刻劃的復雜系統(tǒng).近年來,分數(shù)階擴散方程應(yīng)用于地下水污染遷移、生物學、物理學,甚至在金融學中都起到越來越重要的作用.有關(guān)研究表明： 在地下水污染遷移中,分數(shù)階擴散方程在描述多孔介質(zhì)中具有非費克現(xiàn)象的溶質(zhì)運移過程時,能夠更好地模擬溶質(zhì)提前穿透或嚴重拖尾等非線性遷移行為.若考慮非局域性和空間相關(guān)性,則得到空間分數(shù)階擴散方程.另一方面,一旦建立了一個分數(shù)階擴散的初步模型,求解和分析模型中難以直接測量的參數(shù)就成為一個很重要的問題.這就導致了對于分數(shù)階擴散方程及其參數(shù)反演問題的研究.可以說,分數(shù)階擴散相關(guān)反問題研究對于認識污染物在多孔介質(zhì)中的反常擴散行為具有重要的科學意義和應(yīng)用前景.

近20年以來國內(nèi)外學者對于空間分數(shù)階擴散方程正問題的求解做了大量的研究.對解析解的研究常用的方法有分離變量法[1]、傅里葉變換和拉普拉斯變換法[2-3].數(shù)值方法包括有限差分方法[4-5]、有限元方法[6]和譜方法[7]等.對于給定的分數(shù)階擴散方程,運用一般數(shù)值方法得到的系數(shù)矩陣是稠密的,用高斯消去法計算量為O(N3),存儲量為O(N2).為了節(jié)約計算時間,王宏等[8]提出一種計算空間分數(shù)階擴散方程的快速算法,在保證計算精度的同時,縮短了計算時間.考慮到反問題的計算時間與正問題計算時間的相關(guān)性,本文在文獻[8]的基礎(chǔ)上,對分數(shù)階擴散方程正問題進行數(shù)值求解.

對于分數(shù)階擴散中的反問題的研究,大多是針對時間分數(shù)階擴散方程中反問題的[9-12],由于空間分數(shù)階模型正問題研究的困難性,相應(yīng)的反問題研究較少,如韋徽[13]等利用最佳攝動量算法的思想對一維空間分數(shù)階擴散方程中的源項系數(shù)進行了反演研究.鄭光輝與魏婷[14]提出了求解分數(shù)階擴散方程反問題的兩種正則化方法.池光勝[15]同樣應(yīng)用最佳攝動量算法探討了一維空間分數(shù)階對流彌散方程中依賴空間變量的源項反問題,對于不同形式的源項函數(shù)進行了有效的數(shù)值模擬.賈現(xiàn)正[16]利用同倫正則化算法對時間空間分數(shù)階對流擴散方程中的微分階數(shù)進行了數(shù)值反演.鄧偉華[17]利用正則化算法對空間分數(shù)階擴散方程的源項進行了反問題的研究.

由于微分階數(shù)是刻劃分數(shù)階擴散的首要指標,擴散系數(shù)描繪空間非均質(zhì)性,源項反映外界力屬性,在實際問題中這些量往往都是是未知的,也是難以通過實驗手段來測量的,因此,本文考慮微分階數(shù)、擴散系數(shù)和源項同時反演.由于其高病態(tài)性,使用最佳攝動量算法無法實現(xiàn)問題的求解.考慮到同倫方法在求解線性代數(shù)方程中的廣泛應(yīng)用,將最佳攝動量算法和同倫方法結(jié)合形成新的迭代算法,可以擴大收斂區(qū)域,且成功運用于擴散方程相關(guān)反問題的研究中[16],本文繼續(xù)運用這種同倫正則化算法.

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

對于l>0,T>0,考慮一維(左側(cè))空間分數(shù)階擴散方程

(1)

式中：c(x,t)為狀態(tài)變量; 1<α<2為空間微分階數(shù);D(x)>0為擴散系數(shù);f(x,t)為源/匯項.當α=2時,方程(1)即為經(jīng)典的整數(shù)階擴散方程.方程(1)中的空間分數(shù)階導數(shù)用Riemann-Liouville的定義[3]：

對于方程(1),給定初始條件為

c(x,0)=g(x) 0≤x≤l,

(2)

邊值條件為

c(0,t)=c(l,t)=0 0

(3)

對該類型正問題的數(shù)值求解,本文在文獻[8]的基礎(chǔ)上,構(gòu)造有效的差分格式.

2.1 隱式差分格式

首先給出正問題求解的一般的隱式差分格式.設(shè)xi=ih,h>0,i=0, 1,…,M;tn=nτ,τ>0,n=0, 1,…,N,其中h和τ分別是空間和時間步長. Riemann-Liouville分數(shù)階導數(shù)采用移位的Grünwald公式離散為

(4)

(I+An+1)Cn+1=Cn+τfn+1;

(5)

(6)

注1文獻[4]給出了差分格式(5)的無條件穩(wěn)定性與收斂性的證明.

式(5)即為一般的隱式差分格式.由于差分格式(5)的系數(shù)矩陣為滿秩，因此運算量為O(M3),每個時間步長需要O(M2)的存儲量.下面給出快速求解的差分格式.

(7)

(8)

根據(jù)上述討論,我們將An+1作如下分解：

(9)

即

對Cn+1作近似分解,則

(10)

將式(10)代入式(5)可得新的差分格式為

(11)

由以上的分析可得，格式(11)將式(5)的系數(shù)矩陣計算存儲量由O(M3)降為O(Mlog22M), 但是格式(11)右側(cè)的運算量和存儲需求均為O(M2).下面將給出差分格式(11)右側(cè)的快速算法.

引理2[19]若Bn=Bn(i,j)=b(j-i),i,j=1,2,…,n是一個循環(huán)矩陣,則

其中：b=(b0,bn-1,…,b2,b1)T是矩陣Bn的第一列；Fn是n×n階傅里葉變換矩陣.

1)IMERG日平均降水量的空間分布連續(xù)性較好,能夠較好地反映我國降水量的分布特點,但局部地區(qū)較CGDPA偏低。IMERG與CGDPA日平均降水量的相關(guān)系數(shù)為0.91,均方根誤差為0.66 mm/d,相對偏差為5.32%。IMERG與CGDPA在東部地區(qū)的相關(guān)系數(shù)為0.94,而西部地區(qū)為0.72。兩者在東部地區(qū)的均方根誤差為0.52 mm/d,西部地區(qū)為0.81 mm/d;東部地區(qū)的相對誤差為-1.62%,西部為-18.6%。IMERG在東部地區(qū)的降水估計精度明顯優(yōu)于西部。

步驟2令w2M-2,L=F2M-2C2M-2,F2M-2C2M-2是C2M-2的傅里葉變換;

步驟3令v2M-2,L=F2M-2b2M-2,L,b2M-2,L是B2M-2,L的第一列;

通過以上分析可知，利用快速算法可將格式(11)右側(cè)運算量由O(M2)降為O(Mlog2M),存儲需求由O(M2)降為O(M).

下面就以一個分數(shù)階擴散方程為例,說明如何利用差分格式(11)進行數(shù)值計算.

2.2 數(shù)值試驗

先看離散點數(shù)對數(shù)值計算的影響.不妨取空間微分階數(shù)α=1.9,應(yīng)用差分格式(11)計算獲得數(shù)值解,表1給出了t=1時空間、時間離散點數(shù)與解的相對誤差,其中M,N分別表示空間、時間離散點數(shù).

由表1可以看出,隨著空間、時間離散點數(shù)的增加,解誤差逐步變小.下面取定M=N=27,考察微分階數(shù)變化對數(shù)值解的影響,計算結(jié)果列于表2.

表1 t=1時空間、時間離散點數(shù)對誤差的影響Tab.1 The influence of the time and space discrete points on the error at t=1

表2 t=1時空間微分階數(shù)對誤差的影響Tab.2 The influence of the space differential order on the error at t=1

從表2可以看出,空間微分階數(shù)對誤差基本沒有影響.下面考察快速算法(Fast Finite Difference, FFD)和一般的有限差分方法(Finite Difference, FD)對數(shù)值解的影響,其中TCPU(單位： s)表示計算時間,計算結(jié)果列于表3.

表3 t=1時FFD與FD對數(shù)值解的影響Tab.3 The influence of FFD and FD on numerical solution at t=1

圖1 α=1.9時的精確解與數(shù)值解(t=1)Fig.1 The comparison between exact solution and numerical solution at α=1.9 and t=1

從表3看出,新的差分格式利用快速算法比一般的有限差分方法除了計算精度略高以外還縮短了計算時間.圖1繪制了α=1.9,在t=1時刻的精確解與數(shù)值解(M=N=27).

3 多參數(shù)聯(lián)合反問題與反演算法

對于分數(shù)階擴散初邊值問題泛定方程(1),當擴散系數(shù)D(x)為常數(shù)D,源項具有形式f(x,t)=a(x)λ(t)時,給定附加數(shù)據(jù)為

c(x,T)=θ(x) 0≤x≤l.

(12)

當α,D,a(x)給定時,稱(1)為正問題,可以通過上節(jié)中的差分方法求解;若α,D,a(x)未知,則初邊值問題(1)～(3)與附加數(shù)據(jù)(12)就構(gòu)成了參數(shù)識別的反演問題.本節(jié)就考慮同時確定微分階數(shù)α,擴散系數(shù)D和源項系數(shù)a(x)的多參數(shù)反演問題.雖然反問題解的唯一性、穩(wěn)定性等理論分析很重要,不過本文將關(guān)注反演算法及數(shù)值反演.

Qn+1=Qn+δQn.

(13)

所以Q的確定問題可以轉(zhuǎn)化為δQn的確定問題,并且δQn可以由下列目標函數(shù)的局部極小值來確定：

這里同倫參數(shù)λ參考人工神經(jīng)網(wǎng)絡(luò)的傳遞函數(shù),即

(14)

其中：N0為初值選取參數(shù);β為下降速率參數(shù);λ(N)為第N個迭代步時正則參數(shù)的取值.這里有兩個可調(diào)參數(shù),β和N0.其中β用于調(diào)整正則化參數(shù)的下降速率,一般該參數(shù)取0.5左右的數(shù)值.參數(shù)N0一般可選為0～5之間的整數(shù),主要保證計算初期的迭代穩(wěn)定進行.因為δQn是一個微小的擾動量,把c(x,T;Qn+δQn)利用多元泰勒公式展開得到

于是有

進一步將區(qū)域[0,l]離散為0

求解minF(δQn)相當于求解規(guī)范方程

((1-λ(N))HTH+λ(N)I)δQn=(1-λ(N))HT(V-U),

(15)

其中：

H=(hms)M×Nm=1,2,…,M,s=1,2,…,N;

U=(c(x1,T;Qn),c(x2,T;Qn),…,c(xM,T;Qn))T;

V=(θ(x1),θ(x2),…,θ(xM))T；

其中：z是數(shù)值微分步長;I是單位矩陣.這樣,適當選擇同倫參數(shù)λ(N),可從式(15)中解得對應(yīng)于xm的攝動量δQn,進而再通過迭代過程(13)可逐步求得原來反問題的解.

在給定真解QR=(αR,DR,a(x)R)的情況下,反演算法實施步驟為：

步驟1給定數(shù)值微分步長z及迭代控制收斂精度eps,由真解代入正問題的差分格式進行計算,得到附加數(shù)據(jù)向量V;

步驟2給定初始迭代Qn(n=0,1,…),求得U與矩陣H;

步驟3按式(15)計算最佳攝動量δQn,其中同倫參數(shù)按式(14)選取;

步驟4判斷是否滿足‖δQn‖≤eps,若成立,則Qn+1=Qn+δQn即為所求的反演解,算法終止;否則由Qn+1代替Qn,再轉(zhuǎn)到步驟2繼續(xù)進行.

4 數(shù)值反演

本節(jié)對于微分階數(shù)、擴散系數(shù)和源項多參數(shù)聯(lián)合反問題進行數(shù)值反演.如前所述,每個算例中給定真解,并利用其得到終值時刻的附加數(shù)據(jù),然后應(yīng)用同倫正則化算法進行反演重建.以下若無特別說明,正問題計算中取區(qū)域長度l=1,終值時刻T=1,擴散系數(shù)D=1,且取零邊界條件,初始函數(shù)g(x)=x2(1-x),離散點數(shù)M=N=16且λ(t)=e-t;反演算法中取數(shù)值微分步長z=0.01,控制迭代收斂精度eps=10-6,同倫參數(shù)(14)中取j0=5,β=0.3.

4.1 精確數(shù)據(jù)下的反演

本小節(jié)先考慮精確數(shù)據(jù)下的情形,取方程中的微分階數(shù)α=1.9,源項系數(shù)a(x)=1-x,基向量{φi(x)=xi-1,i=1,2,3,4,則真解為QR=(1.9,1,1,-1,0,0),記Q*為數(shù)值反演解,則解的誤差表示為err=‖Q*-QR‖/‖QR‖,j是迭代次數(shù).

考察不同的數(shù)值微分步長對算法的影響,取初始迭代Q0=(1.1,0,0,0,0,0),其余參數(shù)值不變,計算結(jié)果見表4.

表4 α=1.9時,數(shù)值微分步長對反演算法的影響Tab.4 The influence of numerical differential step on inversion algorithm with α=1.9

由表4可以看出,數(shù)值微分步長對反演結(jié)果有一定的影響,不宜取的過小.再取定z=0.05,考察空間微分階數(shù)變化對反演算法的影響,計算結(jié)果列于表5,其中α表示空間微分階數(shù),Q*,err及j表示意義與表4中相同.

表5 空間微分階數(shù)對反演算法的影響Tab.5 The influence of space differential order on inversion algorithm

由表5可以看出,空間分數(shù)階數(shù)對反演結(jié)果基本沒有影響.考察下降速率參數(shù)對反演結(jié)果的影響,取α=1.9,z=0.05,其余參數(shù)不變,計算結(jié)果見表6.

表6 下降速率參數(shù)對反演算法的影響Tab.6 The influence of descent rate parameters on inversion algorithm

由表6可以看出,下降速率參數(shù)對反演結(jié)果有一定的影響,不宜取得過大,當β≥0.58時,無法計算.考察源項形式對反演結(jié)果的影響,取α=1.9,z=0.1,β=0.5,其余參數(shù)不變,計算結(jié)果見表7.

表7 源項形式對反演算法的影響Tab.7 The influence of source term on inversion algorithm

由表7可以看出,源項形式對反演結(jié)果基本沒有影響.考察離散節(jié)點對反演結(jié)果的影響,取α=1.9,z=0.1,β=0.5,a(x)=1-x,其余參數(shù)不變,計算結(jié)果見表8，TCPU(單位： s)表示一次反演計算所用的CPU時間.

表8 離散節(jié)點數(shù)對反演算法的影響Tab.8 The influence of discrete point number on inversion algorithm

由表8可以看出,離散節(jié)點數(shù)的增加對反演解精度影響不大,但是對反演計算量影響較大.隨著網(wǎng)格的加密,迭代次數(shù)稍微減少,但每次反演計算耗用的時間明顯增加.當離散節(jié)點超過27時,計算失敗.由此可見,離散節(jié)點不宜選取過多.

4.2 擾動數(shù)據(jù)下的反演

實際問題中,附加數(shù)據(jù)往往帶有某種誤差,對于擾動數(shù)據(jù)實施反演算法是反問題數(shù)值方法研究的重要方面.帶有擾動數(shù)據(jù)的附加數(shù)據(jù)表示為

θε(x)=θ(x)(1+ξε),

(28)

從表9可以看出,隨著數(shù)據(jù)擾動水平的減小,反演解與精確解誤差逐漸變小,表明了反演算法的穩(wěn)定性.在ε=0.1%的水平下,考察快速算法(FFD)和有限差分方法(FD)對數(shù)值反演的影響,其中TCPU(單位： s)表示反演計算所用的時間,計算結(jié)果見表10.

表9 不同擾動水平下的反演結(jié)果Tab.9 The inversion results with different noisy levels

表10 給定擾動水平下FFD與FD對反演結(jié)果的影響Tab.10 The influence of FFD and FD on the inversion results with a constant noisy level

從表10可以看出,利用新的差分格式對反演解的誤差并沒有太大的影響,但是可以有效地節(jié)約計算時間.

5 結(jié) 語

本文探討了(左側(cè))空間分數(shù)階擴散方程正問題求解的差分格式,應(yīng)用分數(shù)階導數(shù)離散系數(shù)和Toeplitz矩陣的性質(zhì)及快速傅里葉變換完成了由有限差分方法到快速求解的過渡.進一步,對于空間微分階數(shù)、擴散系數(shù)和源項的多參數(shù)聯(lián)合反演問題,應(yīng)用一種同倫正則化算法給出了精確數(shù)據(jù)與擾動數(shù)據(jù)情形下的數(shù)值反演.結(jié)果表明： 在恰當?shù)倪x取數(shù)值微分步長和下降速率參數(shù)的情況下,同倫正則化算法可以有效地處理對于空間分數(shù)階反常擴散中的多參數(shù)聯(lián)合反演問題.而高維的情形和多參數(shù)聯(lián)合反演的存在唯一性分析,由于問題的高病態(tài)性,將是我們今后的主要工作.

[1] 池光勝.分數(shù)微分對流彌散方程反問題的研究 [D].淄博： 山東理工大學, 2010.

[2] MILL J V. Theory and applications of fractional differential equations [M]. New York,USA： Elsevier Science Ltd, 2006.

[3] PODLUBNY I. Fractional differential equations[M]. New York, USA： Academic Press, 1999.

[4] MEERSCHAERT M M, TADJERAN C. Finite difference approximations for fractional advection-dispersion flow equations [J].JournalofComputationalandAppliedMathematics, 2004,172(1)： 65-77.

[5] MEERSCHAERT M M, TADJERAN C. Finite difference approximations for two-dimensional fractional dispersion equation [J].JournalofComputationalPhysics, 2006,211(1)： 249-261.

[6] DENG W. Finite element method for the space and time fractional Fokker-Planck equation [J].SIAMJournalonNumericalAnalysis, 2008,47(1)： 204-226.

[7] LI X, XU C. Existence and uniqueness of the weak solution of the space-time fractional diffusion equation and a spectral method approximation [J].CommunicationsinComputionalPhysics, 2010,8(5)： 1016-1051.

[8] WANG H, WANG K, SIRCAR T. AdirectO(Nlog22N) finite difference method for fractional diffusion equations [J].JournalofComputationalPhysics, 2010,229(21)： 8095-8104.

[9] CHENG J, NAKAGAWA J, YAMAMOTO M,etal. Uniqueness in an inverse problem for a one-dimensional fractional diffusion equation [J].InverseProblems, 2009,25(11)： 115002.

[10] LI G S, GU W J, JIA X Z. Numerical inversions for space-dependent diffusion coefficient in the time fractional diffusion equation [J].JournalofInverseandIll-PosedProblems, 2012,20(3)： 339-366.

[11] LI G S, ZHANG D L, JIA X Z,etal. Simultaneous inversion for the space-dependent diffusion coefficient and the fractional order in the time-factional diffusion equation [J].Inverseproblems, 2013,29(6)： 065014.

[12] XU X, CHENG J, YAMAMOTO M. Carleman estimate for a fractional diffusion equation with half order and application [J].ApplicableAnalysis, 2011,90(9)： 1355-1371.

[13] WEI H, CHEN W. A coupled method for inverse source problem of spatial fractional anomalous diffusion equations [J].InverseProbleminScienceandEngineering, 2010,18(7)： 945-956.

[14] ZHENG G H, WEI T. Two regularization methods for solving a Riesz-Feller space-fractional backward diffusion equation [J].Inverseproblem, 2010,26(11)： 115071.

[15] CHI G S, LI G S, JIA X Z. Numerical inversions of a source term in the FADE with Dirichlet boundary condition using final ovservations [J].ComputersandMathematicswithapplications, 2011,62(4)： 1619-1626.

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

[17] TIAN W Y, LI C, DENG W H,etal. Regularization methods for unknown source in space fractional diffusion equation [J].MathematicsandComputersinsimulation, 2012,85(3)： 45-56.

[18] CHEN Z Q, MEERSCHAERT M M, NANE E,etal. Space-time fractional diffusion on bounded domains [J].JournalofMathematicalAnalysisandApplications, 2012,393(2)： 479-488.

[19] GRAY R M. Toeplitz and circulant matrices： A review [J].FoundationsandTrendsinCommunicationsandInformationTheory, 2006,2(3)： 155-239.

NumericalInversionsofJointParametersinFractionalDiffusionEquation

CHIGuangsheng1,LIGongsheng1,2

(1.CollegeofScience,InnerMongoliaUniversityofTechnology,Hohhot010051,China;2.SchoolofScience,ShandongUniversityofTechnology,Zibo255049,China)

A joint parameters inversion problem of determining the space fractional diffusion order, diffusion coefficient and source term simultaneously in 1-dimentional space fractional diffusion equation is studied. Using the property of the discrete coefficient of fractional derivative and the fast Fourier transformation, a new finite difference scheme for solving the forward problem is given. Numerical inversions with accurate data and noisy data for the joint parameters inversion problem using homotopy regularization algorithm is discussed. The numerical results show that the homotopy regularization algorithm is efficient for the joint parameters inversion in the space fractional diffusion.

fractional diffusion equation; finite difference scheme; joint parameters inversion; numerical inversion

0427-7104(2017)05-0767-09

2016-03-28

國家自然科學基金(11371231)

池光勝(1985—),男,博士研究生, E-mail： chiguangsheng007@163.com.

O242.1

A