張榮培, 霍俊蓉, 楊程程

(沈陽(yáng)師范大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院, 沈陽(yáng) 110034)

0 引 言

為描述晶體中的反相位邊界運(yùn)動(dòng), Allen和Cahn于1979年提出Allen-Cahn方程[1],可以用于描述鐵合金溶液冷卻過(guò)程中結(jié)晶固體相分離運(yùn)動(dòng)。它是一個(gè)二階非線性拋物方程,對(duì)于不同的物理體系具有不同的自由能形式。Allen-Cahn方程及其各種變形在圖像分析[2]、晶體的生長(zhǎng)[3]、平均曲率運(yùn)動(dòng)[4]和隨機(jī)擾動(dòng)[5]等實(shí)際問(wèn)題中都發(fā)揮著極為重要的作用。由于此類(lèi)相場(chǎng)模型沒(méi)有真解,所以采用有效的數(shù)值方法來(lái)進(jìn)行模擬就顯得尤為重要。

近年來(lái),有很多專(zhuān)家和學(xué)者對(duì)Allen-Cahn 方程的數(shù)值逼近方法進(jìn)行了研究,包括有限元方法[6],有限差分法[7]、譜方法[8-9]、間斷有限元方法[10]和無(wú)網(wǎng)格方法[11]等。另外,在處理高維問(wèn)題時(shí),算子分裂方法[12]也是一種求解復(fù)雜問(wèn)題的有效策略。

本文應(yīng)用二階中心差分方法離散Allen-Cahn方程,利用Kronecker積寫(xiě)出二維拉普拉斯算子的微分矩陣,得到一組非線性常微分方程組。接下來(lái)發(fā)展積分因子方法進(jìn)行時(shí)間離散。在時(shí)間離散過(guò)程中,應(yīng)用Kronecker積的性質(zhì)將微分矩陣進(jìn)行譜分解,結(jié)合快速余弦變換,可以快速地進(jìn)行時(shí)間離散。

1 數(shù)值方法

本文考慮以下形式的非線性Allen-Cahn方程的數(shù)值解:

(1)

邊界和初始條件如下:

1.1 空間離散

設(shè)u在網(wǎng)格節(jié)點(diǎn)(xi,yj)的數(shù)值解為ui,j,由齊次Neumann邊界條件,在網(wǎng)格內(nèi)部點(diǎn)(xi,yj)和邊界處的二階導(dǎo)數(shù)差分格式分別為

定義離散解ui,j,1≤i≤Nx,1≤j≤Ny為Nx×Ny階的矩陣U。由差分格(4)式和(5)式,得到該方程在x,y方向的微分矩陣分別為Bx和By,其中矩陣Bx為Nx×Nx階的矩陣,(Bx)i,i=2,i=2,…,Nx-1,(Bx)i,i-1=-1,i=2,…,Nx,(Bx)i,i+1=-1,i=1,…,Nx-1,(Bx)1,1=(Bx)Nx,Nx=1,矩陣Bx中其他位置的元素均為0。同理有By為Ny×Ny階的矩陣。

設(shè)Ix,Iy分別為Nx,Ny階單位矩陣,將解矩陣向量化后得

U=vec(U)=(u1,1…uNx,1,u1,2…uNx,2,…,u1,Ny…uNx,Ny)T

(6)

利用(4)式～(5)式以及Kronecker積的定義,可以將(1)式的二階中心差分格式寫(xiě)成如下形式:

(7)

證明 首先考慮另外一個(gè)特征值問(wèn)題:

(8)

μ是特征方程的特征值,v為相應(yīng)的特征函數(shù)。將求解區(qū)域離散成Nx個(gè)網(wǎng)格,

(9)

(8)式的特征方程為

λ2+μ=0

(10)

2) 當(dāng)μ=0時(shí),特征方程通解為v=C1+C2ξ,由條件(8),可得v=C1,求得特征值μ=0。

yj=(cosjπξ1,…,cosjπξNx)T,j=0,1,…,Nx-1

(11)

利用三角函數(shù)的和差化積公式可以驗(yàn)證Bxyj=λjyj,j=0,1,…,Nx-1,見(jiàn)附錄。定理1證畢。

(12)

同理,對(duì)By也可得類(lèi)似的結(jié)果。

(13)

(14)

同理有

By?Ix=(Cy?Cx)-1(Λy?Ix)(Cy?Cx)

(16)

結(jié)合(15)式和(16)式可以得到

(17)

(18)

1.2 時(shí)間離散

下面應(yīng)用積分因子法對(duì)(18)式進(jìn)行求解,將(18)式兩端同時(shí)左乘e-At,并從tn到tn+1=tn+Δt進(jìn)行積分得到

(19)

(20)

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

應(yīng)用本文提出的有限差分法和積分因子法求解下面的數(shù)值算例。考慮二維計(jì)算區(qū)域Ω=[-1,1]2,針對(duì)具有齊次Neumann邊界條件的方程(1)進(jìn)行求解。初值選取形狀像啞鈴的函數(shù),見(jiàn)圖1(a)所示。

選取參數(shù)ε=0.05,并在每個(gè)空間方向上取512個(gè)等分點(diǎn),時(shí)間步長(zhǎng)Δt=5×10-5。選取t1=3.53×10-2,t2=7.98×10-2,t3=1.18×10-1這3個(gè)時(shí)間點(diǎn),數(shù)值結(jié)果列在圖1。

圖1 算例在t=0,3.53×10-2,7.98×10-2,1.18×10-1時(shí)的數(shù)值解Fig.1 Numerical solution of the case t=0,3.53×10-2,7.98×10-2,1.18×10-1

圖1顯示了在不同時(shí)間點(diǎn)上Allen-Cahn方程解的數(shù)值結(jié)果。結(jié)果表明,隨著時(shí)間增大,其形狀由初始時(shí)間t0的啞鈴狀不斷向圓形聚攏,數(shù)值結(jié)果與文獻(xiàn)[11]中的結(jié)果吻合。

3 結(jié) 語(yǔ)

本文應(yīng)用快速離散余弦變換結(jié)合積分因子方法求解非線性Allen-Cahn方程。該方法能夠快速地將空間離散得到的非線性常微分方程組進(jìn)行時(shí)間離散,得到方程的數(shù)值解,提高計(jì)算效率。通過(guò)求解齊次Neumann邊界條件下的非線性Allen-Cahn方程,可以發(fā)現(xiàn)隨時(shí)間推移數(shù)值結(jié)果形狀的變化規(guī)律。

附錄驗(yàn)證定理1中,Bxyj=λjyj,j=0,1,…,Nx-1,其中λj=2-2cos(jπ/Nx),yj=(cosjπξ1,…,cosjπξNx)T,j=0,1,…,Nx-1,由網(wǎng)格中心點(diǎn)的定義ξi=(i-1/2)h,i=1,2,…,Nx,在左右各增加一個(gè)虛擬網(wǎng)格,其中心點(diǎn)坐標(biāo)為ξ0=-h/2,ξNx+1=1+h/2。

利用和差化積公式cos(jπξi-1)=cos(jπξi)cos(jπ/Nx)+sin(jπξi)sin(jπ/Nx)和cos(jπξi+1)=cos(jπξi)cos(jπ/Nx)-sin(jπξi)sin(jπ/Nx),可得-cos(jπξi-1)+2cos(jπξi)-cos(jπξi+1)=(2-2cos(jπ/Nx))cos(jπξi)。因?yàn)閏os(jπξ0)=cos(jπξ1),cos(jπξNx)=cos(jπξNx+1),于是cos(jπξ1)-cos(jπξ2)=-cos(jπξ0)+2cos(jπξ1)-cos(jπξ2)=(2-2cos(jπ/Nx))cos(jπξ1),-cos(jπξNx-1)+cos(jπξNx)=-cos(jπξNx-1)+2cos(jπξNx)-cos(jπξNx+1)=(2-2cos(jπ/Nx))cos(jπξNx)。

驗(yàn)證完畢。