童艷蕾 姚曉珍 林夢渟 李怡雯 翁智峰

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

1 引言

1.1 研究背景

Allen-Cahn 方程作為描述相場模型最基本的方程之一,可用于模擬二元合金相位分離的過程.此外,Allen-Cahn 方程也常用于描述一定溫度和氣壓下晶體的生長[1]、生物種群之間的競爭和排斥現(xiàn)象[2]、河床的遷移現(xiàn)象[3],平均曲率-流量[4]等問題.

Allen-Cahn 方程是具有如下形式的方程:

其中,ε為界面寬度,?為拉普拉斯算子,f(u) 為給定的非線性函數(shù),滿足f(u)=F′(u),其中F(u)=是一種雙阱勢函數(shù).

Allen-Cahn 方程的一個重要特征是:它可以被看成能量泛函

的梯度流形式.因此,其能量泛函滿足遞減關(guān)系

由于Allen-Cahn 方程的解析解通常很難得到,因此研究其數(shù)值解變得十分重要.學(xué)者們對其提出了許多數(shù)值解法.比如,Chen 等[5]基于有限差分法求解Allen-Cahn 方程并分析了其數(shù)值格式的收斂性? Zhai 等[6]提出了求解三維Allen-Cahn 方程的高階緊致ADI 方法? Li 等[7]給出了Allen-Cahn 方程的二階無條件能量穩(wěn)定的有限元格式?Feng 和Li[8]構(gòu)造了一種對稱內(nèi)罰間斷有限元格式?Weng 和Tang[9]利用算子分裂格式結(jié)合譜方法求解Allen-Cahn 方程?Wang 等[10]提出了一種線性能量穩(wěn)定和保極大值原理的半隱格式求解Allen-Cahn 方程?Li 等[11]構(gòu)造了Allen-Cahn 方程的降階修正有限差分格式?Li 和Wang[12]給出了一種修正有限體積元方法求解Allen-Cahn 方程?Hou 和Leng[13]構(gòu)造了一個穩(wěn)定的Crank-Nicolson/Adams-Bashforth 格式?Feng 等[14]給出了一個無條件能量穩(wěn)定的線性二階有限差分格式.

上述方法都是基于網(wǎng)格剖分的思想求解Allen-Cahn 方程.文獻(xiàn)[15]提出一種新型無網(wǎng)格的重心插值配點(diǎn)法,該方法具有計(jì)算格式簡單、精度高、程序?qū)嵤┓奖?、?jié)點(diǎn)適應(yīng)性好等特點(diǎn),被廣泛應(yīng)用于求解各類微分方程,比如Cahn-Hilliard 方程[16]、Helmholtz 方程[17]、Black-Scholes 方程[18]、Allen-Cahn 方程[19]等.最近,Yi 和Yao[20]利用重心插值配點(diǎn)法求解時間分?jǐn)?shù)階電報(bào)方程并給出了對應(yīng)的誤差估計(jì)?Li 和Cheng[21]分析了熱傳導(dǎo)方程在重心有理插值配點(diǎn)格式下的誤差估計(jì).

針對非線性微分方程,Ascher 等[22]提出了一種顯隱方法.該方法對方程的非線性項(xiàng)采用顯式處理,對線性項(xiàng)做隱式處理,避免了非線性迭代,可以提高計(jì)算效率.沈和馬[23]提出了一種半隱差分法求解Zakharov 方程?潘[24]給出了Allen-Cahn 方程的一種高階顯隱格式?湯和喬[25]給出了相場方程的一種高效數(shù)值算法.

基于前人的工作,本文利用重心插值配點(diǎn)法結(jié)合半隱格式求解Allen-Cahn 方程,對方程的空間一維半離散格式和二維全離散格式進(jìn)行相容性分析,并通過數(shù)值算例驗(yàn)證算法的高精度性和數(shù)值格式的能量遞減物理性質(zhì).

1.2 重心Lagrange 插值與Chebyshev 配點(diǎn)法

重心Lagrange 插值設(shè)有m+1 個插值節(jié)點(diǎn)xj和相對應(yīng)的一組實(shí)數(shù)fj,j=0,1,···,m.在次數(shù)不超過m的多項(xiàng)式空間中,存在唯一一個多項(xiàng)式p(x),滿足p(xj)=fj(j=0,1,···,m).多項(xiàng)式p(x)即為如下的Larange 插值多項(xiàng)式:

其中,Lagrange 插值基函數(shù)κj(x)滿足性質(zhì):

由于Lagrange 插值公式的計(jì)算量達(dá)到O(n2),且增加新節(jié)點(diǎn)時需要重新計(jì)算,為了減少重復(fù)計(jì)算,我們令l(x)=(x-x0)(x-x1)···(x-xm),定義重心權(quán)

將插值基函數(shù)表示為:

再代入式(1.2),得

由式(1.3),(1.4)可得

再代入式(1.5)可得重心Lagrange 插值公式:

其中

Chebyshev 配點(diǎn)法多項(xiàng)式插值在節(jié)點(diǎn)數(shù)m較大時,若采用等距節(jié)點(diǎn),其插值是病態(tài)的.因此本文將以第二類Chebyshev 節(jié)點(diǎn)作為空間方向的插值節(jié)點(diǎn).

在區(qū)間[-1,1]上的第二類Chebyshev 節(jié)點(diǎn)為:

相應(yīng)的重心Lagrange 插值的權(quán)為:

任意區(qū)間[a,b]上的Chebyshev 節(jié)點(diǎn)可通過坐標(biāo)變換y=x(b-a)/2+(a+b)/2 得到.

重心插值的微分矩陣函數(shù)p(x)在節(jié)點(diǎn)xi處的k階導(dǎo)數(shù)為:

由Lj的定義可得到一階、二階重心插值微分矩陣,其元素如下:

由數(shù)學(xué)歸納法可以得到k階重心插值微分矩陣,其元素如下:

2 一維Allen-Cahn 方程的空間半離散格式及其相容性

2.1 一維Allen-Cahn 方程的空間半離散格式

考慮一維Allen-Cahn 方程的初邊值問題:

式中u(x,t)為未知函數(shù),ε為界面寬度系數(shù),g(x),θ(t),ω(t)為已知函數(shù),且非線性項(xiàng)f(u)=u-u3.

將區(qū)間[a,b] 離散為m+1 個第二類chebyshev 節(jié)點(diǎn):a=x0

將其代入Allen-Cahn 方程,可得m+1 個方程

其中U=[u0(t),u1(t),···,um(t)]T.

2.2 一維Allen-Cahn 方程空間半離散格式的相容性

設(shè)p(x)為u(x)的重心Lagrange 插值近似函數(shù).定義誤差函數(shù)

下面,給出重心插值配點(diǎn)法的逼近性質(zhì).

引理2.1([20])如果u(x)∈Cm+1[a,b],則有

為了討論一維Allen-Cahn 方程的空間半離散格式(2.4)的相容性,定義算子

易知,對于Allen-Cahn 方程的未知函數(shù)的重心插值近似函數(shù)(2.2),有Γu(xm,t)=0 和=0.

定理2.1 如果u(x,t)∈Cm+1[a,b]×C0[0,T],由重心Lagrange 插值配點(diǎn)法得到的方程(2.4)的數(shù)值解為u(xm,t),且方程的非線性項(xiàng)滿足Lipschitz 條件,則有:

證明

其中,

對R1,我們有

再由引理2.1,可得

同樣地,對R2與R3,我們分別有

結(jié)合(2.6),(2.7),(2.8),可得

其中C=.證畢.

3 一維Allen-Cahn 方程的全離散格式

通過對空間方向的離散,一維Allen-Cahn 方程轉(zhuǎn)換為關(guān)于時間變量t的一階常微分方程.對于空間重心插值半離散系統(tǒng)(2.4),我們將在時間方向上分別采用向后差分法和Crank-Nicolson 方法設(shè)計(jì)方程的全離散格式.

3.1 一階半隱全離散格式

將時間域[0,T]等距剖分為n份,得到n+1 個時間節(jié)點(diǎn):tk=kτ,τ=T/n,k=0,1,···,n.令Uk=U(tk).對方程的非線性項(xiàng)采用顯格式處理,時間方向采用向后差分離散,可得

整理可得如下一階半隱全離散格式:

對方程(2.1)的Dirichlet 邊值條件u(a,tk)=θ(tk),u(b,tk)=ω(tk),0≤k≤n,按下述方法離散化.

3.2 二階半隱全離散格式

與一階格式相同,等距剖分時間域,然后對方程(2.4)的非線性項(xiàng)采用顯格式處理,時間方向采用Crank-Nicolson 格式離散,可得:

將式(3.3),(3.4)代入(3.2),整理后即可得一維Allen-Cahn 方程的二階半隱全離散格式:

4 二維Allen-Cahn 方程的數(shù)值格式

考慮二維Allen-Cahn 方程的初邊值問題:

本節(jié)對方程的非線性項(xiàng)f(u)=u-u3仍用顯格式處理,對方程的空間方向基于重心Lagrange 插值配點(diǎn)法進(jìn)行離散,時間方向采用差分格式進(jìn)行離散,給出Allen-Cahn 方程的數(shù)值格式.

4.1 二維重心Lagrange 插值半離散系統(tǒng)

與一維Allen-Cahn 處理方式相同,將方程的矩形空間域?=[a,b]×[c,d]離散為(m+1)×(l+1)個第二類Chebyshev 節(jié)點(diǎn):a=x0

其中,U=,ui(t)=[ui0(t),ui1(t),···,uil(t)]T,i=0,1,···,m,f(U)=U-U3為方程的非線性項(xiàng).

4.2 二維全離散格式

一階半隱全離散格式對(4.2)的非線性項(xiàng)采用顯格式處理,時間方向采用向后差分格式離散,可得到方程(4.1)的一階半隱全離散格式:

二階半隱全離散格式對(4.2)的非線性項(xiàng)采用顯格式處理,時間方向采用CN 格式離散,可得到方程(4.1)的二階半隱全離散格式:

4.3 二維Allen-Cahn 方程全離散格式的相容性

設(shè)p(x,y)為u(x,y)的重心Lagrange 插值近似函數(shù).定義誤差函數(shù):

引理4.1([20]) 如果u∈,則有如下估計(jì):

關(guān)于Allen-Cahn 方程二維全離散格式的相容性,我們有如下的定理.

定理4.1 設(shè)u(x,y,t)∈×C4(0,T],?=[a,b]×[c,d].采用CN 格式、重心Lagrange插值配點(diǎn)法得到的非線性格式的數(shù)值解記為uh(xm,yl,tn).如果非線性項(xiàng)f(u)滿足Lipschitz 條件,則有:

證明設(shè)對方程(1.1)在時間方向采用CN 格式離散得到的數(shù)值解為.則

然后對(4.5)在空間方向利用重心Lagrange 插值配點(diǎn)法進(jìn)行離散,可得

其中,ηm,l是空間截?cái)嗾`差.

由式(4.5)和(4.7),有

由非線性項(xiàng)滿足Lipschitz 條件及引理4.1,類似定理2.1 的推導(dǎo),可得:

再利用(4.6)即可得到定理的結(jié)論.證畢.

5 數(shù)值實(shí)驗(yàn)

下面通過算例驗(yàn)證上文給出的兩種格式的有效性.絕對誤差Erra和相對誤差Errr定義如下:

其中,uc和ue分別是方程(1.1)或(4.4)的數(shù)值解和精確解,∥·∥∞為矩陣的無窮范數(shù).

5.1 算例1

考慮如下二維Allen-Cahn 方程:?={(x,y,t)|0≤x≤1,0≤y≤1,0≤t≤1},

其中h(x,y,t)=-sin(πx)sin(πy)sin(t)+sin3(πx)sin3(πy)cos3(t)-(1-2επ2)sin(πx)sin(πy)cos(t).該方程的精確解為u=sin(πx)sin(πy)cos(t).

為了驗(yàn)證方程(5.1)的BDF/重心插值格式(4.3)和CN/重心插值格式(4.4)的時間收斂階,取參數(shù)ε=0.1,M=m+1=L=l+1=10,及不同的時間步長τ,計(jì)算結(jié)果見表1.表中Rate=log 2(Errr2/Errr1)/log 2(τ2/τ1).

表1 方程(5.1)兩種離散格式的時間收斂階(ε=0.1,M=L=10)

由表1可知,BDF/重心插值格式的時間收斂階為1,CN/重心插值格式的時間收斂階為2,且隨著時間步長τ的減小,相對誤差Errr不斷減小.

表2顯示了在固定參數(shù)ε=0.3,τ=10-4,改變空間節(jié)點(diǎn)數(shù)(M,L)時,時間一階半隱配點(diǎn)格式的絕對誤差和相對誤差結(jié)果.由表2 可知,隨著空間節(jié)點(diǎn)的增加,數(shù)值解的相對誤差越來越小.

表2 方程(5.1)兩種離散格式的誤差(ε=0.3,τ=10-4)

下面,我們將本文給出的CN/重心插值格式與空間方向分別采用二階中心差分離散和緊差分離散的兩種經(jīng)典格式的精度進(jìn)行比較,數(shù)值結(jié)果如表3所示.

表3 不同空間離散格式的精度對比(ε=0.3,τ=10-4)

從表3可知,與空間方向采用中心差分的格式比較,重心插值使用6×6 個節(jié)點(diǎn)達(dá)到的精度比二階中心差分使用32×32 個節(jié)點(diǎn)達(dá)到的精度還高?重心插值使用10×10 個節(jié)點(diǎn)達(dá)到的精度比緊差分使用32×32 個節(jié)點(diǎn)達(dá)到的精度要高得多.

由三種空間離散格式的收斂精度圖(圖1)可以看出,重心插值配點(diǎn)法的誤差從10-2下降到10-8,具有指數(shù)收斂的性質(zhì).

圖1 不同空間離散格式的收斂精度(ε=0.3,τ=10-4,M=L)

5.2 算例2

考慮如下二維Allen-Cahn 方程:?={(x,y,t)|-1≤x≤1,-1≤y≤1,0≤t≤0.5},

其二維能量函數(shù)的離散形式為:

其中hi=xi-xi-1,hj=yj-yj-1.

首先,對BDF/重心插值格式(4.3),取τ=10-4,M=L=10,在不同的參數(shù)值ε下計(jì)算該方程相應(yīng)數(shù)值解的能量函數(shù),結(jié)果如圖2(a) 所示.其次,對CN/重心插值格式(4.4),取τ=10-3,M=L=10,在不同的參數(shù)值ε下計(jì)算該方程相應(yīng)數(shù)值解的能量函數(shù),結(jié)果如圖2(b)所示.

圖2 方程(5.2)兩種格式的能量函數(shù)圖(M=L=10)

由圖2可知,方程(5.2)的兩種數(shù)值格式的能量函數(shù)都隨時間遞減,且ε越大,能量越早達(dá)到穩(wěn)態(tài).

圖3是取參數(shù)ε=0.3,τ=10-3,M=L=10 時,方程(5.2)的兩種不同格式的能量函數(shù)圖,可以看出,CN/重心插值格式的能量函數(shù)線先達(dá)到穩(wěn)定狀態(tài),說明其精度更高.

圖3 方程(5.2)的兩種不同格式能量函數(shù)圖(ε=0.3,τ=10-3,M=L=10)

5.3 算例3

考慮如下二維Allen-Cahn 方程:?={(x,y,t)|-1≤x≤1,-1≤y≤1,0≤t≤50},

取參數(shù)ε=0.01,τ=10-2,M=L=30.方程(5.3)的二階半隱配點(diǎn)格式的數(shù)值解快照圖像如圖4所示.由圖可以觀察到隨機(jī)生成的無序界面逐漸變淡直至到達(dá)穩(wěn)定狀態(tài).

圖4 方程(5.3)的二階半隱格式的數(shù)值解在不同時刻t 的快照圖像(ε=0.01,τ=10-2)

5.4 算例4

考慮如下二維Allen-Cahn 方程:?={(x,y,t)|-2≤x≤2,-2≤y≤4,0≤t≤4},

取參數(shù)ε=0.1,τ=10-3,M=L=40.二階半隱配點(diǎn)格式數(shù)值解的快照圖像如圖5所示.由圖可以觀察到連接的界面逐漸分成兩個曲面,然后慢慢變小變圓直至最后消失.

圖5 方程(5.4)二階半隱格式的數(shù)值解在不同時刻t 的快照圖像(ε=0.1,τ=10-3)

5.5 算例5

在?=[0,1]2×[0,T]上考慮帶初值條件

的Allen-Cahn 方程.

設(shè)置參數(shù)為ε=0.03,τ=10-3,T=2.該方程的二階格式的數(shù)值解在不同時刻的快照圖像如圖6所示.由圖可以看出,圓環(huán)隨著時間的推移逐漸變淡直至消失.

圖6 方程的二階階半隱格式的數(shù)值解在不同時刻t 的快照圖像(ε=0.03,τ=10-3)

6 總結(jié)

本文分別將向后歐拉格式和Crank-Nicolson 格式與重心Lagrange 插值Chebyshev 配點(diǎn)法方法相結(jié)合,建立求解Allen-Cahn 方程的數(shù)值格式,并對方程的空間一維半離散格式和二維全離散格式的相容性進(jìn)行了分析.針對方程的非線性項(xiàng)采用顯格式處理,避免了非線性迭代而帶來的大計(jì)算量,提高了計(jì)算效率.通過與經(jīng)典的兩種差分格式進(jìn)行比較,數(shù)值算例表明我們的數(shù)值格式可以用較少的節(jié)點(diǎn)達(dá)到較高的精度并滿足能量遞減規(guī)律.本文的重心插值半隱配點(diǎn)格式可推廣到其他非線性微分方程,為解決同類問題提供參考.