二維Allen-Cahn方程的有限差分法/配點法求解

鄧楊芳，姚澤豐，汪精英，翁智峰

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

Allen-Cahn方程是一種非齊次半線性泊松方程，應用范圍相當廣泛，例如,晶體生長[1]、圖像分析[2]、平均曲率-流量[3]、生物種群的競爭、河床的遷移過程[4]、排斥現象[5]等許多擴散現象的研究.

Allen-Cahn方程有很多的數值解法，如有限差分法[6-7]、有限元法[8-10]、譜方法[11].Chen等[6]給出求解Allen-Cahn方程差分格式的收斂性的分析.Zhai等[7]利用高階緊致ADI差分格式求解三維Allen-Cahn方程.Feng等[8]提出對稱內罰間斷有限元格式求解帶有平均曲率流Allen-Cahn方程.Li等[9]利用無條件能量穩定的二階有限元方法求解Allen-Cahn方程.吳龍淵等[10]利用交替方向方法求解Allen-Cahn方程.文獻[11-13]利用算子分裂方法結合有限差分法、譜方法、有限元法求解Allen-Cahn方程.上述方法都是利用網格剖分求解微分方程問題.李淑萍等[14]利用重心插值配點法求解微分方程初值問題，很多學者將該方法推廣到求解各類微分方程，如平面彈性問題[15]和分數階Fredholm積分方程[16]等.翁智峰等[17]對時間方向和空間方向都采用重心插值配點法求解一維Allen-Cahn方程.本文將該方法推廣到求解二維Allen-Cahn方程，時間方向采用有限差分法離散，空間方向采用重心插值Chebyshev配點法離散.

1 一般的二維偏微分方程的離散化及推導

考慮一個矩形區域Ω=[a,b]×[c,d]，令a=x0

對于二維的偏微分方程,固定y,考慮變量x的方向，重心Lagrange插值公式為

(1)

可以推得如下表達式，即

(2)

令

(3)

為x方向的重心插值基函數，則式(2)為

(4)

同理，固定x,考慮變量y的方向，式(2)可為

(5)

式(5)中：βj(y)是y方向的重心插值基函數.

由式(4),(5),u(x,y)在節點{(xi,yi),i=0,1,2,…,n，j=0,1,2,…,m}上的重心插值為

(6)

推得u(x,y)的l+k階偏導數為

(7)

偏導數在節點(xp,yq)處的函數值近似為

(8)

所以，二維偏微分方程的重心插值配點法的格式為

(9)

式(9)中：fi,j為f(x,y)在節點(xi,yj)處的函數值.式(9)用微分矩陣表示為

(D(2,0)+D(0,2))ui,j=fi,j.

(10)

2 二維Allen-Cahn方程在重心插值配點法下的計算格式

二維Allen-Cahn方程為

(11)

由于自變量有3個，可以先將時間t離散化.在t固定后，即可利用重心插值配點法求解xy空間上的數值解.對時間向后差分，得到的變形式為

(12)

對非線性項Talor展開后格式變形為

(13)

由式(13)獲得重心插值配點法的計算格式為

(14)

3 數值算例

為便于分析，定義絕對誤差Erra和相對誤差Errr分別為

(15)

式(15)中：yc，ye分別為數值解和解析解的列向量；‖·‖2為向量的二范數.

3.1 算例1

為了驗證算法的精度，只對一維Allen-Cahn方程驗證算法的收斂階.一維Allen-Cahn方程為

(16)

式(16)的精確解為

(17)

首先，對一維Allen-Cahn方程的時間項向后差分，轉化為

(18)

對非線性項Talor展開后，可變形為

(19)

依據上述的格式，可以直接得到重心插值配點法的計算格式為

(20)

利用重心Lagrange插值配點法計算離散式(8)，Δt取不同值時，n=20，n=30的相對誤差,以及時間的收斂階，如表1所示.表1中：時間收斂階為Rate=(log2eΔt1-log2eΔt2)/(log2Δt1-log2Δt2).由表1可知：時間收斂階應為1階.

固定Δt=0.001，x∈[-1,1]對x軸進行不同程度的剖分，可以得到空間方向的相對誤差和絕對誤差數據，計算結果，如表2所示.由表2可知：空間方向誤差具有高精度，而且當剖分的細密程度增加時，誤差也會相對的變得更小.

表1 相對誤差及時間收斂階Tab.1 Relative errors and orders of time convergence

表2 空間方向誤差Tab.2 Error in space

3.2 算例2

對時間向后差分，在時間t固定的情況下，對xy空間網格進行剖分.取n=30,m=30進行剖分后，在t=1時，二維Allen-Cahn方程截面圖像,如圖1所示.圖1中：uc為數值解.

(a) ε=0.3 (b) ε=0.1圖1 二維Allen-Cahn方程的截面圖像Fig.1 Cross section image of two dimension Allen-Cahn equation

3.3 算例3

二維Allen-Cahn方程的能量函數為

(21)

由式(21)，將xy網格取n=40,m=40進行剖分，取不同的ε值，能量遞減圖像,如圖2所示.

(a) ε=0.05 (b) ε=0.10

(c) ε=0.30 (d) ε=0.40圖2 能量遞減圖像Fig.2 Energy decline image

圖2中：E為能量；t為時間.由圖2可知：不論ε取何值，格式滿足能量遞減規律，并且隨著網格剖分細密程度的加深，能量函數E(u)的曲線將會變得更加光滑.

4 結束語

對時間方向向后差分離散，空間方向采用重心插值配點法離散求解二維Allen-Cahn方程，通過數值算例驗證格式中時間方向一階精度，空間方向的高精度,以及滿足方程的能量遞減性質.