二維調和分數階擴散方程的數值模擬

陳景華,陳雪娟,章紅梅

(1.集美大學理學院,福建 廈門 361021;2.福州大學數學與計算機科學學院,福建 福州 350108)

分數階導數可模擬物理領域中的許多問題[1-5]，與離散隨機游走過程聯系緊密[6]. 一個階數α(α<2)的空間分數階導數可模擬具有長尾冪律性態的粒子跳躍,其跳躍步長的分布為P[J>x]≈x-α，而一個階數β<1的時間分數階導數可模擬長尾冪律性態粒子的跳躍,其兩次跳躍的等候時間為P[W>t]≈t-β.本文中我們關注的是調和分數計算，其分數階導數是通常的分數階導數乘上一個指數因子的卷積核.這個調和分數階導數可模擬調和的指數冪律，粒子的跳躍長度為P[J>x]≈x-αe-λx[7-8]. 調和擴散方程已廣泛應用于地球物理[9-12]和金融[13-14]等領域.在金融領域里，調和穩定的過程可模擬具有半長尾性態的價格截斷.近年來關于分數階數值方法已有不少的成果[15-22]，但是關于調和分數階擴散方程的數值方法還比較少.Sabzikar 等[23]給出一個調和分數階差分格式求解調和分數階擴散方程.Baeumer等[24]給出粒子追蹤法和差分法來求解具有漂流項的調和分數階方程.Zhang等[25]發展了一個新的數值算法求解Riesz調和空間分數階擴散方程.以上所考慮的方程都是一維調和分數階方程. 本文中考慮二維調和分數階方程，采用隱式交替方向法(ADI)和Crank-Nicolson(C-N)離散格式求解二維調和分數階方程，并證明所給出的差分格式是無條件穩定和收斂的.

1 數值離散

Podlubny等[26]給出了左側和右側Riemann-Liouville (R-L)分數階導數的定義：

這里u(x)定義在[a,b]區間，1<α≤2,Γ(·)是Gamma函數.

在此基礎上給出α階左側和右側的R-L標準調和分數階導數的定義[23-24].

定義1如果u(x)∈[a,b]，則?λ≥0,α階左側和右側R-L標準調和分數階導數的定義分別用式(1)和(2)表示：

(1)

(2)

這里1<α≤2.

本文中考慮如下具有左側R-L標準調和分數階導數的二維擴散方程：

(x,y)∈Ω=(xL,xR)×(yL,yR),

0≤t≤T.

(3)

其中：1<α≤2,1<β≤2; 擴散系數:d(x,y)>0，e(x,y)>0;q(x,y,t)是源項. 方程(3)可模擬長時間隨機游走極限,其平面跳躍的隨機點(X,Y)坐標的分布為：P[X>x]≈x-αe-λ1x和P[Y>y]≈y-βe-λ2y，且X,Y是相互獨立的. 假設這個分數階擴散方程(3)在以下初始和邊界條件下有光滑的唯一解,并設初始條件為：

u(x,y,0)=f(x,y),(x,y)∈Ω,

(4)

Dirichlet邊界條件為：

u(x,y,t)=B(x,y,t),(x,y)∈?Ω,

0≤t≤T.

(5)

這里加上一個限制條件B(xL,y,t)=B(x,yL,t)=0. 這個條件的設置是有實際意義的：可以設想平面區域足夠大時，溶質在地下水層的擴散，在左邊界和下邊界濃度為0.

注2如果方程(3)中α=β=2,λ=0,則二維調和分數階擴散方程就退化成經典的二維擴散方程.

接下來采用左移位的Grünwald估計和ADI,建立方程(3)的二階C-N格式.

引理1[24]設1<α≤2,f∈Wn+3,1(R)，n是整數且n≥1.設p∈R,h>0,λ≥0.定義左移位的 Grünwald-Letnikov(G-L)調和算子:

這里h是空間步長,權重系數為

則?x∈R存在不依賴于h、f、x和λ的常數aj使得：

O(hn).

αλα-1f′(xj)+O(h),

(6)

(7)

這里

(8)

(9)

方程(10)可寫成算子的形式：

0≤i≤Nx,0≤j≤Ny,

n=0,1,…,N-1.

(11)

(12)

(13)

(14)

(15)

(16)

以下是求解式(15)、(16)的具體步驟. 首先將式(15)展開寫成：

i=1,2，…,Nx-1,

固定j=k(k=1,2,…,Ny-1).

(17)

從式(17)可以得到一個Nx-1個的線性方程，這里第i個方程表示如下

i=1,2，…,Nx-1.

系數Ai,m，(i=1,2,…,Nx-1,m=0,1,…,i+1) 定義如下：

小兒臟腑柔弱,氣血未充,腠理不固,易為外邪侵襲。外感時邪,正邪相搏而致外感發熱。熱極易生風,加之小兒肝常有余,肝風內動易致驚風。風性善行而數變,故小兒病情發展變化迅速,合并癥多,應盡快予以診治,以免貽誤病情。根據中醫“異病同治”的治療原則,由各種不同疾病引起的發熱經中醫辨證論治后都可選用相應的中成藥物治療。而西醫多采用靜脈給予抗感染抗病毒藥物,同時給予口服退熱藥或冰袋物理降溫等治療,雖亦有效,但因小兒畏懼打針吃藥,往往不配合醫護人員的操作,增加了治療難度,而選擇熱毒寧注射液灌腸則為臨床醫護人員提供了方便,減輕患兒的痛苦。

Ai,m=

(18)

系數Bi,m(i=1,2，…,Nx-1,m=0,1，…,k+1)定義如下:

Bi,m=

(19)

類似的,固定i=k(k=1,2，…,Nx-1),j=1,2，…,Ny-1，把方程 (16)重新寫為:

(20)

因此上述的方程組就是求解Ny-1個線性方程,在(xk,yj)處,第j個方程形式如下：

j=1,2,…,Ny-1,

這里系數

m=0,1,…,j+1) 分別定義如下

(21)

(22)

2 穩定性和收斂性證明

這一節，證明方程(3)的差分格式(10)是相容且無條件穩定的，因此是收斂的.

證明這個證明類似文獻[27]. 主要的差別是Grünwald 權重系數由調和加權系數來替換。

定理2方程(3)的差分格式(10)是無條件穩定的.

證明主要采用圓盤定理和矩陣的特征值來證明穩定性.

方程(12)可寫成矩陣的形式

(I-S)(I-T)Un+1=(I+S)(I+T)Un+Rn+1,

si,j=

(23)

注意到

容易計算得到

0，

(24)

3 數值例子

這節考慮平面矩形區域Ω=(0,1)×(0,1)上的二維調和分數階擴散方程

q(x,y,t)=

(25)

初始條件參數設置為λ1=λ2=0.5,α=β=1.7,k1=k2=3，

邊界條件為

u(0,y,t)=0,

u(x,0,t)=0,

精確解為

同時取誤差范數為最大誤差:

表1給出了外推前的最大誤差、誤差率和外推后的最大誤差、誤差率的比較. 可以觀察到外推前的誤差率接近：

外推后的誤差率接近：

表1 外推前后最大誤差、誤差率的比較

4 結 論

本文中發展了一個分數階交替迭代法并結合C-N格式和Richardson外推法數值模擬二維調和分數階擴散方程，采用矩陣分析方法結合圓盤定理證明差分格式的穩定性和收斂性. 此數值方法也可運用于非線性源項的高維調和分數階反常擴散方程.