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

嵇雯蕙

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

1 引言

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

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

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

2 問題的模型

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

(2.1)

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

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

則

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

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

則

由

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

3 有限差分格式的建立

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

將區間[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對應的隱式有限差分格式如下:

為了將數值格式改寫為矩陣形式, 令

其中

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

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

(3.1)

其中系數矩陣A(j)為

(3.2)

其中

(3.3)

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

我們將式 (3.3) 代入 (3.2) 化簡后, 系數矩陣A(j)變為

(3.4)

4 問題的求解

本節主要任務是求解線性方程組 (3.1). 由于系數矩陣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)開始, 構造一組相互正交且范數為1的向量v(1),v(2),…,v(m).

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

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

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

從上述過程可以發現:

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

(1)帶狀矩陣

此時,預處理矩陣M為

(4.2)

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

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

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

圖2 α=0.9

(2)改進的帶狀矩陣

由公式(3.4), 我們注意到系數矩陣是由對稱正定矩陣和低秩矩陣L構成的. 為了構造一個更逼近系數矩陣的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循環預處理矩陣

(4.3)

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

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

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

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

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

圖5 α=0.6

圖6 α=0.9

(4)T. Chan預處理矩陣

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

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

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

圖7 α=0.6

圖8 α=0.9

5 數值實驗

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

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

其中

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

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

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

6 總結

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

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