SAV方法求解Allen-Cahn方程的數(shù)值比較

陳 航,吳 哲,翁智峰

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

0 引言

S.M. Allen等[1]提出反相邊界運動的微觀理論,對應(yīng)的Allen-Cahn方程在描述相場變化上有著非常廣泛的應(yīng)用,如研究生物種群之間的競爭和排斥[2]、在已知溫度氣壓下晶體的生長[3]、河床的遷移[4]、在材料學(xué)中界面擴散動力學(xué)[5]等.

本文考慮能量泛函

(1)

其中F(u)=(u2-1)2/(4ε2)為雙位勢阱函數(shù),則可以寫出式(1)對應(yīng)的梯度流方程(即Allen-Cahn方程)

ut=-Mμ,μ=-Δu+F′(u),(x,t)∈Ω×(0,T],u0=u(x,0),x∈Ω,u|?Ω=0,

(2)

其中M為遷移率系數(shù).

近年來,許多學(xué)者對Allen-Cahn方程的數(shù)值解進行了深入廣泛研究，主要包括有限差分法[6]、有限元法[7]、譜方法[8]等.如D.J. Eyre[9]提出了唯一解且無條件能量穩(wěn)定的凸分裂方法；Chen Longqing等[10]構(gòu)造了滿足能量穩(wěn)定的1階精度和2階精度的半隱格式;Choi Jeongwhan等[11]提出了一種求解Allen-Cahn無條件梯度穩(wěn)定的數(shù)值格式,可以滿足能量耗散;Yang Xiaofeng等[12]提出了不變能量2次化(IEQ)的方法來求解分子束外延生長模型，只要滿足非線性項有下界，就可以保證修正的能量穩(wěn)定;Shen Jie等[13]提出了標(biāo)量輔助變量(Scalar Auxiliary Variable,SAV)方法求解梯度流問題，該方法只要非線性項F(u)的積分項有下界條件，就可以保證修正的能量無條件穩(wěn)定.SAV方法保留了不變能量2次化的優(yōu)點，也不需要每一步求解變系數(shù)方程，只需每一步求解常系數(shù)的線性方程組，并且在差分格式下對應(yīng)的常系數(shù)方程組有快速方法進行求解，數(shù)值格式非常高效.由于SAV數(shù)值格式的諸多優(yōu)點，所以它被廣泛應(yīng)用于求解各類微分方程,如Allen-Cahn方程和Cahn-Hilliard方程[14]、研究細菌趨化性的Keller-Segel方程[15]、在流體力學(xué)中的Navier-Stokes方程[16]等.

眾所周知，在Lagrange插值公式描述經(jīng)過插值點的近似函數(shù)時，若節(jié)點過多則會產(chǎn)生Runge現(xiàn)象，而N.J. Higham提出的重心型Lagrange插值[17]解決了該問題，并且若將節(jié)點改成Chebyshev節(jié)點則可以使插值達到很高的精度.重心插值配點法[18]作為一種新型的無網(wǎng)格計算方法，能有效地避免差分格式帶來的累積誤差，使用Chebyshev節(jié)點有效地克服了Runge現(xiàn)象，且其具備計算格式簡單、精度高、程序?qū)嵤┓奖恪⒐?jié)點適應(yīng)性好等特點.如今，已經(jīng)應(yīng)用在求解各類微分方程(如分數(shù)階Fredholm積分方程[19]、平面彈性問題[20]、Allen-Cahn方程[21]等)中.

本文主要利用SAV格式結(jié)合無網(wǎng)格型的重心Lagrange插值配點法來求解Allen-Cahn方程.首先,基于SAV數(shù)值格式的Allen-Cahn方程，在空間方向上采用了高精度的重心Lagrange插值格式，并與常見的2階中心差分格式做比較，觀察2者誤差和收斂階;然后,再利用差分矩陣的性質(zhì)，運用離散正弦變換(Discrete Sine Transform,DST)求解2階中心差分導(dǎo)出的線性方程組和快速傅里葉變換(Fast Fourier Transform,FFT)計算Toeplitz矩陣與向量的乘積;最后,驗證重心Lagrange插值配點數(shù)值格式指數(shù)收斂且耗時少.

1 Allen-Cahn方程的無條件能量穩(wěn)定格式

1.1 Allen-Cahn方程的SAV形式

ut=-Mμ,

(3)

(4)

(5)

進而,只需對式(3)內(nèi)積μ,式(4)內(nèi)積ut,式(5)乘以2r,整理即可得到

(6)

所以式(6)對應(yīng)修正的能量泛函E(u)=(u,-Δu)/2+r2是無條件能量穩(wěn)定的.

1.2 時間方向的2種離散格式

對方程的時間偏導(dǎo)數(shù)部分,采用2階向后差分(BDF2)格式和Crank-Nicolson(CN)格式進行離散.記在[0,T]上的離散點為0=t0

1.2.1 2階向后差分格式 利用BDF2格式,取正整數(shù)N,記τ=T/N,tn=nτ(0≤n≤N),給出半離散格式:

(3un+1-4un+un-1)/(2τ)=-Mμn+1,

(7)

(8)

(9)

滿足E(un+1,un)≤E(un,un-1).

1.2.2 Crank-Nicolson格式 利用CN格式,取正整數(shù)N,記τ=T/N,tn=nτ(0≤n≤N),給出半離散格式

(un+1-un)/τ=-Mμn+1/2,

(10)

(11)

(12)

滿足E(un+1)≤E(un).

1.3 空間方向的2種離散格式

對方程的空間2階偏導(dǎo)數(shù)部分,采用2階中心差分離散和重心Lagrange插值配點法進行離散.記在[a,b]上m等分后的離散點為a

從而可以得到2階空間離散Δun.其中容易證明Δ是對稱負定矩陣，且該空間離散的收斂階是2階的.由于對稱負定矩陣的特殊性質(zhì)，所以可以用快速算法計算在迭代過程中的線性方程組的解和矩陣向量的乘積部分,具體步驟將在后面的計算過程中給出.

1.3.2 重心Lagrange插值配點法 為了得到更高收斂階的空間離散格式，采用重心Lagrange插值配點法離散空間，并與2階中心差分格式離散空間進行比較.對于在區(qū)間[a,b]上的離散點，函數(shù)u(x)在節(jié)點xj處的函數(shù)值為uj=u(xj),則可以寫出u(x)的重心型插值函數(shù)

(13)

記重心Lagrange插值的微分矩陣為D,利用插值函數(shù)(13)可以寫出在插值節(jié)點處的導(dǎo)數(shù)

通過對Lj(x)求導(dǎo),再利用數(shù)學(xué)歸納法則可以得到D(p)的遞推公式

從而可以得到2階微分矩陣D(2),即Δu=D(2)u.

由于重心Lagrange插值對于Chebyshev節(jié)點有較高的數(shù)值穩(wěn)定性，所以本文將選用第2類Chebyshev節(jié)點及其重心插值權(quán),其計算公式如下:

xj=cos(jπ/n),j=0,1,2,…,m,

從而利用Chebyshev節(jié)點作為Lagrange插值配點，對應(yīng)地就可以給出下面的插值逼近性質(zhì).

定理1[22]若u(x,t) 在[a,b]×(0,T]上連續(xù),且邊界是Lipschitz連續(xù).記h=(b-a)/2,則存在常數(shù)c0、c1有誤差估計如下:

由定理1可知,重心Lagrange插值配點法是指數(shù)收斂的.

1.4 SAV方法計算步驟和技巧

(I/τ-2MΔ/3)un+1+Mβn+1(βn+1,un+1)/3=(4un-un-1)/(3τ)-2Mβn+1(4rn-rn-1-(βn+1,4un-un-1)/2)/9=Qn,

(14)

令A(yù)=I/τ-2MΔ/3,對于中心差分格式離散空間方向,矩陣A顯然是三對角對稱可逆矩陣.則可對式(14)先左乘A-1,再內(nèi)積βn+1,得到

(βn+1,un+1)=(βn+1,A-1Qn)/(1+M(βn+1,A-1·βn+1)/3)=α.

最后計算

un+1=A-1Qn-MA-1βn+1α/3,rn+1=(4rn-rn-1)/3+(βn+1,3un+1-4un+un-1)/6,

(I/τ-MΔ/2)un+1+Mβn+1/2(βn+1/2,un+1)=(I/τ+MΔ/2)un-Mβn+1/2(2rn-(βn+1/2,un))=Qn.

(15)

令A(yù)=I/τ-MΔ/2,因為A是三對角對稱可逆矩陣,所以同理可對式(15)先左乘A-1,再內(nèi)積βn+1/2,得到

(βn+1/2,un+1)=(βn+1/2,A-1Qn)/(1+M(βn+1/2,A-1βn+1/2))=α.

最后計算

un+1=A-1Qn-MA-1βn+1/2α,rn+1=rn+(βn+1/2,un+1-un),

這2種格式在計算的過程中都需要用到前2步的信息，可以先采用1階格式，類似上述步驟利用u0求出u1,再進行迭代.

1.4.2 離散正弦變換和快速傅里葉變換 為了用快速算法求解導(dǎo)出的線性方程組,先介紹特殊矩陣的性質(zhì),有如下2個引理.

引理1[23]對于m×m的三對角矩陣T=tridiag(a,b,c),它的m個特征值為

則可將T分解為T=PΛP-1=Pdiag(λ1,λ2,…,λm)P-1,特征值λk對應(yīng)的特征向量為Pk,其中第j個分量為

Pkj=(a/c)j/2sin(jkπ/(m+1))(j=1,2,…,m).

引理2[23]對于m×m的循環(huán)矩陣C=circle(a0,a1,…,am-1),它的m個特征值為

ωk=cos(2kπ/m)+jsin(2kπ/m)=e2kjπ/m(0≤k≤m-1),

則可將C分解為C=FΛF-1=Fdiag(g(ω0),g(ω1),…,g(ωm-1))F-1,特征值g(ωk)對應(yīng)的特征向量Fk為

由于空間采用2階中心差分離散,且A是三對角對稱Toeplitz可逆矩陣,所以可以利用離散正弦變換(DST)快速計算Ab,A-1b.

先利用引理1計算出A的特征值，令A(yù)=PΛP-1.由DST和IDST的性質(zhì)[23]知

由矩陣A的特征可知,Ab、A-1b可以用DST和IDST快速計算,即

A-1b=PΛ-1P-1b=fDST(Λ-1fIDST(b)),

Ab=PΛP-1b=fDST(ΛfIDST(b)).

將Toeplitz矩陣T擴充成循環(huán)矩陣

再利用離散傅里葉變換(DFT)性質(zhì)知

最后取所得向量的前半部分即可.

2 數(shù)值實驗

2.1 算例1

選取u(x,0)=u0=0.5sin(πx),(x,t)∈[-1,1]×[0,1],取M=1,ε=0.1.由于沒有準(zhǔn)確值，所以取τ=1×10-5,使計算出的un作為參考解u(tn),則計算最后時刻T的誤差e=|un-u(tn)|,收斂階

or=log2(ei+1/ei)/log2(τi+1/τi)

取空間離散點數(shù)m=16,計算時間收斂階如表1所示.

表1 在m=16時的L∞誤差和時間收斂階

從計算結(jié)果可以看出,BDF2格式和CN格式的收斂階都是2階的(見圖1和圖2).在相同空間離散格式下，CN格式的精度比BDF2格式的精度高.

圖1 BDF2格式收斂階

圖2 CN格式收斂階

考慮空間離散方式不一樣,假設(shè)真解u(x,t)=costsin(πx),代入Allen-Cahn方程中計算得到函數(shù)

g(x,t)=sin(πx)(-sint+Mπ2cost+Mcost·((costsin(πx))2-1)/ε2),

取τ=1×10-5,M=1,ε2=0.1,計算結(jié)果如表2所示.

由上述分析結(jié)果可以看出,空間離散使用中心差分的收斂階是2階的.重心Lagrange插值配點法可以在較小的剖分中得到較高的精度，收斂階是指數(shù)遞減的.在相同的時間離散格式下，重心Lagrange插值配點法的精度比2階中心差分格式的精度高很多，前者只需要6個節(jié)點就可以達到后者32個節(jié)點的精度，前者在10個節(jié)點時的誤差就遠小于后者在128個節(jié)點時的誤差.

接下來討論BDF2格式和CN格式的運行時間，運行環(huán)境為8G(intel)i5-7200U.選取N=104,對于差分離散，選取空間剖分數(shù)m=210.因為由上述實驗可知,重心Lagrange插值配點法的計算結(jié)果的誤差較小,所以對于重心Lagrange插值配點法離散選取m=10.由上述實驗可知重心Lagrange插值配點法的誤差較小，下面比較常規(guī)計算和快速計算以及重心Lagrange插值配點法離散的運行時間(見表3).

表3 不同格式的CPU計算時間

由表3可以很明顯看出，快速算法對差分格式的計算時間有較大提升，而重心Lagrange插值配點法能經(jīng)過很少節(jié)點就達到高精度從而縮短了計算時間.

2.2 算例2

取u0=sin(πx),(x,t)∈[-1,1]×[0,1]取M=1,ε=0.1,m=210,τ=1×10-4,1維等距剖分的能量離散形式為

由圖3和圖4可見:在空間是2階中心差分離散空間的情況下，時間離散的BDF2格式和CN格式都滿足能量耗散，且到達穩(wěn)定的時間幾乎一致.

圖3 BDF2/差分格式

圖4 CN/差分格式

取u0=sin(πx),(x,t)∈[-1,1]×[0,1],M=1,ε=0.1,m=12,τ=1×10-4,在Chebyshev節(jié)點下的能量離散形式為

同樣,對于重心Lagrange插值配點法空間的情況，時間離散的BDF2格式和CN格式都滿足能量耗散.如圖5和圖6所示，BDF2格式在t=0.6之前達到穩(wěn)定，而CN格式在t=0.6后達到穩(wěn)定，這說明對于同種情況下，CN格式需要更多時間才達到能量穩(wěn)定.

圖5 BDF2/重心插值格式

圖6 CN/重心插值格式

3 結(jié)束語

本文利用重心Lagrange插值配點法結(jié)合SAV格式來求解Allen-Cahn方程，對比了幾種離散格式的精度和能量耗散情況.通過算例可以發(fā)現(xiàn)：在同種空間離散格式下，CN格式的精度比BDF2格式更高；在同種時間離散格式下，重心Lagrange插值配點格式的精度比2階中心差分的精度更高，且收斂階是指數(shù)收斂.雖然DST和FFT快速算法結(jié)合2階中心差分格式求解所用的CPU時間變少，但是重心Lagrange插值配點法求解所用的CPU時間仍然比快速算法的差分格式求解所用的CPU時間更少.進一步地，數(shù)值算例驗證在SAV方法下的幾種數(shù)值格式都滿足能量耗散規(guī)律.