二維調(diào)和分?jǐn)?shù)階擴(kuò)散方程的數(shù)值模擬

陳景華,陳雪娟,章紅梅

(1.集美大學(xué)理學(xué)院,福建 廈門(mén) 361021;2.福州大學(xué)數(shù)學(xué)與計(jì)算機(jī)科學(xué)學(xué)院,福建 福州 350108)

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

1 數(shù)值離散

Podlubny等[26]給出了左側(cè)和右側(cè)Riemann-Liouville (R-L)分?jǐn)?shù)階導(dǎo)數(shù)的定義：

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

在此基礎(chǔ)上給出α階左側(cè)和右側(cè)的R-L標(biāo)準(zhǔn)調(diào)和分?jǐn)?shù)階導(dǎo)數(shù)的定義[23-24].

定義1如果u(x)∈[a,b]，則?λ≥0,α階左側(cè)和右側(cè)R-L標(biāo)準(zhǔn)調(diào)和分?jǐn)?shù)階導(dǎo)數(shù)的定義分別用式(1)和(2)表示：

(1)

(2)

這里1<α≤2.

本文中考慮如下具有左側(cè)R-L標(biāo)準(zhǔn)調(diào)和分?jǐn)?shù)階導(dǎo)數(shù)的二維擴(kuò)散方程：

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

0≤t≤T.

(3)

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

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)

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

注2如果方程(3)中α=β=2,λ=0,則二維調(diào)和分?jǐn)?shù)階擴(kuò)散方程就退化成經(jīng)典的二維擴(kuò)散方程.

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

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

這里h是空間步長(zhǎng),權(quán)重系數(shù)為

則?x∈R存在不依賴(lài)于h、f、x和λ的常數(shù)aj使得：

O(hn).

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

(6)

(7)

這里

(8)

(9)

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

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

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

(11)

(12)

(13)

(14)

(15)

(16)

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

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

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

(17)

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

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

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

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

Ai,m=

(18)

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

Bi,m=

(19)

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

(20)

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

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

這里系數(shù)

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

(21)

(22)

2 穩(wěn)定性和收斂性證明

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

證明這個(gè)證明類(lèi)似文獻(xiàn)[27]. 主要的差別是Grünwald 權(quán)重系數(shù)由調(diào)和加權(quán)系數(shù)來(lái)替換。

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

證明主要采用圓盤(pán)定理和矩陣的特征值來(lái)證明穩(wěn)定性.

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

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

si,j=

(23)

注意到

容易計(jì)算得到

0，

(24)

3 數(shù)值例子

這節(jié)考慮平面矩形區(qū)域Ω=(0,1)×(0,1)上的二維調(diào)和分?jǐn)?shù)階擴(kuò)散方程

q(x,y,t)=

(25)

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

邊界條件為

u(0,y,t)=0,

u(x,0,t)=0,

精確解為

同時(shí)取誤差范數(shù)為最大誤差:

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

外推后的誤差率接近：

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

4 結(jié) 論

本文中發(fā)展了一個(gè)分?jǐn)?shù)階交替迭代法并結(jié)合C-N格式和Richardson外推法數(shù)值模擬二維調(diào)和分?jǐn)?shù)階擴(kuò)散方程，采用矩陣分析方法結(jié)合圓盤(pán)定理證明差分格式的穩(wěn)定性和收斂性. 此數(shù)值方法也可運(yùn)用于非線(xiàn)性源項(xiàng)的高維調(diào)和分?jǐn)?shù)階反常擴(kuò)散方程.