多項分數(shù)階非線性波動方程的數(shù)值方法及其快速實現(xiàn)

邵林馨，沈卓旸，馬葛沁舟，閔婕，黃健飛

(揚州大學 數(shù)學科學學院，江蘇 揚州 225002)

分數(shù)階偏微分方程在反?，F(xiàn)象和復(fù)雜系統(tǒng)的建模中發(fā)揮著越來越重要的作用[1-3],本文將要研究的帶有空間四階導(dǎo)數(shù)的多項時間分數(shù)階非線性波動方程可以較準確地建模氧氣通過毛細血管的動力學特征[3-4].大量文獻表明：研究人員很難獲得非線性分數(shù)階偏微分方程的解析解[5-6]，有時候即使能獲得，但由于這類解析解中往往含有無窮級數(shù)或無窮限積分，也仍然無法精確使用.因此，研究分數(shù)階偏微分方程的數(shù)值方法越來越受到研究人員的關(guān)注[7-12]，以期為分數(shù)階模型的實際應(yīng)用和模擬提供算法支撐.

到目前為止，對于空間上具有四階導(dǎo)數(shù)的時間分數(shù)階波動方程的數(shù)值計算方法已有較多的研究.例如：Jafari等[13]利用Adomian分解方法得到了定義在有界空間域中的四階分數(shù)階擴散-波動方程的半解析解，并通過實例討論了該方法的收斂性；Hu等[14]提出了一種四階分數(shù)階擴散波系統(tǒng)的有限差分方法，并證明了該方法的唯一性、穩(wěn)定性和收斂性；Li等[15]對四階分數(shù)階擴散波動方程提出了一種帶參數(shù)的5次樣條方法，并用能量法嚴格地建立了數(shù)值方法的可解性、穩(wěn)定性和收斂性；Gao等[16]給出了一類四階時間二項分數(shù)階波動方程的空間緊致差分格式，并證明了該差分格式的無條件穩(wěn)定性和收斂性；劉新龍等[17]對時間分數(shù)階四階擴散方程構(gòu)造了顯-隱和隱-顯2個差分格式，并證明了2個格式的穩(wěn)定性和超收斂性，在精度一致的情況下，顯-隱格式的計算效率要高于隱-顯格式；Elmahdi等[18]對帶空間四階導(dǎo)數(shù)的時間分數(shù)階非線性波動方程，基于等價變形的方法，構(gòu)造了一個線性化差分格式和一個緊差分格式，并進行了相關(guān)的數(shù)值分析.另外，需要指出的是：對于空間上具有四階導(dǎo)數(shù)的時間分數(shù)階偏微分方程的數(shù)值方法也有很多研究，可以參考文獻[19-22].

另一方面，關(guān)于多項時間分數(shù)階偏微分方程數(shù)值方法的研究也受到了廣泛關(guān)注.例如：Ren等[23]利用空間離散的緊致差分法和多項時間分數(shù)階Caputo導(dǎo)數(shù)的L1 逼近，提出了一維和二維多項時間分數(shù)階擴散波方程的有效數(shù)值格式，并加以證明無條件穩(wěn)定性和全局收斂性；Zheng等[24]提出了一種基于時空譜方法求解多項時間分數(shù)擴散方程的高階數(shù)值方法,該時空譜方法在時間和空間方向上都具有指數(shù)衰減；Abdel-Rehim等[25]通過數(shù)值計算得到了時間分數(shù)多項波動方程的近似解，包括時間分數(shù)波、強迫波和阻尼波動方程，并討論了馮-諾依曼穩(wěn)定性條件；Chen等[26]提出了一種具有可變系數(shù)的多項時間分數(shù)階擴散和擴散波動方程的統(tǒng)一格式，該格式在時間方向采用有限差分近似和在空間方向采用 Legendre 譜近似；王芬玲等[27]對具有 Caputo 導(dǎo)數(shù)的二維多項時間分數(shù)階擴散方程，采用線性三角元和改進 L1 格式，建立了一個全離散格式，并進行了高精度分析；Liu等[28]提出了一種二維多項時間分數(shù)混合擴散和擴散波動方程初邊值問題的ADI譜方法，該方法可以處理非光滑問題，并能模擬粘彈性非牛頓流體的擴散和輸運；Huang等[29]對一類二維時空多項分數(shù)階偏微分方程設(shè)計了一個快速的ADI格式，該格式在時間方向只需要相對較低的光滑性.

從現(xiàn)有的文獻資料來看，對帶有空間四階導(dǎo)數(shù)的多項時間分數(shù)階非線性波動方程數(shù)值方法的研究還非常少.因此，本文將對上述多項分數(shù)階非線性波動方程構(gòu)造一個線性化數(shù)值方法，并進行誤差分析.該線性化格式是基于等價的積分-偏微分方程來構(gòu)造的，可以避免求解非線性方程組，能計算初始奇異性問題，并且還能進行快速計算.

1 問題的提出及預(yù)備知識

本文將考慮為以下形式的多項時間分數(shù)階非線性波動方程[3-4]構(gòu)造一個線性化的數(shù)值格式：

(1)

u(x,0)=φ(x),ut(x,0)=φ(x),0≤x≤L,

u(0,t)=u(L,t)=ux(0,t)=ux(L,t)=0,0

給出一個四階導(dǎo)數(shù)的緊差分逼近，用來離散方程(1)中的空間方向四階導(dǎo)數(shù).

引理1[21]設(shè)u(x,·)∈C8[0,L],u(0,·)=u(L,·)=ux(0,·)=ux(L,·)=0,并令

(2)

為了便于計算，式(2)寫成如下的矩陣形式

容易證明矩陣S是1個正定矩陣.

引理2[12]設(shè)u(·,t)∈C1[0,T]∩C2(0,T],且1<αK<αK-1<…<α2<α1<2，則

利用引理2，方程(1)可等價地轉(zhuǎn)化為如下形式的偏微分-積分方程

(3)

引理3[29]設(shè)0<α<1,且u(·,t)∈C1[0,T],則有以下2個對Riemann-Liouville積分的逼近成立,

引理4[29]設(shè)u(·,t)∈C1[0,T]∩C2(0,T],0

在理論分析中將利用引理5，為此需要給出如下引理來說明引理3中定義的系數(shù)滿足引理5的條件和結(jié)論:

2 格式的構(gòu)造及其數(shù)值分析

2.1 格式的構(gòu)造

假設(shè)u(·,t)∈C1[0,T]∩C2(0,T],u(x,·)∈C8[0,L],u(0,·)=u(L,·)=ux(0,·)=ux(L,·)=0.下面考慮在點(xi,tn)處來離散方程(3).具體地，用引理1來離散方程(3)中的空間四階導(dǎo)數(shù)，用引理4來離散方程(3)中的時間一階導(dǎo)數(shù)和時間Caputo導(dǎo)數(shù)，分別用引理3中的第1個逼近和第2個逼近來離散方程(3)中等號左邊和右邊的Riemann-Liouville積分，可得

其中，

上式兩邊同乘以τ，可得

(4)

(5)

從式(5)可以發(fā)現(xiàn)：該式是一個線性式，可以避免求解非線性方程.另外，可以采用文獻[30-31]中的指數(shù)和技術(shù)來對式(5)進行快速實現(xiàn)，以提高計算的效率.值得一提的是：從解函數(shù)u(x,t)在時間方向的光滑性假設(shè)來看，即u(·,t)∈C1[0,T]∩C2(0,T]，格式(5)同樣適用于如下初始奇異性的情況：u(·,t)=1+O(tαK),即當t→0時，u(x,t)在時間方向的二階導(dǎo)數(shù)將趨于無窮.

2.2 數(shù)值分析

定理1式(5)存在唯一解.

證明：由于

因此，式(5)可表示為

(6)

對式(6)進行整理，得

下面，對式(5)進行誤差估計.首先，定義1個網(wǎng)格函數(shù)空間

對于任意νn,wn∈V,引入下列內(nèi)積和范數(shù)

在進行誤差估計時，需要借助于以下結(jié)論.

(7)

其中，Bn-m是式(6)中定義的系數(shù).

證明：用數(shù)學歸納法來證明.當J=2時，左邊=B1(V0-V1)=右邊，此時式(7)成立.現(xiàn)假設(shè)當J=Q時，式(7)成立，即

則當J=Q+1時，結(jié)合上面的假設(shè)，有

從而可得式(7)成立.

‖un-Un‖≤C(τ+h4),

其中，C是與離散參數(shù)無關(guān)的正的常數(shù)，且在不同的地方可能具有不同的值.

證明：式(4)和式(6)相減，得

{Bk}是正的單調(diào)遞減序列，且函數(shù)g滿足利普希茨條件，則對等式右邊使用柯西-施瓦茨不等式，可得

(1+2B0) ‖en‖2-‖en-1‖2≤(B0-B1)(‖en-1‖2+‖en‖2)+

其中，L是利普希茨常數(shù).因為

恒成立，則有

經(jīng)過簡單的代數(shù)運算，上式可變?yōu)?/p>

對上面不等式中的指標n從1～J進行求和，J為正整數(shù)，可得

根據(jù)式(7)和‖e0‖=0，可知

是負的，另外，注意到S是正定矩陣，則存在可逆矩陣H,使得S=HTH,從而上式可變?yōu)?/p>

由引理5可知，上式右端第1項是負的，從而上式變?yōu)?/p>

由引理3中{b(α1-1)n-m}的定義可知，

因而，上式可轉(zhuǎn)化為

根據(jù)Gronwall不等式，可推得

‖eQ‖≤C(τ+h4).

證畢.

3 數(shù)值實驗

該部分將給出2個數(shù)值算例來展現(xiàn)式(5)及其快速算法的計算效率，并驗證上述誤差估計的正確性.使用MATLAB 2012a來編程計算，電腦配置為：Intel(R) Core(TM) i7-8750 H CPU @ 2.20 GHz，內(nèi)存為16 GB.在具體數(shù)值實驗中，將使用max‖eJ‖作為計算的誤差.

算例1考慮以下問題

u(x,0)=sin2(πx),ut(x,0)=0,u(0,t)=u(1,t)=ux(0,t)=ux(1,t)=0,0

其中,T=1.令上述問題的精確解為u(x,t)=(tα1+1+1)sin2(πx).顯然，該解在時間方向的二階導(dǎo)數(shù)連續(xù)，且滿足定理2的光滑性假設(shè).此外，取非線性項為g(u)=u2，則

首先，計算式(5)及其快速算法在時間方向的誤差和數(shù)值收斂階、平均CPU時間結(jié)果見表1.具體地，對α1,α2,α3取2組不同的組合，取一個足夠小的空間步長h=0.001，這樣可以忽略空間方向的計算誤差.把時間步長τ從0.02開始逐步減半，計算時間方向的誤差.從表1可以發(fā)現(xiàn)，式(5)及其快速算法在時間方向的收斂階為1，符合定理2中的結(jié)論.此外，式(5)的快速算法比式(5)本身具有更高的計算效率，特別是當時間步長較小的時候，快速算法優(yōu)勢更大.

表1 在算例1中，式(5)及其快速計算在時間方向的誤差、數(shù)值收斂階和平均CPU時間Tab.1 Errors,temporal numerical convergence orders and average CPU times of scheme (5) for Example 1

其次，來計算空間方向的誤差和數(shù)值收斂階.由于空間方向具有高階精度，因此時間方向的步長需要非常小，故取τ=1/2 000 000.顯然，這個步長太小了，用快速算法來完成表2耗時約4.5 h.如果直接使用式(5)來計算，大概需要幾星期，再次說明了快速算法在時間步長較小時的計算優(yōu)勢.另一方面，從表2可知：構(gòu)造的算法在空間方向具有四階精度，支持了定理2中的理論分析結(jié)果.

表2 算例1中，式(5)的快速計算在空間方向的誤差和數(shù)值收斂階Tab.2 Errors and spatial numerical convergence orders of the fast implementation of Scheme (5) for Example 1

算例2考慮以下問題

u(x,0)=0,ut(x,0)=0,u(0,t)=u(1,t)=ux(0,t)=ux(1,t)=0,0

其中,T=1.令該問題的精確解為u(x,t)=tα3sin2(πx),g(u)=u3，則

8π4tα3cos(2πx)-t3α3sin6(πx).

不難發(fā)現(xiàn)，該精確解在時間方向具有初始奇異性，但仍滿足定理2中的光滑性假設(shè)，因此式(5)同樣適用該問題.從表3可知：式(5)及其快速算法在時間方向的收斂階為1，且快速算法的計算效率要明顯高于式(5)本身.由于空間收斂階的數(shù)值結(jié)果與表2較類似，這里不再列出.

表3 算例2中，格式(5)及其快速計算在時間方向的誤差、數(shù)值收斂階和平均CPU時間Tab.3 Errors,temporal numerical convergence orders and average CPU times of Scheme (5) and its fast implementation for Example 2

4 總結(jié)

本文主要對帶空間四階導(dǎo)數(shù)的多項時間分數(shù)階非線性波動方程構(gòu)造了一個數(shù)值格式，并證明了該格式的可解性和收斂性.該格式是一個線性化的格式，能適用于初始奇異性問題，并且可以進行快速計算.2個數(shù)值算例驗證了理論分析的正確性，并展現(xiàn)了快速算法的計算優(yōu)勢.