二维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方程，通过数值算例验证格式中时间方向一阶精度，空间方向的高精度,以及满足方程的能量递减性质.