田洪材,曹俊英,王自強(qiáng)
(貴州民族大學(xué) 數(shù)據(jù)科學(xué)與信息工程學(xué)院,貴州 貴陽(yáng) 550025)
分?jǐn)?shù)階微積分應(yīng)用范圍廣泛,可以描述一些具有記憶性、非局部性的自然現(xiàn)象。分?jǐn)?shù)階導(dǎo)數(shù)和分?jǐn)?shù)階微分方程被許多科學(xué)和工程領(lǐng)域的研究人員所研究[1]。其中,具有空間Riesz分?jǐn)?shù)階導(dǎo)數(shù)和時(shí)間Caputo分?jǐn)?shù)階導(dǎo)數(shù)的時(shí)空分?jǐn)?shù)擴(kuò)散方程涉及許多應(yīng)用,如金融學(xué)[2]和分?jǐn)?shù)動(dòng)力學(xué)方程[3]。隨著分?jǐn)?shù)階微分方程模型的增多和精確解的不易獲得,人們對(duì)開(kāi)發(fā)快速的數(shù)值方法產(chǎn)生了極大的興趣。
Riesz分?jǐn)?shù)階導(dǎo)數(shù)被視為左右Riemann-Liouville分?jǐn)?shù)階導(dǎo)數(shù)的線性組合。本文研究時(shí)間為Caputo分?jǐn)?shù)階導(dǎo)數(shù),空間為Riesz分?jǐn)?shù)階導(dǎo)數(shù)的時(shí)空分?jǐn)?shù)階擴(kuò)散方程,對(duì)它的一個(gè)高階快速的數(shù)值格式進(jìn)行了構(gòu)造和分析。
考慮如下方程:
(1)
和初邊值條件
u(x,0)=φ(x),x∈[a,b],
u(a,t)=u(b,t)=0,?t∈[0,T],
0<α<1,1<β<2,
這里Γ(·)是伽馬函數(shù),f(x,t)和φ(x)為已知函數(shù)。
(2)
其中:
(3)
其中,kα=Γ(3-α)τα,pα=c1+2-0.5α。
當(dāng)m=3時(shí),
當(dāng)m≥4時(shí),
其中
(j+1)2-α-j2-α
其中U0=[φ1,φ2,…,φN-1]。
利用文獻(xiàn)[6]中的加權(quán)移位Grünwald-Letnikov差分方法在空間上和時(shí)間上利用式(2)與式(3),可得:
當(dāng)m=1,2時(shí),有:
(4)
當(dāng)m≥3時(shí),有:
(5)
其中,
D=
這里
wβ,0=λ1gβ,0,wβ,1=λ1gβ,1+λ2gβ,0,
wβ,k=λ1gβ,k+λ2gβ,k-1+λ3gβ,k-2,k≥2。
當(dāng)m≥3時(shí),數(shù)值格式(5)系數(shù)矩陣是一個(gè)Toeplitz矩陣,我們利用文獻(xiàn)[7]和[8]的思想方法,結(jié)合GMRES迭代方法與FFT方法建立數(shù)值格式(5)的快速迭代方法(FGMRES)。
利用文獻(xiàn)[9]中的思想方法,我們可以證明如下定理。
定理1:數(shù)值格式(5)的收斂階為O(h2+τ3-α)。
考慮一維時(shí)空分?jǐn)?shù)階擴(kuò)散方程模型,我們有[a,b]=[0,1],T=1,該方程的精確解和右端項(xiàng)分別為:u(x,t)=t4x2(1-x)2,
f(x,t)=24x2(1-x)2(t4-α)/Γ(5-α)+
t4((3-β)(4-β)(x2-β+(1-x)2-β)-
6(4-β)(x3-β+(1-x)3-β)+12(x4-β+
(1-x)4-β))/(Γ(5-β)cos(βπ/2))
我們定義數(shù)值解與精確解的最大誤差:
將空間步長(zhǎng)相對(duì)時(shí)間步長(zhǎng)取到很小,時(shí)間收斂階為:
Order1=log2(E(2τ,h)/E(τ,h))
將時(shí)間步長(zhǎng)相對(duì)空間步長(zhǎng)取到很小,空間收斂階為:
Order2=log2(E(τ,2h)/E(τ,h))
表1 誤差與時(shí)間收斂階(h=1/2000)
表2 誤差與空間收斂階(τ=1/800)
從表1、表2可看出,當(dāng)α=0.5,β=1.75時(shí),表1的時(shí)間收斂階為2.5階,表2的空間收斂階是2階。這樣的結(jié)論符合定理1,因此從誤差效果來(lái)看,我們所構(gòu)造的高階數(shù)值格式是有效的。
下面為了說(shuō)明本文所用算法所需時(shí)間的有效性,我們用直接LU分解方法和本文所用算法FGMRES進(jìn)行了比較。從表3可以看出,本文所用算法FGMRES的時(shí)間效率顯然比直接LU分解方法的時(shí)間效率要好。
表3 FGMRES與LU方法所用時(shí)間(單位:s)
以上所有結(jié)果均在MATLAB R2016a在Intel(R)Core(TM)i5-4590 CPU @3.30GHz 3.30GHz和8.00GB內(nèi)存的電腦中實(shí)現(xiàn).