一種基于微分求積的空間分?jǐn)?shù)階擴(kuò)散方程求解方法

朱曉鋼， 聶玉峰

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

1 引 言

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

考慮變系數(shù)空間分?jǐn)?shù)階擴(kuò)散方程

(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的非負(fù)函數(shù)且不同時(shí)恒為0，當(dāng)κ(x) ≡ 0時(shí)g1(t)可不為0，當(dāng)υ(x) ≡ 0時(shí)g2(t)可不為0。在式(1)中,分?jǐn)?shù)階導(dǎo)數(shù)的定義為

(4)

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

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

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

2 預(yù)備知識(shí)

回顧分?jǐn)?shù)階導(dǎo)數(shù)的定義、基本性質(zhì)及與徑向基函數(shù)方法相關(guān)的一些輔助結(jié)論。

2.1 分?jǐn)?shù)階導(dǎo)數(shù)

設(shè)α∈+和Λ= [a,b]，那么

(5)

(6)

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

(7)

(8)

(9)

(10)

關(guān)于分?jǐn)?shù)階導(dǎo)數(shù)的更多性質(zhì)參見(jiàn)文獻(xiàn)[3,27]。

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

(?x∈Λ)(11)

式中

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

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

Reciprocal Multiquadrics:

Thin-Plate Splines:

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

3 分?jǐn)?shù)階導(dǎo)數(shù)的微分求積公式

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

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

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

(16)

(17)

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

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

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

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

(20)

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

(?x∈Λ)(21)

(22)

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

(23)

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

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

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

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

(26)

特別地，當(dāng)k= 0,1，有

(27)

這是因?yàn)?α1/?±xα= 0和?αx/?±xα= 0。然后，針對(duì)每個(gè)節(jié)點(diǎn)xi，整理式(24～27)，最終得

(28)

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

(29)

(30)

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

(32)

(33)

式中P∈+，且當(dāng)k(x)為k(x)時(shí)，有

(34)

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

(x1-x0)

(35)

4 分?jǐn)?shù)階擴(kuò)散方程的微分求積法

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

(36)

(37)

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

(38)

5 高維推廣

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

(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) 是非負(fù)函數(shù)且不同時(shí)恒為0。邊界條件(41)應(yīng)與式(39)相容，如當(dāng)κ+(x) ≠ 0時(shí)，g(x,y,t)必須在x=xL處為0;反之，如果κ-(x) ≠ 0 ，則必須在x=xR處為0。式(39)的分?jǐn)?shù)階導(dǎo)數(shù)的結(jié)構(gòu)類似于式(5,6)，如在y-軸方向的定義為

(42)

為構(gòu)造二維分?jǐn)?shù)階導(dǎo)數(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ù)值算例

通過(guò)5個(gè)數(shù)值算例來(lái)驗(yàn)證所提方法(分別簡(jiǎn)稱RM-DQ和TP-DQ方法)的精確性和可靠性。徑向基函數(shù)的分?jǐn)?shù)階導(dǎo)數(shù)的計(jì)算統(tǒng)一選擇50個(gè)積分節(jié)點(diǎn)和權(quán)值，運(yùn)用公式

(48)

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

(49)

來(lái)計(jì)算，ν=2和∞。固定σ= 2，而對(duì)于Reciprocal Multiquadrics的形狀參數(shù)，如何選取最優(yōu)值(誤差最小)仍然是一個(gè)開(kāi)放性問(wèn)題，有待深入研究，但通過(guò)參考文獻(xiàn)[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ù)值結(jié)果列入表1。 可知，分?jǐn)?shù)階導(dǎo)數(shù)的微分求積公式是有效的，在給定的參數(shù)下除α= 1.6和1.8的最大模誤差外，TP-DQ方法的誤差均小于RM-DQ方法。

算例2考慮空間分?jǐn)?shù)階擴(kuò)散方程

f(x,t)

(51)

計(jì)算區(qū)域?yàn)閇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ù)通過(guò)解析解代入式(51)獲得。由于u(x,t)的一階導(dǎo)數(shù)及其本身在邊界點(diǎn)xL和xR處均為0，根據(jù)性質(zhì)(9，10)，Caputo和Riemann-Liouville導(dǎo)數(shù)等價(jià)。取α= 1.6和c= 0.68，比較有限元法[11]、無(wú)網(wǎng)格點(diǎn)插值法[14]和微分求積法在t= 0.5處的均方誤差，對(duì)比結(jié)果列入表2，其中，有限元法取τ= 1/M2，無(wú)網(wǎng)格點(diǎn)插值法和微分求積法選用τ= 1/M。由

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

表2 當(dāng)α= 1.6時(shí)有限元法、無(wú)網(wǎng)格點(diǎn)插值法和微分求積法在t= 0.5處的均方誤差比較

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

(53)

計(jì)算區(qū)域?yàn)閇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ù)值結(jié)果列入表3?？梢钥闯?，TP-DQ方法在空間方向的收斂速率大于二階；此外，即使離散節(jié)點(diǎn)很少，該方法也能取得令人滿意的計(jì)算精度。

算例4考慮二維空間分?jǐn)?shù)階擴(kuò)散方程

(55)

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

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

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

算例5考慮三維空間分?jǐn)?shù)階擴(kuò)散方程

(56)

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

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

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

7 結(jié) 論

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