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

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

?

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

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

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

摘要：對(duì)于一類帶有三個(gè)時(shí)間分?jǐn)?shù)階的一維反常擴(kuò)散問題,基于Caputo意義下時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的離散,給出了一個(gè)有限差分求解格式,并利用分離變量法及Laplace變換得到該問題的解析解.進(jìn)一步應(yīng)用同倫正則化算法,根據(jù)內(nèi)點(diǎn)處的濃度觀測(cè)數(shù)據(jù)對(duì)確定微分階數(shù)的反問題進(jìn)行數(shù)值反演,并討論時(shí)間-空間步長(zhǎng)及數(shù)據(jù)擾動(dòng)等因素對(duì)反演算法的影響.

關(guān)鍵詞：時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)； 含三個(gè)時(shí)間分?jǐn)?shù)階的擴(kuò)散； 反問題； 同倫正則化算法； 數(shù)值反演

在過去的幾十年,分?jǐn)?shù)階偏微分方程在反常擴(kuò)散現(xiàn)象模擬等領(lǐng)域開始發(fā)揮重要作用,其對(duì)復(fù)雜系統(tǒng)的描述具有建模簡(jiǎn)單、參數(shù)物理意義清楚、描述準(zhǔn)確等優(yōu)勢(shì).陳文案[1]較詳細(xì)地描述了分?jǐn)?shù)階微積分理論在復(fù)雜力學(xué)行為建模和數(shù)值模擬方面的研究成果，強(qiáng)調(diào)分?jǐn)?shù)階微積分建模的物理和力學(xué)背景和概念.分?jǐn)?shù)階反常擴(kuò)散方程是由經(jīng)典擴(kuò)散方程推廣而得到的,按照其含有的時(shí)間-空間導(dǎo)數(shù)的形式可分為三類：分別是時(shí)間分?jǐn)?shù)階擴(kuò)散、空間分?jǐn)?shù)階擴(kuò)散和空間-時(shí)間分?jǐn)?shù)階擴(kuò)散.當(dāng)通常的整數(shù)階擴(kuò)散方程中關(guān)于時(shí)間變量的1階導(dǎo)數(shù)換為α(0<α≤1)時(shí),得到時(shí)間分?jǐn)?shù)階擴(kuò)散方程；當(dāng)整數(shù)階擴(kuò)散方程中關(guān)于空間變量的2階導(dǎo)數(shù)換為α(0<α≤2)時(shí),得到空間分?jǐn)?shù)階擴(kuò)散方程；而當(dāng)整數(shù)階擴(kuò)散方程中關(guān)于時(shí)間與空間變量的整數(shù)階導(dǎo)數(shù)都換為相應(yīng)的分?jǐn)?shù)階導(dǎo)數(shù)時(shí),就得到空間-時(shí)間分?jǐn)?shù)階擴(kuò)散方程.

對(duì)于時(shí)間分?jǐn)?shù)階擴(kuò)散方程正問題的研究,劉發(fā)旺等人做了很多研究.于強(qiáng)與劉發(fā)旺[2]給出了時(shí)間分?jǐn)?shù)階反應(yīng)-擴(kuò)散方程的隱式差分近似；沈淑君與劉發(fā)旺[3]推導(dǎo)出了時(shí)間分?jǐn)?shù)階擴(kuò)散方程初邊值問題的數(shù)值解,并提出了方程的顯式守恒差分格式,證明了此格式的穩(wěn)定性和收斂性,然后將得到的結(jié)果推廣到時(shí)間分?jǐn)?shù)階對(duì)流-擴(kuò)散方程.劉法旺和Meerschaert[4]等人,給出了含有多個(gè)時(shí)間分?jǐn)?shù)階的擴(kuò)散-波動(dòng)方程的數(shù)值計(jì)算方法,并討論了隱式差分求解格式.Luchko[5]與Liu等[6]分別研究了多個(gè)時(shí)間分?jǐn)?shù)階的對(duì)流-擴(kuò)散方程,得到了廣義Mittag-Leffler函數(shù)形式的解析解.

對(duì)于分?jǐn)?shù)階擴(kuò)散中的反問題研究,已引起了大家的關(guān)注.程晉等[7]研究了一維時(shí)間分?jǐn)?shù)階擴(kuò)散中由單個(gè)邊界點(diǎn)處的連續(xù)測(cè)量值同時(shí)確定微分階數(shù)與空間依賴擴(kuò)散系數(shù)的反問題,應(yīng)用特征展開方法證明了該反問題解的唯一性.劉繼軍和Yamamoto[8]應(yīng)用準(zhǔn)可逆性正則化(quasi-reversibilityregularization)方法探討了一個(gè)時(shí)間分?jǐn)?shù)階擴(kuò)散中由終值數(shù)據(jù)恢復(fù)初始狀態(tài)的倒向反問題.鄭光輝與魏婷等[9-10]應(yīng)用Fourier譜正則化方法研究了帶狀區(qū)域一維時(shí)間分?jǐn)?shù)階對(duì)流擴(kuò)散的Cauchy問題與空間Riesz分?jǐn)?shù)階擴(kuò)散的倒向反問題.Sakamoto與Yamamoto[11]利用特征分解方法探討了時(shí)間分?jǐn)?shù)階擴(kuò)散-波動(dòng)方程倒向反問題的適定性.徐翔、程晉與Yamamoto等[12-13]研究了具有1/2階分?jǐn)?shù)階導(dǎo)數(shù)的時(shí)間分?jǐn)?shù)階擴(kuò)散方程的Carleman估計(jì)方法,這是構(gòu)建分?jǐn)?shù)階擴(kuò)散相關(guān)反問題條件穩(wěn)定性的新思路.熊向團(tuán)等[14]構(gòu)建了時(shí)間分?jǐn)?shù)階二維熱傳導(dǎo)方程的反向熱傳導(dǎo)問題的一種條件穩(wěn)定性,李功勝等[15]利用最佳攝動(dòng)量算法對(duì)于一維時(shí)間分?jǐn)?shù)階擴(kuò)散方程的空間依賴的擴(kuò)散系數(shù)進(jìn)行了有效的數(shù)值反演,魏婷等[16]基于邊界元離散提出一種正則化方法,研究了一維時(shí)間分?jǐn)?shù)階擴(kuò)散方程中時(shí)間依賴源項(xiàng)的重建問題.

本文將探討含有三個(gè)時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的一維反常擴(kuò)散方程求解及確定微分階數(shù)的反演問題.首先利用有限差分法給出正問題數(shù)值求解的差分格式,同時(shí)利用分離變量法及拉普拉斯變換得到正問題的解析解表達(dá)式,并且給出正問題的數(shù)值算例.進(jìn)一步,應(yīng)用同倫正則化算法對(duì)反常擴(kuò)散方程中的三個(gè)時(shí)間分?jǐn)?shù)階微分階數(shù)進(jìn)行數(shù)值反演,給出數(shù)值反演算例,并討論時(shí)間空間步長(zhǎng)及數(shù)據(jù)擾動(dòng)水平等因素對(duì)反演算法的影響.

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

考慮區(qū)域Ω=(0,l)×(0,T)上含三個(gè)時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的一維反常擴(kuò)散方程

(1)

其中(x,t)∈Ω;D>0是擴(kuò)散系數(shù),α∈(0,1)和αj(j=1,2)∈(0,1)是時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的微分階數(shù)，且0<α2<α1<α<1.對(duì)于方程(1),給定非零初值條件和第一類零邊值條件:

u(x,0)=f(x),u(0,t)=u(l,t)=0

(2)

這樣,由方程(1)及初邊值條件(2)就構(gòu)成含三個(gè)時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的一維反常擴(kuò)散方程的正問題.下面給出這個(gè)正問題求解的一個(gè)有限差分格式，并運(yùn)用分離變量法推導(dǎo)其解析解的表達(dá)式.

1.1正問題求解的差分格式

(3)

(4)

對(duì)于方程中的整數(shù)階導(dǎo)數(shù)項(xiàng),按照通常的差分離散方法,即有

(5)

(6)

1.2正問題的解析解

本小節(jié)利用分離變量法及Laplace變換推導(dǎo)正問題的解析解.

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

(7)

(8)

(9)

(10)

(11)

(12)

應(yīng)用Laplace變換,得到

(13)

代入(12)式,可得

記c=Dλn為常數(shù),可得

(14)

(15)

于是求得解析解為

(16)

1.3數(shù)值試驗(yàn)

利用上一節(jié)的有限差分法進(jìn)行數(shù)值計(jì)算,取l=π,T=1,且離散點(diǎn)數(shù)M=100,N=100，微分階數(shù)α=0.8,α1=0.5,α2=0.3,擴(kuò)散系數(shù)D=1,系數(shù)r1=0.5,r2=0.5.另外,取初始函數(shù)u(x,0)=f(x)=sinx.此時(shí)解析解可表示為

(17)

對(duì)上述解析解,級(jí)數(shù)取前50項(xiàng)截?cái)嘤?jì)算,記數(shù)值解u*(x,t),相對(duì)誤差表示為Err=‖u(x,t)-u*(x,t)‖2/‖(u,t)‖2.數(shù)值結(jié)果分別列于表1～表3和圖1.

表2 時(shí)間步長(zhǎng)對(duì)解誤差的影響(t=0.5)

表3　不同微分階數(shù)對(duì)正問題求解的影響(t=0.5)

圖1 　t=0.5時(shí)的數(shù)值解與解析解

從表1的計(jì)算結(jié)果可以看出,空間步長(zhǎng)對(duì)誤差的影響不大;由表2的計(jì)算結(jié)果可以看出,隨著時(shí)間步長(zhǎng)的減小,誤差逐漸減小;由表1、表2及圖1的結(jié)果可以看到，正問題的數(shù)值解和解析解吻合的較好.從表3看出,時(shí)間微分階數(shù)對(duì)正問題求解具有一定的影響,微分階數(shù)越接近,解誤差有所減小.

2確定微分階數(shù)的反問題與反演算法

2.1反問題的提出

假設(shè)方程(1)中的時(shí)間微分階數(shù)未知,那么為了確定各個(gè)微分階數(shù)的值,需要補(bǔ)充關(guān)于正問題解的部分條件信息,并聯(lián)合正問題(1)～(2)形成一個(gè)確定微分階數(shù)的反問題.本文給定內(nèi)點(diǎn)x=x0,時(shí)刻tj(j=1,2,…，N)的觀測(cè)值為附加數(shù)據(jù),記為u(x0,tj)=θ(tj),則可定義附加數(shù)據(jù)向量V：

(18)

這樣,由附加數(shù)據(jù)(18)聯(lián)合正問題(1)～(2)構(gòu)成了一個(gè)確定微分階數(shù)(α,α1,α2)的數(shù)值反演問題.下面給出同倫正則化算法,并對(duì)上述確定微分階數(shù)的反問題進(jìn)行數(shù)值反演模擬.

2.2同倫正則化算法

記a=(α,α1,α2),利用上一節(jié)的差分方法求解正問題可得其解,并在x=x0處賦值,記之為u(x0,t;a),這可以看作對(duì)應(yīng)于輸入數(shù)據(jù)a=(α,α1,α2)的一個(gè)計(jì)算輸出.

min‖U(a)-V‖2

(19)

為了克服病態(tài)性,應(yīng)用同倫思想可將上式轉(zhuǎn)化為如下極小值問題：

min{(1-λ)‖U(a-V‖2+λ‖a‖2}

(20)

(21)

根據(jù)同倫正則化思想,上述極小問題(20)的求解又轉(zhuǎn)化為對(duì)于給定的an,通過求解最佳攝動(dòng)量δan進(jìn)而確定an+1的一種迭代算法：an+1=an+δan,n=0,1,2，…

(22)

(23)

將u(x0,tj;an+δan)在an處作泰勒展開,并略去高階項(xiàng),可以得到

則目標(biāo)函數(shù)F(δan)可近似表為

(24)

(25)

式中ω為數(shù)值微分步長(zhǎng).不難驗(yàn)證泛函極小問題等價(jià)于求解規(guī)范方程：

((1-λ)GTG+λI)δan=(1-λ)GT(V-aU)，

(26)

則得到每一步迭代的最優(yōu)攝動(dòng)量δan的計(jì)算公式：

δan=((1-λ)GTG+λI)-1(1-λ)GT(V-U)

(27)

對(duì)于給定的精度eps,判斷‖δan‖≤eps是否成立,若滿足則δan即為所求;否則由式(22)得到an+1,再通過規(guī)范方程繼續(xù)求解.

3數(shù)值反演

本節(jié)應(yīng)用同倫正則化算法對(duì)確定微分階數(shù)的反問題(1)～(2)及(18)進(jìn)行數(shù)值反演.取初始函數(shù)f(x)=sinx,以下若無特殊說明,正問題計(jì)算中均取l=π,擴(kuò)散系數(shù)D=1,終值時(shí)刻T=1,且M=100,N=100.附加數(shù)據(jù)取為在x0=0.5處的觀測(cè)數(shù)據(jù)，數(shù)值微分步長(zhǎng)ω=0.01,停止準(zhǔn)則eps=1e-6,預(yù)估迭代次數(shù)N0=5,β=0.8,設(shè)式(1)中r1=r2=0.5,取α=0.9,α1=0.7,α2=0.5，則微分階數(shù)真值a=(0.9，0.7，0.5)，ainv表示反演解,Err=‖ainv-a‖2/‖a‖2表示反演解與真解的誤差.

3.1附加數(shù)據(jù)不帶擾動(dòng)的情形

(1)時(shí)間步長(zhǎng)對(duì)反演結(jié)果的影響

令初始迭代值a0=(0.3,0.2,0.1),反演計(jì)算結(jié)果列于表4,n表示迭代次數(shù).

表4　時(shí)間步長(zhǎng)對(duì)反演結(jié)果的影響

表5　空間步長(zhǎng)對(duì)反演結(jié)果的影響

由表4可以看出,時(shí)間步長(zhǎng)對(duì)反演結(jié)果影響不大,只是隨時(shí)間步長(zhǎng)的變小迭代的次數(shù)稍微減少.

(2)空間步長(zhǎng)對(duì)反演結(jié)果的影響

取初始迭代值a0=(0.3,0.2,0.1),反演計(jì)算結(jié)果列于表5,n表示迭代次數(shù).

由表5可以看出,空間步長(zhǎng)對(duì)反演結(jié)果的影響很小.

(3)微分階數(shù)對(duì)反演結(jié)果的影響

令初始迭代值a0=(0.3,0.2,0.1),考察時(shí)間分?jǐn)?shù)階α,α1,α2的不同取值對(duì)反演結(jié)果的影響,計(jì)算結(jié)果見表6,n表示迭代次數(shù).

表6　微分階數(shù)對(duì)反演結(jié)果的影響

通過表6可以看出微分階數(shù)對(duì)反演結(jié)果有一定的影響,隨著微分階數(shù)的接近誤差有所增加,并且迭代步數(shù)也略微增加.

3.2附加數(shù)據(jù)有擾動(dòng)的情形

(1)擾動(dòng)水平對(duì)反演結(jié)果的影響仍取微分階數(shù)真值a=(0.9,0.7,0.5),數(shù)值微分步長(zhǎng)τ=0.001,分別取擾動(dòng)水平δ=1%,0.1%,0.05%,0.01%,20次反演計(jì)算的平均值列于表7.

通過表7可以看出隨著數(shù)據(jù)擾動(dòng)水平的減小,反演結(jié)果精度越來越高.

(2)計(jì)算次數(shù)對(duì)反演結(jié)果的影響

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

表8 給定擾動(dòng)水平下計(jì)算次數(shù)對(duì)反演結(jié)果的影響

設(shè)擾動(dòng)水平δ=1%,考察計(jì)算次數(shù)對(duì)反演結(jié)果的影響,結(jié)果見表8,其中p表示計(jì)算次數(shù).

由表8可以看出,對(duì)于給定的擾動(dòng)數(shù)據(jù),計(jì)算次數(shù)對(duì)反演結(jié)果的影響不大.

4結(jié)束語

本文探討了含有三個(gè)時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的一維反常擴(kuò)散方程求解及確定微分階數(shù)的反問題.正問題的數(shù)值模擬結(jié)果表明,數(shù)值解與解析解吻合較好,所建立的有限差分格式是有效的.對(duì)于確定三個(gè)微分階數(shù)的反問題,雖然理論上還沒有解決反演的穩(wěn)定性問題,但應(yīng)用同倫正則化算法可以實(shí)現(xiàn)對(duì)多個(gè)微分階數(shù)的數(shù)值反演.通過反演計(jì)算結(jié)果發(fā)現(xiàn),當(dāng)數(shù)據(jù)擾動(dòng)水平較大時(shí),解誤差變得較大.這說明含多個(gè)時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的反常擴(kuò)散方程的微分階數(shù)反問題具有較高的病態(tài)性,其他更有效的反演算法是需要進(jìn)一步研究的課題.

參考文獻(xiàn):

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

[2] 于強(qiáng),劉發(fā)旺.時(shí)間分?jǐn)?shù)階反應(yīng)-擴(kuò)散方程的隱式差分近似[J],廈門大學(xué)學(xué)報(bào):自然科學(xué)版, 2006,45(3): 315-319.

[3] 沈淑君.分?jǐn)?shù)階對(duì)流-擴(kuò)散方程的基本解和數(shù)值方法[D].廈門：廈門大學(xué), 2008.

[4] Liu F, Meerschaert M M, McGough R J,etal. Numerical methods for solving the multi-term time-fractional wave-diffusion equations[J]. Fractional Calculus and Applied Analysis, 2013, 16: 9-25.

[5] Yury L.Initial-boundary-value problems for the generalized multi-term time-fractional diffusion equation[J].Journal of Mathematical Analysis and Applications, 2011, 374:538-548.

[6] Jiang H, Liu F, Turner I,etal. Burrage,Analytical solutions for the multi-term time-fractional diffusion-wave/diffusion equations in a finite domain [J]. Computers and Mathematics with Applications,2012, 64:3 377-3 388.

[7] Cheng J, Nakagawa J, Yamamoto M,etal. Uniqueness in an inverse problem for a one-dimensional fractional diffusion equation [J].Inverse Problems, 2009,25:115 002-115 017.

[8] Liu J J,Yamamoto M,A backward problem for the time-fractional diffusion equation [J], Applicable Analysis, 2010,89:1 769-1 788.

[9] Zheng G H,Wei T.Spetral regularization method for a Cauchy problem of the time fractional advection-dispersion equation [J]. Journal of Computational and Applied Mathematics, 2010, 233:2 631-2 640.

[10] Zheng G H, Wei T.Two regularization methods for solving a Riesz-Feller space-fractional backward diffusion equation [J]. Inverse Problems, 2010, 26:115 017-115 038.

[11] Sakamoto K,Yamamoto M. Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems [J]. Journal of Mathematical Analysis and Applications, 2011,382:426-447.

[12] Xu X, Cheng J,Yamamoto M.Carleman estimate for a fractional diffusion equation with half order and application [J].Applicable Analysis, 2011,90:1 355-1 371.

[13] Yamamoto M,Zhang Y.Conditional stability in determining a zeroth-order coefficient in a half-order fractional diffusion equation by a Carleman estimate [J].Inverse Problems, 2012,28(10): 105 010-105 019.

[14] Xiong X T, Zhou Q,Hon Y C.An inverse problem for fractional diffusion equation in 2-dimensional case: Stability analysis and regularization [J].Journal of Mathematical Analysis and Applications, 2012,393,185-199.

[15] Li G S, Gu W J, Jia X Z. Numerical inversions for space-dependent diffusion coefficient in the time fractional diffusion equation [J].Journal of Inverse and Ill-Posed Problems, 2012, 20(3): 339-366.

[16] Wei T,Zhang Z Q.Reconstruction of a time-dependent source term in a time-fractional diffusion equation [J]. Engineering Analysis with Boundary Elements, 2013,37:23-31.

[17] Igor Podlubny, Fractional Differential Equations [M].San Diego:Academic press,1999.

(編輯：劉寶江)

Thesolutiontothreetime-fractionalanomalousdiffusion

equationandnumericalinversionofthefractionalorders

SUNChun-long,LIGong-sheng,JIAXian-zheng,DUDian-hu

(SchoolofScience,ShandongUniversityofTechnology,Zibo255049,China)

Abstract:Afinitedifferenceschemeisintroducedtosolvethe1-Dthree-termtime-fractionalanomalousdiffusionequationbasedonCaputo'sdiscretizationtothetimefractionalderivatives.UsingthemethodofseparationofvariablesandLaplacetransform,theanalyticalsolutionoftheforwardproblemisobtained.Furthermore,thehomotopyregularizationalgorithmisappliedtodeterminethethreetime-fractionalordersbytheadditionalmeasurementsataninteriorpointinthedomain.Numericalinversionsareperformedtodemonstrateeffectivenessoftheproposedalgorithm,andinfluencesofthetime-spacemeshgridsandthedatanoisesontheinversionalgorithmarediscussed.

Keywords:timefractionalderivative;threetime-fractionaldiffusion;inverseproblem;homotopyregularizationalgorithm;numericalinversion

中圖分類號(hào)：O175

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

文章編號(hào)：1672-6197(2015)03-0001-07

作者簡(jiǎn)介：孫春龍，男，sunchunlong527@163.com; 通信作者： 李功勝，男，ligs@sdut.edu.cn

基金項(xiàng)目：國(guó)家自然科學(xué)基金資助項(xiàng)目( 11071148, 11371231); 山東省自然科學(xué)基金資助項(xiàng)目(ZR2011AQ014)

收稿日期：2014-09-20