Maximum principle preserving exponential time differencing schemes for the nonlocal Allen-Cahn equation

TL;DR

This paper introduces and analyzes exponential time differencing schemes for the nonlocal Allen-Cahn equation, ensuring maximum principle preservation, stability, and convergence to classical solutions as parameters tend to zero.

Contribution

The paper develops unconditionally maximum principle preserving ETD schemes for the nonlocal Allen-Cahn equation with proven error estimates and asymptotic compatibility.

Findings

Schemes are unconditionally maximum principle preserving.

Error estimates are optimal in maximum norm.

Schemes are energy stable and asymptotically compatible.

Abstract

The nonlocal Allen-Cahn (NAC) equation is a generalization of the classic Allen-Cahn equation by replacing the Laplacian with a parameterized nonlocal diffusion operator, and satisfies the maximum principle as its local counterpart. In this paper, we develop and analyze first and second order exponential time differencing (ETD) schemes for solving the NAC equation, which unconditionally preserve the discrete maximum principle. The fully discrete numerical schemes are obtained by applying the stabilized ETD approximations for time integration with the quadrature-based finite difference discretization in space. We derive their respective optimal maximum-norm error estimates and further show that the proposed schemes are asymptotically compatible, i.e., the approximate solutions always converge to the classic Allen-Cahn solution when the horizon, the spatial mesh size and the time step size go to zero. We also prove that the schemes are energy stable in the discrete sense. Various experiments are performed to verify these theoretical results and to investigate numerically the relation between the discontinuities and the nonlocal parameters.

Tables

Table 1. Table 1: Temporal convergence rates in the maximum-norm sense in Example 6.1 .

$\tau =0.05$		$\alpha =1$ (integrable kernel)				$\alpha =3$ (non-integrable kernel)
		$\delta =0.2$		$\delta =2$		$\delta =0.2$		$\delta =2$
		Error	Rate	Error	Rate	Error	Rate	Error	Rate
ETD1	$\tau$	1.082e-2	$-$	1.090e-2	$-$	1.084e-2	$-$	1.087e-2	$-$
	$\tau /2$	5.535e-3	0.9670	5.580e-3	0.9666	5.545e-3	0.9669	5.561e-3	0.9668
	$\tau /4$	2.800e-3	0.9833	2.823e-3	0.9831	2.805e-3	0.9833	2.813e-3	0.9832
	$\tau /8$	1.408e-3	0.9916	1.420e-3	0.9915	1.410e-3	0.9916	1.415e-3	0.9915
	$\tau /16$	7.060e-4	0.9957	7.121e-4	0.9957	7.074e-4	0.9957	7.095e-4	0.9957
	$\tau /32$	3.536e-4	0.9979	3.566e-4	0.9978	3.542e-4	0.9979	3.553e-4	0.9979
	$\tau /64$	1.769e-4	0.9989	1.784e-4	0.9989	1.772e-4	0.9989	1.778e-4	0.9989
	$\tau /128$	8.849e-5	0.9995	8.924e-5	0.9995	8.865e-5	0.9995	8.892e-5	0.9995
ETDRK2	$\tau$	6.410e-4	$-$	6.464e-4	$-$	6.422e-4	$-$	6.441e-4	$-$
	$\tau /2$	1.676e-4	1.9352	1.690e-4	1.9350	1.679e-4	1.9352	1.684e-4	1.9351
	$\tau /4$	4.287e-5	1.9672	4.323e-5	1.9671	4.294e-5	1.9672	4.308e-5	1.9671
	$\tau /8$	1.084e-5	1.9835	1.093e-5	1.9835	1.086e-5	1.9835	1.089e-5	1.9835
	$\tau /16$	2.726e-6	1.9917	2.749e-6	1.9917	2.730e-6	1.9917	2.739e-6	1.9917
	$\tau /32$	6.834e-7	1.9959	6.892e-7	1.9959	6.846e-7	1.9959	6.867e-7	1.9959
	$\tau /64$	1.711e-7	1.9981	1.725e-7	1.9981	1.714e-7	1.9981	1.719e-7	1.9981
	$\tau /128$	4.278e-8	1.9997	4.314e-8	1.9997	4.285e-8	1.9997	4.299e-8	1.9997

Table 2. Table 2: Spatial convergence rates in the maximum-norm sense in Example 6.1 .

$h=\frac{\pi }{8}$	$\alpha =1$		$\alpha =3$
	Error	Rate	Error	Rate
$h$	1.554e-4	$-$	1.258e-4	$-$
$h/2$	2.430e-5	2.6774	3.328e-5	1.9182
$h/4$	4.491e-6	2.4358	7.441e-6	2.1608
$h/8$	8.679e-7	2.3714	2.017e-6	1.8833
$h/16$	2.068e-7	2.0694	4.701e-7	2.1011
$h/32$	4.944e-8	2.0643	1.344e-7	1.8066
$h/64$	6.590e-9	2.9075	3.483e-8	1.9479

Table 3. Table 3: Rates of convergence to the local limits in the maximum-norm sense in Example 6.1 .

$\delta =0.2$	$\alpha =1$		$\alpha =3$
	Error	Rate	Error	Rate
$\delta$	1.076e-5	$-$	5.371e-6	$-$
$\delta /2$	2.703e-6	1.9927	1.344e-6	1.9991
$\delta /4$	6.250e-7	2.1124	3.153e-7	2.0912
$\delta /8$	1.580e-7	1.9835	6.373e-8	2.3068

Table 4. Table 4: Theoretical and numerical jumps in steady state solutions in Example 6.3 .

	$\delta =0.2$	$\delta =0.8$	$\delta =1.6$	$\delta =3.2$
Theoretical jumps	0	1.802776	1.952562	1.988247
Numerical jumps	0	1.804496	1.952713	1.988242

Equations

$$u_t-{\epsilon }^2{\mathcal{L}}_{\delta }u+u^3-u=0,$$

$$\text{u(⋅,t) is Ω-periodic},$$

$$u(\mathbf{x},0)=u_0(\mathbf{x}),$$

$$\mathcal{L}_{\delta}u(\bm{x})=\frac{1}{2}\int_{B_{\delta}(\bm{0})}\rho_{\delta}(|\bm{s}|)\big{(}u(\bm{x}+\bm{s})+u(\bm{x}-\bm{s})-2u(\bm{x})\big{)}\,\mathrm{d}\bm{s},\qquad\bm{x}\in\Omega$$

$${\int }_0^{\delta }{r}^{1+d}{\rho }_{\delta }(r)\mathrm{d}r=\frac{2d}{S_d}$$

$${\int }_{B_{\delta }(0)}|\mathbf{s}{|}^2{\rho }_{\delta }(|\mathbf{s}|)\mathrm{d}\mathbf{s}=2d.$$

$$\underset{\mathbf{x}\in \Omega }{\max }|{\mathcal{L}}_{\delta }u(\mathbf{x})-{\mathcal{L}}_{0}u(\mathbf{x})|\le C{\delta }^2\Vert u{\Vert }_{C^4},\forall u\in C_{\text{per}}^4(\overline{\Omega }),$$

$${\phi }_t-{\epsilon }^2{\mathcal{L}}_0\phi +{\phi }^3-\phi =0,$$

$$\text{φ(⋅,t) is Ω-periodic},$$

$$\phi (\mathbf{x},0)={\phi }_0(\mathbf{x}),$$

$$\Vert u(\cdot ,t){\Vert }_{L^{\infty }}\le 1,\forall t>0.$$

$$E(u)=\int_{\Omega}\Big{(}\frac{1}{4}(u^{2}(\bm{x})-1)^{2}-\frac{\varepsilon^{2}}{2}u(\bm{x})\mathcal{L}_{\delta}u(\bm{x})\Big{)}\,\mathrm{d}\bm{x},$$

$$E(u(t_2))\le E(u(t_1)),\forall t_2\ge t_1\ge 0,$$

$${\mathcal{L}}_{\delta }u({\mathbf{x}}_{\mathbf{i}})=\frac{1}{2}{\int }_{B_{\delta }(0)}\frac{u({\mathbf{x}}_{\mathbf{i}}+\mathbf{s})+u({\mathbf{x}}_{\mathbf{i}}-\mathbf{s})-2u({\mathbf{x}}_{\mathbf{i}})}{|\mathbf{s}{|}^2}\Vert \mathbf{s}{\Vert }_1\cdot \frac{|\mathbf{s}{|}^2}{\Vert \mathbf{s}{\Vert }_1}{\rho }_{\delta }(|\mathbf{s}|)\mathrm{d}\mathbf{s},$$

$$\mathcal{L}_{\delta,h}u(\bm{x}_{\bm{i}})=\frac{1}{2}\int_{B_{\delta}(\bm{0})}\mathcal{I}_{h}\bigg{(}\frac{u(\bm{x}_{\bm{i}}+\bm{s})+u(\bm{x}_{\bm{i}}-\bm{s})-2u(\bm{x}_{\bm{i}})}{|\bm{s}|^{2}}\|\bm{s}\|_{1}\bigg{)}\frac{|\bm{s}|^{2}}{\|\bm{s}\|_{1}}\rho_{\delta}(|\bm{s}|)\,\mathrm{d}\bm{s},$$

$${\mathcal{I}}_{h}v(\mathbf{s})=\sum _{{\mathbf{s}}_{\mathbf{j}}}v({\mathbf{s}}_{\mathbf{j}}){\psi }_{\mathbf{j}}(\mathbf{s}),$$

$${\mathcal{L}}_{\delta ,h}u({\mathbf{x}}_{\mathbf{i}})=\sum _{0\ne {\mathbf{s}}_{\mathbf{j}}\in B_{\delta }(0)}\frac{u({\mathbf{x}}_{\mathbf{i}}+{\mathbf{s}}_{\mathbf{j}})+u({\mathbf{x}}_{\mathbf{i}}-{\mathbf{s}}_{\mathbf{j}})-2u({\mathbf{x}}_{\mathbf{i}})}{|{\mathbf{s}}_{\mathbf{j}}{|}^2}\Vert {\mathbf{s}}_{\mathbf{j}}{\Vert }_1{\beta }_{\delta }({\mathbf{s}}_{\mathbf{j}}),{\mathbf{x}}_{\mathbf{i}}\in {\Omega }_h,$$

$${\beta }_{\delta }({\mathbf{s}}_{\mathbf{j}})=\frac{1}{2}{\int }_{B_{\delta }(0)}{\psi }_{\mathbf{j}}(\mathbf{s})\frac{|\mathbf{s}{|}^2}{\Vert \mathbf{s}{\Vert }_1}{\rho }_{\delta }(|\mathbf{s}|)\mathrm{d}\mathbf{s}.$$

$$\underset{{\mathbf{x}}_{\mathbf{i}}\in {\Omega }_h}{\max }|{\mathcal{L}}_{\delta ,h}u({\mathbf{x}}_{\mathbf{i}})-{\mathcal{L}}_{\delta }u({\mathbf{x}}_{\mathbf{i}})|\le C{h}^2\Vert u{\Vert }_{C^4},$$

$$\frac{\mathrm{d}U}{\mathrm{d}t}$$

$$U(0)$$

$$L_h:=-{\epsilon }^2{D}_h+\kappa I_{dN},f(U):=(\kappa +1)U-U^{.3},$$

$$\frac{\mathrm{d}U}{\mathrm{d}t}+L_{h}U=f(U),$$

$$U(t+\tau )={\mathrm{e}}^{-L_h\tau }U(t)+{\int }_0^{\tau }{\mathrm{e}}^{-L_h(\tau -s)}f(U(t+s))\mathrm{d}s,\forall t\ge 0,\tau >0.$$

$${\varphi }^{(j)}({\lambda }_i),0\le j\le n_i-1,1\le i\le m$$

$$U(t_{n+1})={\mathrm{e}}^{-L_h\tau }U(t_n)+{\int }_0^{\tau }{\mathrm{e}}^{-L_h(\tau -s)}f(U(t_n+s))\mathrm{d}s.$$

$$U^{n+1}={\mathrm{e}}^{-L_h\tau }{U}^n+{\int }_0^{\tau }{\mathrm{e}}^{-L_h(\tau -s)}f(U^n)\mathrm{d}s,$$

$$U^{n+1}={\varphi }_0(L_h\tau )U^n+\tau {\varphi }_1(L_h\tau )f(U^n),$$

$${\varphi }_0(a):={\mathrm{e}}^{-a},{\varphi }_1(a):=\frac{1-{\mathrm{e}}^{-a}}{a},a\ne 0.$$

$${\tilde{U}}^{n+1}$$

$$U^{n+1}$$

$${\tilde{U}}^{n+1}$$

$$U^{n+1}$$

$${\varphi }_2(a):=\frac{{\mathrm{e}}^{-a}-1+a}{a^2},a\ne 0.$$

Full text

\slugger

mmsxxxxxxxx–x

1 Introduction

In this paper, we consider numerical solutions of the initial-boundary-value problem of the nonlocal Allen-Cahn (NAC) equation as follows:

[TABLE]

where $u(\bm{x},t)$ denotes the unknown function, $\Omega=(0,X)^{d}$ is a hypercube domain in $\mathbb{R}^{d}$, $\varepsilon>0$ is an interfacial parameter, and $\mathcal{L}_{\delta}$ is a nonlocal operator, parameterized by the positive horizon parameter $\delta>0$ measuring the range of nonlocal interactions. Assume that $\mathcal{L}_{\delta}$ is defined by

[TABLE]

with $B_{\delta}(\bm{0})$ denoting the ball in $\mathbb{R}^{d}$ centered at the origin with the radius $\delta$ and $\rho_{\delta}:[0,\delta]\to\mathbb{R}$ being a nonnegative kernel function. To enforce the consistency, as $\delta\to 0$, of the nonlocal operator $\mathcal{L}_{\delta}$ with the standard Laplacian operator $\mathcal{L}_{0}:=\Delta$, we further assume the kernel $\rho_{\delta}$ satisfies

[TABLE]

with $S_{d}$ being the area of the unit sphere in $\mathbb{R}^{d}$, or equivalently,

[TABLE]

Note that (3) also means that the kernel $\rho_{\delta}$ has a finite second order moment. The continuum property of the nonlocal operator $\mathcal{L}_{\delta}$ gives [10, 11]

[TABLE]

where $C>0$ is a constant independent of $\delta$. The local limit of the NAC problem (1) is exactly the classic (local) Allen-Cahn (LAC) equation taking the following form:

[TABLE]

The LAC equation (5) is a well-known phase field model used to describe the motion of anti-phase boundaries in crystalline solids [1].

In recent years, the nonlocal models involving the nonlocal operator (2), such as the NAC equation (1), have appeared in a variety of applications ranging from physics, materials science to finance and image processing, for instance, phase transition [4, 19], peridynamics continuum theory [35, 36], image analyses [20, 21], and nonlocal heat conduction [7]. Rigorous mathematical analysis of nonlocal models can be found in the literatures, e.g., [3, 4, 16], and a more systematic mathematical framework of nonlocal problems was developed in [11, 12] in parallel to the analysis for the classic partial differential equations. Since the exact/analytic solutions of these nonlocal models are usually not available, numerical methods play an important role in studying these models. Bates et al. [5] considered a finite difference discretization of the NAC equation with an integrable kernel and developed an $L^{\infty}$ stable and convergent numerical scheme by treating the nonlinear and nonlocal terms explicitly. A similar technique was applied on the NAC-type problem coupled with a heat equation and an $L^{\infty}$ stable and convergent numerical scheme was obtained [2]. For the nonlocal diffusion models with more general kernels and variable boundary conditions, finite difference and finite element approximations were addressed in [14, 38, 40, 48]. To illustrate the limiting behaviors of the numerical solution of the nonlocal model to the exact solution of the corresponding local counterpart, Tian and Du proposed in [41] the concept of asymptotic compatibility, and the spectral-Galerkin approximation of the NAC equation was then proved to be asymptotically compatible in [15]. The convergence of asymptotically compatible schemes is insensitive to the choices of modeling and discretization parameters so that such schemes provide robust numerical approximations of nonlocal models.

As a nonlocal analogue of the LAC equation (5), the NAC equation (1) possesses some similar properties. First, it can be shown that the NAC equation (1) satisfies a maximum principle: if the initial value and the boundary conditions are bounded by $1$, then the entire solution is also bounded by $1$, i.e.,

[TABLE]

Second, as a phase field type model, the NAC equation (1a) can be viewed as an $L^{2}$ gradient flow with respect to the energy functional

[TABLE]

and thus, the solution $u$ to (1) decreases the energy (6) in time, that is,

[TABLE]

which is often called the energy dissipation law. Such two properties are important in the study of the stability of the solution to (1), and whether they could be inherited in the discrete level is a significant issue in numerical simulations. A major objective of this work is develop maximum principle preserving and energy stable numerical schemes for approximating the NAC equation (1).

Energy stability has been widely investigated for numerical schemes of classic PDE-based phase field models, such as convex splitting schemes [32, 43], stabilized schemes [34, 42], invariant energy quadratization methods [45, 46] and so on. It is interesting to study whether similar analysis can be applied to nonlocal phase field models due to the lack of the high-order diffusion term. Guan et al. [22] constructed a convex splitting scheme for the nonlocal Cahn-Hilliard equation by treating the nonlinear term implicitly and setting the nonlocal term into the explicit part. Their scheme allows one to evaluate the nonlocal term explicitly only once at each time step, but the nonlinear iterations are still inevitable. In order to avoid the nonlinear iterations, a stabilized scheme was the linear stabilization strategy was adopted in [13] to develop the stabilized linear schemes which can be solved efficiently by using the fast Fourier transform. The energy stability of the fully discrete schemes were only shown under the assumption that the stabilizer depends implicitly on the uniform bound of the numerical solution.

The numerical method we will adopt in this work is the so-called exponential time differencing (ETD), which involves exact integration of the governing equations followed by an explicit approximation of a temporal integral involving the nonlinear terms. The ETD schemes was systematically studied in [6] and then further developed by Cox and Matthews with the applications on stiff systems [9], where higher-order multistep and Runge-Kutta versions of these schemes were described. Hochbruck and Ostermann provided several nice reviews on the ETD Runge-Kutta methods [24] and the ETD multistep methods [25]. In addition, the convergence of these methods were analyzed in detail under the analytical framework therein. The linear stabilities of some ETD and modified ETD schemes were investigated by Du and Zhu [17, 18].

A distinctive feature of ETD schemes is the exact evaluation of the contribution of the linear part, which provides satisfactory stability and accuracy even though the linear terms have strong stiffness. Such an advantage leads to some successful applications of ETD schemes on phase field models which usually yield highly stiff ODE systems after suitable spatial discretizations. Ju et al. developed stable and compact ETD schemes and their fast implementations for Allen-Cahn [29, 49], Cahn-Hilliard [28], and elastic bending energy models [44] by utilizing suitable linear splitting techniques. All the proposed ETD schemes are explicit and thus highly efficient for practical implementations. A localized compact ETD algorithm based on the overlapping domain decomposition was firstly used in [47] for extreme-scale phase field simulations of three-dimensional coarsening dynamics in the supercomputer, and the results showed excellent parallel scalability of the method. In [27], the ETD multistep method was applied on the epitaxial growth model without slope selection [31], and the energy stability and the error estimates were established rigorously, which is the first work to analyze the energy stability and convergence of the ETD schemes for phase field models in the theoretical level. To complete the theoretical analysis, there is no need for any assumptions on the numerical solutions due to the specific property of the logarithm term in the no-slope-selection model. However, for other phase field models, such as the Cahn-Hilliard equation, the assumptions on the uniform boundedness of the numerical solutions or the Lipschitz continuity of some nonlinear functions are inevitable to ensure the energy stability. Therefore, for the models whose solutions satisfy the maximum principle essentially, it is highly desired to develop numerical approximations preserving the maximum principle in the discrete sense.

One of the typical phase field models satisfying the maximum principle is the LAC equation (5). Recently, there have been some investigations on the maximum principle preserving numerical schemes for (5). Tang and Yang [37] proved that the first order implicit-explicit schemes, with or without the stabilizing term, preserve the maximum principle under some condition on the time step size. Then, the energy stability and the maximum-norm error estimates are obtained by using the discrete maximum principle. Shen et al. [33] generalized the results presented in [37] to the case of the Allen-Cahn-like equation in a more abstract form with the potential and mobility satisfying some certain conditions. Hou et al. [26] studied the numerical approximation of the fractional Allen-Cahn equation by considering the conventional Crank-Nicolson scheme. They proved that the Crank-Nicolson scheme preserves the maximum principle and this is the first work on the second order schemes preserving the maximum principle. More than ten years ago, Du and Zhu [18] showed that the first order ETD scheme in the space-continuous version for (5) satisfies the maximum principle, where some properties of the heat kernel were used in their proof. However, the fully discrete ETD schemes were never studied.

The organization of this paper is as follows. In Section 2, we construct the first and second order ETD time-stepping schemes for the NAC equation with the quadrature-based finite difference approximation being used for spatial discretization. Efficient implementation issues of the schemes are also briefly discussed. In Section 3, both schemes are shown to satisfy the discrete maximum principle unconditionally. Error estimates and asymptotic compatibility of the schemes are obtained in Section 4 and the discrete energy stability proved in Section 5. Various numerical experiments are carried out in Section 6 to verify the theoretical results and to investigate the effects of the nonlocal parameters. Finally, some concluding remarks are given in Section 7.

2 Fully discrete exponential time differencing schemes

In this section, we present the fully discrete ETD schemes for the NAC equation in general dimensions, where the finite difference method, based on the contribution made in [14], is adopted for the spatial discretization of the nonlocal diffusion operator. In particular, we also give the specific expression of the discrete nonlocal operator in 2D later.

2.1 Quadrature-based finite difference semi-discretization

Given a positive integer $N$, we set $h=X/N$ as the uniform square mesh size and define $\bm{x}_{\bm{i}}=h{\bm{i}}$ as the nodes in the mesh, where ${\bm{i}}\in\mathbb{Z}^{d}$ denotes a multi-index. Let $\Omega_{h}$ be the set of nodes in the domain $\overline{\Omega}$. At any node $\bm{x}_{\bm{i}}$, the nonlocal operator (2) can be rewritten as

[TABLE]

where $\|\cdot\|_{1}$ stands for the vector $1$-norm. Then, a quadrature-based finite difference discretization of the nonlocal operator (8) can be defined as [14]

[TABLE]

where $\mathcal{I}_{h}$ represents the piecewise $d$-multilinear interpolation operator with respect to $\bm{s}$ associated with the mesh. More precisely, for a function $v(\bm{s})$, the interpolation $\mathcal{I}_{h}v(\bm{s})$ is piecewise linear with respect to each component of the spatial variable $\bm{s}$ and

[TABLE]

where $\psi_{\bm{j}}$ is the piecewise $d$-multilinear basis function satisfying $\psi_{\bm{j}}(\bm{s}_{\bm{i}})=0$ when ${\bm{i}}\not={\bm{j}}$ and $\psi_{\bm{j}}(\bm{s}_{\bm{j}})=1$. Therefore, the resulting quadrature-based finite difference discretization of the nonlocal operator (8) reads

[TABLE]

where the periodicity conditions are used for the nodes not in $\Omega_{h}$, and

[TABLE]

It is easy to check that the operator $\mathcal{L}_{\delta,h}$ is self-adjoint and negative semi-definite.

The discretized scheme (9) is proposed in [14] for the problem with a homogeneous Dirichlet-type nonlocal constraint and it has been proved [14, 39] that, for any fixed $\delta>0$, the discrete operators $\mathcal{L}_{\delta,h}$ is consistent to $\mathcal{L}_{\delta}$ with the errors $\mathcal{O}(h^{2})$ as $h\to 0$. For the case of periodic boundary condition considered here, all similar estimates also hold, so we give the following consistency estimates without proof.

Lemma 1**.**

Assume that $u\in C_{\text{\rm per}}^{4}(\overline{\Omega})$, then it holds that

[TABLE]

where $C>0$ is a constant independent of $\delta$ and $h$.

By ordering the nodes in the lexicographical order, we can obtain the nonlocal stiffness matrix, denoted by $D_{h}\in\mathbb{R}^{dN\times dN}$, associated with $\mathcal{L}_{\delta,h}$. It is obvious that $D_{h}$ is symmetric, negative semi-definite, and weakly diagonally dominant with all negative diagonal entries. The space-discrete scheme of (1) is to find a vector-valued function $U:[0,T]\to\mathbb{R}^{dN}$ such that

[TABLE]

where $U^{.3}=(U_{1}^{3},U_{2}^{3},\dots,U_{dN}^{3})^{T}$ and $U_{0}\in\mathbb{R}^{dN}$ is given by the initial data. For the sake of the stability of the time-stepping schemes developed later, we introduce a stabilizing parameter $\kappa>0$ and define

[TABLE]

where $I_{dN}$ is the $dN\times dN$ identity matrix, so $L_{h}$ is symmetric, positive definite, and strictly diagonally dominant with all positive diagonal entries. Then, the ODE system (12a) could be written as

[TABLE]

whose solution satisfies

[TABLE]

In the above we have used a property of the differentiation of matrix exponentials (see Lemma 2 (5)), and we list below some other properties of matrix functions (see [23]) useful to the analysis later.

Lemma 2** (see [23]).**

Let $\phi$ be defined on the spectrum of $A\in\mathbb{C}^{m\times m}$, that is, the values

[TABLE]

exist, where $\{\lambda_{i}\}_{i=1}^{m}$ are the eigenvalues of $A$, and $n_{i}$ is the order of the largest Jordan block where $\lambda_{i}$ appears. Then

(1)* $\phi(A)$ commutes with $A$;*

(2)* $\phi(A^{T})=\phi(A)^{T}$;*

(3)* the eigenvalues of $\phi(A)$ are $\{\phi(\lambda_{i}):1\leq i\leq m\}$;*

(4)* $\phi(P^{-1}AP)=P^{-1}\phi(A)P$ for any nonsingular matrix $P\in\mathbb{C}^{m\times m}$;*

(5)* $\dfrac{\mathrm{d}}{\mathrm{d}s}(\mathrm{e}^{As})=A\mathrm{e}^{As}=\mathrm{e}^{As}A$ for any $s\in\mathbb{R}$.*

2.2 Exponential time differencing schemes for time-stepping

Given a positive integer $K_{t}$, we divide the time interval by $\{t_{n}=n\tau:0\leq n\leq K_{t}\}$ with a uniform time step $\tau=T/K_{t}$. Setting $t=t_{n}$ in (14) gives us

[TABLE]

The first order ETD (ETD1) scheme comes from approximating $f(U(t_{n}+s))$ by $f(U(t_{n}))$ in $s\in[0,\tau]$ and calculating the produced integral exactly [9]. The ETD1 scheme of (1) reads: for $n=0,1,\cdots,K_{t}-1$,

[TABLE]

that is,

[TABLE]

where

[TABLE]

The second order ETD Runge-Kutta (ETDRK2) scheme is obtained by approximating $f(U(t_{n}+s))$ by a linear interpolation based on $f(U(t_{n}))$ and $f(\widetilde{U}^{n+1})$, where $\widetilde{U}^{n+1}$ is an approximation of $U(t_{n+1})$. The ETDRK2 scheme of (1) takes the form: for $n=0,1,\cdots,K_{t}-1$,

[TABLE]

or equivalently,

[TABLE]

where

[TABLE]

We know that $\phi_{0}(a)$, $\phi_{1}(a)$, and $\phi_{2}(a)$ are all positive when $a>0$.

2.3 Efficient implementations of the ETD schemes

We close this section by giving a brief illustration on the practical implementation of the proposed schemes (17) and (19). Using the 2D case as the example, we first give the explicit formula of the discrete operator $\mathcal{L}_{\delta,h}$ (as illustrated in [14]) and then discuss efficient implementation of the actions of the matrix exponentials.

Let $u_{i,j}$ be the nodal value of the numerical solution at the mesh point $(x_{i},y_{j})\in\Omega_{h}\subset\mathbb{R}^{2}$ and $r=[\delta/h]+1$ be the smallest integer larger than $\delta/h$. Then, we have

[TABLE]

where $c_{0,0}=0$ and

[TABLE]

with $\psi_{p,q}$ denoting the bilinear basis function located at the point $(ph,qh)$ and $B_{\delta}^{+}$ the first quadrant of the disc centered at the origin with radius $\delta$. Note that $c_{p,q}=c_{q,p}$ for any $p$ and $q$. One can apply efficient quadrature rules on the double integrals in (21).

We represent $U^{n}\in\mathbb{R}^{N\times N}$ in the matrix form with entries $u_{i,j}^{n}$ and define the operator $\mathcal{L}_{h}=\kappa\mathcal{I}-\varepsilon^{2}\mathcal{L}_{\delta,h}$ whose matrix form is given by $L_{h}$ defined in (13), where $\mathcal{I}$ is the identity mapping. The key process of calculating $U^{n+1}$ from the scheme (17) or (19) is the efficient implementation of the actions of the operator exponentials $\phi_{\gamma}(\mathcal{L}_{h}\tau)$, $\gamma=0,1,2$. Since $\mathcal{L}_{h}$ comes from the discretization of $\mathcal{L}_{\delta}$ with the periodic boundary condition, the exponentials $\phi_{\gamma}(\mathcal{L}_{h}\tau)$ can be implemented by the 2D discrete Fourier transform (DFT). More precisely, if we denote by $\mathcal{F}$ the 2D DFT operator, then, for any $V=(V_{k,l})\in\mathbb{C}^{N\times N}$, the action of the operator $\widehat{\mathcal{L}}_{h}:=\mathcal{F}\mathcal{L}_{h}\mathcal{F}^{-1}$ can be implemented via

[TABLE]

where $\lambda_{k,l}$’s, the eigenvalues of $\mathcal{L}_{h}$, are given by

[TABLE]

According to Lemma 2 (4), we have

[TABLE]

The actions of $\mathcal{F}$ and $\mathcal{F}^{-1}$ can be implemented by the 2D fast Fourier transform (FFT) and its inverse transform, respectively. Such implementation can be naturally generalized to higher-dimensional spaces and the computational complexity is thus $\mathcal{O}(N^{d}\log N)$ per time step.

3 Discrete maximum principle

Denote by $\|\cdot\|_{\infty}$ the standard vector or matrix $\infty$-norm, and by $\|\cdot\|_{2}$ the standard vector or matrix $2$-norm. The following lemma is a special case of Theorem 2 in [30].

Lemma 3**.**

Let $A=(a_{ij})\in\mathbb{R}^{m\times m}$ with $a_{ii}<0$, $1\leq i\leq m$, and there exists $\kappa>0$ such that

[TABLE]

then the nontrivial solution $\theta=\theta(t)$ to the linear differential system

[TABLE]

satisfies

[TABLE]

The following result is a key ingredient to prove the discrete maximum principle.

Lemma 4**.**

For any $\kappa>0$ and $\tau>0$, we always have $\|\mathrm{e}^{-L_{h}\tau}\|_{\infty}\leq\mathrm{e}^{-\kappa\tau}$.

Proof.

Since $L_{h}$ is strictly diagonally dominant with all positive diagonal entries, the matrix $A:=-L_{h}$ satisfies the conditions of Lemma 3 with $\kappa>0$. For any nonzero $\theta_{0}\in\mathbb{R}^{dN}$, we know that the solution to the linear differential system (22) with the initial value $\theta(0)=\theta_{0}$ is given by $\theta(t)=\mathrm{e}^{-L_{h}t}\theta_{0}$ and, by using Lemma 3, satisfies

[TABLE]

Therefore, the result follows from the arbitrariness of $\theta_{0}$. ∎

Remark 3.1**.**

Although our deduction above is restricted to the case of periodic boundary condition, the result of Lemma 4 is also suitable for the case of the Dirichlet boundary condition, since the corresponding nonlocal stiffness matrix $D_{h}$ is still weakly diagonally dominant with all negative diagonal entries. The analysis results in this paper could be obtained similarly for the Dirichlet boundary condition.

Since the nonlinear mapping $f:\mathbb{R}^{dN}\to\mathbb{R}^{dN}$ defined in (13) is actually a set of $dN$ independent one-variable functions, we just need to consider anyone of them.

Lemma 5**.**

Define $f_{0}(\xi)=(\kappa+1)\xi-\xi^{3}$ for any $\xi\in\mathbb{R}$. If $\kappa\geq 2$, then

[TABLE]

Proof.

Obviously, $f_{0}(-1)=-\kappa$ and $f_{0}(1)=-\kappa$. For any $\xi\in[-1,1]$, if $\kappa\geq 2$, we have

[TABLE]

which gives us the result. ∎

Theorem 6**.**

Assume that the initial data satisfies $\|u_{0}\|_{L^{\infty}}\leq 1$. Then, for any time step size $\tau>0$, the ETD1 scheme (16) preserves the discrete maximum principle, i.e.,

[TABLE]

provided the stabilizing parameter $\kappa\geq 2$.

Proof.

We prove this theorem by induction. Obviously, it holds $\|U^{0}\|_{\infty}\leq\|u_{0}\|_{L^{\infty}}\leq 1$. Now assume that the result holds for $n=k$, i.e., $\|U^{k}\|_{\infty}\leq 1$. Next we check this holds for $n=k+1$. According to the scheme (16), we have

[TABLE]

It follows from Lemma 4 that

[TABLE]

Using Lemma 5 and $\|U^{k}\|_{\infty}\leq 1$, we have $\|f(U^{k})\|_{\infty}\leq\kappa$. Consequently,

[TABLE]

which completes the proof. ∎

Theorem 7**.**

Assume that the initial data satisfies $\|u_{0}\|_{L^{\infty}}\leq 1$. Then, for any time step size $\tau>0$, the ETDRK2 scheme (18) preserves the discrete maximum principle, i.e.,

[TABLE]

provided the stabilizing parameter $\kappa\geq 2$.

Proof.

We again prove this by induction. Obviously, it holds $\|U^{0}\|_{\infty}\leq\|u_{0}\|_{L^{\infty}}\leq 1$. Now assume that the result holds for $n=k$, i.e., $\|U^{k}\|_{\infty}\leq 1$. Next we check this for $n=k+1$. According to the formula (18a) and the proof of Theorem 6, we have $\|\widetilde{U}^{k+1}\|_{\infty}\leq 1$. According to the formula (18b), we have

[TABLE]

Since $\|U^{k}\|_{\infty}\leq 1$ and $\|\widetilde{U}^{k+1}\|_{\infty}\leq 1$, using Lemma 5, we have

[TABLE]

and then, for $s\in[0,\tau]$,

[TABLE]

Again, by using (24), we obtain from (25) that

[TABLE]

which completes the proof. ∎

Remark 3.2**.**

We will then always require $\kappa\geq 2$ for the proposed ETD1 and ETDRK2 schemes in the rest of the paper so that they preserve the discrete maximum principle.

4 Error estimates and asymptotic compatibility

We will analyze the two types of convergence behaviors of the numerical solution to the ETD1 scheme (16) and the ETDRK2 scheme (18), respectively. First, for any fixed $\delta>0$, we prove that the numerical solution converges to the exact solution of the NAC equation (1) as the spatial mesh size $h$ and the time step size $\tau$ go to zero. Second, we show that the numerical solution converges to the exact solution of the LAC equation (5) as the horizon parameter $\delta$, the spatial size $h$, and the temporal step $\tau$ approach to zero. The latter convergence behavior of the numerical solution is often called the asymptotic compatibility [41] in nonlocal modeling.

We first establish the $L^{\infty}$ error estimates for the numerical solution produced by the ETD1 scheme (16) for the NAC equation (1) with any fixed $\delta>0$.

Theorem 8**.**

Given a fixed $\delta>0$. Assume that the exact solution $u$ to the NAC equation (1) belongs to $C^{1}([0,T];C_{\text{\rm per}}^{4}(\overline{\Omega}))$ and $\{U^{n}\}_{n=0}^{K_{t}}$ is generated by the ETD1 scheme (16) with $U^{0}=I^{h}u_{0}$. If $\|u_{0}\|_{L^{\infty}}\leq 1$, then we have

[TABLE]

for any $h>0$ and $\tau>0$, where the constant $C>0$ depends on the $C^{1}([0,T];C_{\text{\rm per}}^{4}(\overline{\Omega}))$ norm of $u$, but independent of $\delta$, $h$ and $\tau$.

Proof.

Recalling the construction of the ETD1 scheme (16), we observe that, for a known $U^{n}$, the solution $U^{n+1}$ is actually given by $U^{n+1}=W_{1}(\tau)$ with the function $W_{1}:[0,\tau]\to\mathbb{R}^{dN}$ determined by the following evolution equation

[TABLE]

Then, for the NAC equation (1), we can give a similar illustration as follows: for given $u(\bm{x},t_{n})$, the solution $u(\bm{x},t_{n+1})$ is determined by $u(\bm{x},t_{n+1})=w(\bm{x},\tau)$ with the function $w(\bm{x},s)$ satisfying

[TABLE]

Let $e_{1}(s)=W_{1}(s)-I^{h}w(s)$, where $I^{h}$ is the operator limiting a function on the mesh $\Omega_{h}$. Then, the difference between (27) and (28) yields

[TABLE]

where $R_{h\tau}^{(1)}(s)$ is the truncated error, that is,

[TABLE]

Since $|f^{\prime}(\xi)|\leq\kappa+1$ for any $\xi\in[-1,1]$ and both the exact and numerical solutions satisfy the maximum principles if $\kappa\geq 2$, we have

[TABLE]

According to the consistency result given by Lemma 1, if the exact solution $u$ is sufficiently smooth, at least $u\in C^{1}([0,T];C_{\text{\rm per}}^{4}(\overline{\Omega}))$, then we have

[TABLE]

and for $s\in[0,\tau]$,

[TABLE]

where $C_{1}$ and $C_{2}$ depend on the $C^{1}([0,T];C_{\text{\rm per}}^{4}(\overline{\Omega}))$ norm of $u$, but independent of $\delta$, $h$ and $\tau$. Thus, we obtain

[TABLE]

where $C=\max\{C_{1}\varepsilon^{2},C_{2}(S+1)\}$. Integrating the ODE in (29) leads to

[TABLE]

Setting $t=\tau$ and using (24), we have

[TABLE]

where in the last step we have used the fact that $1-\mathrm{e}^{-s}\leq s$ for any $s>0$. Using the Gronwall’s inequality, we obtain

[TABLE]

Therefore, we obtain (26) since $e_{1}^{0}=0$ and $n\tau=t_{n}$. ∎

Now, we turn to the $L^{\infty}$ error estimates for the ETDRK2 scheme (18) with any fixed $\delta>0$.

Theorem 9**.**

Given a fixed $\delta>0$. Assume that the exact solution $u$ of the NAC equation (1) belongs to $C^{2}([0,T];C_{\text{\rm per}}^{4}(\overline{\Omega}))$ and $\{U^{n}\}_{n=0}^{K_{t}}$ is generated by the ETDRK2 scheme (18) with $U^{0}=I^{h}u_{0}$. If $\|u_{0}\|_{L^{\infty}}\leq 1$, then we have

[TABLE]

for any $h>0$ and $1\geq\tau>0$, where the constant $C>0$ depends on the $C^{2}([0,T];C_{\text{\rm per}}^{4}(\overline{\Omega}))$ norm of $u$, but independent of $\delta$, $h$ and $\tau$.

Proof.

For a known $U^{n}$, the solution $U^{n+1}$ to the ETDRK2 scheme (18) is actually given by $U^{n+1}=W_{2}(\tau)$ with the function $W_{2}:[0,\tau]\to\mathbb{R}^{dN}$ determined by the evolution equation

[TABLE]

where $\widetilde{U}^{n+1}$, defined by (18a), is given by $\widetilde{U}^{n+1}=W_{1}(\tau)$ with $W_{1}(s)$ satisfying (27). Let $e_{2}(s)=W_{2}(s)-I^{h}w(s)$. The difference between (38) and (28) leads to

[TABLE]

where $R_{h\tau}^{(2)}(s)$ is the truncated error given by

[TABLE]

According to error estimates for the linear interpolation, we have, for $s\in[0,\tau]$, that

[TABLE]

where $C_{3}$ depends on the $C^{2}([0,T];C_{\text{\rm per}}^{4}(\overline{\Omega}))$ norm of $u$, but independent of $\delta$, $h$ and $\tau$. Thus, combining with (31), we obtain

[TABLE]

where $C_{4}=\max\{C_{1}\varepsilon^{2},C_{3}\}$. According to (35) in the proof for the ETD1 scheme, we have

[TABLE]

where $C_{5}$ depends on the $C^{1}([0,T];C_{\text{\rm per}}^{4}(\overline{\Omega}))$ norm of $u$. Then, using the Lipschitz continuity of $f$, we obtain

[TABLE]

Combining with (30), we have, for any $s\in(0,\tau)$, that

[TABLE]

Integrating the ODE in (39) leads to

[TABLE]

Setting $t=\tau$ and using (24), we have

[TABLE]

where $C=C_{4}+\frac{1}{2}C_{5}(\kappa+1)$. The condition $\tau\leq 1$ is used in the last two steps of the above derivation. Using the Gronwall’s inequality, we obtain

[TABLE]

which then gives us (37). ∎

Now, let us investigate the asymptotic compatibility of both ETD1 and ETDRK2 schemes. Combining (4) with the uniform estimates (11) of the consistency of $\mathcal{L}_{\delta,h}$, we obtain

[TABLE]

where $C>0$ is a constant independent of $\delta$ and $h$. Then, we can further obtain the asymptotic compatibility of the solution to the ETD1 scheme (16). Denote by $\varphi(\bm{x},t)$ the solution of the LAC equation (5). For given $\varphi(\bm{x},t_{n})$, the solution $\varphi(\bm{x},t_{n+1})$ is determined by $\varphi(\bm{x},t_{n+1})=v(\bm{x},\tau)$ with the function $v(\bm{x},t)$ satisfying

[TABLE]

Let $\widehat{e}(t)=W_{1}(t)-I^{h}v(t)$, where $W_{1}(t)$ is defined by (27). Then, the difference between (27) and (45) yields

[TABLE]

where the remainder $\widehat{R}^{\delta}_{h\tau}(t)$ is given by

[TABLE]

and, according to the estimates (44) and the Lipschitz continuity of $f$ under the condition $\kappa\geq 2$, is bounded by

[TABLE]

where $C>0$ depends on $\varepsilon$, $\kappa$ and the $C^{1}([0,T];C_{\text{\rm per}}^{4}(\overline{\Omega}))$ norm of $\varphi$, but independent of $\delta$, $h$ and $\tau$. By conducting similar analysis as done for error estimates, we can obtain the asymptotic compatibility of the numerical solution to the ETD1 scheme (16). So does the case for the ETDRK2 scheme (18). Therefore, we have the following results.

Theorem 10** (Asymptotic compatibility).**

Assume that the solution $\varphi$ of the local Allen-Cahn equation (5) belongs to $C^{1}([0,T];C_{\text{\rm per}}^{4}(\overline{\Omega}))$ (resp. $C^{2}([0,T];C_{\text{\rm per}}^{4}(\overline{\Omega}))$) and $\{U^{n}\}_{n=0}^{K_{t}}$ is generated by the ETD1 scheme (16) (resp. the ETDRK2 scheme (18)) with $U^{0}=I^{h}\varphi_{0}$. If $\|\varphi_{0}\|_{L^{\infty}}\leq 1$, then we have

[TABLE]

for any $\delta>0$, $h>0$ and $\tau>0$ (resp. $\tau\in(0,1]$), where the constant $C$ depends on the $C^{1}([0,T];C_{\text{\rm per}}^{4}(\overline{\Omega}))$ norm (resp. $C^{2}([0,T];C_{\text{\rm per}}^{4}(\overline{\Omega}))$ norm) of $\varphi$, but independent of $\delta$, $h$ and $\tau$.

5 Discrete energy stability

We first show that the ETD1 scheme (16) inherits the energy decay law (7) in the discrete sense, with respect to the discretized energy $E_{h}$ defined by

[TABLE]

Theorem 11**.**

The approximating solution $\{U^{n}\}_{n=0}^{K_{t}}$ generated by the ETD1 scheme (16) satisfies the energy inequality

[TABLE]

for any $\tau>0$, i.e., the ETD1 scheme is unconditionally energy stable.

Proof.

The difference between the discrete energies at two consecutive time levels yields

[TABLE]

It is easy to verify that

[TABLE]

Since $\kappa\geq 2$, it follows from Theorem 6 that $\|U^{n}\|_{\infty}\leq 1$ and $\|U^{n+1}\|_{\infty}\leq 1$, then we have

[TABLE]

On the other hand, direct calculations lead to

[TABLE]

due to the negative semi-definiteness of the matrix $D_{h}$. Thus, we obtain from (49) that

[TABLE]

We solve $f(U^{n})$ from (17), together with Lemma 2 (1), to get

[TABLE]

and then

[TABLE]

where $B_{1}:=L_{h}-(I-\mathrm{e}^{-L_{h}\tau})^{-1}L_{h}$. Define a function

[TABLE]

then $g_{1}(a)<0$ for any $a>0$ and $B_{1}=\tau^{-1}g_{1}(L_{h}\tau)$. Since $L_{h}$ is symmetric and positive definite, by using Lemma 2 (2) and (3), we know that $B_{1}$ is symmetric and negative definite. Therefore, we obtain

[TABLE]

which completes the proof. ∎

For the ETDRK2 scheme (18), we can prove the uniform boundedness of the discretized energy $E_{h}$.

Theorem 12**.**

Under the assumptions of Theorem 9, the approximate solution $\{U^{n}\}_{n=0}^{K_{t}}$ generated by the ETDRK2 scheme (18) satisfies

[TABLE]

for any $h>0$ and $1\geq\tau>0$, where the constant $\widetilde{C}$ is independent of $h$ and $\tau$. Furthermore, if $h\leq 1$ and $\tau=\lambda\sqrt{h}$ for some constant $\lambda>0$, we have

[TABLE]

where the constant $\widehat{C}$ is independent of $h$ and $\tau$, i.e., the discrete energy is uniformly bounded.

Proof.

The first step is to calculate the increment $E_{h}(U^{n+1})-E_{h}(U^{n})$ directly, which is completely identical to the proof for the ETD1 scheme, and we obtain

[TABLE]

Using (19a) and (19b), we then get

[TABLE]

Acting $(\tau\phi_{1}(L_{h}\tau))^{-1}$ on both sides of (51) gives us

[TABLE]

and then, using the notation $B_{1}$ defined in the proof for the ETD1 scheme, we obtain

[TABLE]

where $B_{2}=g_{2}(L_{h}\tau)$ with

[TABLE]

Since $0<g_{2}(a)<1$ for any $a>0$, we know that $B_{2}$ is symmetric and positive definite and $\|B_{2}\|_{2}<1$. Using the mean-value theorem, we have

[TABLE]

where $G^{n}$ is a diagonal matrix with, according to (23), the diagonal entries between [math] and $\kappa+1$, which implies that $\|G^{n}\|_{\infty}\leq\kappa+1$. Then, by the negative definiteness of $B_{1}$, we obtain

[TABLE]

According to Theorem 9, we can derive

[TABLE]

for any $h>0$ and $\tau\in(0,1]$, where the constants $C_{6}$ and $C_{7}$ depend on the $C^{2}([0,T];C_{\text{\rm per}}^{4}(\overline{\Omega}))$ norm of $u$ and $C_{8}=2C_{6}\mathrm{e}^{T}+C_{7}$. Similarly, using Theorems 8, we can obtain

[TABLE]

with the constant $C_{9}>0$ depending on $T$ and the $C^{1}([0,T];C_{\text{\rm per}}^{4}(\overline{\Omega}))$ norm of $u$. In addition, we have

[TABLE]

Therefore, we obtain

[TABLE]

where the constant $\widetilde{C}=C_{8}C_{9}\sqrt{dX}(\kappa+1)$ is independent of $h$ and $\tau$. By induction, we further obtain

[TABLE]

When $h\leq 1$ and $\tau=\lambda\sqrt{h}$, it holds $h^{-\frac{1}{2}}\tau^{-1}(h^{2}+\tau)^{2}\leq(\lambda+1)^{2}/\lambda$. The proof is then completed by setting $\widehat{C}=\widetilde{C}T(\lambda+1)^{2}/\lambda$. ∎

6 Numerical experiments

In this section, we will carry out some numerical experiments in the 2D space to demonstrate the effectiveness and efficiency of the ETD schemes (16) and (18) for solving the NAC equation (1). The fractional power kernels

[TABLE]

is chosen and satisfies the finite second order moment condition (3). When $\alpha\in[0,2)$, the kernel satisfies $\rho_{\delta}(|\bm{s}|)\in L^{1}(\mathbb{R}^{2})$, which means that the nonlocal diffusion operator $\mathcal{L}_{\delta}$ is a bounded linear operator in this case. The kernel is non-integrable when $\alpha\in[2,4)$. We first verify the temporal and spatial convergence rates of the fully discrete schemes with a smooth initial data, and then check the discrete maximum principle and energy stability of the evolutions beginning with a random initial state. Next, we present a further numerical investigation on the steady state solutions to the model with integrable kernels. The ETDRK2 scheme (18) is adopted in all the simulations while the ETD1 scheme (16) is only considered for the temporal convergence tests due to the lack of high accuracy. The domain $\Omega=(0,2\pi)\times(0,2\pi)$ is used in all examples. We also take the stabilizing parameter $\kappa=2$ for the numerical schemes in all experiments.

6.1 Convergence tests

Example 6.1**.**

We consider the NAC equation (1) with a smooth initial data $u_{0}(x,y)=0.5\sin x\sin y$. We set the interfacial parameter $\varepsilon=0.1$ and the terminal time $T=0.5$. For the kernel (52), $\alpha=1$ (integrable) and $\alpha=3$ (non-integrable) are adopted, respectively.

First, by setting $N=256$, we tested the convergence in time for the cases $\delta=0.2$ and $\delta=2$. We calculated the numerical solutions of the NAC equation using the ETD1 scheme (16) and the ETDRK2 scheme (18) with various time step sizes $\tau=0.05\times 2^{-k}$ with $k=0,1,\dots,7$. To compute the errors, we treated the solution obtained by the ETDRK2 scheme with $\tau=10^{-6}$ as the benchmark. The maximum-norms of the numerical errors and corresponding convergence rates are given in Table 1, where the expected temporal convergence rates ($1$ for ETD1 and $2$ for ETDRK2) are obviously observed in both cases of integrable and non-integrable kernels. It is also easy to see that the numerical errors are almost independent of the choices of $\delta$ and $\alpha$.

Next, we tested the convergence with respect to the spatial size $h$ by fixing $\delta=2$ and $\tau=T$. The numerical solution of the NAC equation obtained by the ETDRK2 scheme with $N=4096$ is treated as the benchmark for computing the errors of the numerical solutions obtained with $N=2^{k}$ with $k=4,5,\dots,10$. The numerical errors in the maximum-norm sense are presented in Table 2, and it is observed that the convergence rates with respect to $h$ are almost of second order in both cases of integrable and non-integrable kernels, which is again consistent with the theoretical results.

We also investigated the limit behaviors of the numerical solutions of (1) as $\delta\to 0$. By fixing $N=4096$ and $\tau=T$, we calculated the numerical solutions of the NAC equation obtained by the ETDRK2 scheme (18) with various $\delta$’s and compared them with the numerical solution of the LAC equation. Table 3 collects the errors between the nonlocal and local numerical solutions in the maximum-norm sense and the second order convergence with respect to $\delta$ is obviously observed.

6.2 Stability tests

For the case $\rho_{\delta}(|\bm{s}|)\in L^{1}(\mathbb{R}^{2})$, i.e., $\alpha\in[0,2)$, it has been proved in [15] that the steady state solution $u^{*}$ to the NAC equation (1) is continuous if $\varepsilon^{2}C_{\delta}\geq 1$, where

[TABLE]

Under certain assumptions, if $\varepsilon^{2}C_{\delta}<1$, the locally increasing $u^{*}$ has a discontinuity at $x_{*}$ with the jump

[TABLE]

Example 6.2**.**

We simulate the NAC equation (1) with a random initial data ranging from $-0.9$ to $0.9$ uniformly generated on the $512\times 512$ mesh. We set the interfacial parameter $\varepsilon=0.1$ and adopt the kernel (52) with $\alpha=1$ and various $\delta$’s. For the comparison, we also simulate the LAC equation (5) with the same settings. The time step is set to be $\tau=0.01$ for all cases.

Under these settings, the critical value of $\delta$ to satisfy $\varepsilon^{2}C_{\delta}=1$ is $\delta_{0}=2\sqrt{3}\varepsilon$. The three rows in Fig. 1 correspond to the evolutions of phase structures governed by the LAC equation and the NAC equation with $\delta=3\varepsilon$ and $\delta=4\varepsilon$ at times $t=6$, $14$, $50$, and $180$, respectively. Fig. 2 presents the evolutions of the corresponding maximum-norms and the energies of the numerical solutions, respectively. It is observed in all cases that the discrete maximum principle is preserved perfectly and the discrete energy decays monotonically. It is easy to see that the dynamics of the NAC equation with $\delta=3\varepsilon$ is quite similar to that of the LAC equation. The evolution processes of these two cases reach the steady states at about $t=190$ and $t=370$, respectively, while the evolution of the NAC equation with $\delta=4\varepsilon$ lasts much longer time. In addition, the NAC equation with $\delta=3\varepsilon$ has thinner and sharper interface than the LAC equation but wider interface than the NAC equation equation with $\delta=4\varepsilon$. Actually, the interface in the case $\delta=4\varepsilon$ is discontinuous since the condition $\varepsilon^{2}C_{\delta}<1$ holds. The discontinuities in the solutions will be investigated further in the next example.

6.3 Discontinuity in the steady state solution

Example 6.3**.**

We simulate the evolution of a bubble governed by the NAC equation (1) starting with a smooth initial configuration (See Fig. 3). Again, we set the interfacial parameter $\varepsilon=0.1$ and adopt the kernel (52) with $\alpha=1$ and various $\delta$’s. The parameters of space-time mesh are set to be $\tau=0.01$ and $N=2048$ for all cases.

This example is devoted to the relationship between the discontinuities in the steady state solutions and the horizon parameter $\delta$. Under the settings of the parameters given above, it is known from (53) that the theoretical values of the jumps occurring at the discontinuity points can be formulated as

[TABLE]

We chose several $\delta$’s ($\delta=0.8,1.6,3.2$) larger than $\delta_{0}$ to observe the discontinuities and the jumps in the numerical results, and for the comparison, we also considered one case ($\delta=0.2$) with $\delta$ smaller than the critical value. Table 4 collects the theoretical and numerically computed jumps occurring at the discontinuity points in the steady state solutions with various $\delta$’s. It is observed that the numerical jumps match the theoretical values very well.

Fig. 4 presents the evolutions of the bubble governed by the NAC equation with $\delta=0.2$ ($<\delta_{0}$), $\delta=0.8$ and $\delta=3.2$ (both $>\delta_{0}$), respectively. In each row, the first three graphs give the surface-projection views of the numerical solutions at several times and the last graph cross-section views with $y=\pi$ by zooming-in around the interface. For the case $\delta=0.2$, the bubble shrinks quickly and disappears finally, which is similar to the process of the shrinkage occurring in the case of the LAC equation (see [8]). The evolutions for the cases $\delta=0.8$ and $\delta=3.2$ are similar: the bubble does not shrink and the interface turns sharper and sharper so that the solution preforms discontinuity on the interface after some times and reaches the steady state with the expected jump. It is seen from this example that the NAC equation with small $\delta$ has more similar dynamics with the local model, which is consistent with the observations in Example 6.2, while the NAC equation with large $\delta$, especially larger than $\delta_{0}$, leads to the steady state solution within the discontinuity even though the initial state is smooth.

7 Conclusions

We designed and analyzed maximum principle preserving numerical schemes of for solving nonlocal Allen-Cahn equation by using the quadrature-based finite difference method for spatial discretizations and the exponential time differencing method for temporal integrations. Especially, we developed the first order ETD and second order ETD Runge-Kutta schemes, derive for both schemes the error estimates, and prove their energy stability as well as the asymptotic compatibility, a special convergence considered for the numerical approximations of nonlocal models. Numerical experiments are carried out to verify the theoretical results and to study some more interesting properties of the solutions caused by the nonlocality. The maximum principle preserving schemes studied here are up to the second order in time. Whether higher order numerical schemes can preserve the maximum principle still remains open and is one of our future works. In addition, for some other models, for instance, the nonlocal Cahn-Hilliard equation [13, 22], the solution does not possess the maximum principle but is $L^{\infty}$ stable instead. Numerical schemes naturally inheriting the $L^{\infty}$ stability, weaker than the maximum principle, are also worthy of study.