秦 瀟，呂 蓬，楊曉忠

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

§1 引言

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

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

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

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

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

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

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

方程(1)的隱式格式

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

§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)時(shí)間分?jǐn)?shù)階對(duì)流擴(kuò)散方程(1)的E-I差分格式(9)的解是存在唯一的.

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

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

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

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

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

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

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

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

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

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

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

3.3 格式的收斂性

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

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

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

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

C是一正常數(shù).

證其證明過(guò)程與定理1-3類似.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖1 三種格式的計(jì)算時(shí)間比較

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

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

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

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

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

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

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

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

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

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

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