一類空間分數(shù)階擴散方程的預處理方法*

嵇雯蕙

(昆明理工大學理學院,昆明,650000)

1 引言

隨著自然科學和生產技術的不斷發(fā)展, 微分方程逐漸成為現(xiàn)代科學技術中分析問題和解決問題的一個強有力的工具. 近幾十年來, 分數(shù)階微分方程在水文、生物、物理、化學、金融等學科有著廣泛的應用[1-5]. 分數(shù)階微分方程之所以能夠在諸多領域得到廣泛的應用, 主要是因為分數(shù)階微分算子具有非局部特征, 從而能更精確地反映一些物理現(xiàn)象.

在空間擴散模型中, 用分數(shù)階導數(shù)代替空間擴散二階導數(shù), 將導致更強的擴散. 最近, Mao and Shen[6]考慮了一類基于階數(shù)α∈(1/2,1)左右導數(shù)的分數(shù)階微分算子的偏微分方程, Xu, Sun, and Sheng[7]從變分的角度考慮了平衡中心分數(shù)階導數(shù), 它們保持了所需的變分性質以及L2框架上的對稱雙線性形式. 由于分數(shù)階偏微分方程的解析解往往很難用簡單的函數(shù)表示, 因此研究分數(shù)階偏微分方程的數(shù)值方法和理論分析就顯得尤為重要.

本文的主要目的是研究求解空間分數(shù)階擴散方程的有限差分格式離散系統(tǒng)的快速迭代方法.

2 問題的模型

我們考慮的空間分數(shù)階擴散方程(FDEs)如下[9]：

(2.1)

為了對方程 (2.1) 建立差分格式, 我們先給出如下幾個引理.

引理1[12]假設u(x)∈C2[xL,xR], 令

則

利用同樣的數(shù)值微分方法, 可以得到右Caputo分數(shù)導數(shù)的以下推論.

推論1假設u(x)∈C2[xL,xR], 令

則

由

以及引理1和推論1, 我們可以得到Riemann-Liouville分數(shù)導數(shù)的差分逼近.

3 有限差分格式的建立

類似文獻[14]的離散方法, 我們把方程 (2.1) 離散如下.

將區(qū)間[xL,xR]均勻分為n+1個等分, 空間步長Δx=(xR-xL)/(n+1), 則xi=xL+iΔx, 其中i=0,1,…,n+1. 同理, 將[0,T]均勻分為m個等分, 時間步長Δt=T/m, 則tj=jΔt, 其中j=0,1,…,m.

令

通過引理1和2以及推論1, (2.1) 中FDEs對應的隱式有限差分格式如下:

為了將數(shù)值格式改寫為矩陣形式, 令

其中

于是, 我們得到有限差分格式的矩陣表達式

A(j)u(j)=ηu(j-1)+Δx2α(f(j)-σ(j)),j=1,2,…,m,

(3.1)

其中系數(shù)矩陣A(j)為

(3.2)

其中

(3.3)

關于 (2.1) 離散的詳細情況, 見文獻[14].

我們將式 (3.3) 代入 (3.2) 化簡后, 系數(shù)矩陣A(j)變?yōu)?/p>

(3.4)

4 問題的求解

本節(jié)主要任務是求解線性方程組 (3.1). 由于系數(shù)矩陣A(j)是非對稱的, 我們考慮應用Krylov子空間方法的GMRES法來求解.

下面先簡單地介紹一下GMRES方法[15].

4.1 GMRES法

設要求解的線性方程組為

Ax=b.

任取一n維實向量x(0), 令x=x(0)+z, 則上述方程組等價于

Az=r(0),

(4.1)

其中r(0)=b-Ax(0). 故下面直接討論(4.1)的求解問題.

從r(0)開始, 構造一組相互正交且范數(shù)為1的向量v(1),v(2),…,v(m).

類似地, 計算v(3)的過程如下:

(1)計算u=Av(2),h12=(u,v(1)),h22=(u,v(2));

繼續(xù)這一過程, 經(jīng)過m步, 即可得到矩陣Vm=[v(1),v(2),…,v(m)].

從上述過程可以發(fā)現(xiàn):

Av(1)=h11v(1)+h21v(2),

Av(2)=h12v(1)+h22v(2)+h32v(3),

…

Av(m)=h1mv(1)+h2mv(2)+h3mv(3)+…+hm+1,mv(m+1).

寫成矩陣形式, 即

其中

假設存在一個y∈Rn滿足z=Vmy, 則有

分解得到.

具體算法如下[15]：

GMRES算法1. 計算r0=b-Ax0, β∶=‖r0‖2, v1∶=r0/β2. For j=1,2,…,m, Dowj∶=AvjFor i=1,…,j, Dohij∶=(wj,vi)wj∶=wj-hijviEnd Dohj+1,j=‖wj‖2. If hj+1,j=0 set m∶=j and go to 3vj+1=wj/hj+1,jEnd Do3. 定義(m+1)×m的Hessenberg矩陣H-m=hij 1≤i≤m+1,1≤j≤m4. 計算ym, 極小化‖βe1-H-my‖2, 且xm=x0+Vmym.

4.2 預條件子的建立

為了提高收斂速度, 我們建立預處理的GMRES方法. 針對系數(shù)矩陣的結構,我們分別構造帶狀矩陣、改進的帶狀矩陣、Strang循環(huán)矩陣以及T. Chan循環(huán)矩陣作為預處理矩陣M來逼近矩陣A(j), 使得M-1A(j)的特征值聚集在1附近, 從而達到收斂速度提高的效果. 四種不同的預處理矩陣的構造如下:

(1)帶狀矩陣

此時,預處理矩陣M為

(4.2)

易證此矩陣為三對角陣. 而我們知道三對角陣的線性方程組容易求解, 可以利用“追趕法”求出來, 所以M的構造合理.

為了更直觀地說明預處理后矩陣譜的聚集情況, 基于不同的α, 我們分別給出預處理前后系數(shù)矩陣特征值分布的對比圖如下(其中n=1024).

由圖1和圖2可見, 對于不同的α, 預處理后矩陣譜的普遍聚集在1附近, 只有個別點距1較遠.

圖2 α=0.9

(2)改進的帶狀矩陣

由公式(3.4), 我們注意到系數(shù)矩陣是由對稱正定矩陣和低秩矩陣L構成的. 為了構造一個更逼近系數(shù)矩陣的M, 對公式(4.2)作如下改進:

因為一個可逆矩陣加上一個低秩矩陣的逆也是容易求的, 若記

令

L=XGY,

其中X是n×r,Y是r×n,G是r×r的非奇異矩陣,r是L的秩，則

M-1=C-1-C-1X(G-1+YC-1X)-1YC-1.

故這里的M構造合理.

當預條件子為改進的帶狀矩陣時, 預處理前后的特征值分布圖如下.

由圖3和圖4可見, 采用改進后的帶狀矩陣, 個別遠點消失, 所有的特征值都聚集在1附近.

圖3 α=0.6

圖4 α=0.9

(3)Strang循環(huán)預處理矩陣

(4.3)

對于實Toeplitz矩陣B=[bj-k]0≤j,k

s(B)=[sj-k]0≤j,k

其中,sj是Strang預條件子對角線上的元素, 且有：

相對于文獻[14]中的預處理矩陣, 此處的預處理矩陣更逼近A(j). 當預條件子為Strang預處理矩陣 (4.3) 時, 預處理前后系數(shù)矩陣的特征值分布圖如下.

通過圖5和圖6, 我們可以清晰地看到基于該預條件子預處理后系數(shù)矩陣的特征值更集中且在1附近.

圖5 α=0.6

圖6 α=0.9

(4)T. Chan預處理矩陣

T. Chan預條件子c(B)=[ck-j]0≤k,j

當預條件子為T. Chan預處理矩陣時, 預處理前后系數(shù)矩陣的特征值分布圖如下.

由圖7和圖8可見, 預處理后系數(shù)矩陣的特征值顯然在1周圍集中, 故該預條件子選取合理.

圖7 α=0.6

圖8 α=0.9

5 數(shù)值實驗

為了驗證算法的有效性, 我們考慮如下的方程[14]

其中擴散系數(shù)為d+(x,t)=(1-x)α,d-(x,t)=xα,源項為

其中

方程的精確解為u(x,t)=e-tx2(1-x)2.

表1 四種不同的預條件子的比較

結果表明, 在不進行預處理的情況下, GMRES方法收斂需要多次迭代且隨著n變大迭代次數(shù)增多, 而預處理的結果迭代次數(shù)更少更穩(wěn)定且條件數(shù)大大減少. 以上的預條件子中, Strang和T. Chan循環(huán)預處理條件子處理效果最好.

6 總結

本文采用預處理的GMRES方法, 對一類FDEs(2.1)產生的離散線性系統(tǒng)進行求解. 數(shù)值試驗證明了所提出的預處理方法的有效性.

但是, 如果在由FDEs(2.1)產生的離散化方案中考慮其變分性質或有限體積法, 那么得到的線性系統(tǒng)有可能是對稱正定的. 然后, 可以采用共軛梯度法來有效地解決這些問題. 因此, 在今后的研究中, 我們應進一步研究這類FDEs的數(shù)值格式以及相關的預條件.