時空分數階擴散方程的高階快速數值方法分析*

田洪材,曹俊英,王自強

(貴州民族大學 數據科學與信息工程學院,貴州 貴陽 550025)

1 背景

分數階微積分應用范圍廣泛,可以描述一些具有記憶性、非局部性的自然現象。分數階導數和分數階微分方程被許多科學和工程領域的研究人員所研究[1]。其中,具有空間Riesz分數階導數和時間Caputo分數階導數的時空分數擴散方程涉及許多應用,如金融學[2]和分數動力學方程[3]。隨著分數階微分方程模型的增多和精確解的不易獲得,人們對開發快速的數值方法產生了極大的興趣。

Riesz分數階導數被視為左右Riemann-Liouville分數階導數的線性組合。本文研究時間為Caputo分數階導數,空間為Riesz分數階導數的時空分數階擴散方程,對它的一個高階快速的數值格式進行了構造和分析。

2 數值格式構造

考慮如下方程:

(1)

和初邊值條件

u(x,0)=φ(x),x∈[a,b],

u(a,t)=u(b,t)=0,?t∈[0,T],

0<α<1,1<β<2,

這里Γ(·)是伽馬函數,f(x,t)和φ(x)為已知函數。

(2)

其中:

(3)

其中,kα=Γ(3-α)τα,pα=c1+2-0.5α。

當m=3時,

當m≥4時,

其中

(j+1)2-α-j2-α

其中U0=[φ1,φ2,…,φN-1]。

利用文獻[6]中的加權移位Grünwald-Letnikov差分方法在空間上和時間上利用式(2)與式(3),可得:

當m=1,2時,有:

(4)

當m≥3時,有:

(5)

其中,

D=

這里

wβ,0=λ1gβ,0,wβ,1=λ1gβ,1+λ2gβ,0,

wβ,k=λ1gβ,k+λ2gβ,k-1+λ3gβ,k-2,k≥2。

當m≥3時,數值格式(5)系數矩陣是一個Toeplitz矩陣,我們利用文獻[7]和[8]的思想方法,結合GMRES迭代方法與FFT方法建立數值格式(5)的快速迭代方法(FGMRES)。

利用文獻[9]中的思想方法,我們可以證明如下定理。

定理1:數值格式(5)的收斂階為O(h2+τ3-α)。

3 數值實驗

考慮一維時空分數階擴散方程模型,我們有[a,b]=[0,1],T=1,該方程的精確解和右端項分別為: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))

我們定義數值解與精確解的最大誤差:

將空間步長相對時間步長取到很小,時間收斂階為:

Order1=log2(E(2τ,h)/E(τ,h))

將時間步長相對空間步長取到很小,空間收斂階為:

Order2=log2(E(τ,2h)/E(τ,h))

表1 誤差與時間收斂階(h=1/2000)

表2 誤差與空間收斂階(τ=1/800)

從表1、表2可看出,當α=0.5,β=1.75時,表1的時間收斂階為2.5階,表2的空間收斂階是2階。這樣的結論符合定理1,因此從誤差效果來看,我們所構造的高階數值格式是有效的。

下面為了說明本文所用算法所需時間的有效性,我們用直接LU分解方法和本文所用算法FGMRES進行了比較。從表3可以看出,本文所用算法FGMRES的時間效率顯然比直接LU分解方法的時間效率要好。

表3 FGMRES與LU方法所用時間(單位：s)

以上所有結果均在MATLAB R2016a在Intel(R)Core(TM)i5-4590 CPU @3.30GHz 3.30GHz和8.00GB內存的電腦中實現.