A new study on the Newell-Whitehead-Segel equation with Caputo-Fabrizio fractional derivative

1.   Introduction

In this article, we present an explicit spatially fourth-order accurate compact method for the Allen-Cahn (AC) equation ^[1,2] in one-, two-, and three-dimensional spaces. The AC equation is as follows:

where the order parameter $\phi({\bf x}, t)$ is the difference of the local concentrations of the two components in a domain $\Omega$, $F(\phi) = 0.25(\phi^2-1)^2$ is a double well free energy density, and $\epsilon$ is a small parameter. The homogeneous Neumann boundary condition is assumed:

where $\partial/ \partial { \rm \bf n}$ is the normal derivative on the boundary of domain. The AC equation (1.1) is $L^2$-gradient flow of the Ginzburg–Landau functional:

Developing high-order accurate numerical schemes is important, therefore, it has been widely studied. In studying high-order accurate numerical schemes in space, wide stencils are required. In this case, we should eliminate the inefficiency caused by calculating high-order derivatives with respect to space, such as the Laplace operator. Therefore, the compact scheme has been actively developed for high-order accurate schemes with relatively less wide stencils.

Abide ^[3] presented finite difference preconditioning for compact scheme discretizations of the Poisson equation with variable coefficients. Li and Liao ^[4] developed an explicit high-order compact finite difference method (FDM) to solve three-dimensional (3D) acoustic wave equations with spatially variable acoustic velocity. In ^[5], based on FDM, a high-order accurate compact method for the heat equation was studied on the inverse Lax-Wendroff boundary treatment. In ^[6], the authors developed a highly accurate local one-dimensional (1D) explicit compact method for the acoustic wave equation in two-dimensional (2D) space. Patel and Mehra ^[7] studied the fourth-order compact method for space fractional advection-diffusion reaction equations using the alternating direction implicit scheme. An explicit fourth-order compact scheme was developed for heat transfer of boundary layer flow ^[8], which has the advantage of finding solutions of nonlinear and linear convection-diffusion type equations that cannot be solved with the existing implicit compact scheme. Qiu et al. ^[9] developed a general conservative eighth-order compact finite difference method to solve the coupled Schrödinger-KdV equations, this method is decoupled and preserves several physical invariants in a discrete sense. Elmahdi and Huang ^[10] presented linearized finite difference and fourth-order compact finite difference schemes for the time fractional nonlinear diffusion-wave equations based on their equivalent partial integro-differential equations. Abdi et al. ^[11] presented high-order compact finite difference schemes for the fractional Black-Scholes equation in the case of European option pricing. Using implicit-type schemes, spatially fourth-order accurate compact schemes were presented in ^[12,13,14]. There are some studies for solving partial differential equation (PDE), including the Laplace terms using a fourth-order accurate method on hexagonal grids ^{[15,16,17,18]}. Zhai et al. ^[12] presented a linearized, temporally second-order, spatially fourth-order accurate compact method for the 3D AC equations with different boundary conditions. Long et al. ^[13] described the unconditional stable linear FDM for the 3D AC equation, combined with backward differentiation and fourth-order compact schemes. Bo et al. ^[14] theoretically proved the discrete maximum principle and energy stability of the compact difference scheme for AC equation in the 2D space and presented numerical examples. In ^[19], the authors presented a fourth-order FDM for the AC equation with a stabilized term which has an effect on structure-preserving.

Moreover, the compact finite difference scheme has been applied for solving the Cahn-Hilliard (CH) equation ^[20,21]:

which is a phase-field model, $H^{-1}$-gradient flow of the energy functional (1.2), modeling the phase separation of spinodal decomposition in a binary alloy. Li et al. ^[22] analyzed a three-level linearized compact method for the CH equation as proving the unique solvability and unconditional convergence of the numerical solution. Ju et al. ^[23] proposed a fast and stabilized compact exponential time differencing multi-step method for the CH equation. Lee ^[24] presented a fourth-order spatial and second-order temporal accurate, unconditionally stable compact FDM for the CH equation. Lee and Shin ^[25] presented a compact scheme for solving the CH equation with a periodic boundary condition. In ^[26], a fourth-order compact method for the pure stream function formulation of 3D unsteady incompressible Navier-Stokes equation was presented using the Crank-Nicolson method for time discretization.

Many studies for the compact scheme have adopted implicit-type schemes to numerically solve various equations. Implicit-type schemes have the advantage that they can solve the equations efficiently by using large time-step sizes compared to explicit-type schemes. On the other hand, in the case of a nonlinear partial differential equation, we need to use efficient numerical solvers such as multigrid methods. In addition, highly accurate results require small time-step sizes. For example, when we solve the CH equation which is a fourth-order nonlinear PDE, an explicit method imposes a severe constraint on time-step size for stability. A very small time-step size makes it almost impossible to compute the fine-grid solution. In this case, the multigrid method is useful. However, in this study, we focus on the fourth-order accurate compact method for the AC equation, which is a second-order nonlinear parabolic PDE. When we solve the AC equation using an explicit scheme, the linear term restricts the time-step size as $\Delta t < 0.5h^2/d$, where $h$ is a mesh size and $d$ is the spatial dimension ^[27,28,29]. The constraint on the time-step size is not too severe, therefore, the fourth-order compact explicit scheme increases the accuracy without greatly losing efficiency.

The main purpose of this study is to present a fully explicit spatially fourth-order accurate compact scheme for the AC equation in 1D, 2D, and 3D spaces. To the authors' knowledge, this is the first study of a fully explicit spatially fourth-order accurate compact scheme for the AC equation. The advantages of this study focusing on space are the accuracy and efficiency of the numerical solution, and the simplicity of implementation. In the follow-up study, we will improve the proposed scheme to be higher-order accuracy in time.

The layout of this paper is as follows. In Section 2, it is presented that the description of the fourth-order compact FDM for the AC equation. In Section 3, we perform various computational experiments to demonstrate the basic properties and advantages of the proposed method. Conclusions are made in Section 4.

2.   Numerical method

2.1.   Discretization

Let $N_x$, $N_y$, and $N_z$ be the number of grid points in the $x$, $y$, and $z$ direction, respectively, $\Delta t = T/N_t$ be the time-step size, $T$ be the final time, and $N_t$ be the total number of time steps.

2.1.1.   One-dimensional space

In a 1D domain $\Omega = (a_x, b_x)$, let $h = (b_x-a_x)/N_x$ be the uniform mesh size. A discrete numerical domain is denoted by $\Omega_h = \{x_i\ :\ x_i = a_x+(i-0.5)h, \ 1\leq i\leq N_x\}$. Let $\phi_{i}^n$ be the approximation of $\phi(x_i, n\Delta t)$. The homogeneous Neumann boundary conditions for $\phi$ are given as follows: $\phi_{0} = \phi_{1}, \; \phi_{-1} = \phi_{2}, \; \phi_{N_x+1} = \phi_{N_x}, \; \phi_{N_x+2} = \phi_{N_x-1}.$ The discrete differentiation operator is $D_x \phi_{i+ \frac{1}{2}} = (\phi_{i+1} - \phi_i)/h.$

2.1.2.   Two-dimensional space

We discretize the AC equation in a 2D domain $\Omega = (a_x, b_x)\times(a_y, b_y)$. Let $h = (b_x-a_x)/N_x = (b_y-a_y)/N_y$ be the uniform mesh size. Let $\Omega_h = \{(x_i, y_j)\ :\ x_i = a_x+(i-0.5)h, \ y_j = a_y+(j-0.5)h, \ 1\leq i\leq N_x, \ 1\leq j\leq N_y\}$ be a discrete numerical domain. Let $\phi_{ij}^n$ be the approximation of $\phi(x_i, y_j, n\Delta t)$. The homogeneous Neumann boundary conditions for $\phi$ are given as follows:

The discrete differentiation, gradient, and divergence operators are

$\nabla_c \phi_{ij} = \left(D_x \phi_{i+ \frac{1}{2}, j}, D_y \phi_{i, j+ \frac{1}{2}}\right)$, and $\nabla_d \cdot (u, v)_{ij} = (u_{i+ \frac{1}{2}, j}-u_{i- \frac{1}{2}, j})/h+ (v_{i, j+ \frac{1}{2}} -v_{i, j- \frac{1}{2}})/h$, respectively. The discrete $l_2$-inner products and norms are defined as

We define the discrete total energy functional by

where ${\bf 1}$ is a vector with all entries are 1.

2.1.3.   Three-dimensional space

We discretize the AC equation in a 3D domain $\Omega = (a_x, b_x)\times(a_y, b_y)\times(a_z, b_z)$. Let $h = (b_x-a_x)/N_x = (b_y-a_y)/N_y = (b_z-a_z)/N_z$ and $\Omega_h = \{(x_i, y_j, z_k)\ :\ x_i = a_x+(i-0.5)h, \ y_j = a_y+(j-0.5)h, \ z_k = a_z+(k-0.5)h, \ 1\leq i\leq N_x, \ 1\leq j\leq N_y, \ 1\leq k\leq N_z\}$. Let $\phi_{ijk}^n$ be the approximation of $\phi(x_i, y_j, z_k, n\Delta t)$. The homogeneous Neumann boundary conditions for $\phi$ are given as follows:

The discrete differentiation operators are

and then $D_y \phi_{i, j+ \frac{1}{2}, k}$ and $D_z \phi_{ij, k+ \frac{1}{2}}$ are similarly defined above. We then define the discrete $l_2$-inner products as

2.2.   Compact finite difference scheme

2.2.1.   One-dimensional space

In a 1D space, a compact Laplace operator $\Delta_c$ is defined as follows:

By the Taylor expansion, we can obtain

Then, we can derive the following equation from Eq (2.8) by replacing $\Delta x$ with $\pm h$.

By dividing both sides of Eq (2.9) by $h^2$, we obtain

where $\Delta^2 \phi = \Delta(\Delta \phi)$ is the biharmonic operator. Note that the standard fourth-order 5-point Laplace operator $\Delta_s$ is defined as

From Eq (2.8), we can derive $\Delta_s \phi_{i} = \Delta \phi(x_i) + O(h^4)$ in the similar manner. Using Eq (2.10), we discretize the governing Eq (1.1) in a 1D domain as follows:

Therefore, up to $O(\Delta t)$ and $O(h^4)$, we have

for $n \geq 2$. Given $\phi^0$, $\phi^1$ is computed by using the standard fourth-order Laplacian $\Delta_s$.

2.2.2.   Two-dimensional space

In a 2D space, the compact Laplace operator $\Delta_c$ ^[30] is defined as

By the Taylor series in two variables, we can obtain

By replacing $\Delta x$ and $\Delta y$ with different values $\pm h$, we get

From Eqs (2.13) and (2.14), we have

Finally, we have

where $\Delta^2 \phi = \Delta (\Delta \phi)$ is the biharmonic operator. Note that another standard fourth-order 9-point Laplace operator $\Delta_s$ is defined as

In a similar way, $\Delta_s \phi_{ij} = \Delta \phi(x_i, y_j) +O(h^4)$. Figure 1 shows the nodes used by the compact and standard Laplace operators at a point, respectively.

Using Eq (2.15), similar to Eq (2.11), we discretize the governing Eq (1.1). Therefore, up to $O(\Delta t)$ and $O(h^4)$, we have

for $n \geq 2$. Given $\phi^0$, $\phi^1$ is computed by using the standard fourth-order Laplacian $\Delta_s$.

2.2.3.   Three-dimensional space

In a 3D space, a compact Laplace operator $\Delta_c$ ^[31] is defined as

By the Taylor expansion in three variables, we can obtain

Then, we can derive the following equations from Eq (2.17) by replacing $\Delta x$, $\Delta y$, and $\Delta z$ with different values $\pm h$.

Multiplying both sides of Eqs (2.18)–(2.20) by $14$, $3$, and $1$, respectively, and adding the equations, we have

By dividing both the sides of Eq (2.21) by $30h^2$, we obtain

Note that the 3D standard fourth-order 13-point Laplace operator $\Delta_s$ is defined as

In a similar way, $\Delta_s \phi_{ijk} = \Delta \phi(x_i, y_j, z_k) +O(h^4)$. Figure 2 shows the grid point usages of standard $\Delta_h$ and compact $\Delta_h$ in a 3D space.

Using Eq (2.15), we discretize the governing Eq (1.1). Therefore, up to $O(\Delta t)$ and $O(h^4)$, we have

for $n \geq 2$. Given $\phi^0$, $\phi^1$ is computed by using the standard fourth-order Laplacian $\Delta_s$.

3.   Numerical results

We demonstrate the representative properties of the explicit compact AC equation: phase separation, fourth-order accuracy in space, non-increase of total energy, and motion by mean curvature. We show the efficiency of our method using numerical comparison and measuring computational time. In addition, we perform the computational experiments for a non-convex initial condition with and without an obstacle to highlight the advantage of the proposed scheme. It can easily handle the computation on arbitrary complex domains.

3.1.   Phase separation

We consider one of the most basic computational simulations, that is, the dynamics of phase separation. First, we conduct the numerical experiment in a 1D computational domain $\Omega = (0, \; 1)$ with the initial condition given as

where ${\rm rand}$ is a random number in $[-1, \; 1]$. To observe the equilibrium profile, we cease the iterative computation when the following criterion is met:

where $n$ is a positive integer, and $\phi^n$ is called a numerical equilibrium solution. The parameters are taken as $N_x = 100, \; h = 1/N_x, \; \Delta t = 0.1h^2,$ and $\epsilon = 4h/(2\sqrt{2}\tanh^{-1}(0.9))$. Figure 3 indicates the initial condition ($t = 0$), numerical solution ($t = 30\Delta t$), and numerical equilibrium solution ($t = 2242\Delta t$). From this result, we can also numerically confirm that the numerical solutions are bounded between $-1$ and $1$.

We next observe the phase separation in a 2D domain $\Omega = (0, \; 1)^2$. The initial condition is given as

The parameters are taken as $N_x = N_y = 100, \; h = 1/N_x, \; \Delta t = 0.1h^2,$ and $\epsilon = 4h/(2\sqrt{2}\tanh^{-1}(0.9))$. Figure 4 indicates the initial condition ($t = 0$) and numerical solutions ($t = 30\Delta t, \; 100\Delta t$). From this result, we can also numerically confirm that the numerical solutions are bounded between $-1$ and $1$.

3.2.   Convergence test

We investigate a quantitative estimate of the rate of convergence. Note that the discrete $l_2$-norms of numerical solutions is defined as

in 1D and 2D domains, respectively. First, in a 1D domain $\Omega = (0, \; 2)$, we consider the traveling wave solution of the AC equation ^[32,33]:

where $s = 3\epsilon/\sqrt{2}$. At $t = 0$, the initial profile is

Figure 5 shows the initial, numerical, and analytic profiles. Here, the parameters used are $N_x = 100, \; h = 2/N_x, \; \epsilon = 8h/(2\sqrt{2}\tanh^{-1}(0.9)), \; \Delta t = 0.1h^2,$ and the final $T = 120\Delta t$.

We perform a convergence test with the initial condition (3.4) on a set of increasingly finer mesh grid sizes. The computational solutions are computed on the uniform grids and time steps, $h = 2/N_x$ for $N_x = 40, \; 80, \mbox{ and } 160$. The fixed time-step size is $\Delta t = 5.0e$-$9$, the final time is $T = 50000\Delta t$, and $\epsilon = 0.025$. The convergence rate is defined as the ratio of successive errors, that is, $\log_2 (||e_{N_x}||_2 / ||e_{2N_x}||_2)$, where the discrete $l_2$-norm error is defined as $||e_{N_x}||_2 = \sqrt{ \frac{1}{N_x}\sum_{i = 1}^{N_x} (\phi(x_i, n\Delta t)-\phi_{i}^n) }$ and $||e_{2N_x}||_2 = \sqrt{ \frac{1}{2N_x}\sum_{i = 1}^{2N_x} (\phi(x_i, n\Delta t)-\phi_{i}^n) }$. Table 1 lists the errors and rates of convergence and demonstrates that the proposed scheme is spatially fourth-order accurate.

In a 2D domain $\Omega = (0, \; 2)^2$, we consider the following solution:

The parameters are taken as $N = N_x = N_y = 100, \; h = 2/N, \; \Delta t = 0.1h^2,$ $\epsilon = 8h/(2\sqrt{2}\tanh^{-1}(0.9)),$ and the final time $T = 90\Delta t$. In this test, the Dirichlet boundary condition is used and the values of order parameter on the boundary are set as the analytic solutions. Figures 6(a) and 6(b) show the initial profile at $t = 0$ and the numerical solution at $t = 90\Delta t$, respectively. Figure 6(c) shows the contours of the initial, numerical, and analytic solutions at the $\phi = 0.5$ level.

We conduct a convergence test with the initial condition (3.5) on a set of increasingly finer mesh grid sizes. The computational solutions are computed using $h = 2/N$ for $N = 40, \; 80,160$, the fixed time-step size $\Delta t = 1.0e$-$9$, the final time $T = 250000\Delta t$, and $\epsilon = 0.025$. Table 2 lists the errors and rates of convergence. Here, the discrete $l_2$-norm error in a 2D domain is defined as $||e_{N}||_2 = \sqrt{ \frac{1}{N^2}\sum_{i = 1}^{N}\sum_{j = 1}^{N} (\phi(x_i, y_j, n\Delta t)-\phi_{ij}^n) }$. The results in the 2D domain suggest that the scheme is indeed fourth-order accurate in space.

3.3.   Equilibrium profile

We verify that a numerical equilibrium solution in which the profile is a hyperbolic tangent is consistent with a theoretical solution when the time is large enough. The initial condition on a domain $\Omega = (-1, 1)$ is

The parameters are taken as $N_x = 200, \; h = 2/N_x, \; \Delta t = 0.1h^2, \; t = 5,$ and $\epsilon = 4h/(2\sqrt{2}\tanh^{-1}(0.9))$. In Figure 7, the numerical equilibrium profile agrees well with the analytic one

Here, the $l_2$-norm error between two equilibrium profiles is less than 1.0e-16.

3.4.   Total energy decrease

We numerically validate that the total energy functional (1.2) decreases over time. In a 2D domain, the discrete total energy is defined as

The initial condition is defined as Eq (3.2) on $\Omega = (0, \; 1)^2$. The parameters are used as follows: $N = 100, \; h = 1/N, \; \Delta t = 0.1h^2, \; T = 2000\Delta t,$ and $\epsilon = 4h/(2\sqrt{2}\tanh^{-1}(0.9))$. Figure 8 shows the time evolution of the discrete total energy and numerical solution. The normalized discrete total energy $\mathcal{E}(\phi^n)/\mathcal{E}(\phi^0)$ is non-increasing over time $t$. The snapshots of the numerical solutions are at $t = 0, \; 0.002, \; 0.006,$ and $0.02$, respectively.

3.5.   Efficiency

In this section, we verify the efficiency of the proposed method. First, we present numerical comparisons between the proposed scheme and previous work ^[13], and compare the performance using their computational time. In ^[13], the authors used the linear multigrid method to solve the AC equation using the fourth-order compact FDM. One V-cycle in the linear multigrid procedure consists of three parts: a pre-smoothing, a coarse grid correction, and a post-smoothing ^[34]. For the sake of convenience, we define the cost of performing one smoothing sweep on the finest grid as a work unit (WU). We assume that each number of pre- and post-smoothing sweeps is equal to the positive integer $\nu$ larger than 1. The smoothing steps mainly adopt the Gauss-Seidel type method, however, we apply the explicit method whose computational cost is less than that. The computational cost of the proposed method is estimated to be less than 2 WU, whereas that of the multigrid method in a 3D domain is

where $N = N_x = N_y = N_z$. In general, one V-cyle is not enough even though it depends on the given tolerance ^[35,36]. Hence, the computational cost of our method is estimated to be less than that of the multigrid method.

We use the initial condition as $\phi(x, y, z, 0) = \sin(2\pi x)\sin(2\pi y)\sin(2\pi z)$ on a domain $\Omega = (0, 1)^3$, the mesh size $N = 32, \; 64, \; 128, \; 256$, and the other parameters as $h = 1/N, \; \Delta t = 0.1h^2, \; T = 10\Delta t,$ and $\epsilon = 0.1$. Our computation is performed using MATLAB 2022b on an Intel Core i9-12900H CPU @ 2.50 GHz with 32 GB of RAM. We measure total CPU time (sec) for 10 time iterations and calculate it divided by the number of time iterations. Table 3 lists the average CPU times. Our method is more than 10 times faster than the implicit-type method.

3.6.   Motion by mean curvature

Motion by mean curvature, the representative property of the AC equation, means that the interface $\Gamma(t)$ of the solution evolves in the direction of the normal velocity proportional to the magnitude of the mean curvature ^[37]. The interface $\Gamma(t)$ is defined by the zero-level set of the solution $\phi$, i.e., $\Gamma(t) = \left\{ {\rm \bf x}\in\Omega \; |\; \phi({\rm \bf x}, t) = 0 \right\}$. The normal velocity $V$ of $\Gamma(t)$ at each point ${\rm \bf x}$ is

where $\kappa_1$ and $\kappa_2$ are the principal curvatures of the interface, and $R_1$ and $R_2$ are the principal radii. In this section, we only consider the case where the two principal radii are the same, therefore, we set the radius as $R(t)$. In a 2D space, there is one principal radius, therefore, $V = -1/R(t)$. In a 3D space, $V = -2/R(t)$. Thus, for the spatial dimension $d$,

Because the normal velocity $V$ is the change in radius with time $t$,

From Eqs (3.7) and (3.8) and the initial radius $R_0 = R(0)$, the analytic radius is

For the details, refer to the references ^[37,38]. We demonstrate the motion by mean curvature in 2D and 3D domains. First, we perform a 2D simulation with the initial condition

in which the shape of zero-level contour is a circle with center $(0.5, \; 0.5)$ and radius $R_0 = 0.35$, on $\Omega = (0, \; 1)^2$. The parameters used are $N = 100, \; h = 1/N, \; \Delta t = 0.1h^2, \; \epsilon = 4h/(2\sqrt{2}\tanh^{-1}{(0.9)}),$ and $T = 6000\Delta t$. Figure 9(a) shows that the given circle shrinks in the direction of the normal vector according to motion by mean curvature. The zero-level contours are drawn every $600\Delta t$ from $t = 0$ to $t = 6000\Delta t$. From the contour spacing, it can be observed that the circle shrinks faster as the radius of the circle becomes smaller. Figure 9(b) indicates that the analytic and numerical radii agree with each other. Figure 9(c)–9(e) show the snapshots of the computational solution at time $t = 0, \; 3600\Delta t,$ and $6000\Delta t$, respectively.

We consider three circles with different centers and radii as shown in Figure 10(a). The initial condition is

Here, the parameters, $N, \; h, \; \epsilon, \; \Delta t$, are used as described above. Figures 10(b) and 10(c) are the snapshots of the computational solutions at $t = 1200\Delta t$ and $t = 3000\Delta t$, respectively. Figure 10(d) shows the zero-level contours every $300\Delta t$ from $t = 0$ to $t = 3000\Delta t$. The interface of the computational solution shrinks with the different velocities depending on the mean curvature of each circle. In other words, because the curvature of a circle with a smaller radius is large, the smaller-radius circle disappears first, and the larger-radius circle shrinks slowly.

Figure 11 shows the snapshots of the initial and numerical solution with the torus-type initial condition:

on $\Omega = (0, \; 1)^2$. Here, the parameters, $N, \; h, \; \epsilon, \; \Delta t$, are used as described above. In Figure 11, the top and bottom rows show mesh plots and zero-level contours of the solution, respectively, where the time $t$ is $0, \; 3000\Delta t,$ and $5000\Delta t$. It is observed that the radius of the inner contour is smaller than that of the outer contour as shown in Figure 11. Hence, with motion by mean curvature, the inner contour decreases faster than the outer contour.

In addition, we observe the results of 2D simulation with the star-shaped initial condition on $\Omega = (0, \; 1)^2$ as shown in Figure 12(a). We use $N = 150, \; h = 1/N, \; \Delta t = 0.1h^2,$ and $\epsilon = 8h/(2\sqrt{2}\tanh^{-1}{(0.9)})$. Even though an initially complex shape is given, the part with the largest mean curvature is smoothed first, and then becomes circular. Figure 12(b) and 12(c) show mesh and contour plots of the computational solution at $t = 0, \; 100\Delta t,$ and $3000\Delta t$.

Next, we consider several computational simulations in a 3D domain to demonstrate motion by mean curvature. First of all, we verify the decreasing trend in the radius of a sphere. We conduct the 3D simulation with the following initial condition:

i.e., a sphere with center $(0.5, \; 0.5, \; 0.5)$ and radius $R_0 = 0.35$, on the computational domain $\Omega = (0, \; 1)^3$. The parameters used are $N = N_x = N_y = N_z = 80, \; h = 1/N, \; \Delta t = 0.1h^2,$ and $\epsilon = 4h/(2\sqrt{2}\tanh^{-1}{(0.9)})$. Figure 13 shows that the given sphere shrinks in the direction of the normal vector according to motion by mean curvature.

Figure 14 shows the temporal evolution of a torus in a domain $\Omega = (0, \; 1)\times(0, \; 1)\times(0, 0.5)$. The initial torus with center (0.5, 0.5, 0.25), large radius $R_1 = 0.25,$ and small radius $R_2 = 0.12$ is defined as:

The parameters are used as $N_x = N_y = 100, \; N_z = 50, \; h = 0.01,$ $\epsilon = 8h/(2\sqrt{2}\tanh^{-1}{(0.9)}), \; \Delta t = 0.1h^2,$ and $T = 700\Delta t$.

On a computational domain $\Omega = (-1, \; 1)^3$, we consider a sphere of center $(0, \; 0, \; 0)$ and radius $R_0$ perturbed by a spherical harmonic $Y_{l, m}$ ^[39,40]:

where $A$ denotes amplitude, $\theta$ is the polar angle (measured from $z$-axis), and $\varphi$ is the azimuthal angle (measured counterclockwise from the $x$-axis to the $y$-axis). Here, $l = 10, \; m = 7, \; R_0 = 0.7,$ and $A = 0.7$. The parameters are used as $N = 100, \; h = 2/N,$ $\epsilon = 4h/(2\sqrt{2}\tanh^{-1}{(0.9)}), \; \Delta t = 0.1h^2,$ and $T = 1000\Delta t$. Figure 15 shows the snapshots at $t = 0, \; 200\Delta t,$ and $1000\Delta t$. From the results, the temporal evolution works according to motion by mean curvature.

3.7.   Temporal evolution with obstacles

The implicit-type methods, compared to the explicit method, can use a relatively large $\Delta t$. However, it is difficult to implement their algorithms for the computational simulation. For example, the multigrid method is one of the most famous algorithms for solving implicitly discretized PDEs. If the computational domain is complex or there are obstacles that interfere with the evolution of order parameter $\phi$, it gets more complicated. In contrast, our proposed explicit compact method can obtain the spatially fourth-order accurate solution for the AC equation efficiently even in complex domains or in the presence of obstacles. To demonstrate these advantages, we perform the computational simulations considering obstacles that are solid and fixed structures in 2D and 3D domains.

We use the non-convex initial condition and the rectangular bar-shaped obstacle in the $2$D domain as shown in Figure 16(a). In Figure 16, from top to bottom, each row represents without the obstacle, with the obstacle, and overlapping the two cases. Here, the solid line and the purple area are the zero-level set of the solutions with and without the obstacle, respectively, and the green area inside the dotted line is the obstacle. The parameters used for the test are $N = 256$, $h = 1/N$, $\Delta t = 0.1h^2$, $T = 39000\Delta t$, and $\epsilon = 4h/(2\sqrt{2}\tanh^{-1}(0.9))$ in the domain $\Omega = (0, \; 1)^2$. The first row of Figure 16 shows the temporal evolution by mean curvature flow in the case without the obstacle as the numerical simulations presented in previous subsections. To perform the simulations with the obstacle, we perform iterative calculations for numerical solutions only at grid points in the whole domain except for obstacles. For the numerical solution near the obstacle, we always use the fixed solution, $\phi = -1$, at points inside the obstacle. The second row of Figure 16 shows the temporal evolution of the zero-level set of numerical solutions with the obstacle. The mean curvature of the inward curves near the obstacle becomes larger than in the case without the obstacle. Therefore, the part with larger mean curvature, such as near the obstacle, shrinks faster according to motion by mean curvature. It is observed that the different mean curvature causes different results of temporal evolution depending on whether the obstacle exists as shown in the bottom row of Figure 16.

Next, we conduct the 3D experiments using the non-convex initial condition, where the values inside the surface are $\phi = 1$ and the outside is $\phi = -1$ as shown in Figure 17(a). Figures 17(a)–17(d) illustrate the temporal evolutions of zero-level isosurface of the numerical solutions without and with the obstacle. In the second row of Figure 17, the green object represents the obstacle where the inside and boundary of the obstacle are fixed to $\phi = -1$. Here, we use the numerical parameters as $N = 150$, $h = 1/N$, $\Delta t = 0.1h^2$, and $\epsilon = 4h/(2\sqrt{2}\tanh^{-1}(0.9))$ in the domain $\Omega = (0, 1)^3$. Even in 3D space, the mean curvature of the area with the obstacle is maintained to be large, and the solution contracts faster than when there is no obstacle. The difference between 2D and 3D simulations is that there are two principal curvatures in the 3D domain, which complicates the analysis of the numerical results more than in the 2D cases. However, the basic mechanism is the same as the 2D simulations.

4.   Conclusions

We presented the explicit spatially fourth-order accurate compact method for the AC equation. We use a fourth-order compact FDM, not only being more accurate but needing relatively less wide stencils for the discrete Laplace operator. The time step restriction that is the typical problem of the explicit Euler scheme is not severe because the governing equation is a second-order parabolic PDEs, and small temporal step sizes should be used for an accurate numerical solution. Moreover, it is simple to implement. Therefore, the proposed numerical algorithm is fast and efficient. We observed the phase separation phenomena using random number initial conditions, calculated the rate of convergence to verify the fourth-order accuracy in space, verified the hyperbolic tangent equilibrium profile, and numerically demonstrated the dissipation of discrete total energy. We verified the efficiency of the proposed method by comparing the computational cost between our method and the existing method. We performed various numerical simulations for motion by mean curvature, the representative property of the AC equation. Moreover, experiments comparing the numerical solutions with and without solid obstacles were performed to indicate that the proposed numerical algorithm is efficient and simple to implement. In previous research ^[41], the maximum bound principle and energy dissipation were proven for a one-step explicit method to the AC equation. On the other hand, this is a non-trivial and challenging task for the current proposed scheme, which involves a two-step method using the ($n-1$)-, $n$-, and ($n+1$)-th time step solutions. In this study, we focused on proposing a fully explicit fourth-order compact scheme and presenting numerical experiments of the proposed method for the AC equation. In future work, theoretical proof for the maximum bound principle, energy dissipation, and error estimates will be considered, drawing insights from the studies conducted by ^[42,43].

Use of AI tools declaration

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

Acknowledgments

The first author (C. Lee) was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. 2022R1C1C2003896). The authors extend their thanks to the reviewers for the valuable and constructive input they provided during the revision of the article.

Conflict of interest

The authors declare there is no conflicts of interest.