時(shí)間分?jǐn)?shù)階次擴(kuò)散方程的多層擴(kuò)充算法

陳 劍， 曾泰山

(1. 佛山科學(xué)技術(shù)學(xué)院數(shù)學(xué)與大數(shù)據(jù)學(xué)院, 佛山 528000； 2. 華南師范大學(xué)數(shù)學(xué)科學(xué)學(xué)院, 廣州 510631)

分?jǐn)?shù)階導(dǎo)數(shù)具有非局部性,相比整數(shù)階導(dǎo)數(shù)，它能更精確地描述現(xiàn)實(shí)世界中具有記憶和遺傳性質(zhì)的問題. 因此，分?jǐn)?shù)階微分方程已被廣泛應(yīng)用于眾多領(lǐng)域，如半導(dǎo)體物理、彈性力學(xué)、生物學(xué)異速生長、金融學(xué)期權(quán)定價(jià)、信號(hào)處理和分形理論等[1-2]. 但分?jǐn)?shù)階導(dǎo)數(shù)的非局部性使得這類方程難以獲得解析解,或者解析解需要用特殊函數(shù)來表達(dá),這不利于實(shí)際應(yīng)用. 所以,研究其數(shù)值解顯得尤為重要.

本文研究如下時(shí)間分?jǐn)?shù)階次擴(kuò)散方程的初邊值問題:

(1)

這里Γ(·)表示Gamma函數(shù).

時(shí)間分?jǐn)?shù)階次擴(kuò)散方程是一類重要的數(shù)學(xué)模型，常被用于描述運(yùn)輸中的反常擴(kuò)散行為，能很好地刻畫長時(shí)記憶的流通過程. 已有很多學(xué)者研究了有關(guān)時(shí)間分?jǐn)?shù)階次擴(kuò)散方程的數(shù)值方法[3-12]，如：給出了一種隱式差分格式,討論了格式的收斂階，證明了該格式無條件穩(wěn)定[3]；研究了二維問題的一種高階緊致差分格式[5]；給出了一種有限元方法，并證明了所得全離散格式具有超收斂性[4]；給出了一種時(shí)間方向具有二階收斂階的有限元計(jì)算格式，證明了該格式無條件穩(wěn)定且空間方向具有最優(yōu)收斂階[6]；利用經(jīng)典有限元方法討論了時(shí)間分?jǐn)?shù)階對(duì)流擴(kuò)散方程和電報(bào)方程的數(shù)值解[7-8]；給出了一種最小二乘譜方法[9]；運(yùn)用時(shí)空譜方法討論了方程(1)的數(shù)值解，并獲得了先驗(yàn)誤差估計(jì)[10]. 但上述方法并未涉及全離散所得線性方程組的快速求解問題. 事實(shí)上，對(duì)于大尺度、高精度的實(shí)際問題而言，全離散格式將導(dǎo)出一個(gè)大規(guī)模的線性方程組，而求解該方程組需要消耗巨大的計(jì)算量.

為了高效求解全離散所得線性方程組，本文研究了方程(1)的一種多尺度快速方法. 首先，基于時(shí)間方向采用L1逼近和空間方向采用多尺度Galerkin逼近建立全離散格式，給出了全離散格式的誤差估計(jì)；然后，利用矩陣分裂策略設(shè)計(jì)快速求解該全離散格式的多層擴(kuò)充算法,并證明了該算法具有最優(yōu)收斂階.

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

本節(jié)介紹一些相關(guān)記號(hào)和引理,以方便后文引用. 沒有特別說明,后文中出現(xiàn)的字母c表示與時(shí)間和空間步長無關(guān)的常數(shù),不同的位置可能表示不一樣的常數(shù).

Xn=Xn-1⊕⊥Wn=X0⊕⊥W1⊕⊥W2⊕⊥…⊕⊥Wn,

〈u-nu,v〉1=(?x(u-nu),?xv)=0(uXn).

引理1[5]設(shè)m,r,r≥1是正整數(shù),u如果1≤m≤r+1,則

‖u-nu‖p≤c2-(m-p)n‖u‖m(p=0,1).

2 L1-多尺度Galerkin全離散格式

本節(jié)給出時(shí)間分?jǐn)?shù)階次擴(kuò)散方程(1)的L1-多尺度Galerkin全離散格式,并推導(dǎo)其誤差估計(jì).

其中,aj=(j+1)1-α-j1-α(j≥0). 根據(jù)文獻(xiàn)[2],如果uC2([0,T];L2(I)),則用逼近分?jǐn)?shù)階導(dǎo)數(shù)的截?cái)嗾`差為：

(2)

(3)

其中,fk=f(x,tk).

(4)

(5)

接下來論證誤差估計(jì). 由方程(1)可知,真解uk滿足方程

(6)

(7)

由于1=a0>a1>a2>…>an>0,則由Cauchy-Schwarz不等式、三角不等式及引理2,有

再注意到式(7)左端第二項(xiàng)非負(fù),則有

從而由三角不等式及引理2，有

證畢.

3 求解全離散格式的多層擴(kuò)充算法

基底的多尺度特性使得全離散格式(3)對(duì)應(yīng)的線性方程組的系數(shù)矩陣具有高低頻層次結(jié)構(gòu),因此,本文利用多層擴(kuò)充法進(jìn)行高效求解. 為了方便敘述,先將全離散格式改寫成:

(8)

(9)

其中，En=[(w′ij,w′i′j′):(i,j),(i′,j′)Jn],Kn=[[(wij,]wi′j′):(i,j),(i′,j′)由基底的正交性可知En是單位矩陣.

為使用多層擴(kuò)充算法,先將Kn進(jìn)行分塊. 對(duì)m0,p,p′+,記

其中

記m0和m是2個(gè)固定的正整數(shù),m0<

算法1求解線性方程組(8)的多層擴(kuò)充算法

4 數(shù)值算例

本節(jié)以數(shù)值算例來驗(yàn)證本文的理論估計(jì)和算法1的計(jì)算效果. 考慮方程(1)，其中

u0(x)=2sin(2πx),

其真解u=(tα+2+t+2)sin(2πx).

為驗(yàn)證算法1的時(shí)間收斂階,選取二次多尺度正交基(基底具體構(gòu)造可參看文獻(xiàn)[14]),并在多層擴(kuò)充求解時(shí)取(m0,m)=(2,5),時(shí)間步長分別取為1/4、1/8、1/16、1/32、1/64、1/128,對(duì)不同的α進(jìn)行測(cè)試. 由所得結(jié)果(圖1)可以看到：數(shù)值結(jié)果與理論的時(shí)間收斂階2-α相吻合.

圖1 u2,5(x,t)在t=1時(shí)刻的時(shí)間收斂階

在空間方向,分別取線性基底和二次基底進(jìn)行數(shù)值計(jì)算,用算法1時(shí)，初始層m0分別取為4和2,m表示擴(kuò)充的層數(shù),x(n)表示相應(yīng)逼近子空間的維數(shù),n=m0+m. 收斂階為：

這里p可取1和0,分別對(duì)應(yīng)H1和L2范數(shù). 對(duì)于線性和二次基底,其H1范數(shù)的理論收斂階分別為1和2;其L2范數(shù)的理論收斂階分別為2和3. 由表1和表2可以看到：算法1所得的數(shù)值結(jié)果與理論收斂階非常吻合.

利用Matlab中的Cond命令,計(jì)算了系數(shù)矩陣En+Kn的條件數(shù),從表1和表2可知條件數(shù)一致有界.

表1 線性基底的數(shù)值結(jié)果Table 1 The numerical results for linear basis

表 2 二次基底的數(shù)值結(jié)果Table 2 The numerical results for quadratic basis

在計(jì)算效率方面,對(duì)比算法1與Gauss迭代法在線性基底下的數(shù)值結(jié)果. 由對(duì)比結(jié)果(圖2)可知：算法1具有更高的效率,而且問題規(guī)模越大,其計(jì)算優(yōu)勢(shì)越明顯.

圖2 算法1與Gauss迭代法的計(jì)算時(shí)間對(duì)比

Figure 2 The comparison of computational time between algorithm 1 and the Gauss iteration method

5 結(jié)束語