朱曉鋼， 聶玉峰

(1.邵陽學院 理學院，邵陽 422000；2.西北工業(yè)大學 應用數(shù)學系，西安 710129)

1 引 言

分數(shù)階微積分是微積分理論的重要分支，也是經典微積分的任意階推廣，彌補了經典微積分理論在應用方面的不足。分數(shù)階微積分的一個最為突出的特點在于能捕捉動力學行為的記憶效應，更好地反映系統(tǒng)函數(shù)發(fā)展的歷史依賴性。近十幾年，利用帶分數(shù)階導數(shù)的微分方程刻畫自然界中的反常現(xiàn)象一直是學術界關注的焦點之一，潛在的應用已囊括各個學科領域，如電化學、金融、生物工程、材料科學和統(tǒng)計物理等[1-3]。

考慮變系數(shù)空間分數(shù)階擴散方程

(1)

式中1 <α≤2，x∈Λ，Λ= [xL,xR]，0

u(x,0) =u0(x)

(x∈Λ) (2)

u(xL,t) =g1(t),u(xR,t) =g2(t)

(0

此處，κ(x)和υ(x)是x的非負函數(shù)且不同時恒為0，當κ(x) ≡ 0時g1(t)可不為0，當υ(x) ≡ 0時g2(t)可不為0。在式(1)中,分數(shù)階導數(shù)的定義為

(4)

式中Γ(·)為歐拉Gamma函數(shù)。

在數(shù)學物理中，空間分數(shù)階偏微分方程(1～3)可以描述諸多現(xiàn)象，如反常擴散、彈性振動和波的傳播等[4，5]。與整數(shù)階偏微分方程相比，其在反映物理系統(tǒng)隨時間演化上更具優(yōu)勢。由于解析技術的局限性，數(shù)值方法越來越受到人們的重視，并且在研究中獲得了廣泛的應用，典型的方法包括有限差分法[6-9]、有限元法FEM[10-12]、無網格點插值法MPIM[13,14]和有限體積法[15]。Xu等[16]給出了分數(shù)階對流擴散方程的間斷Galerkin方法。Ilic等[17，18]基于導數(shù)的矩陣表示提出了矩陣轉換法。文獻[19]討論了一種穩(wěn)定的方向分裂譜Galerkin方法；文獻[20]構造了求解變系數(shù)空間分數(shù)階擴散方程的徑向基配置法。Ford等[21]將Riemann-Liouville 導數(shù)轉化成等價的Hadamard有限部分積分，并運用分段二次多項式對被積函數(shù)插值建立了一種新的有效方法。Li等[22]研究了二維非線性空間分數(shù)階擴散方程的時空勒讓德譜配置法。

由于分數(shù)階導數(shù)的非局部性質，采用常規(guī)數(shù)值方法求解空間分數(shù)階問題通常面臨著算法復雜、計算量大和效率低等困難，因此尋找一種簡便高效的算法非常必要。微分求積法(DQ Method)是求解微分方程的高效數(shù)值方法，其基本思想是利用函數(shù)值的加權線性組合替換導數(shù)，將微分方程轉化成方程組[23]，試函數(shù)可選擇正交多項式、徑向基函數(shù)、樣條和拉格朗日插值多項式[24-26]。 該方法具有計算量小、對維數(shù)增長不敏感和易于編程實現(xiàn)等特性。鑒于此，本文將微分求積法推廣到方程組(1～3)及與之對應的高維空間分數(shù)階問題，介紹分數(shù)階導數(shù)的微分求積公式，利用Reciprocal Multiquadrics和Thin-Plate Splines兩種徑向基函數(shù)作為試函數(shù)確定加權系數(shù)，形成的常微分方程組采用Crank-Nicolson格式離散。最后，數(shù)值算例的數(shù)值結果及與現(xiàn)有算法的比較展示了其在計算精度和計算效率上的優(yōu)勢。

2 預備知識

回顧分數(shù)階導數(shù)的定義、基本性質及與徑向基函數(shù)方法相關的一些輔助結論。

2.1 分數(shù)階導數(shù)

設α∈+和Λ= [a,b]，那么

(5)

(6)

分別是α-階左Caputo和右Caputo分數(shù)階導數(shù)的定義，此處，f(x) ∈Cm(Λ)，當α?時，m=[α]+1，反之，當α∈時，m=α，[ · ]是向下取整函數(shù)。α-階左Riemann-Liouville和右Riemann-Liouville分數(shù)階導數(shù)的定義分別是

(7)

(8)

(9)

(10)

關于分數(shù)階導數(shù)的更多性質參見文獻[3,27]。

2.2 徑向基函數(shù)方法

(?x∈Λ)(11)

式中

(i= 1,2,…,mp)(13)

考慮兩種類型的徑向基函數(shù)，即

Reciprocal Multiquadrics:

Thin-Plate Splines:

式中k= 0,1,…,M，ε為形狀參數(shù)。對于Reciprocal Multiquadric徑向基函數(shù)，由式(11～13)生成的插值矩陣是嚴格正定的，即使在無多項式項P(x)時也一樣;而對于Thin-Plate Spline徑向基函數(shù)，插值矩陣卻是條件正定的，需滿足mp≥2。所以，在采用徑向基插值構造數(shù)值方法時，P(x)對于Thin-Plate Spline徑向基函數(shù)是不可或缺[28]的。

3 分數(shù)階導數(shù)的微分求積公式

介紹分數(shù)階導數(shù)的微分求積公式，基于Reciprocal Multiquadric和Thin-Plate Spline徑向基函數(shù)，給出計算加權系數(shù)的方法。在[xL,xR]上定義xL=x0

(i,k= 0,1,…,M)(14)

(i,k= 0,1,…,M)(15)

(16)

(17)

情況I Reciprocal Multiquadrics當mp= 0時，由于式(11～13)生成的代數(shù)問題是唯一可解的，因此，該情況的加權系數(shù)的計算公式可直接通過替換式(14，15)的基函數(shù)獲得，即

(i,k= 0,1,…，M)(18)

(i,k= 0,1,…，M)(19)

針對節(jié)點xi，將式(18，19)整理成矩陣形式為

(20)

情況II Thin-Plate Splines由于Thin-Plate Spline徑向基函數(shù)是條件正定，為了保證代數(shù)問題的適定性，需加上多項式項P(x)。為簡便起見，本文取mp= 2，根據(jù)式(11～13)，有

(?x∈Λ)(21)

(22)

利用式(22)消去式(21)的變量λ0(t) 和λ1(t)，可知

(23)

再通過式(14～17)，得到

(i= 0,1,…，M)(24)

(i= 0,1,…，M)(25)

此處，k= 2,3,…,M，并且

(26)

特別地，當k= 0,1，有

(27)

這是因為?α1/?±xα= 0和?αx/?±xα= 0。然后，針對每個節(jié)點xi，整理式(24～27)，最終得

(28)

下面介紹Dα(· )的計算過程。一般來說，函數(shù)分數(shù)階導數(shù)的顯式表達式很難推導，這是由其弱奇異性質決定的。但通過適當?shù)姆e分變換，可以通過數(shù)值積分公式來逼近。分別對左右Caputo分數(shù)階導數(shù)作積分變換ξ=x-(x-xL)(1+ζ)/2和ξ=x+(xR-x)(1+ζ)/2，可以驗證

(29)

(30)

可以看出，式(5,6)的積分已演變成

(32)

(33)

式中P∈+，且當k(x)為k(x)時，有

(34)

Q(x) = (x-xi)2 σ - 2[(2σ-1)σlog(x-xi)2+

(x1-x0)

(35)

4 分數(shù)階擴散方程的微分求積法

通過上述微分求積公式建立求解方程組(1～3)的全離散微分求積格式。在區(qū)間[0,T]上定義等距網格節(jié)點tn=nτ(τ=T/N，N∈+)。將加權線性組合(16,17)代入式(1)，得

(36)

(37)

式中i= 1,2,…,M-1，以及

(38)

5 高維推廣

求解高維空間分數(shù)階問題是一項富有挑戰(zhàn)性的課題。將上述微分求積法推廣到高維情形，由于對維數(shù)增長不敏感，采用其處理高維問題具備較大的便利性。以二維情形為例，三維情形以此類推，考慮二維空間分數(shù)階擴散方程

(39)

式中(x,y)∈Λ,0

u(x,y,0) =u0(x,y)((x,y) ∈Λ)

(40)

u(x,y,t) =g(x,y,t)((x,y) ∈?Λ,0

(41)

此處， 1 <α,β≤2，Λ= [xL,xR]×[yL,yR]，?Λ為邊界，κ+(x) ，κ-(x) ，υ+(x) 和υ-(x) 是非負函數(shù)且不同時恒為0。邊界條件(41)應與式(39)相容，如當κ+(x) ≠ 0時，g(x,y,t)必須在x=xL處為0;反之，如果κ-(x) ≠ 0 ，則必須在x=xR處為0。式(39)的分數(shù)階導數(shù)的結構類似于式(5,6)，如在y-軸方向的定義為

(42)

為構造二維分數(shù)階導數(shù)的微分求積公式，在區(qū)域Λ上定義xL=x0

(43)

(44)

(45)

(46)

(k= 0,1,…,Mx)

而在y-方向，有

(k= 0,1,…,My)

(47)

式中i= 1,2,…,Mx-1，j= 1,2,…,My-1，以及

6 數(shù)值算例

通過5個數(shù)值算例來驗證所提方法(分別簡稱RM-DQ和TP-DQ方法)的精確性和可靠性。徑向基函數(shù)的分數(shù)階導數(shù)的計算統(tǒng)一選擇50個積分節(jié)點和權值，運用公式

(48)

來衡量誤差，Mtotal= (Mx+1)(My+1)，再令0<ε<1，相應的收斂速率利用

(49)

來計算，ν=2和∞。固定σ= 2，而對于Reciprocal Multiquadrics的形狀參數(shù)，如何選取最優(yōu)值(誤差最小)仍然是一個開放性問題，有待深入研究，但通過參考文獻[30-32]，可選擇=cdim/(M+1)0.5(c> 0，dim=，x或者y)。

(50)

(x)k表示階乘，即(x)k=x(x+1)…(x+k-1)。取c= 1.45和M= 20，不同α下的數(shù)值結果列入表1。 可知，分數(shù)階導數(shù)的微分求積公式是有效的，在給定的參數(shù)下除α= 1.6和1.8的最大模誤差外，TP-DQ方法的誤差均小于RM-DQ方法。

算例2考慮空間分數(shù)階擴散方程

f(x,t)

(51)

計算區(qū)域為[0,1]，施加初邊值條件

u(x,0) = 10x2(1-x)2

u(xL,t) =u(xR,t) = 0

(52)

解析解為u(x,t) = 10e-tx2(1-x)2。源函數(shù)通過解析解代入式(51)獲得。由于u(x,t)的一階導數(shù)及其本身在邊界點xL和xR處均為0，根據(jù)性質(9，10)，Caputo和Riemann-Liouville導數(shù)等價。取α= 1.6和c= 0.68，比較有限元法[11]、無網格點插值法[14]和微分求積法在t= 0.5處的均方誤差，對比結果列入表2，其中，有限元法取τ= 1/M2，無網格點插值法和微分求積法選用τ= 1/M。由

表1 算例1中M = 20時的數(shù)值結果

表2 當α= 1.6時有限元法、無網格點插值法和微分求積法在t= 0.5處的均方誤差比較

表2 可知，微分求積法的精度要高于前兩種方法。算例3考慮空間分數(shù)階擴散方程

(53)

計算區(qū)域為[0,1]，施加初邊值條件

u(x,0) = 0,u(xL,t) = 0,u(xR,t) =t2

(54)

解析解為u(x,t) =t2x4。取α= 1.2和τ= 2.0×10-3，TP-DQ方法在t= 1處的數(shù)值結果列入表3?？梢钥闯?，TP-DQ方法在空間方向的收斂速率大于二階；此外，即使離散節(jié)點很少，該方法也能取得令人滿意的計算精度。

算例4考慮二維空間分數(shù)階擴散方程

(55)

設Λ= [0,1]×[0,1]，解析解為u(x,y,t) =e-tx2(1-x)2y2(1-y)2，初邊值條件均由解析解給出，源函數(shù)通過解析解代入式(55)獲得。取α=1.5，β= 1.9和τ= 2.5×10-3，同時采用文獻[12]的有限元法和TP-DQ方法求解上述問題，并比較其在t= 1處的誤差及計算時間，為此，有限元法使用線性元與結構化網格，對比結果列入表4。 顯然，TP -DQ 方法在計算精度和收斂速率上都要高于有限元法；且在相同的參數(shù)下，微分求積法的運行時間遠小于有限元法。

表3 當α= 1.2時TP-DQ方法在t= 1處的數(shù)值結果

表4 當α= 1.5和β= 1.9時有限元法和TP-DQ方法的比較Tab.4 Comparison of FEM and TP-DQ method when α= 1.5 and β= 1.9

算例5考慮三維空間分數(shù)階擴散方程

(56)

設Λ= [0,1]×[0,1]×[0,1]和初邊界條件

u(x,y,z,0) = 0,u(1,y,z,t) =t2y3z3

u(x,1,z,t) =t2x3z3,u(x,y,1,t) =t2x3y3

u(0,y,z,t) =u(x,0,z,t) =u(x,y,0,t) = 0

(57)

κ+，υ+和ν+為間斷系數(shù)，分別定義如下，

解析解為u(x,y,z,t) =t2x3y3z3，源函數(shù)通過解析解代入式(56)獲得。取c= 1.8，α= 1.2，β= 1.3，γ= 1.5和τ= 5.0×10-3，設x-，y-和z-坐標方向的離散節(jié)點數(shù)分別為Mx,My和Mz，誤差同樣采用式(48)來衡量并將其分別表示成e2(tn,Mx,My,Mz)和e∞(tn,Mx,My,Mz)，RM-DQ方法在t= 1處的數(shù)值結果列入表5?？梢钥闯?，RM-DQ方法的誤差隨著節(jié)點增多而迅速衰減且收斂速率足夠高，證明本文方法適用于三維空間分數(shù)階問題及帶間斷系數(shù)的空間分數(shù)階問題。

注：上述算法，除算例2的有限元法和無網格點插值法之外，均采用Matlab R2012a編程實現(xiàn)，并在微機Lenovo H520上運算，配置：Intel(R) Pentium(R) G2030 3.00 GHz處理器和4 GB內存。

表5 當α= 1.2，β= 1.3和γ= 1.5時RM-DQ方法在t= 1處的數(shù)值結果Tab.5 Numerical results of RM-DQ method at t= 1 when α= 1.2，β= 1.3 and γ= 1.5

7 結 論

研究變系數(shù)空間分數(shù)階擴散方程的一種全離散微分求積方法，基于Reciprocal Multiquadric和Thin-Plate Spline徑向基函數(shù)，介紹了兩種近似分數(shù)階導數(shù)的微分求積公式。提出的方法具有精度高、靈活性強、無網格依賴和易于編程實現(xiàn)等特點。求解了5個典型問題，結果表明,在形狀參數(shù)選擇恰當?shù)那闆r下，該方法在計算精度與效率上要優(yōu)于一些現(xiàn)有算法。由于對維數(shù)增長不敏感，本文方法在處理高維分數(shù)階問題上有著較大的優(yōu)勢。