高廣花， 湯 銳， 楊 倩

(南京郵電大學(xué) 理學(xué)院，江蘇 南京 210023)

分?jǐn)?shù)階偏微分方程是指包含未知函數(shù)分?jǐn)?shù)階偏導(dǎo)數(shù)的方程.天文學(xué)[1]、金融學(xué)[2]、醫(yī)藥學(xué)[3]、物理學(xué)[4]等諸多領(lǐng)域中的許多數(shù)學(xué)模型都可以用分?jǐn)?shù)階偏微分方程描述.相較于整數(shù)階偏微分方程，分?jǐn)?shù)階偏微分方程在模擬某些物理以及化學(xué)過程中的精確性更高，所以分?jǐn)?shù)階偏微分方程的理論研究及應(yīng)用正成為備受關(guān)注的熱點(diǎn)問題之一.分?jǐn)?shù)階導(dǎo)數(shù)有許多不同的定義，如Riemann-Liouville導(dǎo)數(shù)、Riesz導(dǎo)數(shù)、Caputo導(dǎo)數(shù)等.一般來說，時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)用Caputo導(dǎo)數(shù)描述，而空間分?jǐn)?shù)階導(dǎo)數(shù)用Riemann-Liouville導(dǎo)數(shù)或Riesz導(dǎo)數(shù)定義.然而，大多數(shù)分?jǐn)?shù)階微分方程的解析解很難給出，尤其是變系數(shù)和非線性情形.因此，對分?jǐn)?shù)階微分方程尋找有效的數(shù)值模擬方法成為當(dāng)前研究的重要課題之一.

目前，對于空間四階時(shí)間分?jǐn)?shù)階擴(kuò)散方程，已有諸多學(xué)者進(jìn)行了研究.Jafari[5]等利用Adomian分解法得到了定義在有界空間域上的空間四階分?jǐn)?shù)階擴(kuò)散波方程的解.Hu和Zhang[6]利用降階法研究了空間四階分?jǐn)?shù)階擴(kuò)散波方程的緊致差分格式，并證明了該格式的穩(wěn)定性和收斂性.結(jié)合降階法，通過定義平均值算子，Ji等[7]給出了一種處理第一類Dirichlet邊界條件的方法.Vong和Wang[8]通過定義平均算子和離散的四階差商算子，直接離散空間四階導(dǎo)數(shù)，建立了求解第一類Dirichlet邊界下空間四階時(shí)間分?jǐn)?shù)階慢擴(kuò)散方程的高階緊致差分格式.Yao和Wang[9]考慮了Neumann邊界下求解空間四階時(shí)間分?jǐn)?shù)階慢擴(kuò)散方程的有限差分方法.對于第一類Dirichlet邊界下的空間四階時(shí)間分?jǐn)?shù)階慢擴(kuò)散方程，Cui[10]構(gòu)造了一種緊致Stephson差分格式進(jìn)行離散，將未知函數(shù)及其一階導(dǎo)函數(shù)同時(shí)作為未知量耦合求解，但是所得格式的理論分析較為復(fù)雜.文獻(xiàn)[11]使用樣條法求解一類空間四階時(shí)間分?jǐn)?shù)階偏微分方程.Xu等[12]提出了求解一類具有弱奇性核的空間四階時(shí)間分?jǐn)?shù)階微積分方程的數(shù)值算法，然后采用離散能量法、Cholesky分解法和降階法證明了該算法的穩(wěn)定性和收斂性.

上述成果多是對單項(xiàng)時(shí)間分?jǐn)?shù)階問題進(jìn)行的研究.然而，某些現(xiàn)象僅用單項(xiàng)分?jǐn)?shù)階方程描述是不夠的，許多過程需要借助多項(xiàng)時(shí)間分?jǐn)?shù)階偏微分方程來描述，如具有損耗的物理過程[13]、黏彈阻尼現(xiàn)象[14]、氧氣通過毛細(xì)血管輸送到組織的過程[15]、在多種復(fù)雜黏彈性材料中發(fā)生的介質(zhì)不規(guī)則擴(kuò)散現(xiàn)象[16]等.特別是，多項(xiàng)時(shí)間分?jǐn)?shù)階擴(kuò)散波方程可以有效地描述連續(xù)時(shí)間隨機(jī)游動模型中的冪律頻率依賴性[17].

查閱文獻(xiàn)可知關(guān)于多項(xiàng)時(shí)間分?jǐn)?shù)階方程已有不少研究.Jin等[16]基于Galerkin有限元方法，對含有多個(gè)Caputo分?jǐn)?shù)階導(dǎo)數(shù)的時(shí)間分?jǐn)?shù)階擴(kuò)散波方程建立了一種簡單的數(shù)值格式.Gao等[18]在超收斂點(diǎn)處構(gòu)造了逼近多項(xiàng)Caputo分?jǐn)?shù)階導(dǎo)數(shù)線性組合的數(shù)值公式，并證明了該公式至少可以達(dá)到二階精度.Gao和Liu[19]采用降階法對第二類Dirichlet邊界下的空間四階時(shí)間多項(xiàng)分?jǐn)?shù)階波方程建立了一種空間緊致差分格式，其中時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)利用L1公式進(jìn)行了逼近.Abdel-Rehim等[20]分別對時(shí)間分?jǐn)?shù)波方程、剪切波方程、阻尼波方程的近似解進(jìn)行了仿真,用G-L公式離散了Caputo時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)，并討論了這些數(shù)值方法的Von-Neumann穩(wěn)定性條件.Wei[21]構(gòu)造了求解一類多項(xiàng)時(shí)間分?jǐn)?shù)階擴(kuò)散方程的全離散局部間斷的Galerkin方法，并證明了該方法的無條件穩(wěn)定性和收斂性.利用降階法，Sun等[22]在超收斂點(diǎn)處建立了求解時(shí)間多項(xiàng)分?jǐn)?shù)階波方程的兩種差分格式，并證明了格式的唯一可解性、穩(wěn)定性和收斂性.Zhang等[23]考慮了多項(xiàng)時(shí)間分?jǐn)?shù)階Burgers流體模型，采用分離變量法得到了解析解，然后基于時(shí)間加權(quán)位移Grünwald差分算子和空間Legendre譜方法建立了數(shù)值方法，并提出了一種改進(jìn)格式以提高收斂精度.

以上問題都是通過常數(shù)階分?jǐn)?shù)階微分方程來描述的.然而，近年來，越來越多的學(xué)者發(fā)現(xiàn)，各種重要的動力學(xué)行為表現(xiàn)出變階分?jǐn)?shù)階行為.這些分?jǐn)?shù)階行為可能隨時(shí)間、空間或其他條件而變化，這表明變階分?jǐn)?shù)階微積分是描述復(fù)雜問題的一種有效的數(shù)學(xué)工具.變階分?jǐn)?shù)階導(dǎo)數(shù)是常數(shù)階分?jǐn)?shù)階導(dǎo)數(shù)的推廣，后者無法表征的一些現(xiàn)象，如擴(kuò)散過程在多孔介質(zhì)中的演化，外場隨時(shí)間的變化等，都可以通過變階分?jǐn)?shù)階導(dǎo)數(shù)進(jìn)行描述.Samko和Ross[24—25]早在1993年和1995年就率先提出了變階Riemann-Liouville和Marchaud分?jǐn)?shù)階微積分及它們的一些基本屬性.Lorenzo和Hartley[26—27]總結(jié)了變階算子的概念，其中變階算子可隨時(shí)間或空間變化.此外，他們對變階積分和微分的概念，以及數(shù)學(xué)概念與物理過程之間的關(guān)系進(jìn)行了深入研究.在此基礎(chǔ)上，他們還進(jìn)一步討論了變階分?jǐn)?shù)階微分模型的一些新的推廣.Coimbra[28]在考慮分?jǐn)?shù)階導(dǎo)數(shù)的拉普拉斯變換的基礎(chǔ)上，給出了變階微分算子的一個(gè)新的定義.Ingman等[29—30]采用時(shí)間變階算子對黏彈性變形過程進(jìn)行了建模.Pedro等[31]利用變階微積分研究了懸浮在帶有阻力的黏性流體中的顆粒的運(yùn)動.Du等[32]提出了求解多維變階時(shí)間分?jǐn)?shù)階慢擴(kuò)散方程的兩種差分格式.二者均具有時(shí)間二階精度，在空間上分別具有二階和四階精度.Wei和Yang[33]提出了求解一類時(shí)間變階分?jǐn)?shù)階擴(kuò)散問題的一種數(shù)值方法，時(shí)間方向用有限差分方法進(jìn)行了離散，空間方向采用局部間斷的 Galerkin方法進(jìn)行逼近，并通過理論分析和數(shù)值算例證明了該方法的高精度.文獻(xiàn)[34]基于位移二進(jìn)制塊區(qū)間劃分和多項(xiàng)式近似，建立了求解變階時(shí)間分?jǐn)?shù)階擴(kuò)散方程的快速緊致差分格式，并從理論上證明了格式的無條件穩(wěn)定性和收斂性.Ganji等[35]考慮了一類非線性的變階微積分方程.Tian等[36]提出了一種求解時(shí)間多項(xiàng)變階分?jǐn)?shù)階偏微分方程的無網(wǎng)格配置方法.

注意到，Ji等[7]建立了求解第一類Dirichlet邊界下空間四階時(shí)間分?jǐn)?shù)階方程的數(shù)值格式，但是該工作在邊界點(diǎn)處引入了對空間二階導(dǎo)數(shù)四點(diǎn)處加權(quán)平均的逼近，且只研究了單項(xiàng)時(shí)間常數(shù)階問題.Cui[10]以空間一階導(dǎo)數(shù)作為中間變量，建立了求解第一類Dirichlet邊界下的空間四階時(shí)間常數(shù)階分?jǐn)?shù)階慢擴(kuò)散方程的有限差分格式，但是其理論分析極具挑戰(zhàn)性.Du等[32]考慮了求解空間二階時(shí)間單項(xiàng)變階分?jǐn)?shù)階慢擴(kuò)散方程的有限差分方法，但是對于時(shí)間多項(xiàng)變階分?jǐn)?shù)階問題沒有進(jìn)行相關(guān)研究.Tian等[36]考慮了求解時(shí)間多項(xiàng)變階分?jǐn)?shù)階方程，但僅考慮了空間二階問題，且未對所提算法進(jìn)行理論分析.本文中，我們考慮空間四階時(shí)間多項(xiàng)變階分?jǐn)?shù)階問題.通過對邊界條件進(jìn)行巧妙處理，建立了求解第一類Dirichlet邊界條件下空間四階時(shí)間多項(xiàng)變階分?jǐn)?shù)階慢擴(kuò)散方程的高階穩(wěn)定的數(shù)值格式，并用離散能量分析方法證明了所得格式的穩(wěn)定性和收斂性.主要貢獻(xiàn)包括:

(ⅰ)在超收斂點(diǎn)處構(gòu)造了逼近時(shí)間多項(xiàng)變階分?jǐn)?shù)階導(dǎo)數(shù)線性組合的數(shù)值公式，并證明了該公式至少可以達(dá)到時(shí)間二階精度.

(ⅱ)對邊界進(jìn)行巧妙處理.本文中在邊界點(diǎn)處只引入了對空間二階導(dǎo)數(shù)兩點(diǎn)處加權(quán)平均的逼近，與文獻(xiàn)[7]相比，算法更為簡潔.應(yīng)用降階法，將空間二階導(dǎo)數(shù)作為中間變量，并用離散能量分析方法對所建立的數(shù)值格式進(jìn)行了穩(wěn)定性和收斂性分析.與文獻(xiàn)[10]相比，本文的理論分析更為直觀易行.

1 預(yù)備知識

考慮如下第一類Dirichlet邊界下時(shí)間多項(xiàng)變階分?jǐn)?shù)階慢擴(kuò)散方程：

f(x,t),x∈(0,L),t∈(0,T],

(1)

u(0,t)=g1(t),u(L,t)=g2(t),t∈(0,T],

(2)

ux(0,t)=γ1(t),ux(L,t)=γ2(t),t∈(0,T],

(3)

u(x,0)=φ(x),x∈[0,L],

(4)

假設(shè)u∈C3[t0,tn],1≤n≤N.記(L2,ku)(s)為函數(shù)u基于(tk-1,u(tk-1)),(tk,u(tk))和(tk+1,u(tk+1))3點(diǎn)插值得到的二次多項(xiàng)式，則有

且

(5)

余項(xiàng)為

u(s)-(L2,ku)(s)=

s∈[tk-1,tk+1],ξk∈(tk-1,tk+1),

1≤k≤n-1,1≤n≤N.

(6)

在區(qū)間[tn-1,tn]上，利用(tn-1,u(tn-1))和(tn,u(tn))兩點(diǎn)對函數(shù)u做線性插值，記為(L1,nu)(s)，則

且

此外，有

(8)

定義

σ>0,n≥1.

引理2設(shè)m≥1， 由

注1引理1和引理2的證明可參考文獻(xiàn)[18]和文獻(xiàn)[32]得出，此處從略.

定理1設(shè)n≥1，u∈C3[t0,tn]，記

則

|rn|≤

證明易知

由

(tn-1+σn-s)-αr(tn-1+σn)-1ds]=

(tn-1+σn-s)-αr(tn-1+σn)-1ds

及

得到

(tn-1+σn-t0)-αr(tn-1+σn)]≤

對于(10)式中的第二項(xiàng)，由(8)式可得

(tn-1+σn-s)-αr(tn-1+σn)ds.

(12)

注意到Fn(σn)=0，從而

因此

(tn-1+σn-s)-αr(tn-1+σn)ds.

進(jìn)一步可得

將(11)式和(13)式代入(10)式，即得(9)式. 定理證畢.

根據(jù)(5)式和(7)式可將數(shù)值微分公式改寫為

[(k+1+σn)2-α(tn-1+σn)-2(k+σn)2-α(tn-1+σn)+

2(k+σn)1-α(tn-1+σn)+(k-1+σn)1-α(tn-1+σn)],

1≤k≤n-2,

(n-2+σn)2-α(tn-1+σn)]-

(n-2+σn)2-α(tn-1+σn)].

引理3[37]記θ(s)=(1-s)3[10-3(1-s)2],ξ(s)=(1-s)3[5-3(1-s)2].

(ⅰ)若函數(shù)g∈C6[x0,x1]，則有

(ⅱ)若函數(shù)g∈C6[xM-1,xM]，則有

(ⅲ)若函數(shù)g∈C6[xi-1,xi+1]，則有

下列引理給出了函數(shù)u(t)在t=tn-1+σn處的值的一個(gè)二階逼近公式.

引理4[9]若u∈C3[tn-1,tn],σn滿足0<σn<1，則

u(tn-1+σn)=σnu(tn)+(1-σn)u(tn-1)+O(τ2).

2 差分格式的建立

對任意網(wǎng)格函數(shù)u∈Uh，引入如下記號：

算子H通常被稱為平均值算子或緊算子.

引入中間變量v(x,t)=uxx(x,t)，則方程(1)～(4)可等價(jià)寫成

f(x,t),x∈(0,L),t∈(0,T],

(18)

v(x,t)=uxx(x,t),x∈(0,L],t∈(0,T],

(19)

u(0,t)=g1(t),u(L,t)=g2(t),t∈(0,T],

(20)

ux(0,t)=γ1(t),ux(L,t)=γ2(t),t∈(0,T],

(21)

u(x,0)=φ(x),x∈[0,L].

(22)

分別在點(diǎn)(xi,tn-1+σn)和(xi,tn)處考慮方程(18)～(19)，得到

vxx(xi,tn-1+σn)+qu(xi,tn-1+σn)=

f(xi,tn-1+σn), 0≤i≤M,1≤n≤N,

(23)

v(xi,tn)=uxx(xi,tn), 0≤i≤M,0≤n≤N.

(24)

對(23)～(24)式的兩端同時(shí)作用算子H可得

Hvxx(xi,tn-1+σn)+qHu(xi,tn-1+σn)=

Hf(xi,tn-1+σn)，1≤i≤M-1， 1≤n≤N，

(25)

Hv(xi,tn)=Huxx(xi,tn),0≤i≤M,0≤n≤N.

(26)

根據(jù)引理3、引理4和定理1，可得

1≤i≤M-1,1≤n≤N,

(27)

(28)

其中存在正常數(shù)c1，滿足

(29)

(30)

在方程(1)中，令x→0+，并注意到邊界條件(2)，可得

(31)

方程(1)兩端同時(shí)關(guān)于x求偏導(dǎo)一次，再令x→0+，并注意到邊界條件(3)，可得

0≤n≤N.

(32)

同理，令x→L-，可以得出類似的2個(gè)方程.

當(dāng)i=0時(shí)，(26)式為Hv(x0,tn)=Huxx(x0,tn),0≤n≤N,由引理3中(15)式及(31)～(32)式可得

其中

并且存在一個(gè)正常數(shù)c2,使得

(34)

同理可得

(35)

其中

并且存在一個(gè)正常數(shù)c3,使得

(36)

注意到初邊值條件(20)和(22)，可得

定理2差分格式(39)～(44)等價(jià)于

(45)

2≤i≤M-2,1≤n≤N,

(46)

1≤n≤N,

(47)

且

(50)

其中pn-1+σn=σnp(tn)+(1-σn)p(tn-1),qn-1+σn類似定義.

證明仿照文獻(xiàn)[37]中定理3.1的證明易得，此處從略.

3 差分格式的穩(wěn)定性和收斂性

先給出3個(gè)引理.

引理5[37]對任意網(wǎng)格函數(shù)v∈Uh，有

引理7[18]對任意網(wǎng)格函數(shù)un∈Uh,0≤n≤N，uσn=σnun+(1-σn)un-1，有

1≤n≤N.

對于單項(xiàng)變階問題，文獻(xiàn)[32]已有結(jié)果. 然而對于多項(xiàng)變階問題，條件一非負(fù)性易得，條件二的證明不再顯而易見. 但是參考對多項(xiàng)常數(shù)階問題[18]和單項(xiàng)變階問題[32]的研究技巧，條件二的證明亦可得出，此處從略.

則有

(59)

證明由(54)式易得

由Cauchy-Schwarz不等式可得

1≤n≤N,

(62)

1≤n≤N.

(63)

將(62)～(63)式代入(61)式，整理得

1≤n≤N.

(64)

將引理5、引理6及引理7應(yīng)用到(64)式，可得

由(55)～(56)式易得

從而有

改寫(67)式可得

注意到

進(jìn)一步可得

1≤n≤N.

遞推可得

注意到(57)式并應(yīng)用引理5，即可證明(59)式．定理證畢．

其中

證明將(27)～(28)式、(33)式、(35)式以及(37)～(38)式分別與(39)～(44)式對應(yīng)相減，可得誤差方程組為

1≤i≤M-1,1≤n≤N,

對上述誤差方程組應(yīng)用定理3，可得

將上式兩邊開方即得定理結(jié)果.定理4證畢．

4 數(shù)值實(shí)驗(yàn)

本節(jié)旨在應(yīng)用一些數(shù)值算例驗(yàn)證本工作中所得算法的計(jì)算效果及數(shù)值精度. 記

-πt4-sin 1.該問題的精確解為u(x,t)=t4sin(πx)+cosx.

選取不同的λr和αr(t)(r=0,1,…,m)，應(yīng)用差分格式(45)～(49)計(jì)算該問題．首先，取固定且足夠小的空間步長h，時(shí)間步長τ從1/10細(xì)化到1/160，記錄最大計(jì)算誤差及收斂階. 由表1可知，差分格式(45)～(49)在時(shí)間方向能達(dá)到二階收斂，計(jì)算結(jié)果與理論分析吻合．

類似地，表2列出了空間步長h取不同值時(shí)的最大誤差及收斂階．由表2可知，差分格式(45)～(49)的空間方向收斂階可達(dá)到理論上的四階，與預(yù)期結(jié)果相同．

5 結(jié)論

本文在超收斂點(diǎn)處建立了逼近時(shí)間多項(xiàng)變階Caputo分?jǐn)?shù)階導(dǎo)數(shù)線性組合的數(shù)值微分公式，并研究了所得微分公式的數(shù)值精度.首先，通過求解非線性方程，在每個(gè)時(shí)間層上找到了一個(gè)特殊的點(diǎn)，使得所得數(shù)值微分公式逼近時(shí)間多項(xiàng)變階Caputo分?jǐn)?shù)階導(dǎo)數(shù)的線性組合在該點(diǎn)的值至少可以達(dá)到二階精度. 其次，利用該數(shù)值微分公式求解第一類Dirichlet邊界下空間四階時(shí)間多項(xiàng)變階分?jǐn)?shù)階慢擴(kuò)散方程，建立了有效的差分格式，并利用能量分析法證明了所得格式的穩(wěn)定性和收斂性. 數(shù)值算例證明了所得格式的有效性及理論結(jié)果的正確性.

另外，本文僅討論了第一類Dirichlet邊界下空間四階時(shí)間多項(xiàng)變階分?jǐn)?shù)階慢擴(kuò)散方程的差分方法，很容易推廣到求解第二類Dirichlet邊界下空間四階時(shí)間多項(xiàng)變階分?jǐn)?shù)階慢擴(kuò)散方程的初邊值問題，此處不再贅述.

表1 M=100時(shí)差分格式(45)～(49)的數(shù)值誤差及時(shí)間方向收斂階

表2 N=10 000時(shí)差分格式(45)～(49)的數(shù)值誤差及空間方向收斂階