Convergence analysis of exponential time differencing scheme for the nonlocal Cahn-Hilliard equation

The Cahn-Hilliard (CH) equation was originally proposed in [6] as one of the most typical phase field models which provides a macroscopic description of phase separation and microstructure evolution in binary alloy systems. As a nonlocal variant of the classic CH equation, the Nonlocal Cahn-Hilliard (NCH) equation has attracted increasing attention and has been widely applied in diverse fields ranging from chemistry, material science and image processing [20, 14, 1, 16, 2]. In this paper, we consider the following NCH equation with the periodic boundary conditions [2, 9]:

where $\Omega=\prod_{i=1}^{d}(-X_{i},X_{i})$ is a rectangular domain in $\mathbb{R}^{d}(d=1,2,3)$ and $u=u(x,t)$ is an order parameter subject to the initial condition $u(x,0)=u_{0}(x)$, $T>0$ is the terminal time, $\varepsilon>0$ is an interfacial parameter. The function $f(u)=F^{\prime}(u)$ and $F(u)=\frac{1}{4}(u^{2}-1)^{2}$ is a double well function. The nonlocal linear operator $\mathcal{L}$ is defined by:

where $J$ is a nonnegative, $\Omega$-periodic radial kernel function and has a finite second moment in $\Omega$.

The linear operator $\mathcal{L}$ with the kernel function $J$ is self-adjoint and positive semi-definite. Further, if $J$ is integrable, then $\mathcal{L}v=(J\ast 1)v-J\ast v$, where

We assume $\gamma_{0}:=\epsilon^{2}(J\ast 1)-1>0$ to ensure the NCH equation (1.1) is positive diffusive [4, 3].

The NCH equation (1.1) can be view as the $H^{-1}$ gradient flow with respect to the free energy functional

where $(\cdot,\cdot)$ denotes the standard $L^{2}$ inner product on $\Omega$. For the smooth solutions of (1.1), the total mass is conserved:

In this paper, we will only consider initial data $u_{0}$ with mean zero. Then, the fractional Laplacian operator $|\nabla|^{s}=(-\Delta)^{\frac{s}{2}}$ for $s<0$ is well defined [25]. Due to the gradient structure of (1.1), the energy dissipation laws hold:

There have been many works on both theoretical analysis and numerical methods for the NCH equation. In mathematical analysis, the well-posedness of the NCH equation equipped with a Neumann or Dirichlet boundary conditions were studied in [3, 4] with an integrable kernel function and in [19] claimed that the existence and uniqueness of the solution to the NCH equation subject to the periodic boundary condition may be established by using a similar technique.

For numerical methods, because of the energetic variational structure inherent in the nonlocal phase field model, an important fact is that the exact solution decrease the energy in time to NCH equation. Therefore, it is essential to develop the so called energy stable numerical algorithms share this key property at the discrete time level. Energy stable numerical schemes of the NCH equation including convex splitting schemes [19, 18], exponential time differencing (ETD) schemes [30, 31, 15], stabilized schemes [9, 27, 28], and so on.

The convergence of numerical solution to NCH equation are established in [19, 17, 26, 27, 28]. In [19, 17], convex splitting schemes were introduced for the NCH equation, demonstrating energy stability and convergence properties. These schemes handled the nonlinear term implicitly and the nonlocal term explicitly, leading to effective numerical solutions but nonlinear iterations were required. Du et al. proposed energy-stable linear semi-implicit schemes in [9], utilizing stabilization techniques to avoid nonlinear iterations. The work in [26] presented the convergence in the discrete $H^{-1}$ norm and established the $\ell^{\infty}$ bound of the numerical solutions for the first-order stabilized linear semi-implicit scheme proposed in [9]. Furthermore, Li et al. [27, 28] developed two energy-stable and convergent second-order linear numerical schemes for the NCH equation. These schemes involved combining a modified Crank-Nicolson (CN) approximation with the Adams-Bashforth extrapolation and stability of the second-order backward differentiation formula (BDF2) for the time discretization, resulting in accurate and stable numerical solutions for the NCH equation.

However, the convergence analysis of ETD schemes for the NCH equation is more challenging and there are currently no results, one of the reasons being that the lack of higher-order diffusion term and the Laplacian of nonlinear term. The ETD schemes [5, 7, 12, 13, 21], which involve exact integration of the target equation followed by an explicit approximation of the temporal integral of the nonlinear term. Exact evaluation of the linear terms makes the ETD schemes achieve high accuracy and satisfactory stability when dealing with stiff differential equations. The first- and second-order ETD schemes were applied to the Nonlocal Allen-Cahn equation and have been analyzed energy stable and $\ell^{\infty}$ convergence [10], and further extended to a class of semi-linear parabolic equations [11]. Convergence and the $\ell^{\infty}$ bound of numerical solutions were proved of ETD schemes to classic CH equation [24].

Compare with the classic CH equation or Nonlocal Allen-Cahn equation, the NCH equation equipped with nonlocal diffusion operator and has more complicated nonlinear term, so that the analytical techniques mentioned above are hardly applicable to the NCH equation. In this paper, we establish the convergence of the fully discrete first-order ETD scheme (ETD1) and second-order ETD multistep scheme (ETD2) for the NCH equation which are developed our recent work [31]. Our analysis is mainly based on two important observations:

• 1.

  Unlike the standard error analysis based on the error equation (3.5), our analysis based on two new error decomposition formulas

  		$\displaystyle e^{n}:=\tilde{e}^{n}+\tau\mathcal{P}_{N}u_{\tau,1}+\tau^{2}\mathcal{P}_{N}u_{\tau,2}\quad(\text{ETD1 scheme});$
  		$\displaystyle e^{n}:=\tilde{e}^{n}+\tau^{2}\mathcal{P}_{N}u_{\tau,3}\hskip 52.63777pt(\text{ETD2 scheme}).$

  See the details in (3.9). This formula allows us to establish the higher order consistency and obtain the $\ell^{\infty}$ bound of the numerical solutions.

• 2.

  We adopt $(-\Delta_{N})^{-1}\tilde{e}^{n+1}$ to test the error functions in (3.23) and (3.2.2), instead of $\tilde{e}^{n+1}$ as in a previous work [24], where $(-\Delta_{N})^{-1}$ is a spatial discrete operator defined in the Appendix A. This method can be used to overcome the possible instability caused by the absence of high-order diffusion terms in the NCH equation.

Based on the above two points, we prove the convergence in the discrete $H^{-1}$ norm and obtain the $\ell^{\infty}$ bound of numerical solutions under some natural constraints on the space-time step sizes.

The rest of this paper is organized as follows. The fully discrete ETD1 and ETD2 schemes are constructed in Section 2, and some notations, definitions and useful lemmas are also introduced. In Section 3, we prove the convergence by the higher-order consistency estimate. In addition, the $\ell^{\infty}$ bound of numerical solutions are also obtained. In Section 4, some numerical experiments are carried out to test the convergence in time level. The coarsening dynamics are simulated to show the long time behaviors of ETD2. Some concluding remarks are given in Section 5.

In this section, we present the fully discrete exponential time differencing schemes for the NCH equation, where the Fourier spectral collocation method is adopted for the spatial discretization.

We consider the square domain $\Omega=(-X,X)^{2}$ and define the spatial grid $\Omega_{h}=\{(x_{i},y_{j})=(-X+ih,-X+jh),\;1\leq i,j\leq N\},$ where the space step size $h=\frac{2X}{N}$ ($N$ is even). We need two spaces. Let $\mathcal{M}_{h}$ be the space of all the $\Omega_{h}$-periodic grid functions. That means if $f\in\mathcal{M}_{h}$, then $f$ is a discrete function defined on the discrete grid $\Omega_{h}$. Let $\mathcal{M}_{h}^{0}$ be the space of the zero-mean grid functions. That means if $f\in\mathcal{M}_{h}^{0}$, then $\overline{f}:=\frac{h^{2}}{4X^{2}}\sum_{(i,j)\in S_{h}}f_{i,j}=0$, where $S_{h}$ is the index set defined in the Appendix A.

We denote $\Delta_{N}$ and $\mathcal{L}_{N}$ are Fourier spectral collocation space approximation operators of $\Delta$ and $\mathcal{L}$. These detailed definitions and some related properties are given in the Appendix A. For any $f,g\in\mathcal{M}_{h}$, the discrete inner product $\langle\cdot,\cdot\rangle$, $\langle\cdot,\cdot\rangle_{-1,N}$ and norm $\|\cdot\|_{2}$, $\|\cdot\|_{\infty}$, $\|\cdot\|_{-1,N}$, and the discrete convolution $f\circledast\phi$ are defined in the Appendix A.

It is well known [12] that a suitable linear operator splitting can improve the stability. So we can introduce a stabilizing parameter $\kappa>0$ and define

Now $L_{h}$ is self-adjoint and positive definite on $\mathcal{M}_{h}^{0}$. Then the space semi-discrete equation for (1.1) is to find $U:[0,T]\to\mathcal{M}_{h}^{0}$ such that

where the initial value $U^{0}\in\mathcal{M}_{h}^{0}$ is given.

Let $\{t_{n}=n\tau:0\leq n\leq N_{t}\}$, where $\tau=\frac{T}{N_{t}}$ is the time step size and $N_{t}$ is a given positive integer. By using the property of the differentiation of matrix exponentials, the solution of the equation (2.2) satisfies

The key to construct the ETD schemes is to approximate the nonlinear term $f_{\kappa}(U(t_{n}+s))$ in (2.3) by polynomial interpolation, and then precisely integrate the interpolation.

The ETD1 scheme comes from approximating $f_{\kappa}(U(t_{n}+s))$ by the constant $f_{\kappa}(U(t_{n}))$, given by

where $\phi_{-1}(a)=e^{-a},\phi_{0}(a)=\frac{1-e^{-a}}{a},a\neq 0$. The ETD2 scheme is obtained by approximating $f_{\kappa}(U(t_{n}+s))$ by a linear extrapolation based on $f_{\kappa}(U(t_{n}))$ and $f_{\kappa}(U(t_{n-1}))$, given by

where $\phi_{1}(a)=\frac{a-1+e^{-a}}{a^{2}},a\neq 0.$

The ETD2 scheme (2.5) is a two-step algorithm, an accurate approximation for the value at $t_{1}$ is needed. Usually, $U^{1}$ can be computed by the ETD1. But in this paper, we choose a higher-order approximation for $U^{1}$, to facilitate the higher-order consistency analysis presented in the later sections. For example, we can apply the second-order Runge-Kutta (RK2) method in first step, which yields a third-order temporal accuracy at $t_{1}$ for exact initial data $U^{0}$.

The following lemmas are useful in our following convergence analysis.

In this section, we analyze the convergence of the ETD1 scheme and ETD2 scheme, and give an optimal rate error estimates in $\ell^{\infty}(0,T;H_{h}^{-1})\cap\ell^{2}(0,T,\ell^{2})$. To perform the convergence analysis, let’s first introduce the concept and preliminary.

For a linear symmetric operator $Q:\mathcal{M}_{h}\to\mathcal{M}_{h}$, we define the norm $\interleave Q\interleave$ by the spectrum radius of $Q$, i.e., $\interleave Q\interleave=\max\{|\lambda|:\lambda\in\sigma(Q)\}$, where $\sigma(Q)$ is the set of all the eigenvalues of $Q$. It holds that $\|Qv\|_{2}\leq\interleave Q\interleave\cdot\|v\|_{2},\;\forall v\in\mathcal{M}_{h}.$

Denote by $u_{e}$ the exact solution of (1.1). The existence and uniqueness of $u_{e}$ may be obtained by using the techniques adopted in [3, 4], from which one can get the following estimate

Moreover, if the initial data is sufficient regularity, we can assume that the exact solution has regularity as

Let $\mathcal{B}^{K}$ be the space of trigonometric polynomials of degree up to $K:=\frac{N}{2}$. Let $\mathcal{P}_{N}:L^{2}(\Omega)\to\mathcal{B}^{K}$ be the $L^{2}$ orthogonal projection operator. We define $u_{N}(\cdot,t):=\mathcal{P}_{N}u_{e}(\cdot,t)$. If $u_{e}\in L^{\infty}(0,T;H^{\ell}_{per}(\Omega))$ for some $\ell\in\mathbb{N}$, the following estimate is standard [29]

By the orthogonality of Fourier projection $\mathcal{P}_{N}$, we have $\int_{\Omega}u_{N}(\cdot,t_{n})\,\mathrm{d}x=\int_{\Omega}u_{e}(\cdot,t_{n})\,\mathrm{d}x$. Noting that the exact solution $u_{e}$ is mass conservative, i.e., $\int_{\Omega}u_{e}(\cdot,t_{n})\,\mathrm{d}x=\int_{\Omega}u_{e}(\cdot,t_{n-1})\,\mathrm{d}x$. Thus, we obtain

Since $u_{N}\in\mathcal{B}^{K}$, we have $\int_{\Omega}u_{N}(\cdot,t_{n})\,\mathrm{d}x=h^{2}\sum_{(i,j)\in S_{h}}u_{N}(x_{i},y_{j},t_{n})$, that means the rectangular quadrature rule holds exactly for all $u_{N}\in\mathcal{B}^{K}$. The values of $u_{N}(t_{n})$ at discrete grid points are still denoted as $u_{N}(t_{n})$, i.e., $u_{N}(x_{i},y_{j},t_{n})=\mathcal{P}_{N}u_{e}(t_{n})|_{i,j}$. Recall  $\overline{u_{N}}(t_{n})=\frac{h^{2}}{4X^{2}}\sum_{(i,j)\in S_{h}}u_{N}(x_{i},y_{j},t_{n}),$ then from (3.4) we can get the mass conservative property of $u_{N}$ in the discrete average sense: $\overline{u_{N}}(t_{n})=\overline{u_{N}}(t_{n-1}).$

On the other hand, the solutions of the numerical schemes (2.4)-(2.5) are also mass conservative [30]: $\overline{U^{n}}=\overline{U^{n-1}},\;\forall n\in\mathbb{N}.$ We use the mass conservative for the initial value $U^{0}=u_{N}(t_{0})$, and define the error grid function

We have $\overline{e^{n}}=0$, so $\|e^{n}\|_{-1,N}$ is well defined. Under the regularity assumption (3.1), for the projection functions $u_{N}(t_{k}):=\mathcal{P}_{N}u_{e}(t_{k})$, we have

The main result of this work is the following theorem.

We provide a heuristic explanation of the main idea behind the proof. The usual convergence analysis starts directly from the error equation $e^{n}:=U^{n}-u_{N}(t_{n})$ defined in (3.5). However, our here adopts a new strategy. We construct two new error decomposition representation formulas

where $\tilde{e}^{n}:=U^{n}-\tilde{u}(t_{n})$ and $\tilde{u}(t_{n})$ is the constructed approximate solution defined in (3.14) for the ETD1 scheme and in (3.37) for the ETD2 scheme, $u_{\tau,1}$, $u_{\tau,2}$ and $u_{\tau,3}$ are temporal correction functions. By doing so, we can obtain an $\mathcal{O}(\tau^{3}+h^{m})$ truncation error. This higher consistency allows us to derive a higher order convergence estimate of the modified error function $\tilde{e}^{n+1}$ as $\|\tilde{e}^{n+1}\|_{-1,N}=\|\tilde{u}(t_{n+1})-U^{n+1}\|_{-1,N}=\mathcal{O}(\tau^{3}+h^{m})$ (see (3.31) and (3.49)), which in turn lead to $\|U^{n+1}\|_{\infty}<\infty$ under reasonable condition $\tau\leq Ch$. Due to the lack of higher-order diffusion term, we adopt $(-\Delta_{N})^{-1}\tilde{e}^{n+1}$ to test the error function in (3.23) and (3.2.2), instead of $\tilde{e}^{n+1}$ as that in [24]. Then, we estimate both sides of the error equation by using the result of the higher-order consistency analysis and induction argument, we can obtain an optimal convergence rate in $\ell^{\infty}(0,T;H_{h}^{-1})\cap\ell^{2}(0,T,\ell^{2})$ for ${e}^{n+1}$ of the proposed two numerical schemes.