SAV/重心插值配點(diǎn)法求解Allen-Cahn方程*

黃 蓉, 鄧楊芳, 翁智峰

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

0 引 言

Allen-Cahn方程是一類求解界面問題的數(shù)學(xué)模型,最早由Allen和Cahn[1]在1979年引入,用來描述晶體反相位邊界的運(yùn)動.此后,Allen-Cahn方程在描述各類復(fù)雜現(xiàn)象中占據(jù)著重要地位,如圖像處理[2]、晶體生長[3]、平均曲率-流量[4]等.本文考慮如下的Allen-Cahn方程:

(1)

式中,u(x,t)是相場函數(shù),F(u)=(u2-1)2/(4ε2),ε>0是界面寬度參數(shù),Ω∈2是有界區(qū)域.

Allen-Cahn方程可以視為自由能泛函的L2梯度流[5]:

(2)

將能量泛函E(u)關(guān)于t求導(dǎo),可知Allen-Cahn方程滿足能量遞減規(guī)律,即

(3)

Allen-Cahn方程具有廣泛的應(yīng)用前景,但其解析解往往難以求出,因此許多學(xué)者致力于其數(shù)值模擬的研究,提出了各種穩(wěn)定高效的數(shù)值模擬方法.Zhai等[6]采用緊致差分交替方向隱式(ADI)方法求解高維Allen-Cahn方程;Weng等[7]基于算子分裂法求解Allen-Cahn 方程,并給出了格式的穩(wěn)定性分析;Li 等[8]采用有限元法建立Allen-Cahn方程的無條件能量穩(wěn)定格式,并對數(shù)值格式進(jìn)行了誤差分析;Liao等在文獻(xiàn)[9]中構(gòu)建了變步長的二階向后差分格式(BDF2)來求解Allen-Cahn方程; Li和Song等[10]采用非線性內(nèi)罰間斷Galerkin有限元(IPDGFE)方法求解Allen-Cahn方程,建立了時間二階精度、能量穩(wěn)定的高效數(shù)值格式; Li和Ju等[11]提出了守恒型Allen-Cahn方程的極大值原理分析;汪精英等[12]基于Laplace變換,結(jié)合譜方法和有限差分法求解了分?jǐn)?shù)階Cahn-Hilliard相場方程.

目前,無網(wǎng)格方法[13-14]逐漸引起學(xué)者們的關(guān)注,重心Lagrange插值配點(diǎn)法是一種新型的無網(wǎng)格方法,由Lagrange插值公式引入重心權(quán)推導(dǎo)而來,可以克服后者數(shù)值不穩(wěn)定的缺陷.當(dāng)節(jié)點(diǎn)數(shù)量增多時,采用Lagrange插值公式來逼近函數(shù)時易出現(xiàn)Runge現(xiàn)象,節(jié)點(diǎn)適應(yīng)性差;而重心Lagrange插值公式在選取特殊的節(jié)點(diǎn)時,如Chebyshev點(diǎn)族,具有優(yōu)良的數(shù)值穩(wěn)定性,解決了Lagrange插值公式逼近函數(shù)時出現(xiàn)的震蕩現(xiàn)象,提高了節(jié)點(diǎn)的適應(yīng)性.此外,重心Lagrange插值配點(diǎn)法兼具精度高、編程簡單、易于計算等優(yōu)勢.近年來,重心Lagrange插值配點(diǎn)法已被應(yīng)用于求解各類問題,如平面彈性問題[15]、Burgers方程[16]、最優(yōu)控制問題[17]、Allen-Cahn方程[18]、高維Fredholm積分方程[19]等.盡管關(guān)于重心Lagrange插值配點(diǎn)法數(shù)值應(yīng)用的研究很多,但目前該方法的收斂性分析、誤差分析的理論工作還較少.最近,Yi 等[20]提出了求解時間分?jǐn)?shù)階電報方程的重心Lagrange插值方法,并給出了離散系統(tǒng)的誤差分析.

近年來提出的標(biāo)量輔助變量(SAV)方法[21-25]是一種新型高效的處理非線性問題的策略.該方法基于IEQ格式改進(jìn),保留了IEQ 格式的優(yōu)勢并克服其缺陷,為梯度流模型構(gòu)建時間離散化方案提供了一種高效而穩(wěn)定的策略.SAV方法具備兩大優(yōu)勢:一方面,SAV格式擴(kuò)大了IEQ格式的使用范圍,不再局限于非線性勢函數(shù)有下界的情形下使用;另一方面,SAV方法借助引入不依賴空間變量的輔助變量,在時間步只涉及到常系數(shù),構(gòu)建線性無條件能量穩(wěn)定的離散系統(tǒng).目前,SAV方法已被應(yīng)用于求解各種問題.Huang等[22]為梯度流模型提出了一種高效且準(zhǔn)確的新型SAV方法.Cheng等[23]提出了廣義SAV(G-SAV)方法,該方法可保持原始能量而非修改能量的無條件能量穩(wěn)定,同時也能保持原始SAV方法的基本優(yōu)勢.Li和Shen[24]基于SAV Fourier譜格式求解相場晶體方程,并給出了無條件能量穩(wěn)定格式的誤差分析.

基于前人的工作,本文著重于空間離散,采用新型的重心插值配點(diǎn)格式來構(gòu)建二維Allen-Cahn方程的SAV無條件能量穩(wěn)定數(shù)值格式.在時間方向上采用的離散方法是Crank-Nicolson(CN2)和BDF2.此外,給出了重心Lagrange插值配點(diǎn)格式的逼近性質(zhì).為驗證所提格式的高精度,將重心插值配點(diǎn)格式與有限差分格式進(jìn)行對比,數(shù)值算例驗證了基于SAV方法的重心Lagrange配點(diǎn)格式具有指數(shù)收斂的特性.

1 重心Lagrange插值配點(diǎn)法

本文采用一種新型的無網(wǎng)格方法——重心Lagrange插值配點(diǎn)法,對空間方向進(jìn)行離散.對于M+1個互異節(jié)點(diǎn)xk(k=0,1,…,M),令函數(shù)V(x)在節(jié)點(diǎn)xk處的函數(shù)值是Vk,插值多項式P(x)在離散節(jié)點(diǎn)xk處成立P(xk)=Vk.

將多項式P(x)改寫成Lagrange插值形式:

(4)

式中, Lagrange插值基函數(shù)為

記L(x)=(x-x0)(x-x1)…(x-xM),則基函數(shù)Lk(x)的另一表達(dá)式為

(5)

將式(5)代入式(4),并令P(x)=1,可得

(6)

據(jù)上式,可推出重心Lagrange插值公式為

(7)

令P(x)關(guān)于x求導(dǎo),并在離散節(jié)點(diǎn)xi處成立,即

(8)

(9)

重心插值配點(diǎn)法在選取Chebyshev點(diǎn)族時具有優(yōu)良的數(shù)值穩(wěn)定性,本文選取的節(jié)點(diǎn)和對應(yīng)的重心權(quán)表達(dá)式為

(10)

(11)

針對選取的第二類Chebyshev點(diǎn)族,令空間x,y方向的節(jié)點(diǎn)數(shù)分別為M,N.類似上述推導(dǎo)過程,則二元函數(shù)u(x,y)在張量型節(jié)點(diǎn)(xk,yj)的重心Lagrange插值公式為

(12)

考慮u(x,y)分別對變量x,y求二階偏導(dǎo)數(shù),即

(13)

2 Allen-Cahn方程的SAV配點(diǎn)格式

2.1 Allen-Cahn方程的離散格式推導(dǎo)

在本小節(jié)中,構(gòu)造二維Allen-Cahn方程的全離散格式.應(yīng)用SAV策略,在時間方向采用CN2格式離散,空間方向采用重心Lagrange插值配點(diǎn)格式離散,建立時間二階精度的SAV數(shù)值格式,記為CN2-BLI.為便于分析,定義(·,·)和‖·‖分別表示L2內(nèi)積和L2范數(shù).

為構(gòu)建SAV方案,在Allen-Cahn方程中引入變量ω,將式(1)改寫為

ω=-Δu+F′(u).

(14)

(15a)

(15b)

(15c)

分別將式(15a)內(nèi)積ω、式(15b)內(nèi)積?u/?t、式(15c)乘以2R,可得

(16)

式(16)表明修正能量E(u)=(u,-Δu)/2+R2是無條件能量穩(wěn)定的.

在時間方向采用CN2格式進(jìn)行離散,則Allen-Cahn方程的半離散格式如下:

(17a)

(17b)

(17c)

將式(17a)、式(17b)分別與ωs+1/2,us+1-us作內(nèi)積,式(17c)乘以2rs+1/2,可得

(us+1-us,ωs+1/2)+τ(ωs+1/2,ωs+1/2)=0,

(18)

(19)

(20)

將式(18)—(20)結(jié)果化簡,可得

(21)

整理如下:

(22)

綜上,Allen-Cahn方程的半離散格式(17a)—(17c)是關(guān)于時間無條件能量穩(wěn)定的.

(23)

式中,空間離散算子Dh=C(2)?IN+IM?D(2),C(2)和D(2)分別是空間x和y方向的二階微分矩陣.

(24)

記A=I-(τ/2)Dh,上式方程兩邊同時乘以A-1,即

(25)

將式(25)與σs作內(nèi)積,即

(26)

經(jīng)整理,可得

(27)

故

(28)

(29)

類似地,若采用BDF2格式結(jié)合重心Lagrange配點(diǎn)法來離散Allen-Cahn方程,構(gòu)造基于SAV方案的另一種二階離散格式,記為BDF2-BLI,算法格式如下:

(30)

2.2 重心插值公式誤差估計

基于插值逼近誤差結(jié)論,分析重心Lagrange配點(diǎn)格式求解Allen-Cahn方程的相容性誤差.

針對函數(shù)u(x,y),對應(yīng)的重心Lagrange插值函數(shù)為uh(x,y),誤差是η(x,y)=u(x,y)-uh(x,y).引用文獻(xiàn)[20]定理內(nèi)容,即有如下引理.

引理1 如果u(x,y)∈C(n+1)(Ω),Ω=[a,b]×[c,d]是具有Lipschitz連續(xù)邊界的非空區(qū)域,則如下的誤差估計結(jié)論成立:

(31)

由引理1結(jié)論,設(shè)函數(shù)u(x,y,t)的重心Lagrange插值逼近函數(shù)為uh(x,y,t),則以下結(jié)論成立:

(32)

根據(jù)重心插值逼近的誤差結(jié)果可知,該配點(diǎn)格式是在空間上具有指數(shù)收斂的特性,并且微分算子的階數(shù)決定了算法求解方程的空間收斂階.

3 數(shù) 值 實 驗

在本節(jié)中,選取Dirichlet邊界條件,給出4個數(shù)值算例來驗證配點(diǎn)格式的有效性與高精度.針對二維Allen-Cahn方程,分別驗證兩種SAV配點(diǎn)格式均具有高精度,且滿足無條件能量穩(wěn)定特性;此外,分別采用重心插值配點(diǎn)法和五點(diǎn)差分法對空間方向離散,將數(shù)值結(jié)果對比分析.

3.1 算例1

為驗證CN2-BLI、BDF2-BLI格式的時間收斂階,且滿足能量遞減規(guī)律,選取真解如下:

u=cos(t)sin(πx)sin(πy).

(33)

令Ω=[-1,1]2×(0,T],ε=0.2,M=N=20,T=1,則兩種數(shù)值格式求解Allen-Cahn方程的誤差、收斂階結(jié)果見表1.由表1可知:CN2-BLI、BDF2-BLI兩種數(shù)值格式的時間收斂階均為二階;在選取相同的剖分節(jié)點(diǎn)時,CN2-BLI格式求解得到的L2誤差小于BDF2-BLI格式的L2誤差,表明前者的數(shù)值求解效果略優(yōu)于后者.

表1 采用CN2-BLI、BDF2-BLI格式求解u的L2誤差

基于SAV策略,分別選擇不同的空間離散技術(shù)求解Allen-Cahn方程,即重心Lagrange插值配點(diǎn)法、五點(diǎn)差分法來離散,得到的L2誤差結(jié)果見表2.由表2可知:CN2-BLI、BDF2-BLI格式在空間上采用重心Lagrange插值配點(diǎn)法離散,當(dāng)選取節(jié)點(diǎn)M=N=15時,L2誤差均可達(dá)到10-8量級;CN2-FD格式空間離散方式是五點(diǎn)差分法,在M=N=120時,誤差僅達(dá)到10-4量級.通過對比誤差結(jié)果可知,前兩種格式在選取少量節(jié)點(diǎn)時即可達(dá)到高精度,即重心Lagrange插值配點(diǎn)法具有高精度的優(yōu)勢.

表2 不同空間離散方案的精度對比結(jié)果

對于上述剖分,驗證不同數(shù)值格式求解Allen-Cahn方程時滿足能量遞減規(guī)律.

根據(jù)圖1所示,令τ=0.001,0.000 1,兩種SAV配點(diǎn)格式的能量耗散曲線均呈現(xiàn)下降趨勢,表明應(yīng)用SAV策略的重心Lagrange配點(diǎn)格式滿足能量遞減的性質(zhì).

圖1 不同數(shù)值格式的能量演化圖(ε=0.1)

令參數(shù)ε=0.3,τ=0.001,T=1,選取不同的空間剖分節(jié)點(diǎn)數(shù)(M=N),則AC方程的空間收斂階結(jié)果如圖2所示.圖2表明當(dāng)空間采用重心Lagrange插值配點(diǎn)法離散時,CN2-BLI、BDF2-BLI格式在空間方向上均具有指數(shù)收斂的特性.

圖2 不同數(shù)值格式的空間收斂階

3.2 算例2

對于問題(1),令Ω=[-1,1]2×(0,1],給定初值如下:

u0(x,y)=(x2-1)(y2-1)sin(πx)sin(πy).

(34)

選取空間節(jié)點(diǎn)數(shù)M=N=40,設(shè)定參數(shù)ε=0.06,Δt=0.002,則數(shù)值解隨時間的演化如圖3所示.圖3是在給定初值的條件下,數(shù)值解在t=0 s,0.2 s,0.7 s,1 s的演化情況.隨著時間的累積,數(shù)值解隨之出現(xiàn)變化層,然后達(dá)到亞穩(wěn)態(tài),最后達(dá)到穩(wěn)態(tài).

(a)t=0 s (b)t=0.2 s (c)t=0.7 s (d)t=1 s

3.3 算例3

針對Allen-Cahn方程, 考慮選取隨機(jī)初值如下:

u0(x,y)=0.95rand(x,y)+0.05.

(35)

令區(qū)域為Ω=[-1,1]2×(0,1],選取參數(shù)ε=0.05,τ=0.000 1,M=N=40, 則u在不同時刻的數(shù)值解如圖4所示.圖4展示了在初值隨機(jī)分布的情形下,SAV配點(diǎn)格式求解Allen-Cahn方程在t=0 s,0.03 s,0.08 s,0.5 s,1 s時的數(shù)值解快照,隨著時間變化,相函數(shù)從分離狀態(tài)逐漸達(dá)到聚合狀態(tài).

(a)t=0 s (b)t=0.03 s (c)t=0.08 s (d)t=0.5 s (e)t=1 s

3.4 算例4

在本算例中,研究數(shù)值解的粗化現(xiàn)象.給定參數(shù)ε=0.05,τ=0.002,T=0.2,初值如下:

u0(x,y)=h1(x,y)h2(x,y),

(36)

式中

(37)

(38)

令Ω=[-1,1]2×(0,T],空間剖分為M=N=60,則數(shù)值解在不同時刻的快照如圖5所示.圖5展示了相場函數(shù)u的數(shù)值解隨著時間變化而逐漸融合的過程:初始時刻呈現(xiàn)“啞鈴狀”圖像;在t=0.1 s演變成橢圓狀;其后在t=0.15 s融合成小型圓點(diǎn);最后呈現(xiàn)完全融合狀態(tài).

(a)t=0 s (b)t=0.05 s (c)t=0.1 s (d)t=0.15 s (e)t=0.2 s

4 結(jié) 束 語

本文主要構(gòu)建了基于SAV方案的二維Allen-Cahn方程全離散格式.在時間方向上采用的離散方法是CN2格式、BDF2格式, 空間上采用重心Lagrange插值配點(diǎn)法離散,應(yīng)用SAV策略,得到線性無條件能量穩(wěn)定的數(shù)值格式.數(shù)值實驗結(jié)果顯示:兩種SAV 配點(diǎn)格式均具有高精度的特性,在選取第二類Chebyshev點(diǎn)時具有指數(shù)收斂的特性.與經(jīng)典有限差分格式的誤差結(jié)果比較可知,重心Lagrange 配點(diǎn)格式采用非常少的節(jié)點(diǎn)數(shù)即可達(dá)到更高精度.