1. 引言
Allen-Cahn方程是一类非齐次半线性泊松方程。1979年，为了描述晶体中反相位边界运动，Allen和Cahn [1] 引入了Allen-Cahn方程。该方程是材料科学中描述流体动力学问题和反应扩散问题的一类重要方程，且在描述生物种群的竞争与排斥现象 [2] 、河床的迁移过程 [3] 等许多扩散现象的研究中也提出同样的数学模型。对Allen-Cahn方程的研究是在上个世纪七十年代以后开始的。
设基本能泛函为E(u)，是由u离开物相的惩罚项F(u)和扩散项∇u决定，可表示为：
(1)
其中
关于u求导得到
。作图得F(u)的图像如图1所示。从而
为图像的稳定点，也是极小值点。
通过对式(1)进行变分、Taylor展开从而可得到如下Allen-Cahn方程：
(2)
Allen-Cahn方程广泛运用于处理各种问题，例如图像分析 [4] ， [5] ，平均曲率-流量 [6] ，和晶体生长 [7] 。人们在对Allen-Cahn方程进行数值计算时，采用算子分裂算法 [8] 进行方程的求解计算，能够将一个复杂的算子分裂成几个较简单的子算子之积，从而把一个复杂的数学物理问题分裂成一些简单的问题来求解。它既适用于典型的双曲型方程和抛物型方程，也适用于更为复杂的方程的初边值问题，且分裂后的方程更加容易求解、格式灵活、稳定性好。
经典的算子分裂格式有Marchuk-Strang分裂格式、Trotter-Lie分裂格式 [9] ， [10] ， [11] ， [12] 以及对称加权分裂格式 [13] ， [14] 等。本文采用二阶精度的Marchuk-Strang算子分裂格式 [15] ，将原始方程分裂为线性部分和非线性部分。其中线性部分的数值解算子记为
，非线性部分的精确解算子记为
，则上述问题可通过以下三步格式进行求解，即
(3)
该算子分裂算法的主要优点是算子
与算子
可以分别通过不同的数值计算方法进行求解。这里我们对线性部分
采用构造差分的方法，实验证明Crank-Nicolson格式 [16] 是最为稳定，且误差值是最小的。本文采用的高精度有限差分方法是紧致差分法。最后，通过数值实验比较了两种格式数值解差异及运行效率。
2. 算子分裂法求解Allen-Cahn方程
2.1. 两种求解Allen-Cahn方程的高效算子分裂方法
运用算子分裂方法 [17] 求解Allen-Cahn方程数值解具体过程如下：
1) 计算前半步算子
，其中
；
2) 计算完整一步算子
，其中
，且
；
3) 计算后半步算子
，其中
，且
；
在Allen-Cahn方程
中，根据上述算子分裂算法，可将方程分解为：
热传导方程
(4)
和非线性方程
(5)
因此，Allen-Cahn方程的求解可通过3个子步实现 [18] ，即：
(6)
或者
(7)
其中热传导方程
可通过下文2.3及2.4节的格式计算得到。非线性方程
可通过如下解析式 [19] 得到：
(8)
2.2. 二阶中心差分格式计算热传导方程
取空间节点数为N，空间步长为
，空间节点可表示为
。取时间节点数为M，时间步长为
，时间节点可表示为
。同时记
。
令：
(9)
则方程(4)有如下CN格式：
(10)
容易验证上式精度为
，进一步化简可得：
(11)
为进一步验证格式的有效性，将对其做收敛性分析。
定理1. 差分方程按谱范数稳定的充分必要条件是对于任意均成立
，
，均成立：
证明：
令
，式(11)可化为：
利用Fourier方法，令
，得：
由此得到增长因子：
从而
，此格式无条件稳定。
由于
是精确求解，所以此定理表明整个分裂算法的格式都是稳定的。
2.3. 四阶紧致差分格式计算热传导方程
上文2.2节所用格式得到的精度为
，为了提高计算精度，本节引入如下四阶帕德逼近格式：
(12)
结合(12)及CN格式，可得求解
的如下高阶紧致差分格式：
(13)
容易验证上式精度为
，进一步可化简得：
(14)
为进一步验证格式的有效性，将对其做收敛性分析。
定理1. 差分方程按谱范数稳定的充分必要条件是对于任意均成立
，
，均成立：
证明：
令
，式(14)可化为：
利用Fourier方法，令
，得：
由此得到增长因子：
从而
，此格式无条件稳定。
由于
是精确求解，所以此定理表明整个分裂算法的格式都是稳定的。
3. 数值算例
在本节，将通过数值算例来验证算子分裂格式的收敛阶和稳定性。为方便分析，对如下符号进行解释：
CI-ABA：
CI-ABA：
其中
通过二阶中心差分法求得；
CII-ABA：
CII-ABA：
其中
通过二阶中心差分法求得；
3.1. 算例一
为验证收敛阶，本文选取文献 [20] 作为数值算例。Allen-Cahn方程的精确解为：
(15)
其中
。且其区域的取值范围为
，
。
现欲验证所用算法格式ABA与BAB的高效性和准确性。
假设
，当h取充分小值时，得到
，且
(16)
当
取充分小值时，可得到
，且
(17)
由上得到收敛阶的表达式，以下记为Rate。
在2.2节二阶中心差分格式求解热传导方程的基础上，我们将网格以N = 100，200，400，800，1600；M对应取M = 20，40，80，160，320进行剖分，满足
，其中c = 0.2。
在迭代格式CI-ABA和CI-BAB中，通过算子分裂算法分别进行计算得到表1、表2。
Table 1. Calculation results of the CI-ABA scheme
表1. 格式CI-ABA的计算结果
Table 2. Calculation results of the CI-BAB scheme
表2. 格式CI-BAB的计算结果
对比表1和表2，在用二阶中心差分格式求解Allen-Cahn中，随着网格剖分变细，最大误差值
和2-范数
的数量级从
下降到
，且收敛精度逐渐接近预期的算子分裂方法的2阶收敛率。再对比它们CPU时间，在M、N取值较小时，格式CI-BAB计算所用时间比格式CI-ABA长，而随着网格的剖分变细，格式CI-BAB的CPU远小于格式CI-ABA。
同理，在2.4节四阶紧致差分格式求解热传导方程的基础上，我们将网格以N = 20，40，80，160，320；M对应取M = 100，400，1600，6400，25600进行剖分。满足
，其中c = 0.25。
在迭代格式CII-ABA和CII-BAB中，分别计算得到表3、表4。
通过对比表3和表4，在四阶紧致差分格式下求解Allen-Cahn方程，随着网格剖分变细，最大误差值
和2-范数
的数量级从
下降到
，收敛精度逐渐接近预期的算子分裂方法的4阶收敛率。再对比两者的CPU时间，在M、N取较小值时，迭代格式CII-BAB计算所用时间比格式CII-ABA长，而随着网格的剖分变细，格式CII-BAB的CPU远小于格式CII-ABA。
综上，结合实际情况，在求解热传导方程时，由于迭代过程中产生的误差不断累积，如果采用格式CII-BAB进行计算，误差能够被大大减小。
3.2. 算例二
在下面的数值实验中，我们定义能量函数E(u)：
(18)
初值
，取Dirichlet边界条件，左边界值
，右边界值
。
在2.2节的二阶中心差分格式情况下，取
，
。
对此，我们将网格取N = 40，M = 750进行剖分，此时
，
。其中图2和图3为
的运行结果；图4和图5为
的运行结果；图6和图7为
的运行结果。
Table 3. Calculation results of the CII-ABA scheme
表3. 格式CII-ABA的计算结果
Table 4. Calculation results of the CII-BAB scheme
表4. 格式CII-BAB的计算结果
Figure 2. Curve: Energy variation with time in the second-order center difference scheme
图2. ξ = 0.1.二阶中心差分格式能量随时间变化
Figure 3. Numerical solution of the second-order center difference scheme
图3. ξ = 0.1.二阶中心差分格式数值解图像
Figure 4. Curve: Energy variation with time in the second-order center difference scheme
图4. ξ = 0.075.二阶中心差分格式能量随时间变化
Figure 5. Numerical solution of the second-order center difference scheme
图5. ξ = 0.075.二阶中心差分格式数值解图像
Figure 6. Curve: Energy variation with time in the second-order center difference scheme
图6. ξ = 0.125.二阶中心差分格式能量随时间变化
Figure 7. Numerical solution of the second-order center difference scheme
图7. ξ = 0.125.二阶中心差分格式数值解图像
Figure 8. Curve: Energy variation with time in the fourth-order compact difference scheme
图8. ξ = 0.1.四阶紧致差分格式能量随时间变化
Figure 9. Numerical solution of the fourth –order compact difference scheme
图9. ξ = 0.1.四阶紧致差分格式数值解图像
Figure 10. Curve: Energy variation with time in the fourth-order compact difference scheme
图10. ξ = 0.075.四阶紧致差分格式能量随时间变化
Figure 11. Numerical solution of the fourth –order compact difference scheme
图11. ξ = 0. 075.四阶紧致差分格式数值解图像
Figure 12. Curve: Energy variation with time in the fourth-order compact difference scheme
图12. ξ = 0.125.四阶紧致差分格式能量随时间变化
Figure 13. Numerical solution of the fourth –order compact difference scheme
图13. ξ = 0.125.四阶紧致差分格式数值解图像
首先，观察
时，图2表示的能量函数E(u)随着时间t的增加在不断减少，图3的数值解图像表明首先形成边界层，然后到达亚稳态，最后到达稳态。从中我们可以观察到从初始值到亚稳态这一个快速的动力学行为。当
时，比较图2与图4，随着
的变小，要出现图3的行为，所需的网格剖分需越细密；当
时，比较图2与图6，随着
的变大，所需的网格剖分越稀疏，即M、N需取较小的值。
在2.3节的四阶紧致差分格式情况下，我们将网格取N=20，M=750进行剖分，此时
，
，此时
，
。其中图8和图9为
的运行结果；图10和图11为
的运行结果；图12和图13为
的运行结果。
同理，可得到类似结论。
综上，不论二阶中心差分格式或者四阶紧致差分格式，
的取值越小，所取的网格剖分则需越细密，才能减少误差，满足能量函数E(u)呈现如上图所示的递减规律。对比两种差分格式，发现出现同样结果下，四阶紧致差分格式的剖分可以更稀疏，从而可见高阶情况下所需的网格剖分要求不那么高。
4. 结论
本文利用算子分裂算法将Allen-Cahn方程分解为两个子问题，然后用二阶中心差分格式和四阶紧致差分格式求数值解，并利用傅里叶方法给出了稳定性的分析。最后，通过数值算例比较两种格式的误差和CPU时间及验证其数值解满足能量递减规律。在本文中对算子分裂算法列出了详细的计算步骤，通过第一个算例对两种算子分裂格式的优缺点进行了细致的比较，对比分析之后，我们发现采用格式CII-BAB的误差最小，且能有效减少计算时间。通过第二个算例，说明不论采用哪种格式进行迭代计算，得到的能量函数都满足递减的规律，其数值解也呈现对应的变化。本文对Allen-Cahn方程的求解方法，也能为我们今后解决同类问题提供了简单的途径和有力的支持。
基金项目
国家自然科学基金(11526094)和福建省青年基金(2016J05007)资助项目。