Stability analysis of an explicit numerical scheme for the Allen-Cahn equation with high-order polynomial potentials

1.   Introduction

We investigate a stability analysis of an explicit finite scheme to numerically solve the Allen-Cahn (AC) model with a high-order polynomial potential [1]:

where $\phi({\bf x}, t)$ is the order parameter at two-dimensional space position ${\bf x} = (x, y)$ and time $t$, $\epsilon$ is positive, $F_\alpha(\phi)$ is a double-well potential, and ${\bf n}$ is normal to the boundary $\partial \Omega$. We consider the following high-order polynomial free energy potential [1,2,3]:

where $\alpha$ is a positive integer. Equation (1.1) can be rewritten as follows:

We note that if $\alpha = 1$, then Eq (1.4) becomes the classical AC equation [4], which is a partial differential equation describing phase separation in materials, such as solidification or liquid-vapor transitions. It models the evolution of a scalar order parameter to minimize interfacial energy and has been used in volume repair [5], oil pollution dynamics [6], 3D volume reconstruction from slice data [7,8], shape transformation [9], data assimilation [10], image segmentation [11], structural topology optimization [12], mathematical physics for phase transition studies [13], and numerous applications related to interfacial evolution.

Due to the fact that analytical solutions for the AC equation are often impractical or unattainable, numerical methods are necessary to obtain approximate solutions for the equation. Recently, many numerical methods for the AC equation have been proposed. Lee [14] introduced a linear convex splitting scheme and an energy-minimizing method on the mass-conserving space for solving the AC equation. Numerical simulations verified phase separation, dissipation of energy, and mass conservation. Zhang et al. [15] developed a class of up to third-order explicit structure-preserving schemes for solving two modified conservative AC equations. Furthermore, the authors developed and analyzed an unconditionally structural-preserving parametric single-step method of up to fourth-order accuracy for the conservative AC equation in [16]. Feng et al. [17] developed a maximum-principle preserving, unconditionally stable, and second-order FDM for the AC equation. Zhang et al. [18] developed the Runge-Kutta method for temporal integration and FDM for spatial discretization to propose high-order maximum principle-preserving schemes. Characteristics such as maximum principle preservation and energy dissipation were demonstrated through numerical experiments. Lee [19] presented a mathematical model and numerical method of area-minimizing surfaces using the AC equation. Deng and Weng [20] developed the barycentric interpolation collocation scheme based on the Crank-Nicolson method for the AC model. Park et al. [21] presented an unconditionally stable computational algorithm for the AC equation with logarithmic free energy. Zhang et al. [22] proposed and analyzed a class of up to eighth-order inequality-preserving two-step integrating factors using the Runge-Kutta method to solve the semilinear parabolic equation. Typically, fourth-order equations such as the Cahn-Hilliard (CH) type equations [23,24] require implicit numerical schemes for stability due to their severe time step restrictions with explicit schemes.

However, the AC equation is a second-order equation, and its explicit numerical scheme requires only a moderate time step condition for stability. Furthermore, as demonstrated in a series of previous research papers [25], the fully explicit numerical schemes of the AC equation have performed well in terms of both efficiency and accuracy. The primary objective of this paper is to study the stability condition of an explicit FDM to solve the AC model with high-order polynomial potential.

The contents of this article are as follows: In Section 2, a fully explicit computational method is described and stability analysis is conducted for the AC equation with the high-order polynomial free energy. In Section 3, numerical tests are performed to validate the derived theorem. In Section 4, we conclude.

2.   Numerical scheme and stability analysis

Now, we present a fully explicit numerical scheme of the AC equation with the high-order (higher than fourth) polynomial free energies and derive its temporal step restriction, which satisfies the maximum principle. Let $\Omega = (L_x, U_x) \times (L_y, U_y)$ be the domain, and $\Omega_{h} = \{ (x_i, y_j)|x_i = L_x+(i-0.5)h, \; y_j = L_y+(i-0.5)h, \; 1\le i \le N_x, \; 1\le j \le N_y \}$ be its discrete domain, where $h$ is spatial step size and $N_x$ and $N_y$ are integers. Let $\phi_{ij}^n = \phi(x_i, y_j, n \Delta t)$, where $\Delta t$ is the time step size. We define the discrete maximum norm by $\| \phi^n \|_\infty = \max\limits_{1\leq i \leq N_x, \; 1\leq j \leq N_y}|\phi_{ij}^n|$. The AC model (1.1) is discretized by using an explicit method as follows:

where we use a 5-point stencil for the discrete Laplacian operator. We can use the 9-stencil isotropic finite difference scheme [26] for the isotropic solutions. Here, we use the homogeneous Neumann boundary condition:

Now, we shall prove the following maximum principle theorem of the fully explicit scheme for the AC equation with high-order polynomial potentials.

Theorem 1. Assume that the initial condition satisfies $\| \phi^0 \|_\infty \leq 1$. Then, the fully explicit numerical method (2.1) maintains the boundedness of the numerical solutions at any discrete time:

if the temporal step size satisfies

Proof. We have $\Vert \phi^0 \Vert_{\infty} \leq 1$ by the assumption. By mathematical induction, suppose that $\Vert \phi^n \Vert_{\infty} \leq 1$. The explicit numerical method for the AC model (2.1) can be reformulated as follows:

Let us first consider the following case:

Case 1.

In Eq (2.3), we have $\phi_{ij}^{n+1} \geq -1$. From Eq (2.3), we can derive the following inequality:

We want to determine the condition for $\Delta t$ such that

Then, it can be reformulated as

If $\phi_{ij}^n = 1$, the time step is free from any constraints; and otherwise, i.e., $\phi_{ij}^n < 1$, by dividing both sides of (2.6) by $(1-\phi_{ij}^n) > 0$, we have

Hence, $\Delta t$ is subject to the following restriction:

Here, because $\left| \phi_{ij}^n \right| \leq 1$, we derive the following condition:

Hence, when $\Delta t \leq 0.5\epsilon^2 h^2 /\left(h^2 +2\epsilon^2 \right)$, we have $\| \phi^{n+1} \|_\infty \leq 1$. Next, let us consider the following case:

Case 2.

From Eq (2.3), $\phi_{ij}^{n+1}\leq1$ always holds true. From Eq (2.3), we get the following inequality:

From (2.8), to satisfy that $-1\leq\phi_{ij}^{n+1}$, we consider

It can be reformulated as

If $\phi_{ij}^n = -1$, the time step is free from any constraints; and otherwise, i.e., $\phi_{ij}^n > -1$, by dividing both sides of Eq (2.10) by $-(1+\phi_{ij}^n) < 0$, we have

If $-1 \le \phi_{ij}^n \le 1$, since $0 \le \left(1 -(\phi_{ij}^n)^{\alpha} \right)\left(1 -\phi_{ij}^n +\cdots +(-1)^{\alpha-1}(\phi_{ij}^n)^{\alpha-1} \right) \le 2\alpha$, we obtain the following inequality as

Thus, we derive the following condition:

Therefore, $\Delta t \leq 0.5h^2\epsilon^2/(h^2+2\epsilon^2)$ results in $\| \phi^{n+1} \|_\infty \leq 1$. □

We note that the case of $\alpha = 1$ was proved in [27]. Here, we generalized the stability of the numerical methods of the AC model with high-order polynomial potentials for all $\alpha \ge 1$ by mathematical induction. Additionally, this time step constraint can be evaluated by using Lemma 3.7 in [28].

Next, we verify that the time step size condition derived from Theorem 1 is optimal; namely, the bound is the maximum time step size that satisfies the maximum principle. Let

be the least upper bound of the time steps that satisfy the maximum principle for the AC equation. To verify $\Delta t_{\max}$ is optimal, let us assume $\Delta t > \Delta t_{\max}$. Then, we have $\Delta t > \Delta t_{\max}$, which implies $0 < \Delta t_{\max}/\Delta t < 1$. Let us consider the following polynomial function on $(0, 1)$:

which is a strictly monotonically increasing function because $f'(\psi) > 0$ on $(0, 1)$. Because $f(0) = 0$, $f(1) = 1$, and $f(\psi)$ is continuous, there exists a unique $0 < \psi < 1$ satisfying $\Delta t_{\max}/\Delta t = f(\psi)$. Using this obtained value of $\psi$, let us consider an initial condition on $\Omega_h$ as follows:

where $1 < p < N_x$ and $1 < q < N_y$. From Eq (2.3), we have

which implies that the maximum principle is not satisfied for any time step sizes larger than the maximum time step size $\Delta t_{\max}$. Therefore, we can conclude that $\Delta t_{\max}$ is the optimal bound.

3.   Computational experiments

In this section, we perform several characteristic computational tests to validate the maximum principle property, motion by mean curvature, and the effect of parameter $\alpha$ on the dynamics of the equation.

3.1.   Convergence test

Before conducting computational experiments, we investigate the accuracy of the numerical scheme. Consider the following initial condition:

in the computational domain $\Omega = (-1, \; 1)\times(-1, \; 1)$. The parameters used to investigate spatial accuracy are $\epsilon = 1/32$, $\Delta t = \epsilon^2(1/512)^2/(2((1/512)^2+2\epsilon^2))$, and $T = 2000\Delta t$. We use the relative error. Table 1 shows the spatial error and convergence rate. When $h = 1/32$, the error is the relative error between using $h = 1/32$ and $h = 1/16$.

The temporal accuracy is also investigated using the same initial condition (3.1). The other parameters are used as follows: $N_x = N_y = 2^{10}$, $\epsilon = h$, $\Delta t_{\max} = \epsilon^2h^2/(2(h^2+2\epsilon^2))$, and $T = 100\Delta t_{\max}$. Table 2 provides the temporal error and convergence rate.

Consequently, it was found that the method has second-order accuracy in space and first-order accuracy in time.

3.2.   Maximum principle

We consider the maximum principle of the fully explicit method with $\Delta t_{\max}$. The random initial condition is given as follows:

where $\mbox{rand}(x, y)$ is a random value between $-1$ and $1$. Here, the parameters used are $\alpha = 1$, $h = 0.02$, $\epsilon = h$, and a final time $T = 2000\Delta t_{\max}$. Figure 1 displays the temporal progressions of the highest and lowest values of $\phi$, using time step sizes $\Delta t = \Delta t_{\max}$ and $\Delta t = 1.21\Delta t_{\max}$. When employing a time step greater than $\Delta t_{\max}$, it is confirmed that the $\phi$ value exceeds the bounds of $-1$ and $1$.

3.3.   Motion by mean curvature

Motion by mean curvature in the AC equation means the evolution of interfaces between different phases. It describes how the interface between two distinct phases evolves over time due to the minimization of the total interface energy. The motion is governed by the mean curvature of the interface, which tends to smooth out irregularities and minimize the interface area. This phenomenon plays a very important role in various physical processes, including phase separation, pattern formation, and material microstructure evolution. We consider the motion by mean curvature using the fully explicit scheme of the AC equation with a high-order polynomial potential. The initial condition on the computational domain $\Omega = (-1, 1)\times (-1, 1)$ is defined as

We used the following parameters: $N_x = N_y = 128$, $\alpha = 3$, $r = 0.5$, $\epsilon = 0.013$, $\Delta t = 0.5h^2\epsilon^2/(h^2+2\epsilon^2)$, and $N_t = 3400$. The analytic radius is defined as $R(t) = \sqrt{r^2-2t}$. Figure 2 shows the temporal evolution of the theoretical radius $R(t)$ and numerical radius with the zero-level filled contours of the numerical solutions at $t = 560\Delta t$, $1600\Delta t$, and $3100\Delta t$. From the computational experiment, it becomes evident that the numerical results demonstrate a strong correspondence with the theoretical radius, indicating a high level of agreement between the two results.

3.4.   Effect of $\alpha$ on the dynamics of phase-field

Next, we consider the effect of $\alpha$ on the dynamics of the phase-field of separated four squares as an initial condition on the computational domain $\phi = (0, 2\pi)\times(0, 2\pi)$ with the following numerical parameters: $N_x = N_y = 256$, $h = 2\pi/256$, and $\Delta t = 2.0 \times 10^{-5}$. The initial condition is given as follows:

Figures 3(a),(b) show the temporal evolution of the zero-level filled contour representing the numerical solution with $\alpha = 1$ and $\alpha = 10$, respectively. We note that when $\alpha = 1$, the proposed model reduces to the conventional AC equation with fourth-order polynomial free energy. In Figure 3(a), the separated four squares merge into one structure. Conversely, for $\alpha = 10$, the four squares shrink individually without merging, as shown in Figure 3(b). Here, $\epsilon = 0.0512$ is used for $\alpha = 1$, and $\epsilon = 0.0085$ is used for $\alpha = 10$.

4.   Conclusions

In this paper, we presented a detailed stability analysis of the fully explicit finite difference scheme to solve the AC equation with a high-order polynomial free energy, which was recently proposed to preserve a more intricate structure of interfaces. To demonstrate the accuracy and optimal estimation of our stability limitations, we conducted a series of numerical experiments using the derived time step formulas that ensure the maximum principle. The obtained time step estimation, which guarantees stability, can be useful for applying the AC equation in various practical applications, such as data classification, image inpainting, and image segmentation. In data classification, for instance, the equation's stability ensures robust and accurate categorization of complex datasets. In the context of image inpainting, where missing or damaged portions of images are reconstructed, the stability provided by our estimation enhances the fidelity of the inpainting process. Similarly, in image segmentation tasks aimed at partitioning images into meaningful regions, the reliability of the AC equation enables precise description of boundaries. In conclusion, the stability contributions outlined in this paper not only advance the theoretical understanding of stability in finite difference schemes but also provide the methodology for practical implementations across a range of applications.

Appendix

The following MATLAB code is used for the numerical simulation as shown in Figure 2 and can be accessed from the corresponding author's webpage: https://mathematicians.korea.ac.kr/cfdkim/open-source-codes/.

Author contributions

Jaeyong Choi: Formal Analysis, Validation, Investigation, Writing-original draft, Writing-review & editing; Seokjun Ham: Software; Formal Analysis, Writing-original draft, Writing-review & editing; Soobin Kwak: Data curation, Software, Visualization, Writing-original draft; YoungJin Hwang: Formal Analysis, validation, Writing-original draft, Writing-review & editing; Junseok kim: Conceptualization, Supervision, Validation, Methodology, Writing-original draft, Writing-review & editing. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by BK21 FOUR.

The authors are grateful to the reviewers for their constructive feedback on this article.

Conflict of interest

The authors declare that there is no conflict of interest.