陳雪娟，陳景華,2*，章紅梅

(1.集美大學(xué)理學(xué)院,福建 廈門 361021 ；2.集美大學(xué)數(shù)字福建大數(shù)據(jù)建模與智能計算研究所,福建 廈門 361021；3.福州大學(xué)數(shù)學(xué)與計算機(jī)科學(xué)學(xué)院,福建 福州 350108)

近年來，研究者們在湍流、核磁共振、多孔介質(zhì)中的滲透等系統(tǒng)中都觀察到偏離布朗運(yùn)動特性的擴(kuò)散行為，而分?jǐn)?shù)階微分方程能夠較精確地描述具有記憶和遺傳、路徑依賴性質(zhì)的物理過程[1]，因此分?jǐn)?shù)階微分方程目前已被用來描述諸多領(lǐng)域[2-4]中的反常擴(kuò)散問題.從現(xiàn)實(shí)問題中抽象出分?jǐn)?shù)階微分方程之后，如何對分?jǐn)?shù)階微分方程進(jìn)行求解成為一個迫切的研究課題.目前已有很多文獻(xiàn)討論了分?jǐn)?shù)階微分方程的解析解，但是這些解大多含有特殊函數(shù)(比如多變量的Mittag-Leffler函數(shù)[5-9])，要計算這些特殊函數(shù)相當(dāng)困難，而且并非所有的分?jǐn)?shù)階微分方程都能得到其解析解，因此研究分?jǐn)?shù)階方程的數(shù)值模擬有著重要的理論與應(yīng)用價值.國內(nèi)外對分?jǐn)?shù)階微分方程數(shù)值解的研究也有一定的進(jìn)展，主要采用的數(shù)值方法是有限差分[10-14]、有限元、譜方法[15-17]等；但是采用二次多項(xiàng)式樣條函數(shù)進(jìn)行數(shù)值逼近的研究文獻(xiàn)卻較缺乏.Sousa[18]研究線性分?jǐn)?shù)階擴(kuò)散方程，借助樣條函數(shù)近似離散分?jǐn)?shù)階導(dǎo)數(shù)，并在時間上采用Crank-Nicolson離散.Rashidinia等[19]提出一種新的非多項(xiàng)式樣條方法數(shù)值求解二階雙曲線方程.Siraj等[20]構(gòu)造了三次非多項(xiàng)式樣條函數(shù)來求解非線性薛定諤方程.陳雪娟等[21]采用基于二次樣條函數(shù)的數(shù)值方法求解分?jǐn)?shù)階擴(kuò)散方程，但是處理的是單側(cè)分?jǐn)?shù)階導(dǎo)數(shù).較少文獻(xiàn)涉及用樣條方法數(shù)值求解雙側(cè)分?jǐn)?shù)階導(dǎo)數(shù).由于雙側(cè)分?jǐn)?shù)階導(dǎo)數(shù)涉及到區(qū)間兩個邊界的函數(shù)值，比單側(cè)分?jǐn)?shù)階導(dǎo)數(shù)能更好地模擬反常擴(kuò)散現(xiàn)象，所以多用于定義空間變量的分?jǐn)?shù)階導(dǎo)數(shù).

本文提出一種基于多項(xiàng)式樣條函數(shù)的數(shù)值方法，用來求解具有非線性源的雙側(cè)空間分?jǐn)?shù)階擴(kuò)散方程，同時證明了該數(shù)值方法是無條件穩(wěn)定和收斂的.

考慮具有非線性源項(xiàng)的雙側(cè)空間分?jǐn)?shù)階擴(kuò)散方程：

1<α≤2,a

(1)

初始條件：

u(x,0)=g(x),a

(2)

具有Dirichlet零邊界條件：

u(a,t)=u(b,t)=0, 0≤t≤T.

(3)

這里擴(kuò)散系數(shù)λ1(x,t)>0和λ2(x,t)>0.非線性源項(xiàng)f(u,x,t)滿足局部Lipschitz連續(xù)條件[22-23], 即存在某個常數(shù)L>0，使得

‖f(u,x,t)-f(v,x,t)‖≤L‖u-v‖.

(4)

α階的左、右側(cè)Caputo導(dǎo)數(shù)分別定義如下[24]：

(5)

(6)

1 樣條配置法

1.1 二次多項(xiàng)式樣條函數(shù)

考慮如下二次多項(xiàng)式樣條函數(shù)Pi(x,tj)(xi≤x

Pi(x,tj)=ai(tj)(x-xi)2+bi(tj)(x-xi)+

ci(tj),i=0,1,2,…,m-1.

(7)

為了確定函數(shù)Pi(x,tj)的系數(shù)表達(dá)式，定義：

(8)

(9)

(10)

由方程(7)～(9)可以計算出：

(11)

再由方程(7)和(10)，及Caputo分?jǐn)?shù)階導(dǎo)數(shù)的定義可以推出：

(12)

從而可以確定方程(7)的系數(shù)為：

(13)

1.2 構(gòu)造數(shù)值方法

為了使二次函數(shù)pi(x,t)在x=xi,i=0,1,2,…，m-1處滿足連續(xù)性條件

利用方程(13), 得到如下等式：

(14)

進(jìn)一步，可以計算出等式(14)的局部截斷誤差.

O(h2).

(15)

(16)

這里,

(17)

以下建立分?jǐn)?shù)階擴(kuò)散方程(1)～(3)的差分格式.

由方程(1)可得：

(18)

對時間導(dǎo)數(shù)采用向后差分格式，得到：

o(τ)，

(19)

及

(20)

從方程(14)和(18)可得以下方程：

(21)

(22)

(23)

(24)

因此得到分?jǐn)?shù)階擴(kuò)散方程(1)～(3)的數(shù)值逼近格式(22)～(24)，且此數(shù)值格式的收斂階為O(hα+τ).

2 數(shù)值方法的穩(wěn)定性和收斂性分析

本節(jié)討論數(shù)值方法的收斂性和穩(wěn)定性.將方程組(22)寫成更簡單的形式：

(25)

為了證明數(shù)值方法的穩(wěn)定性和收斂性，引入引理2.

引理2離散的Gronwall不等式[25].假設(shè)fk≥0，ζk≥0，k=0,1,2,…,并且ζk+1≤μζk+τfk,μ=1+C0τ,k=0,1,2,…,ζ0=0,這里C0≥0是常數(shù), 則

2.1 穩(wěn)定性分析

由迭代格式(25), 得到誤差方程：

(26)

對j=1,2,…,n, 分別定義網(wǎng)格函數(shù)[26]

則εj(x)和ρj(x) 可分別展開成Fourier級數(shù)：

其中系數(shù)為

及

有

這里j=0,1,…,n.

(27)

定理1數(shù)值離散格式(22)～(24)是無條件穩(wěn)定的.

證明將式(27)代入迭代格式(26), 得到

則有

經(jīng)過簡單計算得到下列不等式：

因此有

|ξj|≤(|ξj-1|+τ|ηj-1|).

(28)

則|ηj|≤L|ξj|.因此

|ξj|≤(1+τL)|ξj-1|，j=1,2,…,n.

(29)

運(yùn)用上式歸納得到

‖εn‖2≤(1+Lτ)n‖ε0‖2≤eLT‖ε0‖2.

因此，證明了數(shù)值離散格式(25)是無條件穩(wěn)定的.

2.2 收斂性分析

根據(jù)迭代格式(25)的相容性，得到以下誤差方程：

(30)

j=0,1,…,n.

(31)

定理2數(shù)值離散格式(22)～(24)是無條件收斂的，并且存在正常數(shù)C>0使得：

j=1,2,…,n.

證明將式(31)代入到迭代格式(30)中，得到：

經(jīng)過簡單的計算得到

2,…,n.

(32)

這里C1是正常數(shù). 從而

‖ej‖2≤(1+τL)‖ej-1‖2+C1(hα+τ).

運(yùn)用離散的Gronwall不等式(引理2)得到：

C1TeLT(hα+τ)=C(hα+τ)，

(33)

這里C=C1TeLT.

因此，證明了數(shù)值離散格式(22)～(24)是無條件收斂的，且收斂為o(hα+τ).

3 分?jǐn)?shù)階行方法

一般分?jǐn)?shù)階微分方程的精確解是很難求得，因此為了說明數(shù)值方法的有效性，通常將所得的數(shù)值解結(jié)果與分?jǐn)?shù)階行方法(method of line,MOL)的解進(jìn)行比較.分?jǐn)?shù)階MOL其本質(zhì)是只離散空間變量，將分?jǐn)?shù)階偏微分方程轉(zhuǎn)換成一組常微分方程加以求解.繼而可以采用解決微分代數(shù)系統(tǒng)的工具DASSL來求解.2004年Liu等[27]首次運(yùn)用分?jǐn)?shù)階MOL求解空間分?jǐn)?shù)階Fokker-Planck偏微分方程.

2u(xi-k+1,tj)+u(xi-k,tj)]}+

2u(xk,tj)+u(xk-1,tj)]}+

u(xk-1,tj)]}+o(h2).

由此得到左、右側(cè)Caputo分?jǐn)?shù)階導(dǎo)數(shù)的逼近格式：

2u(xi-k+1,tj)+u(xi-k,tj)]},

2u(xk,tj)+u(xk-1,tj)]}.

因此，具有非線性源項(xiàng)的雙側(cè)空間分?jǐn)?shù)擴(kuò)散方程(1)～(3)的MOL可以寫成以下形式：

(34)

4 數(shù)值例子

本節(jié)給出兩個數(shù)值例子來證明所提出的理論分析的有效性.

例1考慮具有非線性源項(xiàng)的分?jǐn)?shù)階擴(kuò)散方程：

(35)

這里系數(shù)

λ1(x,t)=0.5Γ(5-α)x2+α(1-x)4et,

λ2(x,t)=0.5Γ(5-α)x4(1-x)2+αet,

非線性源項(xiàng)為f(x,t,u)=u(x,t)+u2(x,t)·[-24x2+24x-12+2α(4-α)].

圖1是方程的精確解、數(shù)值格式(22)～(24)的數(shù)值解及MOL的數(shù)值解在t=1.0時刻的比較.取空間步長h=0.05，時間步長τ=0.01，α=1.5.

圖1 t=1.0時刻，取α=1.5，h=0.05，τ=0.01，數(shù)值解與 分?jǐn)?shù)階MOL的解及精確解的比較Fig.1 Comparison of numerical solutions using our method, the MOL, and the exact solution at time t=1.0，α=1.5,h=0.05，τ=0.01

表1顯示了在t=1.0時刻，α=1.5，不同的時間步長和空間步長(空間步長和時間步長的關(guān)系為h≈τ-α)下精確解與數(shù)值解的最大誤差和誤差率.可以看出

表1 在t=1.0時，取α=1.5，數(shù)值解與 精確解的最大誤差及誤差率Tab.1 The maximum errors and error rates for the numerical method relative to the exact solution at time t=1.0，α=1.5

說明數(shù)值解的收斂階為o(hα+τ)，這與本文中理論分析結(jié)果是相一致的.

從圖1和表1中可以看出，基于二次多項(xiàng)式樣條函數(shù)的數(shù)值方法與精確解很好地吻合.

例2考慮以下具有非線性源項(xiàng)的分?jǐn)?shù)階擴(kuò)散方程：

(36)

這里f(x,t,u)=u(1-u)+x2，此題精確解難以求得.因此考慮將數(shù)值解與分?jǐn)?shù)階MOL的數(shù)值解進(jìn)行比較.

圖2顯示在t=1.0時刻, 取α=1.5，h=0.05，τ=0.01時數(shù)值解與分?jǐn)?shù)階MOL數(shù)值解的比較.圖3顯示了當(dāng)α=1.5，h=0.05，τ=0.01時， 在不同時刻t=0.3，0.5,1.0的數(shù)值解的圖形.圖4顯示了在t=1.0時刻取不同的分?jǐn)?shù)階導(dǎo)數(shù)值α=1.5,1.8，2.0時數(shù)值解的圖形.

圖2 t=1.0時刻，取α=1.5，h=0.05，τ=0.01，數(shù)值解與 分?jǐn)?shù)階MOL數(shù)值解的比較Fig.2 Comparison of numerical solutions using our method and the fractional MOL at time t=1.0,α=1.5,h=0.05,τ=0.01

圖3 取α=1.5，h=0.05，τ=0.01，本文方法 在不同時刻的數(shù)值解比較Fig.3 Comparison of numerical solutions using our method at different times when α=1.5，h=0.05，τ=0.01

5 結(jié) 論

本文提出了一種基于多項(xiàng)式樣條函數(shù)的數(shù)值離散格式，用于求解具有非線性源項(xiàng)的雙側(cè)空間分?jǐn)?shù)階擴(kuò)散方程，其中的分?jǐn)?shù)階導(dǎo)數(shù)采用Caputo分?jǐn)?shù)階導(dǎo)數(shù)；同時，提出一種計算分?jǐn)?shù)階微分方程數(shù)值解的計算技巧，即分?jǐn)?shù)階MOL.證明了提出的數(shù)值離散格式是無條件穩(wěn)定和收斂的，且收斂階為o(hα+τ)(1<α≤2).最后，給出數(shù)值例子證明了方法的有效性.這種樣條方法和分析技術(shù)可用于求解和分析其他類型的分?jǐn)?shù)階偏微分方程.