李玉山（甘肅政法學院 網絡空間安全學院，甘肅 蘭州，730070）

1 引 言

進入21世紀以來,隨著分數階微積分在反常擴散、非牛頓流體力學、多孔介質力學、粘彈性材料力學、地球物理、生物醫(yī)學工程、經濟學等諸多領域的建模及廣泛應用[1-6],分數階微積分的理論及應用開始受到越來越多的關注,相關的結果也陸續(xù)出現,參見Podlubny, Kilbas等人的專著[7、8]以及中文專著[9、10]。分數階反常擴散方程是由經典擴散方程推廣而得到的。相比整數階擴散方程，分數階擴散方程更適合描述反常擴散現象，因為分數階導數能夠描述具有記憶和遺傳性質的非均勻物質[12、13]。

時間-空間分數階擴散方程

02α＜＜被運用在反常擴散模型中，通常,若考慮時間相關性或擴散的記憶性，就得到時間分數階擴散方程；若考慮空間相關性或非局域性，則得到空間分數階擴散方程；如果既考慮時間相關性，又考慮空間整體性，則得到所謂的時間-空間分數階擴散方程。對于時間分數階擴散方程的數值求解，已經有了很多的研究成果[15-20]。對于時間-空間擴散方程的數值求解，相對來說，研究的比較少，劉發(fā)旺等[21,22]給出了空間分數階拉普拉斯算子的處理辦法，但是時間項依然是整數階的。

本文將探討一維時間-空間擴散方程的數值求解，利用矩陣轉換技術處理空間分數階拉普拉斯算子和有限差分法處理時間分數階項，得到代數方程組，同時利用分離變量法得到解析解表達式，并且給出數值算例作比較。

2 問題的提出

考慮如下的齊次Dirichlet邊界條件的時間-空間分數階擴散方程的初邊值問題：

其中Γ為伽馬函數。

為階數為α(0 ＜α≤2)的空間分數階拉普拉斯算子，由拉普拉斯算子-（）Δ的譜分解來定義，定義如下：

設

那么對于任意的fγ∈F，定義分數階拉普拉斯算子為：

本文考慮給定f(x),p(t),g(x)和方程(1)-(3)，數值求解u(x,t)。

3 數值方法

利用矩陣轉換技巧求解齊次Dirichlet邊界條件的時間-空間分數階擴散方程的初邊值問題。先考慮標準的一維擴散方程：

取 正 整 數 1N＞ ,空 間 步 長 1/hN= ,ixih=(0iN≤≤)，引進標準的中心差分離散二階空間導數：

利用文獻[21,22]的結果，方程（14）-（16）可以離散為如下矩陣形式

取正整數1M＞,時間步長(0jM≤≤)，用隱式差分格式離散時間分數階項[20]：

將（19）寫為矩陣形式，代入（17）得：

利用Matlab軟件求解代數方程組（20），即可得到問題（1）-（3）的近似解。

為了驗證本文的方法，利用分離變量法，給出問題（1）-（3）的解析解如下：

其中,Eβ β為廣義的Mittag-Leffler函數，定義如下：

在計算（21）的過程中，只計算無窮級數的前50項來近似得到u(x,t)，利用Matlab程序[23]計算Mittag-Leffler函數時取精度為10-6。

4 數值例子

例:取T=1,初始值u(x,0)=g(x)=x(1-x),f(x)=sin(x),p(t)=t,u(0,t)=u(1,t)=0,γ=1,M=N=100。

圖1 α =1.2, β =0.3時的數值解、解析解以及誤差圖

圖2 α =1.8, β =0.3時的數值解、解析解以及誤差圖

圖3 α =1.2, β =0.7時的數值解、解析解以及誤差圖

圖4 α =1.8, β =0.7時的數值解、解析解以及誤差圖

表1 數值解和精確解CPU占用時間對比 單位（s）

圖1-圖4給出了例1中α=1.2,1.8和β=0.3,0.7的數值結果，并與解析表達式（3.15）做了比較，可以看出，本文中的方法很有效，并且具有較高的精度。由表1可以看出，隨著時間和空間離散點的增加，解析解的CPU占用時間遠高于本文提出的數值方法的CPU占用時間，這是由于利用解析表達式（3.15）計算時，牽涉到計算兩個無窮級數、數值積分以及Mittag-Leffler函數，在時效上有很大的缺陷，而本文提出的方法避免了這一缺陷，只需要在每一時間層解一個線性方程組，有較高的精度，且占用極少的CPU時間。

5 結 語

本文研究了一維具有齊次Dirichlet邊界條件的時間-空間分數階擴散方程的數值求解，采用矩陣轉換技巧求解且給出了數值例子，并和解析解做了比較。數值例子表明，1）本文的方法是很有效的，具有較高的計算精度，可以處理實際問題；2）矩陣轉換技巧求解只需要求解線性方程組，而利用解析表達式求解計算比較復雜，且占用CPU時間較長；3）利用矩陣轉換技巧可以推廣到第二類、第三類邊界條件以及非齊次邊界條件的問題。