求解Allen-Cahn 方程的半隱配點格式

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

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

1 引言

1.1 研究背景

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

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

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

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

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

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

上述方法都是基于網格剖分的思想求解Allen-Cahn 方程.文獻[15]提出一種新型無網格的重心插值配點法,該方法具有計算格式簡單、精度高、程序實施方便、節點適應性好等特點,被廣泛應用于求解各類微分方程,比如Cahn-Hilliard 方程[16]、Helmholtz 方程[17]、Black-Scholes 方程[18]、Allen-Cahn 方程[19]等.最近,Yi 和Yao[20]利用重心插值配點法求解時間分數階電報方程并給出了對應的誤差估計?Li 和Cheng[21]分析了熱傳導方程在重心有理插值配點格式下的誤差估計.

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

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

1.2 重心Lagrange 插值與Chebyshev 配點法

重心Lagrange 插值設有m+1 個插值節點xj和相對應的一組實數fj,j=0,1,···,m.在次數不超過m的多項式空間中,存在唯一一個多項式p(x),滿足p(xj)=fj(j=0,1,···,m).多項式p(x)即為如下的Larange 插值多項式:

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

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

將插值基函數表示為:

再代入式(1.2),得

由式(1.3),(1.4)可得

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

其中

Chebyshev 配點法多項式插值在節點數m較大時,若采用等距節點,其插值是病態的.因此本文將以第二類Chebyshev 節點作為空間方向的插值節點.

在區間[-1,1]上的第二類Chebyshev 節點為:

相應的重心Lagrange 插值的權為:

任意區間[a,b]上的Chebyshev 節點可通過坐標變換y=x(b-a)/2+(a+b)/2 得到.

重心插值的微分矩陣函數p(x)在節點xi處的k階導數為:

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

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

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

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

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

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

將區間[a,b] 離散為m+1 個第二類chebyshev 節點:a=x0

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

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

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

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

下面,給出重心插值配點法的逼近性質.

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

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

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

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

證明

其中,

對R1,我們有

再由引理2.1,可得

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

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

其中C=.證畢.

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

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

3.1 一階半隱全離散格式

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

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

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

3.2 二階半隱全離散格式

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

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

4 二維Allen-Cahn 方程的數值格式

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

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

4.1 二維重心Lagrange 插值半離散系統

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

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

4.2 二維全離散格式

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

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

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

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

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

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

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

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

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

其中,ηm,l是空間截斷誤差.

由式(4.5)和(4.7),有

由非線性項滿足Lipschitz 條件及引理4.1,類似定理2.1 的推導,可得:

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

5 數值實驗

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

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

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).

為了驗證方程(5.1)的BDF/重心插值格式(4.3)和CN/重心插值格式(4.4)的時間收斂階,取參數ε=0.1,M=m+1=L=l+1=10,及不同的時間步長τ,計算結果見表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顯示了在固定參數ε=0.3,τ=10-4,改變空間節點數(M,L)時,時間一階半隱配點格式的絕對誤差和相對誤差結果.由表2 可知,隨著空間節點的增加,數值解的相對誤差越來越小.

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

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

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

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

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

圖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},

其二維能量函數的離散形式為:

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

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

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

由圖2可知,方程(5.2)的兩種數值格式的能量函數都隨時間遞減,且ε越大,能量越早達到穩態.

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

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

5.3 算例3

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

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

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

5.4 算例4

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

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

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

5.5 算例5

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

的Allen-Cahn 方程.

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

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

6 總結

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