時間分?jǐn)?shù)階變系數(shù)對流擴(kuò)散方程的數(shù)值解法

郭非凡，張新東，王碩

(新疆師范大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，新疆 烏魯木齊 830017)

近年來，自然界發(fā)現(xiàn)越來越多的現(xiàn)象無法用傳統(tǒng)的整數(shù)階方程來描述。為了解決現(xiàn)實問題，分?jǐn)?shù)階微積分的研究逐漸引起了人們的廣泛關(guān)注。在物理、化學(xué)以及生物等學(xué)科領(lǐng)域，分?jǐn)?shù)階微積分應(yīng)用十分廣泛[1]，例如擴(kuò)散和輸運理論、混沌與湍流、生物組織、高分子材料的解鏈等?？梢哉f，分?jǐn)?shù)階導(dǎo)數(shù)早已成為了描述各種復(fù)雜力學(xué)、生物及物理行為等的重要工具[2]。目前，國內(nèi)外很多學(xué)者都已開始進(jìn)行分?jǐn)?shù)階微積分方程理論以及近似算法的研究。大量研究表明，在分?jǐn)?shù)階偏微分方程的近似計算中，有限差分方法和級數(shù)逼近法占據(jù)主要地位。與此同時，將多種方法結(jié)合運用來得到精度高并且足夠穩(wěn)定的算法也是一個研究熱點。

對流擴(kuò)散方程的研究大多在常系數(shù)或者整數(shù)階的范圍之內(nèi)，為了能更加精確地描述溶質(zhì)的運動特征，將其推廣到變系數(shù)的情形。文獻(xiàn)[3-5]研究了分?jǐn)?shù)階對流擴(kuò)散方程的數(shù)值解法，而對于變系數(shù)的時間分?jǐn)?shù)階對流擴(kuò)散方程的研究較少。然而，這類問題的解析解大多都是由較為復(fù)雜的函數(shù)構(gòu)成[6-8]，求解這類復(fù)雜函數(shù)較為困難。因此，研究分?jǐn)?shù)階微分方程的穩(wěn)定而高效的數(shù)值方法[9-13]備受國內(nèi)外學(xué)者的關(guān)注。

1 預(yù)備知識

定義1[14]Caputo分?jǐn)?shù)階導(dǎo)數(shù)為

用有限差分方法來求解式(1)的時間分?jǐn)?shù)階變系數(shù)對流擴(kuò)散方程，

(1)

其中，α是Caputo分?jǐn)?shù)階導(dǎo)數(shù)的階數(shù)，滿足0<α<1；f(x,t)為源項；a(x)，p(x)為有界函數(shù)。

2 剖分網(wǎng)格和構(gòu)造離散格式

首先對方程(1)進(jìn)行剖分網(wǎng)格和構(gòu)造離散格式：

(2)

(3)

(4)

式(4)減式(3)并除以h，得

(5)

對于空間一階偏導(dǎo)采用如下的離散方式：

。

(6)

第三步，針對時間方向的Caputo導(dǎo)數(shù)，進(jìn)行如下方式的離散：

根據(jù)定義1，得：

(7)

其中bk=(k+1)1-α-(k)1-α。

最終可得到方程(1)的差分格式：

(8)

初邊值條件為

(9)

將式(8)兩邊同乘h2，整理得

當(dāng)k=1時，

當(dāng)k>1時，

(10)

式(10)的矩陣形式為

，

(11)

其中

根據(jù)式(11)可以得到，對于式(10)的解存在唯一性，有下面的定理：

定理1差分方程(8)～(9)的解是存在且唯一的。

3 穩(wěn)定性分析和誤差估計

針對式(10)有如下的穩(wěn)定性分析和誤差估計：

3.1 穩(wěn)定性分析

定義2[12]對于任意初始誤差E0，若存在一個正常數(shù)C且與h和τ無關(guān)，使|Ek|∞≤C|E0|，那么該差分格式穩(wěn)定。

利用數(shù)學(xué)歸納法和最大模方法來證明算法是穩(wěn)定的。針對上述的離散格式，有定理2穩(wěn)定性結(jié)論。

定理2差分方程(8)～(9)為無條件穩(wěn)定。

當(dāng)k=1時，

≤G|E0|∞。

設(shè)|Em|∞≤C|E0|∞,m=2,3,…,k，又0

則

≤G|E0|∞。

得|Ek+1|∞≤|E0|∞,證畢。

3.2 誤差估計

式(8)關(guān)于u(x,t)，有如下截斷誤差：

≤C1(τ+h)。

當(dāng)k=1時，

當(dāng)k=2，3，…，N-1時，

因此，可以得到定理3差分格式的收斂性結(jié)論。

當(dāng)k=1時

從而有

則有

G|ek+1|

由kτ≤T，得

4 數(shù)值算例

通過運用一個一維時間分?jǐn)?shù)階變系數(shù)對流擴(kuò)散方程的數(shù)值算例來驗證其收斂階。

其中系數(shù)為a(x)=x2+1，p(x)=x+1，源項為

-20(t+1)2(1-3x+3x2-6x3)+10(t+1)2(2x-x2-3x3)。

上述方程的解析解為u(x,t)=10x2(1-x)(t+1)2。

本文取不同的空間步長及時間步長來證明算法的收斂精度，得到如表1和圖1～2所示結(jié)果。在表1中，τ=h取不同步長，當(dāng)α=0.3，0.5，0.8時，估計出τ和h的收斂精度；圖1～2得到了當(dāng)τ=h=1/50時，α=0.3的數(shù)值解圖像和函數(shù)圖像。

表1 當(dāng)τ=h時，α=0.3，0.5， 0.8分別對應(yīng)的誤差和收斂階

圖1 當(dāng)α=0.3，M=N=50時的數(shù)值解Fig.1 Numerical solutions of example for α=0.3 and M=N=50

圖2 當(dāng)α=0.3，t=1，M=N=50時的函數(shù)圖像Fig.2 Numerical solutions of example for α=0.3,t=1 and M=N=50

數(shù)值實驗結(jié)果表明：時間步長和空間步長當(dāng)取相同值而α取不同值時，收斂精度和理論分析基本一致；真解曲線與數(shù)值解曲線能夠較好地吻合。

5 結(jié)論

在本文中，首先利用有限差分方法將空間二階偏導(dǎo)數(shù)以及一階偏導(dǎo)數(shù)離散，再通過運用Caputo導(dǎo)數(shù)的定義，近似代替了時間分?jǐn)?shù)階偏導(dǎo)數(shù)，得到了該方程的有限差分格式。 理論分析表明，所提出的離散格式，其解是存在并且唯一的，收斂精度為ο(τ+h)，一維數(shù)值算例驗證出理論分析的準(zhǔn)確性。