陳 劍， 曾泰山

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

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

本文研究如下時間分數(shù)階次擴散方程的初邊值問題:

(1)

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

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

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

1 預備知識

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

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

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

引理1[5]設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ù)階次擴散方程(1)的L1-多尺度Galerkin全離散格式,并推導其誤差估計.

其中,aj=(j+1)1-α-j1-α(j≥0). 根據(jù)文獻[2],如果uC2([0,T];L2(I)),則用逼近分數(shù)階導數(shù)的截斷誤差為：

(2)

(3)

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

(4)

(5)

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

(6)

(7)

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

再注意到式(7)左端第二項非負,則有

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

證畢.

3 求解全離散格式的多層擴充算法

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

(8)

(9)

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

為使用多層擴充算法,先將Kn進行分塊. 對m0,p,p′+,記

其中

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

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

4 數(shù)值算例

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

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

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

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

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

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

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

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

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

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

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

圖2 算法1與Gauss迭代法的計算時間對比

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

5 結束語