時間分數(shù)階Allen-Cahn方程的重心插值配點法

黃蓉， 翁智峰

(華僑大學(xué) 數(shù)學(xué)科學(xué)學(xué)院， 福建 泉州 362021)

基于自由能的概念來描述界面且已經(jīng)被應(yīng)用于許多帶有界面問題的多相流系統(tǒng)中.最初的Allen-Cahn方程是以自由能的積分表示，是一種求解界面問題的數(shù)學(xué)模型，被廣泛應(yīng)用于平均曲率-流量[1]、晶體生長[2]、圖像處理[3]、人群擴散現(xiàn)象[4]和材料科學(xué)[5]等研究中.近年來，分數(shù)階算子的發(fā)展促進了非局部問題的迅速發(fā)展，將整數(shù)階Allen-Cahn方程擴展到分數(shù)階Allen-Cahn方程，已經(jīng)引起眾多學(xué)者的關(guān)注.湯濤等[7]首次證明了時間分數(shù)階梯度流的能量耗散性和數(shù)值穩(wěn)定性的分析；王宏等[8]采用快速Caputo算法結(jié)合傅里葉譜方法求解時間分數(shù)階Allen-Cahn方程；杜強等[9]提出分數(shù)階Allen-Cahn方程的凸分裂數(shù)值格式；劉歡等[10]基于SOE思想近似Caputo分數(shù)階導(dǎo)數(shù)，提出時間分數(shù)階Allen-Cahn的一種有限差分高效格式；張楠等[11]分別運用穩(wěn)定化方法和SAV策略結(jié)合向后微分公式(backward differentiation formula，BDF)構(gòu)造了時間分數(shù)階Allen-Cahn方程和Cahn-Hilliard方程的高階算法.任金城等[12]基于Laplace變換近似Caputo型分數(shù)階導(dǎo)數(shù)， 從而求解偏微分方程； 汪精英等[15]利用Laplace變換與算子分裂法、差分法求解分數(shù)階Allen-Cahn方程.采用Laplace變換近似分數(shù)階導(dǎo)數(shù)，可有效地減少由于分數(shù)階導(dǎo)數(shù)的歷史記憶性所引起的儲存量.

近年來，重心插值配點法作為一種新型的無網(wǎng)格計算方法，被廣泛應(yīng)用于求解各類微分方程.重心插值配點法成功地應(yīng)用于求解平面彈性問題[17]、Fredholm方程[18]、整數(shù)階Allen-Cahn方程[19-20].重心插值公式具有良好的數(shù)值穩(wěn)定性，它能以機器精度任意逼近光滑函數(shù)，將Lagrange插值公式改進，并定義重心權(quán)，可得到重心Lagrange公式.重心插值公式具有計算量小、操作方便、易于編程、高精度、穩(wěn)定性好等優(yōu)點.重心插值配點格式求解微分方程的理論分析相對少.最近，文獻[13-14]用重點插值配點法分別求解分數(shù)階電報方程和熱傳導(dǎo)方程，并給出理論分析.

基于此，本文將分數(shù)階Allen-Cahn方程通過Laplace變換近似為整數(shù)階方程，在時間與空間方向均采用重心插值配點法求解整數(shù)階方程，非線性項選取Newton迭代格式求解，并給出相容性誤差分析.

1 分數(shù)階Allen-Cahn方程的數(shù)值解

1.1 一維分數(shù)階Allen-Cahn方程

考慮Caputo類型的時間分數(shù)階Allen-Cahn方程，即

(1)

(2)

式(2)中，Γ(·)是Gamma函數(shù).

Allen-Cahn方程的基本能量泛函E(u)為

(3)

E(u)關(guān)于時間t求Caputo型分數(shù)階導(dǎo)數(shù)，即

(4)

顯然，E(u)關(guān)于時間t具有非增長性，該方程滿足能量耗散定律.

1.2 將分數(shù)階方程通過Laplace變換轉(zhuǎn)變?yōu)檎麛?shù)階方程

一維Allen-Cahn方程中的分數(shù)階導(dǎo)數(shù)可利用Laplace變換逼近，即

(5)

將qα在區(qū)間(0，1]作線性插值，并代入式(5)，即

(6)

對式(6)采用Laplace逆變換，并代入式(1)，推出整數(shù)階方程，即

(7)

(8)

1.3 重心插值配點法求解Allen-Cahn方程

(9)

(10)

(11)

1.3.2 重心有理插值 考慮重心有理插值配點法的基函數(shù).給定n+1個插值節(jié)點a=x0

(12)

將Zk(x)以Lagrange公式改寫，即有

(13)

(14)

由式(12)～(14)可推出重心有理插值公式，即

(15)

1.3.3 重心插值配點法求解Allen-Cahn方程 考慮整數(shù)階Allen-Cahn方程的重心插值配點法計算格式，其方程變形為

(16)

式(16)中：(x，t)∈[a，b]×[0，T]，G(u)=u(u2-1).

區(qū)域Ω=[a，b]×[0，T]，在空間[a，b]上取m+1個節(jié)點a=x0

考慮上述熱傳導(dǎo)方程，先固定t，考慮變量x，則u(x，t)在節(jié)點x1，x2，…，xm上的重心插值公式為

(17)

類似，固定x，考慮變量t，則u(x，t)在節(jié)點t1，t2，…，ts上的重心插值公式為

(18)

式(18)中，φj(t)為t方向上的插值基函數(shù).

由式(17),(18)可知，u(x，t)在節(jié)點{(xi，tj)，i=0，1，…，m，j=0，1，…，s}上的重心插值為

(19)

考慮u(x，t)對變量x，t求l+k階偏導(dǎo)數(shù)，即有

(20)

上述偏導(dǎo)數(shù)在節(jié)點(xr，tp)處的函數(shù)近似值為

(21)

因此，式(16)的重心插值配點法的計算格式為

上式可改寫成微分矩陣形式，即

(23)

1.4 相容性分析

設(shè)函數(shù)u(x，t)運用重心Lagrange插值法逼近的數(shù)值解為pm，s(x，t)，定義誤差函數(shù)e(x，t)=u(x，t)-pm，s(x，t).引用文獻[13]中的重要定理.

引理1若u(x，t)∈C(n+1)(Ξ)，Ξ=[a，b]×[0，T]，n=max{m，s}，定義的e(x，t)成立.即

(24)

(25)

根據(jù)引理1可推出定理1.

(26)

證明非線性項G(u)滿足利普希茨條件，則存在常數(shù)K>0，成立

|G(u(x，t))-G(u(xm，ts))|≤K|u(x，t)-u(xm，ts)|.

(27)

由式(1)，(8)可知

(28)

由引理1可知:R1=ut(x，t)-ut(xm，ts)=et(x，ts)+et(xm，ts)，可推出

(29)

同理，可得

(30)

阿玉奇是土爾扈特汗國開創(chuàng)者書庫岱青之孫、朋楚克之子，生于1642年，卒于1724年2月19日，享年81歲，執(zhí)政50余年。阿玉奇汗是土爾扈特歷史上著名的汗王，他執(zhí)政期間，由于他卓越的統(tǒng)治才能、高超的外交策略，以及輝煌的武功戰(zhàn)績，不僅造就他本人一生的豐功偉業(yè)，使他成為當時遐邇聞名的游牧汗國領(lǐng)袖，同時也將土爾扈特國推向前所未有的鼎盛階段。托忒文歷史文獻《卡爾梅克諸簡史》中曾對阿玉奇汗作過如下評述：他一生“幫助了許多國家和部落，沒有讓卡爾梅克人衰弱和受欺。讓他強盛者尊重他，與他相衡者懼怕他。名義上是俄羅斯臣民，但一切事情均由自己做主，所以，他是伏爾加河卡爾梅克汗中最有威望的一位”。

注1據(jù)文獻[14]的重要定理，類似定理1的推導(dǎo)，易知采用重心有理配點法求解整數(shù)階Allen-Cahn方程的誤差分析如下.

定理2設(shè)u(x，t)∈Cμ1+2[a，b]×Cμ2+2[0，T]，u(xm，ts)是u(x，t)運用重心有理配點法求解的數(shù)值解，成立Du(xm，ts)=0，非線性項G(u)滿足利普希茨條件，則

|u(x，t)-u(xm，ts)|≤C(hμ1-1+τμ2).

(31)

1.5 基于牛頓迭代的Allen-Cahn方程

考慮Allen-Cahn方程的非線性項G(u)，在u0處泰勒展開，可得

(32)

從而，式(23)的Newton迭代格式為

(33)

1.6 一維Allen-Cahn方程在重心插值配點法下的計算格式

由節(jié)1.3和節(jié)1.5可知，整數(shù)階Allen-Cahn方程在重心插值配點法下的計算格式為

(34)

2 數(shù)值算例

2.1 算例一

為了驗證數(shù)值格式的準確性，選取如下一維Allen-Cahn方程.即

(35)

式(35)中：精確解u(x，t)=(1+t3)sin(1.5πx)， (x，t)的取值范圍為[0，1]×[0，1].則有

取空間節(jié)點M=30，時間節(jié)點N=20，在ε=0.3，α=0.5下方程的精確解與誤差分布，分別如圖1,2所示.由圖1,2可知：兩種重心插值配點法的數(shù)值解圖像均與真解圖像逼近，且具有較高的精度.

(a) 精確解 (b) 誤差分布圖1 重心Lagrange插值配點法的數(shù)值解與誤差分布圖像(α=0.5)Fig.1 Numerical solution and error distribution images of barycentric lagrange collocation method (α=0.5)

(a) 精確解 (b) 誤差分布圖2 重心有理插值的數(shù)值解與誤差分布圖圖像(α=0.5)Fig.2 Numerical solution and error distribution images of barycentric rational collocation method (α=0.5)

在ε=0.3下，分別選取不同的剖分節(jié)點數(shù)和不同的α，利用兩種重心插值配點法計算Er(M，N)，結(jié)果如表1所示.

表1 兩種重心插值配點法求解方程的最大相對誤差(ε=0.3)Tab.1 Maximum relative error of solving equations by two barycentric interpolation collocation methods (ε=0.3)

由表1可知：當選取不同的α?xí)r，該配點格式采用較少的節(jié)點數(shù)，數(shù)值解可達到高精度；網(wǎng)格剖分越細，數(shù)值精度更高.與文獻[15]中的經(jīng)典差分法比較，當剖分節(jié)點數(shù)為M=50，N=2 500，最大相對誤差可達到10-3量級；文中提出的重心Lagrange配點法格式選取節(jié)點M=6，N=6，最大相對誤差就達到10-4量級，重心有理配點法格式選取節(jié)點M=4，N=4，最大相對誤差可達到10-3量級.數(shù)值算例表明，兩種重心插值配點格式剖分較少的節(jié)點就能達到很高的精度.

2.2 算例二

一維Allen-Cahn方程的離散能量函數(shù)如下

(37)

選定初值u(x，0)=ε·sin(1.5πx)，取Dirichlet的邊界條件，左邊界為u0=1， 右邊界為uM=-1，t∈[0，2];結(jié)合式(37)，令網(wǎng)格剖分為M=30，N=20，固定ε=0.3，α分別取0.2，0.5，0.9，則兩種重心插值配點法的能量遞減，如圖3，4所示.

(a) 重心Lagrange插值配點法 (b) 重心有理插值配點法圖3 不同插值配點法的能量遞減圖Fig.3 Energy decline images of different interpolation collocation method

(a) 重心Lagrange插值配點法 (b) 重心有理插值配點圖4 不同插值配點法的數(shù)值解圖像(α=0.2)Fig.4 Numerical solution images of different interpolation collocation method (α=0.2)

從圖3可知：當α取不同值時，時間分數(shù)階Allen-Cahn方程的能量耗散特性；當t增加時，能量泛函E(u)隨之遞減，最終趨于穩(wěn)定狀態(tài).系統(tǒng)的能量耗散速度隨著分數(shù)階階數(shù)α不同而變化，α越小，能量衰減越快.此外，E(u)的曲線受剖分細密程度影響，剖分變細，圖像更光滑.從圖4可知：兩種重心插值配點法的數(shù)值解圖像均與真解圖像逼近.

3 結(jié)束語

利用Laplace變換近似Caputo型分數(shù)階導(dǎo)數(shù)，將分數(shù)階Allen-Cahn方程轉(zhuǎn)化為整數(shù)階問題；在時-空方向均利用重心插值配點法離散求解整數(shù)階Allen-Cahn方程，并給出配點法格式的相容性誤差分析.數(shù)值算例結(jié)果表明，與文獻[15]中的有限差分格式比較，該重心配點格式剖分少量節(jié)點數(shù)即可達到格式的高精度，并滿足能量耗散規(guī)律.該方法可以廣泛推廣到求解其他時間分數(shù)階微分方程.