Maximum principle preserving the unconditionally stable method for the Allen–Cahn equation with a high-order potential

1.   Introduction

We consider a simple and unconditionally stable numerical method for the Allen–Cahn (AC) equation with a high-order polynomial potential ^[1]:

where $\phi({\bf x}, t)$ is the phase-field function in a two-dimensional (2D) domain $\Omega \subset {\mathbb R}^2$,

is the high-order free energy potential ($\alpha$ is an even integer, see Figure 1, and $\epsilon$ is an interfacial thickness-related parameter). When $\alpha = 2$, $F_2(\phi) = 0.25(\phi^2 - 1)^2$ ^[2]. The primary purpose of using the high-order polynomial potential is to resolve the fine structure of phase-field interfaces ^[3].

The AC equation has been used as a fundamental equation in the mathematical modeling of phase transitions across various scientific disciplines such as materials science ^[4], image processing ^[5], data assimilation ^[6], phase-field modeling without an artificial curvature effect for the dendritic growth simulation ^[7], topology optimization ^[8], and fluid dynamics ^[9,10]. Its ability to describe the evolution of phase boundaries makes it a valuable tool for investigating scientific processes such as coarsening, phase separation, and domain wall dynamics. A key mathematical property of the AC equation is the maximum principle ^[11], which guarantees that the solution remains bounded within a physically meaningful range. This property is very important for maintaining the accuracy and physical relevance of simulations. However, preserving the maximum principle in numerical approximations remains a challenging task, especially with large time steps. Numerous computational schemes have been developed to solve the AC equation, aiming to achieve stability, efficiency, and accuracy. These methods include the reduced-order modified finite difference method (FDM) ^[12], a fourth-order compact FDM ^[13], Lie–Trotter's and Strang's splitting methods ^[14], spectral methods ^[15], energy-stable methods ^[16], and implicit-explicit (IMEX) schemes ^[17]. Yang et al. ^[18] introduced a discontinuous Galerkin method for the stochastic AC equation with multiplicative noise, using a symmetric interior penalty approach for spatial discretization and an implicit Euler method for temporal discretization, and demonstrated optimal convergence rates in both space and time through rigorous analysis and numerical experiments. Among these, operator splitting methods simplify implementation and improve computational efficiency. Despite these advantages, it is difficult to preserve the maximum principle with high-order polynomial potentials and arbitrarily large time steps. This limitation can lead to numerical efficiency, which requires small time steps and more computational times.

In this work, we propose a novel numerical scheme that overcomes these limitations by guaranteeing the preservation of the maximum principle while achieving unconditional stability. Our approach is specifically designed to resolve the challenges associated with modeling phase transitions using the AC equation with high-order polynomial potentials. By preserving the maximum principle, the proposed method improves both the physical consistency of simulations and the numerical stability under large time steps. The core of our method is an operator splitting technique, which decomposes the AC equation into simpler subproblems that can be solved sequentially. This decomposition results in one nonlinear equation and several linear diffusion equations. The nonlinear equation is solved using the frozen coefficient method, a strategy that approximates variable coefficients by treating them as constants within small regions. This simplification significantly reduces computational complexity while keeping the stability and accuracy of the computational solution. For the linear diffusion equations, we use a fully implicit FDM. The implicit nature of this scheme guarantees unconditional stability, even with large time steps, which makes it highly efficient for long time simulations. To validate the performance of the developed method, we conduct a series of computational tests. These experiments demonstrate the scheme's ability to keep unconditional stability and preserve the maximum principle for various test cases. The computational results highlight the efficiency and reliability of our method, which show its potential for simulating complex phase transitions with high accuracy. Our computational tests indicate that the proposed scheme effectively handles noisy data, which provides a robust methodology for data classification applications.

The potential applications of the proposed numerical method to other problems are as follows: two-phase porous media flow ^[19], the shape optimization problem ^[20], the image classification method combining phase-field and deep neural networks (DNNs) ^[21], thermal-fluid topology optimization with unconditional energy stability and high-order accuracy using the phase-field equation ^[22], and the modified diffuse interface fluid model and its consistent energy-stable computation in arbitrary domains ^[23], which will demonstrate the versatility and practical relevance of the method.

The rest of this paper is structured as follows. In Section 2, we present the computational method for the AC equation. Section 3 presents computational tests. In Section 4, we provide the conclusions.

2.   Computational method

An unconditionally stable computational scheme is proposed for the AC equation with a high-order potential. Let $\Omega = (a, b) \times (c, d)$ and $\Omega_h = \{ (x_i, y_j) | x_i = a+h(i-0.5), \; 1\le i \le N_x \; \mbox{ and }\; y_j = c+h(j-0.5), \; 1 \le j \le N_y \}$ is its discrete domain. Let $\phi^n_{ij} = \phi(x_i, y_j, t_n)$, where $t_n = n \Delta t$ and $\Delta t$ is the time step. We use the operator splitting technique, which is a computational scheme for solving complex differential equations by decomposing them into simpler subproblems. Each subproblem, which is easier to solve, is solved separately; see ^[24] and the references therein. The solutions to these subproblems become the approximate solution to the full problem:

First, we consider the following nonlinear equation in ${Step \; 1}$:

with the initial conditions $\phi_{ij}^n$ for $1 \le i \le N_x$ and $1 \le j \le N_y$. We rewrite Eq (2.4) using the definition of $F_\alpha(\phi({\bf x}, t))$ as follows:

where

To solve Eq (2.5) efficiently, we apply the frozen coefficient method ^[25]. This approximation technique is used to simplify the solution of differential equations with variable coefficients. By freezing or holding the coefficients constant at specific points, the method transforms a complex equation into one with constant coefficients, which is easier to solve. This approach is widely used in numerical analysis to approximate solutions for nonlinear systems and allows for more efficient computations. Then, by applying a frozen coefficient method, we obtain a closed-form solution, for $1 \le i \le N_x$ and $1 \le j \le N_y$:

Here,

where we included a small positive number in the calculation to prevent division by zero. In this study, $tol = 1.0e$-$15$ is used. Note that if $\alpha = 2$, then we have $A_{ij}^n = 1/\epsilon^2$ and Eq (2.7) becomes the conventional form for the fourth-order polynomial potential. The solution for ${Step \; 2}$ is obtained using a fully implicit FDM with the initial condition $\psi_{ij}^*$. For a fixed index $j = 1, \ldots, N_y$, we solve the following system of discrete equations for $i = 1, \ldots, N_x$ by using the Thomas algorithm:

where we use

Similarly, the solution for ${Step \; 3}$ is obtained using a fully implicit FDM with the initial condition $\psi_{ij}^{**}$. For a fixed index $i = 1, \ldots, N_x$, we solve the following system of discrete equations for $j = 1, \ldots, N_y$ by using the Thomas algorithm:

where we use

An efficient Fourier spectral method can be used for both ${Steps \; 2}$ and ${3}$, which was used in studying a maximum principle of the Fourier spectral method for diffusion equations ^[26]. The AC equation is widely recognized for its fundamental property of satisfying the maximum principle, which is essential for evaluating the stability of its solutions. A key challenge in numerical simulations is determining whether these properties can be preserved at the discrete level ^[27]. The proposed numerical method satisfies the maximum principle: if $\|\phi^0 \|_\infty \le 1$, then $\|\phi^k \|_\infty \le 1$ for any $k\ge1$. More specifically, assume $\|\phi^n \|_\infty \le 1$ for any $n$, and then from Eq (2.7), we have

It is widely recognized that the fully implicit FDMs for Eqs (2.9) and (2.11) also satisfy the maximum principle. In this study, we focused on a square domain. However, the proposed algorithm can be applied to arbitrary domains by using the finite element method, which is effective for performing computations in complex domains ^[28,29]. Another possible approach is to use an alternating direction explicit (ADE) scheme on a complex Cartesian grid ^[30], which was used in studying the semi-automatic fingerprint image restoration technique using a partial differential equation.

Note that we can straightforwardly extend the proposed method to a three-dimensional numerical scheme. Specifically, we have the following steps:

We solve Eq (2.14) as

where $\phi^n_{ijk} = \phi(x_i, y_j, z_k, t_n)$ and

Furthermore, Eqs (2.15)–(2.17) can be solved using a method similar to that used in the two-dimensional case.

Although the proposed numerical algorithm can be solved using other established numerical approaches, such as Fourier-spectral methods and multigrid techniques ^[31], our focus was on using the alternating implicit method. Each method has its own advantages and disadvantages, depending on factors such as the shape of the application domain, the required accuracy, and the computational times. For example, while the proposed method is less effective at handing complex domains, it preserves the maximal principle property.

3.   Numerical experiments

3.1.   Stability test

To evaluate the robustness and reliability of the proposed numerical method, a stability test is conducted. This test examines the behavior of the numerical solution under varying conditions, including random initial conditions and different parameter values, to confirm the method's stability over time. We consider a random initial condition on $\Omega = (-1, 1)\times(-1, 1)$:

where ${\rm{rand}}(x, y)$ is a random number between $-1$ and $1$.

Figure 2 shows the temporal evolution of the solution $\phi(x, y, t)$ at $t = 3\Delta t$, $10\Delta t$, and $50\Delta t$, with $h = 0.02$, $\Delta t = h^2$, and $\epsilon = h/\alpha$. The top and bottom rows correspond to $\alpha = 2$ and $\alpha = 10$, respectively. For $\alpha = 2$, the solution retains its randomness, with small-scale oscillations persisting over time. In contrast, for $\alpha = 10$, the solution evolves into larger-scale structures with noticeable changes, particularly at $t = 50\Delta t$, which indicates possible instability. The computational results demonstrate that the numerical solutions are stable for large values of $\alpha$, even with a relatively large time step.

3.2.   Effects of $\alpha$ and $\epsilon$

Next, we investigate the influence of both $\alpha$ and $\epsilon$ on the thickness of the interfacial transition layers. Let us consider an initial condition on a one-dimensional domain $\Omega = (-1, 1)$:

Figure 3 illustrates the effects of $\alpha$ and $\epsilon$ on the equilibrium interfacial transition profile. The equilibrium state is defined as where the numerical solution in which the $l_2$-norm of the difference between consecutive numerical solutions is less than a specified tolerance, $tol = 10^{-10}$. In this test, we used $N_x = 50$, $h = 2/N_x$, and $\Delta t = 0.1h^2$. Figure 3(a) presents the equilibrium profiles for $\alpha = 2$, $4$, and $10$ with $\epsilon$ fixed at $5h$. Larger values of $\alpha$ result in sharper transitions and show its influence on the steepness control. Figure 3(b) shows the equilibrium profiles for $\epsilon = h$, $5h$, and $10h$ with $\alpha$ fixed at $10$. Smaller values of $\epsilon$ result in narrower transitions, whereas larger values of $\epsilon$ produce smoother and broader profiles.

3.3.   Classification with noisy data

Next, we consider a classification with noisy data and an initial condition on $\Omega = (0, 1) \times (0, 1)$ is as follows:

As shown in Figure 4(a), the initial condition, Eq (3.3), shows a noisy field with substantial fluctuations around $\phi = 0$. Two prominent "square seeds" are visible as sharp pillars at specific locations in the domain, where $\phi$ takes the values of $1$ and $-1$. These square seeds are expected to grow or stabilize into distinct regions over time while the surrounding noise diminishes. This configuration highlights the competition between the seeds and the noisy background, which should result in well-separated stable regions as the simulation progresses.

Figure 4(a)–(e) show the temporal evolution of the solution $\phi(x, y, t)$ at $t = 3\Delta t$, $10\Delta t$, and $50\Delta t$, respectively, with $h = 0.02$, $\Delta t = h^2$, and $\epsilon = h/\alpha$. Initially, as shown in Figure 4(a), the solution exhibits high noise from the initial condition. Over time, as seen in Figure 4(b)–(e), the noise diminishes, and distinct regions of $\phi = 1$ and $\phi = -1$ emerge. By $t = 120\Delta t$, as shown in Figure 4(e), the solution stabilizes, which reflects a clear separation between regions. Furthermore, Figure 4(f)–(h) display the computational results at times $t = 20 \Delta t,$ $40\Delta t$, and $120 \Delta t$, respectively, with a larger time step, $\Delta t = 10h^2$. This test highlights the smoothing and classification capability of the high-order AC equation when applied to noisy data, which effectively produces distinct regions over time. The numerical results demonstrated that the computational solutions satisfy the maximum principle. Specifically, if the initial condition satisfies $\|\phi^0 \|_\infty \le 1$, then the numerical solutions ensure that $\|\phi^k \|_\infty \le 1$ for any time step $k \ge 1$, which maintains consistency with the theoretical bounds throughout the computation.

Table 1 lists the computational cost measured in seconds for different the grid resolutions $N_x \times N_y$, while keeping the other parameter values unchanged. The computational cost increases linearly with the number of grid points, as shown in the table. Doubling the grid resolution in each direction results in an approximately fourfold increase in computational cost, consistent with the quadratic relationship between grid resolution ($N_x \times N_y$) and grid points, confirming the algorithm's efficiency.

As a final test, let us consider the effect of the parameter $\alpha$ on the dynamics of the phase-field using the same initial condition specified in Eq (3.3) and other parameters except $\alpha$.

Figure 5 illustrates the evolution of a scalar field $\phi$ over a 2D domain at $t = 40\Delta t$ for different parameter values of $\alpha = 2, 4, 6, 8, 10,$ and $12$. The results in Figure 5(a), (b) indicate that the initial noise has developed from the initial condition. This highlights the sensitivity to the initial conditions, where small perturbations evolve into more complex structures over time. In contrast, the results in Figure 5(c)–(f) demonstrate the formation of clear separation. As $\alpha$ increases ($\alpha = 6, 8, 10, 12$), the system exhibits more stable dynamics, with sharp interfaces emerging and high-frequency noise being suppressed.

4.   Conclusions and future works

We have developed a maximum principle-preserving and unconditionally stable computational method for the high-order AC equation. By using an operator splitting method, the numerical algorithm effectively decomposes the problem into simpler subproblems, which improves both computational efficiency and stability for time-dependent partial differential equations. A frozen coefficient approach is used to handle the nonlinear component, and a fully implicit FDM guarantees unconditional stability for the linear diffusion equations. The proposed method successfully satisfies the maximum principle, an important property for preserving the physical fidelity of phase transition simulations, even under challenging conditions such as high-order potentials. Numerical experiments have confirmed the robustness and efficiency of the approach, and demonstrated its ability to accurately model phase separation phenomena, particularly under noisy data conditions. These computational results highlighted the potential of the method for advancing the numerical analysis of complex phase transition problems. Future work will focus on extending the method to three-dimensional simulations ^[32], the time-fractional AC equation ^[33], image processing applications ^[34], and the AC equation with the logarithmic Flory–Huggins potential ^[35], which will broaden its applicability to modeling complex phase transitions. Furthermore, the proposed algorithm can be integrated with the finite element method, which is highly versatile for solving the high-order AC equation on arbitrary domains ^[36].

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 the National Research Foundation(NRF), Korea, under project BK21 FOUR. The author expresses sincere gratitude to the reviewers for their insightful comments and constructive suggestions, which greatly contributed to improving the quality of this manuscript.

Conflict of interest

The author declares there is no conflict of interest.