SAV/重心插值配點法求解Allen-Cahn方程*

黃 蓉, 鄧楊芳, 翁智峰

(華僑大學 數學科學學院, 福建 泉州 362021)

0 引 言

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

(1)

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

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

(2)

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

(3)

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

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

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

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

1 重心Lagrange插值配點法

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

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

(4)

式中, Lagrange插值基函數為

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

(5)

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

(6)

據上式,可推出重心Lagrange插值公式為

(7)

令P(x)關于x求導,并在離散節點xi處成立,即

(8)

(9)

重心插值配點法在選取Chebyshev點族時具有優良的數值穩定性,本文選取的節點和對應的重心權表達式為

(10)

(11)

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

(12)

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

(13)

2 Allen-Cahn方程的SAV配點格式

2.1 Allen-Cahn方程的離散格式推導

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

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

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

(14)

(15a)

(15b)

(15c)

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

(16)

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

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

(17a)

(17b)

(17c)

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

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

(18)

(19)

(20)

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

(21)

整理如下:

(22)

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

(23)

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

(24)

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

(25)

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

(26)

經整理,可得

(27)

故

(28)

(29)

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

(30)

2.2 重心插值公式誤差估計

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

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

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

(31)

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

(32)

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

3 數 值 實 驗

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

3.1 算例1

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

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

(33)

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

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

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

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

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

根據圖1所示,令τ=0.001,0.000 1,兩種SAV配點格式的能量耗散曲線均呈現下降趨勢,表明應用SAV策略的重心Lagrange配點格式滿足能量遞減的性質.

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

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

圖2 不同數值格式的空間收斂階

3.2 算例2

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

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

(34)

選取空間節點數M=N=40,設定參數ε=0.06,Δt=0.002,則數值解隨時間的演化如圖3所示.圖3是在給定初值的條件下,數值解在t=0 s,0.2 s,0.7 s,1 s的演化情況.隨著時間的累積,數值解隨之出現變化層,然后達到亞穩態,最后達到穩態.

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

3.3 算例3

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

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

(35)

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

(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

在本算例中,研究數值解的粗化現象.給定參數ε=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,則數值解在不同時刻的快照如圖5所示.圖5展示了相場函數u的數值解隨著時間變化而逐漸融合的過程:初始時刻呈現“啞鈴狀”圖像;在t=0.1 s演變成橢圓狀;其后在t=0.15 s融合成小型圓點;最后呈現完全融合狀態.

(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 結 束 語

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