彭 鑫，周曉軍

(貴州師范大學 數(shù)學科學學院，貴州 貴陽 550025)

0 引言

與整數(shù)階導數(shù)相比，分數(shù)階微積分具有描述物質(zhì)記憶功能和遺傳效應的特征。分數(shù)階微積分在電化學過程、色噪聲、控制理論、流體力學、混沌、生物、工程等領(lǐng)域有諸多應用。例如粘彈性材料和波的傳播[1-3]、無序材料中的非布朗輸運現(xiàn)象[4]、非牛頓復雜流體和流變性[5-9]、生物組織中的多尺度模式[10-14]、混沌分形和自動控制[15-16]。

近些年來，多孔、流變、粘彈性和分形介質(zhì)中的反常擴散現(xiàn)象引起了人們的廣泛關(guān)注。對于粒子擴散過程中的隨機運動，經(jīng)典的Fick定律認為粒子只能在單位時間步長內(nèi)跳到相鄰的位置。然而，在異質(zhì)介質(zhì)中的許多實驗觀測中，粒子不僅可以跳躍到相鄰的位置，而且能以較低的概率跳到更遠的位置，這一非局部特征可以利用分數(shù)階導數(shù)來刻畫。分數(shù)階微分算子憑借其特有的全局性，可以更好地描述這些反常擴散現(xiàn)象。例如，針對磁流體動力學模型，Tassaddiq[17]利用分數(shù)階導數(shù)構(gòu)造了一個帶非奇異核的非局部算子來進行研究，實驗數(shù)據(jù)表明，所得到的結(jié)果更符合實際。Sun等[18]給出了粘彈性復雜介質(zhì)的分數(shù)階Langevin方程及其反常擴散。在聲波傳播實驗中，分數(shù)波模型與實驗結(jié)果得到了很好的吻合。Xu 等[19]研究了混凝土狀顆粒材料的輸運特性。在微觀尺度上，這種材料固有的各向異性和異質(zhì)性不符合經(jīng)典的菲克定律。這表明了分數(shù)階階數(shù)與材料結(jié)構(gòu)的異質(zhì)性之間有著重要聯(lián)系。

在研究各種反常擴散現(xiàn)象過程中，分數(shù)階微分算子的導數(shù)是未知的，只有一些實驗觀測數(shù)據(jù)是可用的。對于構(gòu)建分數(shù)階擴散模型而言，關(guān)鍵在于確定其分數(shù)階微分算子的階數(shù)。因此，確定模型的階數(shù)這一反問題是很有意義的。例如，文獻[20]利用實際觀測數(shù)據(jù)來估計簡單反常擴散模型中分數(shù)階導數(shù)的階數(shù)。

本文研究了時間分數(shù)階擴散方程中分數(shù)導數(shù)階的估計問題。首先，定義了一個帶自然對數(shù)核的Caputo分數(shù)階導數(shù)算子，推導出了時間分數(shù)階擴散方程的分數(shù)階導數(shù)所滿足的方程，稱之為時間分數(shù)階擴散的伴隨方程。其次，我們對時間離散構(gòu)造有限差分格式和弱形式，再對弱形式中的半離散解進行Legendre多項式逼近得到全離散格式。然后，在源項已知這一情況下，我們利用序列最小優(yōu)化法對分數(shù)階導數(shù)進行了估計。最后，數(shù)值實驗結(jié)果表明，推導出的伴隨方程是正確的，迭代序列是收斂的。

1 預備知識

定義1[21]α階Caputo分數(shù)階導數(shù)：

定義2 帶自然對數(shù)核的積分導數(shù)算子：

令T>0,Λ∈(0,1),考慮如下時間分數(shù)階擴散方程

(1)

在接下來的討論中，我們假定方程(1)存在唯一光滑解。

定理1 令u是(1)的解，則uα是如下時間分數(shù)階擴散伴隨方程的解

(2)

證明設u是方程(1)的解，則

對上式關(guān)于α求導，即

先對等式左端第一項進行求導，得

同理，對方程(1)左端第二項和右端源項關(guān)于α求導，得

同理，對方程(1)的初、邊值條件關(guān)于α求導變?yōu)?/p>

uα(x,0)=gα(x),x∈Λ，

uα(-1,t)=uα(1,t)=0, 0

我們考慮源項f已知, 但分數(shù)階導數(shù)的階數(shù)α未知。為了估計未知參數(shù)α，設問題(1)的解u的測量值為z。定義二次誤差函數(shù)如下

于是，最小化問題可寫為: 尋找合適參數(shù)α*滿足E(α*)=inf{E(α)，0<α<1}，眾所周知其必要條件是:E′(α)=0，我們構(gòu)造一個最小優(yōu)化序列：

(a)初值估計: 選擇α0∈(0,1)；

(b)問題求解: 通過求解方程(1)和(2)，得出u和uα；

(d)參數(shù)更新:αk+1=αk-ρΔαk，其中ρ>0是一個足夠小的數(shù)。

2 時間離散:有限差分格式

2.1 方程(1)的有限差分格式及弱形式

為了簡化記號且不失一般性，在格式建立過程中考慮f≡0。首先考慮均勻離散(tk)0≤k≤K，時間步長h=T/K。將分數(shù)階導數(shù)寫成如下的式子: 對所有的0≤k≤K-1，

則(1)的有限差分格式為：

其中uk+1(x)是u(x,tk+1)的近似，即

(3)

對于系數(shù)wj，容易發(fā)現(xiàn)：

(i)wj>0,j=0,1,…,k；

(ii)1=w0≥w1≥w2≥…≥wk→0,k→∞。

設λ=Γ(2-α)hα，于是得到(3)的等價形式：

(4)

特別地，當k=0時，格式(4)變?yōu)椋?/p>

(5)

于是，式(4)和(5)，以及初、邊值條件

u0(x)=g(x)，x∈Λ，

uk+1(-1)=uk+1(1)=0,k≥0

組成了一個半離散問題。

(6)

2.2 方程(2)的有限差分格式及弱形式

同樣為了方便起見，在格式建立過程中考慮f≡0。使用如下的式子：對所有的0≤k≤K-1，

則(2)的有限差分格式為：

其中uk+1(x)是u(x,tk+1)的近似，即

(7)

特別地，當k=0時，格式(7)變?yōu)椋?/p>

(8)

于是，式(7)和(8),以及初、邊值條件

組成了一個半離散問題。

(9)

3 全離散

3.1 方程(6)的Legendre 譜配置法

令LN(x)為N次Legendre多項式。ξj,ωj(j=0,1,2,…,N)分別是 Legendre-Gauss-Lobatto(LGL)點和權(quán)。Legendre-Gauss型求積如下給出：

對任意的連續(xù)函數(shù)φ和ψ，其離散內(nèi)積及范數(shù)定義為

(10)

其中雙線性型AN(·,·)和泛函FN(·)分別定義為

3.2 方程(9)的Legendre 譜配置法

(11)

其中其中雙線性型BN(·,·)和泛函SN(·)分別定義為

4 數(shù)值有效性

4.1 數(shù)值實施

(12)

(13)

hj∈PN(Λ)，hj(ξi)=δi,j,?i,j∈{0,1,…,N}，其中δi,j是Kronecker函數(shù)。

對上述方程組使用離散內(nèi)積定義，得

于是問題(1)和(2)的矩陣表述為：

其中Fi=FN(hi),Si=SN(hi),并且?i,j∈{0,1,…,N}，

Hij=Aij+λBij,Aij=ωiδij，

4.2 數(shù)值算例

例1 考慮方程(1)其初值為g(x)=0，精確解u的觀測值z為

z(x,t)=tx(1-x)，

源項f定義為

其中α*=0.5給定的導數(shù)精確階。

我們的目標是，根據(jù)已知的觀測數(shù)據(jù)z，去估計精確階α*的值。為此，采用所介紹的優(yōu)化算法，給α取定初值進行迭代求解。我們?nèi)》謹?shù)階導數(shù)的初始值α0=0.7，根據(jù)所給算法利用MATLAB的編寫程序并運行，得出迭代過程中導數(shù)階的變化情況(見表1)和相應的精確階與迭代階的誤差曲線(見圖1)。

圖1 分數(shù)階導數(shù)的精確階數(shù)和迭代階數(shù)的誤差曲線(例1)

表1 不同迭代次數(shù)下的分數(shù)階導數(shù)的階數(shù)(例1)

例2 考慮方程(1)，假設初值g(x)=0，觀測值z(x,t)=t2sin(2πx)，其源項f為

同樣地，本例中取初始導數(shù)階α0=0.8，編寫程序并運行，得出迭代過程中導數(shù)階的變化情況(見表2)和相應的精確階與迭代階的誤差曲線(見圖2)。

圖2 分數(shù)階導數(shù)的精確階數(shù)和迭代階數(shù)的誤差曲線(例2)

表2 不同迭代次數(shù)下的分數(shù)階導數(shù)的階數(shù)(例2)

5 總結(jié)

在本文中，我們推導出了時間分數(shù)階擴散方程的解關(guān)于分數(shù)階導數(shù)α所滿足的方程，稱之為分數(shù)階擴散伴隨方程。首先，我們對兩個方程進行時間離散，得出了相應的有限差分格式和弱形式。然后，我們使用了Legendre譜配置法對弱形式中的半離散解進行多項式逼近，得到全離散格式。最后，在源項f已知這一情況下，我們利用序列最小優(yōu)化法對分數(shù)階α進行了估計。數(shù)值結(jié)果表明，我們推導出的伴隨方程是正確的，迭代序列是收斂的，且精度良好。