多項(xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程的一類顯-隱和隱-顯差分格式

秦 瀟，呂 蓬，楊曉忠

(華北電力大學(xué) 數(shù)理學(xué)院，北京 102206)

§1 引言

分?jǐn)?shù)階微積分已經(jīng)成功地運(yùn)用于反常擴(kuò)散，控制理論，粘彈性材料等領(lǐng)域[1-3].在這些問題中，分?jǐn)?shù)階對流擴(kuò)散模型能更準(zhǔn)確地模擬物理學(xué)中的傳輸過程[4-5].眾所周知，在復(fù)雜的污染物遷移過程中，它的短時間行為特征與長時間行為特征有很大的區(qū)別，擴(kuò)散運(yùn)動在許多情況下會跨越多個時間尺度.因此，與單項(xiàng)相比，引入多項(xiàng)時間分?jǐn)?shù)階導(dǎo)數(shù)能夠考慮時間尺度的變化對整個過程的影響，更準(zhǔn)確地建立數(shù)學(xué)物理模型.本文考慮多項(xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程[2，6]

多項(xiàng)時間分?jǐn)?shù)階導(dǎo)數(shù)定義為

多項(xiàng)分?jǐn)?shù)階微分方程的定性分析已經(jīng)完成了許多工作.Luchko[7]基于最大值原理，研究了多項(xiàng)時間分?jǐn)?shù)階擴(kuò)散方程的廣義解.Liu等[8]采用積分變換的方法研究了多項(xiàng)時間分?jǐn)?shù)階擴(kuò)散波方程的精確解.Li等[9]研究了多項(xiàng)時間分?jǐn)?shù)階擴(kuò)散方程初邊值問題的適定性和長時間漸近行為.

對于多項(xiàng)時間分?jǐn)?shù)階擴(kuò)散方程的數(shù)值解法已經(jīng)有了許多研究:Ren和Sun[10]對多項(xiàng)時間分?jǐn)?shù)階擴(kuò)散方程構(gòu)造了一種高階緊有限差分格式，證明了格式是無條件穩(wěn)定和收斂的.Jin等[11]采用有限差分方法和Galerkin有限元法構(gòu)造了多項(xiàng)時間分?jǐn)?shù)階擴(kuò)散方程的全離散格式，并給出了完整的誤差分析.基于時空頻譜法，Zheng等[12]對多項(xiàng)時間分?jǐn)?shù)階擴(kuò)散方程提出了一種高精度的數(shù)值格式.Zhao等[13]基于雙線性有限元法和有限差分法，構(gòu)造了多項(xiàng)時間分?jǐn)?shù)階擴(kuò)散方程的數(shù)值格式，討論了其收斂性.同時，采用插值后處理技術(shù)分析了該格式的超收斂性.Li等[14]采用混合有限元和有限差分方法得到了多項(xiàng)時間分?jǐn)?shù)階擴(kuò)散方程的全離散格式.楊曉忠等[15]構(gòu)造了一類差分格式求解雙項(xiàng)時間分?jǐn)?shù)階慢擴(kuò)散方程，證明了格式的無條件穩(wěn)定性，數(shù)值試驗(yàn)表明格式具有較高的計算效率.Kanth和Garg[16]基于時間上的有限差分方法和空間上的指數(shù)B-樣條方法構(gòu)造了多項(xiàng)時間分?jǐn)?shù)階擴(kuò)散方程的數(shù)值格式并證明了格式的無條件穩(wěn)定性.目前，關(guān)于多項(xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程數(shù)值求解的相關(guān)研究還很有限.Liu等[17]采用隱式差分方法對一類兩項(xiàng)時間分?jǐn)?shù)階移動/靜止對流擴(kuò)散方程進(jìn)行了數(shù)值求解，系統(tǒng)地分析了隱式數(shù)值方法的穩(wěn)定性和收斂性.Hussain和Haq[18]采用加權(quán)無網(wǎng)格譜徑向點(diǎn)插值方法對多項(xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程進(jìn)行了數(shù)值求解，并討論數(shù)值格式的穩(wěn)定性.

針對多項(xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程，基于經(jīng)典的顯式和隱式差分格式，本文構(gòu)造了一類顯式-隱式(E-I)差分格式和隱式-顯式(I-E)差分格式.這類格式體現(xiàn)了對稱離散的數(shù)值優(yōu)勢:將離散格式互補(bǔ)的兩個古典格式結(jié)合起來，在確保良好精度的基礎(chǔ)上提高計算效率.將顯式格式和隱式格式相結(jié)合成為一個雙步格式，一半步長用顯式格式計算，一半步長用隱式格式計算，這樣在每一個雙步只需求解一個含有三對角矩陣的方程組，計算復(fù)雜度得到下降，計算效率會明顯改善.這也是本文設(shè)計顯隱交替格式的動機(jī).

§2 多項(xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程的E-I格式

為了構(gòu)造方程(1)的E-I格式，首先給出方程(1)的經(jīng)典顯式格式和經(jīng)典隱式格式.方程(1)的顯式格式

方程(1)的隱式格式

多項(xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程的E-I格式構(gòu)造如下:求奇數(shù)層時，采用顯式格式計算;求偶數(shù)層時，采用隱式格式計算

§3 E-I格式的數(shù)值分析

3.1 格式解的存在唯一性

當(dāng)4-h2≤0，矩陣A的特征值可以寫為p+iq，其中0;

當(dāng)4-h2＞0，有λi≥p0+

因此，在任何情況下，A的特征值都不為零，系數(shù)矩陣A是可逆的.綜上分析，有如下定理1.

定理1多項(xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程(1)的E-I差分格式(9)的解是存在唯一的.

3.2 格式的穩(wěn)定性

εn(x)的傅里葉展開為

引理1存在一正常數(shù)C，使得|vk|≤C|v0|，k1，2，···，N.

證當(dāng)時間層k1時，用顯式格式計算

其中C為一正常數(shù).

當(dāng)時間層k2時，第一層用顯式格式計算，第二層用隱式格式計算，即

其中C為一正常數(shù).

假設(shè)當(dāng)時間層k ≤2n+1時，有|vk|≤C|v0|.則當(dāng)k2n+2時，由(10)可得誤差方程為

因此，存在一正常數(shù)C，使得|vk|≤C|v0|，k1，2，···，N.

定理2多項(xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程(1)的E-I差分格式(9)是無條件穩(wěn)定的.

證由(13)和引理1，可以得出||εn||2≤C||ε0||2，n1，2，···，N.即多項(xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程(1)的E-I差分格式(9)是無條件穩(wěn)定的.

3.3 格式的收斂性

T1(τ，h)中的τuxxt，τuxt和T2(τ，h)中的τuxxt，τuxt形式相同但符號相反.因此，交替使用顯格式和隱格式時這兩項(xiàng)的誤差會被抵消掉.故E-I格式的精度在空間上為2階，時間上為2-α階(αmax{α0，α1，···，αm}).

§4 多項(xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程的I-E格式

類似多項(xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程E-I差分格式的構(gòu)造方法，構(gòu)造隱-顯(I-E)差分格式:求奇數(shù)層時，采用隱格式計算;求偶數(shù)層時，采用顯格式計算.得到多項(xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程的I-E差分格式

定理4假設(shè)u(x，t)滿足光滑條件(16)，多項(xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程(1)的I-E差分格式(21)是唯一可解的，無條件穩(wěn)定且收斂的，滿足

C是一正常數(shù).

證其證明過程與定理1-3類似.

E-I格式和I-E格式均為雙步格式，在時間層上依次采用顯式格式和隱式格式，兩者的區(qū)別僅在于使用顯式格式和隱式格式的順序，因此E-I和I-E格式的計算量是相當(dāng)?shù)?

注本文是在(1)的解滿足光滑條件(16)時，得到格式(9)和格式(21)的收斂階為O(τ2-α+h2).一般來說，解u(x，t)及其導(dǎo)數(shù)在t0附近具有奇異性.因此，當(dāng)解u(x，t)非光滑時，格式(9)和格式(21)的時間收斂階明顯低于2-α階(αmax{α0，α1，···，αm}).實(shí)際上，多項(xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程(1)的E-I格式和I-E格式的時間收斂階在這種情況下為α階(參見文獻(xiàn)[22-23]).也在數(shù)值試驗(yàn)中證實(shí)了這個結(jié)論.

§5 數(shù)值試驗(yàn)

數(shù)值試驗(yàn)基于Inter Core i5-8265 CPU，在Matlab R2017b環(huán)境下運(yùn)行.

例1多項(xiàng)(雙項(xiàng))時間分?jǐn)?shù)階對流擴(kuò)散方程[6]

方程(22)的精確解為u(x，t)t2x2(1-x)2.

取α10.15，α20.95，為保證顯格式的穩(wěn)定性，時間步長τ0.0002，空間步長h0.025，在t0.1處，將本文的精確解，顯格式解，隱格式解與E-I，I-E格式解進(jìn)行比較，計算結(jié)果如表1所示.可以看到E-I格式解和I-E格式解比隱格式解更好地逼近了精確解.

表1 數(shù)值解與精確解的比較

顯隱交替格式展現(xiàn)出對稱離散的數(shù)值優(yōu)勢:由表1看出，古典隱式格式給出精確解的上方近似值，而古典顯式格式給出精確解的下方近似值.顯式格式和隱式格式的數(shù)值誤差具有相反的符號，正負(fù)誤差的相互抵消提高了數(shù)值解的精度[24].

因?yàn)槎囗?xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程的E-I格式和I-E格式十分相近，下面只考慮多項(xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程的E-I格式.

取M1000，使得相較于時間離散誤差空間離散誤差可以忽略不計.時間收斂階的計算結(jié)果如表2所示.這些結(jié)果表明E-I格式在時間上具有2-α階精度(αmax{α1，α2})，在數(shù)值上支持了定理3的結(jié)論.表3可以看出E-I格式具有空間2階精度，與隱格式精度一致.同時表2，表3可以看出相比較于隱格式，E-I格式可以節(jié)省計算時間.

表2 E-I格式和隱格式的數(shù)值誤差，時間收斂階和計算時間(h)

表2 E-I格式和隱格式的數(shù)值誤差，時間收斂階和計算時間(h)

表3 E-I格式和隱格式的數(shù)值誤差，空間收斂階和計算時間(τ h2)

表3 E-I格式和隱格式的數(shù)值誤差，空間收斂階和計算時間(τ h2)

為了更好的比較隱格式，E-I格式和I-E格式的計算效率，令α10.2，α20.8，將時間層固定為N1000，選取空間網(wǎng)格數(shù)M500，1000，2000，3000，4000，5000， 計算結(jié)果如圖1.可以看到對于多項(xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程，當(dāng)M ≤500時，E-I格式和隱格式的計算時間相近，這和表3中數(shù)值結(jié)果相符.隨著空間網(wǎng)格數(shù)的增加，三種差分格式的計算時間都在增加.但顯然，隱格式的計算時間比E-I格式和I-E格式增加速度更快.當(dāng)M5000時，E-I格式和I-E格式的計算時間比隱格式節(jié)省約30%，表明E-I格式和I-E格式的計算效率高于隱格式.

圖1 三種格式的計算時間比較

為了驗(yàn)證E-I格式和I-E格式求解具有初始奇性的多項(xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程的有效性.考慮如下方程.

例2具有初始奇性的多項(xiàng)(雙項(xiàng))時間分?jǐn)?shù)階對流擴(kuò)散方程

取MN300，α10.4，α20.5，在t0.1處，將本文的精確解，隱格式解與E-I，I-E格式解進(jìn)行比較，計算結(jié)果如表4.可以看到三種數(shù)值解與精確解的誤差是很小的.表明E-I和I-E方法同樣適用于求解多項(xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程的解具有初始奇性的情況.

表4 數(shù)值解與精確解的比較

因?yàn)槎囗?xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程的E-I格式和I-E格式十分相近，下面只考慮多項(xiàng)時間分?jǐn)?shù)階對流擴(kuò)散方程的E-I格式.將討論E-I格式和隱格式的時間收斂階和空間收斂階.

表5表明E-I格式的時間精度達(dá)到了α2階，比隱格式的時間精度要高.可以看到，在解非光滑的前提下，E-I格式的時間收斂階要低于定理3證明的2-α階(αmax{α1，α2}).E-I差分格式和隱格式的數(shù)值誤差和空間收斂階如表6所示.從表6可以看出E-I格式和隱格式具有2階空間精度.

表5 E-I格式和隱格式的數(shù)值誤差和時間收斂階(h)

表5 E-I格式和隱格式的數(shù)值誤差和時間收斂階(h)

表6 E-I格式和隱格式的數(shù)值誤差和空間收斂階(τh4)

表6 E-I格式和隱格式的數(shù)值誤差和空間收斂階(τh4)

基于對表1至表6中數(shù)據(jù)的分析和對圖1的描述，可以看到數(shù)值理論的正確性和數(shù)值算法的有效性.