Barycentric interpolation collocation method based on Crank-Nicolson scheme for the Allen-Cahn equation

1.   Introduction

In this paper, we consider the following Allen-Cahn equation

with the Dirichlet boundary condition

where $\Omega$ is a bounded domain in $R^d$ $(d = 1, 2)$. The positive parameter $\varepsilon$ representing generally the interfacial width is a very small constant and the reaction term $f(u) = F^{'}(u)$ with $F(u) = \dfrac{1}{4}(u^{2}-1)^{2}$. The function $u(\textbf{x}, t)$ can be viewed as the difference between the concentrations of the two components of the mixture or be defined as the order parameter that indicates the local state of the whole system. For instance, $u = 1$ and $u = -1$ stand for two different phases. It is well known that the Allen-Cahn equation possesses energy-decay property associating with the following total free energy functional

In fact, taking the inner product for equation (1.1) with $u_t$, we can obtain

The derivative of the energy $E(u(t))$ with respect to $t$ gives

which implys the total energy is decreasing in time, namely

The Allen-Cahn equation introduced by Allen and Cahn ^[1] in 1979 is a kind of non-homogeneous semi-linear poisson equation, which often occurs in convection diffusion equations in computational fluid dynamics or reaction-diffusion problems in material science. Hitherto, the Allen-Cahn equation has been widely used to model mathematically in crystal growth ^[2], image analysis ^[3] and average curvature-flow rate ^[4]. In addition, Allen-Cahn equation can also be used to describe competition of biological populations and phenomenon of exclusion ^[5]. It has become a basic model to study interfacial dynamics and phase transitions in material science ^[6]. However, the analytical solutions of the phase-field models can not be obtained so that it is greatly necessary and significant to develope stable, efficient and highly accurate numerical methods.

During the past two decades, enormous works have been devoted to numerical simulations and numerical analysis for the Allen-Cahn equation such as finite difference ^[7,8,9], finite element ^[10,11], spectral method ^[12]. Zhai et al. in ^[8] proposed a novel linearized high-order compact difference method for Allen-Cahn equation with different boundary conditions. Hou et al. presented numerical analysis of fully discretized Crank-Nicolson scheme for fractional-in-space Allen-Cahn equations and proved that the second order temporal discretization scheme preserves the discrete maximum principle in ^[9]. Feng and Li ^[10] developed two fully discrete interior penalty discontinuous Galerkin scheme for Allen-Cahn equation. Li et al. ^[11] introduced an unconditionally energy stable finite element method for Allen-Cahn equation. But the scheme is only second-order accurate in both space and time. Weng and Tang ^[12] presented two unconditionally stable second-order operator splitting approaches based on Fourier spectral for solving the Allen-Cahn equation. In addition, the operator splitting method combined with the finite difference method or finite element method for Allen-Cahn equation have also been considered in ^[13,14].

Most of above-mentioned works depend on the mesh generation to solve differential equations, which causes some difficulties especially for high dimension problems with irregular domains. In recent years, meshless methods have captured attention of scholars due to it's advantage of meshless and ability of treating irregular domains. As a novel meshless method, barycentric interpolation method is a good choice for dealing with polynomial interpolations since it is accurate highly, fast, stable and easy on program implementation. Barycentric interpolation method include barycentric Lagrange interpolation and barycentric rational interpolation, both of which can be written the unified formula of barycentric interpolation. The weight functions become extremely big when the nodes are of uniform distribution, which leads to the Runge phenomenon and ruins the merits of the barycentric Lagrange interpolation. Luckily, if the nodes obey the density proportion $(1-x^2)^{-\frac{1}{2}}$ such as the families of Chebyshev points that is the simplest clustered point sets, the barycentric Lagrange interpolation has a good numerical stability. Barycentric interpolation also can effecively avoid the accumulated errors caused by difference scheme. Readers can refer to ^[15,16] for detailed introduction and derivation of the barycentric interpolation. The collocation method based on barycentric interpolation was lately extended to solve various integral equations and partial differential equations including high-dimensional Fredholm integral equation of the second kind ^[17], Volterra integral equations with weakly singular kernels^[18], nonlinear parabolic equations ^[19] and 2D viscoelastic wave equation ^[20].

To our knowledge, the works of the error analysis for barycentric interpolation collocation method are comparatively sparse. Recently, Yi and Yao in ^[21] put forward a steady barycentric Lagrange interpolation method and presented error analysis of system for solving the time-fractional telegraph equation. Based on the above works, we focus on a fully discrete scheme for the Allen-Cahn equation, which is a second-order Crank-Nicolson scheme in time combined with the barycentric interpolation collocation method in space. Moreover, we will give consistency analysis of the semi-discretized scheme in space.

The remainder of the paper is structured as follows: In section 2, the barycentric interpolation collocation method is introduced. Semi-discretized scheme in space and corresponding consistency analysis are carried out in detail for the Allen-Cahn equation in section 3. In section 4, fully discrete scheme based on Crank-Nicolson scheme is presented. Nmerical examples test the accuracy and efficiency of proposed algorithm in section 5 and some conclusions are given in section 6.

2.   Barycentric interpolation collocation method

2.1.   Barycentric Lagrange interpolation with the Chebyshev points

Suppose $n+1$ distinct interpolation nodes $x_j$ be given, together with corresponding a set of real numbers $y_j$ $(j = 0, 1, \ldots, n)$. Let $p(x)$ denotes the polynomial of degree at most $n$, satisfying $p(x_j) = y_j$ $(j = 0, 1, \ldots, n)$. As we all known, such polynomial $p(x)$ is unique and can be written in Lagrange form as

where $L_{j}(x)$ is the basis function in Lagrange interpolation, which possesses the following properties

and

Let

Defining the barycentric weights by

the $L_{j}(x)$ can be expressed as

Inserting Eq (2.6) into Eq (2.1), we get

Combining Eq (2.3) with Eq (2.6), we have the following identical equation

Dividing Eq (2.8) by Eq (2.7) and cancelling the common factor $l(x)$, the barycentric Lagrange interpolation formula for $p(x)$ can be obtained

In the paper, we choose the Chebyshev points given by

which will ensure that above polynomial interpolation has good numerical stability.

2.2.   Barycentric rational interpolation

The rational function interpolation based on the idea of mixed function can effectively overcome the instability of interpolation.

Give $n+1$ distinct interpolation nodes $x_j$ equipped with corresponding numbers $y_j$ $(j = 0, 1, \ldots, n)$. Choose a integer $d$ $(0 \leq d \leq n)$ and let $p_{i}(x)$ denotes the polynomial of degree at most $d$ interpolating $d+1$ point pairs $(x_i, y_i), (x_{i+1}, y_{i+1}), \ldots, (x_{i+d}, y_{i+d})$ for each $i = 0, 1, \ldots, n-d$. Then, we set

where

Obviously, rational function $r(x)$ satisfy $r(x_i) = y_i$ $(i = 0, 1, \ldots, n)$. In order to obtain barycentric interpolation form of Eq (2.11), we write $p_{i}(x)$ in Lagrange interpolation

Insert Eq (2.13) into numerator of Eq (2.11), we have

with interpolation weight $\omega_{k} = \sum \limits_{i \in J_{k}}(-1)^{i} \prod \limits_{j = i, j \neq k}^{i+d} \dfrac{1}{\left(x_{k}-x_{j}\right)}$, where $J_k = \left \{i \in I: k-d \leq i \leq k\right \}$ represents index set, $I = \left \{0, 1, \ldots, n \right \}$.

Noticing the identical equation

we thus get

Combining Eqs (2.11), (2.14) and (2.16), the barycentric rational interpolation formula for $r(x)$ can be obtained

2.3.   Differential matrix

The derivative of $p(x)$ defined as Eq (2.9) with respect to $x$ as

The first and second order differentiation matrices can be obtained by the following formula ^[15]

By mathematical induction, we have the following $v$-order differential matrix

3.   Semi-discretized system and consistency analysis

3.1.   Semi-discretized scheme based on barycentric interpolation collocation method

Choose a rectangular domain $\Omega = [a, b]\times[c, d]$. Partition respectively the interval $[a, b]$, $[c, d]$ into $M+1$, $N+1$ distinct Chebyshev nodes: $a = x_{0} < x_{1} < \ldots < x_{M} = b$, $c = y_{0} < y_{1} < \ldots < y_{N} = d$. Denoting $u(x_i, y, t) = u_i(y, t), i = 0, 1, \ldots, M$, and fixing variable $y$, the unknown function $u(x, y, t)$ can be written in barycentric interpolation form

where $\gamma_{l}(x)$ is the basis function in barycentric interpolation on the direction of $x$.

Substituting Eq (3.1) into Eq (1.1) and making Eq (1.1) be identical at the nodes $x_{i}, i = 0, \ldots, M$, we have

where $\gamma''_{l}(x_i) = \dfrac {{\mathrm{d}}^2 \gamma_{l}(x_i)}{\mathrm{d} x^2} = C_{il}^{(2)}$.

Equation (3.2) can be written in matrix form, namely

Denote $u_i(y_j, t) = u_{ij}(t)$. Similarly, $u_i(y, t)$ can be written in barycentric interpolation form

where $\beta_{q}(y)$ is the basis function in barycentric interpolation on the direction of $y$.

Inserting Eq (3.4) into Eq (3.3) and making Eq (3.3) be identical at the nodes $y_{j}, j = 0, 1, \ldots, N$, we get the following ODE system

where ${\beta_q ''}(y_{j}) = \dfrac {{\mathrm{d}}^2 \beta_{q}(y_j)}{\mathrm{d} y^2} = D_{j q}^{(2)}$.

Introduce the following notations

Equation (3.5) can be rewritten in the following matrix form

where the sign $\otimes$ represents the Kronecker product of the matrix, $\boldsymbol{C}^{(2)}$ and $\boldsymbol{D}^{(2)}$ stand for second order differential matrix on nodes $x_{0}, x_{1}, \ldots, x_{M}$ and on nodes $y_{0}, y_{1}, \ldots, y_{N}$, severally. $\boldsymbol{I}_{M}$ and $\boldsymbol{I}_{N}$ being the identity matrix of $M+1$, $N+1$ order, respectively.

3.2.   Consistency analysis for the semi-discretized scheme

In this part, we present consistency estimates of the semi-discretized system (3.6) with the collocation method. Let $p(x)$ is the Lagrange interpolation function approximating $u(x)$. According to interpolation remainder theorem, we have

Estimates of (3.7) are considered by the following Lemma that have been proved in Yi ^[21].

Lemma 1. (see Yi ^[21] ) If $u(x) \in {C^{m + 1}}([a, b])$, then the following estimates for functional $e(x)$ defined as (3.7) hold

where $C_1$, $C_1^{*}$ and $C_1^{**}$ are constant independent of $x$, $e$ is natural logarithm and $l_x$ represents the length of the interval $[a, b]$.

Let $p(x, y)$ is Lagrange interpolation function of $u(x, y)$ satisfying $p({x_m}, {y_n}) = u({x_m}, {y_n})$. Defining the error function as

we have the following results based on the Lemma 1.

Theorem 1. If $u\in C^{\bar{m}+1}([a, b]) \times [c, d])$, where $\bar{m} = max\{m, n\}$. Then the following estimates for functional $e(x, y)$ defined as (3.9) hold

where $l_y$ represents the length of the interval $[c, d]$. Similarly, we get

Set $u(x_m, y_n, t)$ be the numerical solution of $u(x, y, t)$, then we have

and

where $\Pi$ is differential operator.

Theorem 2. If $u\in C^0([0, T]), C^{\bar{m}+1}([a, b) \times [c, d])$, where $\bar{m} = max\{m, n\}$, let $u(x_m, y_n, t):\Pi u({x_m}, {y_n}, t) = 0$ and assume that $f(u)$ satisfies the Lipschitz condition, we have

Proof. Following (1.1) and (3.12), we get

where

For $R_1$, we have

Applying Theorem 1, we obtain

Analogously, we estimate $R_2$ and $R_3$ as

For the term $R_4$, we have

Substituting (3.17)-(3.20) into (3.15), this completes the proof.

From Theorem 2, we can find that the order of the difference operator $\Pi$ determines the consistency rate of the semi-discrete scheme.

4.   Fully discretized scheme based on Crank-Nicolson scheme

In the section, we will solve the ODE system (3.6) by the Crank-Nicolson scheme. Partition the interval (0, T] into a uniform mesh with the time step $\tau = \dfrac{T}{l}$: $0 = t_{0} < t_{1} < \ldots < t_{l} = T$. Let ${\bf U}^k = {\bf U}(t_k), k = 0, 1, \ldots, l.$

the nonlinear term of (3.6) is expanded using Taylor formula at the node vector ${{\bf{U}}^k}$ $(k = 0, 1, \ldots, l)$, we get

Crank-Nicolson scheme is used for time discretization, we have

Inserting ${\bf U}^{k+\frac{1}{2}} = \dfrac{1}{2} \left({\bf U}^{k+1}+{\bf U}^{k}\right)$ into Eq (4.2), we obtain fully discretized scheme as

5.   Numerical experiments

In this section, we will perform several numerical examples to test the accuracy and energy decline property of the proposed scheme. we consider convergence tests in the first example and three 2D problems with different initial conditions in the next examples.

For convenience, introducing the following error notations

where $u_{h}$ and $u_{e}$ denote the numerical solution and the exact solution of the problem, respectively, $\parallel \cdot \parallel_{\infty}$ is the $L^{\infty}$ norm.

5.1.   Convergence tests

We test numerical convergence order in time for 1D Allen-Cahn equation at first, the initial and boundary condition are given the following analytical solution

where

Fix $M = 30$, and time step $\tau$ is varied. The errors in maximum norm and convergence rates are displayed in Table 1, from which we see that the scheme is second order in time for 1D Allen-Cahn equation.

Next, we verify the accuracy and convergence rate in time of the algorithm when applying to the 2D Allen-Cahn equation with the inhomogeneous term $g(x, y, t)$. We consider the following problem on $\Omega = {[-1, 1]^2} \times (0, T]$ as

where

Take the following simulation parameters

and vary numebers of mesh node $M, N$. The error comparison results of different discretization schemes in space can be obtained, as shown in Table 2. One can see that both barycentric interpolation collocation methods have high precision, as they only use $8\times 8$ mesh nodes can achieve the accuracy of second-order central difference approach with $40 \times 40$ mesh nodes. In addition, it is also observed that the accuracy of barycentric Lagrange interpolation is slightly higher than that of barycentric rational interpolation especially taking fewer mesh nodes from Table 2. In addition, exponential convergence property of two methods also can be observed from Figure 1.

Choose $M = N = 20$ and vary the temporal step $\tau$. Table 3 indicates the derived scheme is indeed of second order in time for 2D Allen-Cahn equation.

After that, we will show the energy decay property of the 2D Allen-Cahn equation. Define the discrete energy function as

where $u_{ij}^{k}$ and $E^{h}(u^{k})$ represent the numerical solution and the energy value at the $k$-th time separately, $h_{ij}$ denotes spatial step size between adjacent Chebyshev nodes. The space diagram of 2D Allen-Cahn equation is shown in Figure 2. And from Figure 3, it is obvious that the energy is decreasing with time.

5.2.   Example 2

Consider 2D Allen-Cahn equation on $\Omega = {[-2, 2]^2} \times (0, T]$ with the following initial condition

where

We take simulation parameters as

From Figure 4, we can see clearly that the connected interface splits into two curves at first. Then, the two components of the interface develop circular shapes. Eventually, the diameters of the two particles decrease to zero until they collapse.

5.3.   Example 3

Consider 2D Allen-Cahn equation on $\Omega = {(- 1, 1)^2} \times (0, T]$. Taking the following initial condition

and parameters

Figure 5 shows process of the evolution of the numerical solution from forming the transition layer, metastable state and finally reaching the steady state, respectively.

5.4.   Example 4

Consider 2D Allen-Cahn equation in $\Omega = {(- 1, 1)^2} \times (0, T]$ with simulation parameters

The initial condition is taken as the randomly perturbed condition fields as follows

where ${\rm {rand}}(x, y)$ represents a pair of random numbers generated between -1 and 1.

The coarsening phenomena of Allen-cahn equation with the evolution of time are presented in Figure 6, from which we see that the solution reaches the steady state at $t = 0.1s$.

6.   Conclusions

In this work, by a combination of the Crank-Nicolson scheme and barycentric interpolation collocation method, an efficient numerical scheme for Allen-Cahn equation is developed. Besides, the consistency analysis for the semi-discretized scheme in space is derived. The numerical examples show that our scheme is second order in time and convergent exponentially in space. Energy dissipation law is also checked numerically. In future work, we plan to give the error estimate of fully discretized scheme and generalize this method to other models.

Acknowledgments

This work is in part supported by the NSF of China (Nos 11701197), the Fundamental Research Funds for the Central Universities (No. ZQN-702), the Key Laboratory of Intelligent Computing and Information Processing of Ministry of Education (No. 2020ICIP03) and the Natural Science Foundation of Fujian Province (No. 2020J01074).

Conflict of interest

The authors declare that there are no conflicts of interest regarding the publication of this manuscript.