Parents’ perception and willingness to maintain provider care continuity for their children under universal health coverage

1.   Introduction

In this study, we present a novel explicit numerical method for solving the Allen–Cahn equation on triangular meshes. The method is characterized by its simplicity and efficiency. The AC equation, a partial differential equation (PDE), is employed in mathematical physics and materials science to represent the dynamics of phase transitions and pattern formation in binary systems ^[1]. The equation incorporates a free energy density function that characterizes the potential energy of the system and a linear diffusion term that drives the system toward a state of minimum energy ^[2]. In two-dimensional (2D) space, the AC equation with Dirichlet boundary condition on the domain $ \Omega $ as follows

where $ u(x, y, t) $ is the phase-field function, $ \epsilon $ is an interface transition layer related parameter, $ F(u) = 0.25(u^2-1)^2 $ is a double-well potential energy functional, and $ \Delta $ is a Laplace operator.

Choi et al. ^[1] developed an unconditionally gradient-stable numerical system for the AC equation on uniformly square grid and showed the same property for the corresponding discrete problem by using the eigenvalues of the Hessian matrix of the energy functional. In ^[3], the authors developed a simple and practical adaptive finite difference scheme for AC equations in 2D. The developed method uses a temporally adaptive narrow band domain in a uniform discrete rectangular domain. Xiao et al. ^[4] introduced the lumped mass surface finite element method on a triangular mesh to solve the surface AC equation. To maintain mass conservation, Rubinstein and Sternberg ^[5] introduced a conservative AC equation by incorporating a Lagrange multiplier into the AC equation. Construction of the conservative AC equation addresses the issue of the AC equation inability to conserve total mass. Sun and Zhang ^[6] developed a meshless radial basis function approach to solve the conservative AC equation on smooth compact surfaces embedded in $ \mathbb{R}^3 $. The developed method ensures mass conservation while employing an operator splitting scheme. The Laplace–Beltrami operator (LBO) is approximated iteratively by using radial basis functions. In addition, the time discretization of the diffusion equation is achieved using the Euler method. Xia et al. ^[7] developed a second-order accurate method in time and space for the AC and conservative AC equation using triangular discretization. The AC and conservative AC equation can also be solved using the explicit hybrid numerical method, as studied by Choi et al. ^[8]. Joshi and Jaiman ^[9] proposed an adaptive variational procedure for the AC equation with a two-phase field. This scheme is stable, robust, and general for the AC equation on a triangular mesh. Kwak et al. ^[10] proposed the conservative AC equation with a curvature-dependent Lagrange multiplier. Hong et al. ^[11] developed numerical algorithms with fully discrete structure-preserving of arbitrarily high order for the AC model with a nonlocal constraint subject to the Neumann boundary condition. Chai et al. ^[12] proposed a simple multiple-relaxation-time lattice-Boltzmann method for the local and nonlocal AC equation, presenting the property of mass conservation through numerical experiments. Inc et al. ^[13] performed Lie symmetry analysis, explicit solutions, and convergence analysis for the time-fractional Cahn–Allen equation. Kim et al. ^[14] proposed a temporally second-order unconditionally energy-stable computational method for the AC equation with high-order polynomial free energy potential. Bhatt et al. ^[15] devised an attractive and easy-to-implement alternative for integrating the multi-dimensional AC equation with no-flux boundary conditions by combining the Fourier spectral method with the strongly stable exponential time difference method. The method they developed has several key advantages, including the use of discrete fast Fourier transform for efficiency, the ability to extend to two and three spatial dimensions in a manner similar to 1D problems, and the capability to handle various boundary conditions. Sun and Gao ^[16] proposed a high order multiquadric trigonometric quasi-interpolation method for function approximation and derivative approximation based on periodic sampling data. Furthermore, the method was applied to solve different types of PDEs including the AC equation. Sun and Ling ^[17] developed a meshless conservative Galerkin method for solving Hamiltonian wave equations, which is a classical example of Hamiltionian PDEs. The equation was first discretized in space using radial basis functions in a Galerkin-type formulation. Their method used appropriate projection operators for the construction of the Galerkin equation and was shown to conserve global energies.

The AC equation without a nonlinear term is a diffusion equation that is a linear differential equation:

where $ D $ is the diffusion coefficient parameter and we set $ D = 1 $ in this paper for simplicity of exposition.

There exist meshes that do not satisfy the properties of the continuous Laplace operator for continuous functions using the discrete Laplace operator, and it is more difficult to satisfy the properties of the continuous Laplace operator for triangular meshes ^[18]. The method we propose requires defining a discrete Laplace operator to solve the diffusion term. Thus, a simple and efficient algorithm is included that numerically satisfies the simple properties of the Laplace operator.

Triangular meshes are widely used for representing complex domains or curved surfaces. Until recently, the discrete Laplace operator for triangular meshes has been a subject of numerous studies. Xu ^[19] proposed discretization schemes of the Laplace-Beltrami operator on triangulated surfaces with curvatures. Meanwhile, Caissard et al. ^[20] presented a new discretization of the Laplace-Beltrami operator over digital surfaces. Thampi et al. ^[21] demonstrated schemes for isotropic discrete Laplacian operators based on lattice hydrodynamics. McCartin ^[22] studied the discrete Laplacian on an equilateral triangle with both Dirichlet and Neumann boundary conditions using eigenvalues and eigenfunctions of the continuous Laplacian. A finite difference algorithm for solving the generalized Laplace equation on unstructured triangular meshes was developed by Ganzha et al. ^[23] using a support operator method combining discrete divergence and gradient operators. Yoon et al. ^[24] conducted numerical experiments on triangular meshes to investigate dendritic pattern formation in an isotropic crystal growth model on surfaces, and subsequently presented their findings. Paquet and Korikache ^[25] established an a priori error estimate for the fully discretized problem of the dual mixed method for the non-stationary heat equation in 2D polygonal domains. The Cahn–Hilliard (CH) equation is one of the famous phase-field models along with the AC equation ^[26]. Tian et al. ^[27] developed an adaptive finite element method based on superconvergent cluster recovery to solve the Laplace operator of the CH equation on a triangular mesh. Therefore, the construction of simple and efficient algorithms for solving Laplace operators on triangular meshes is important. We propose a numerical method involving the novel discrete Laplace operator and verify the proposed discrete Laplace operator through numerical experiments, including comparisons to general discrete Laplace operators on rectangular meshes or the LBO on triangular meshes.

The contents of this paper are as follows. In Section 2, the proposed numerical scheme is described. In Section 3, we present the numerical experiments. In Section 4, the conclusion is drawn.

2.   Numerical solution algorithm

In this section, we present a numerical solution algorithm. The proposed algorithm solves the AC equation by using an explicit Euler method. The numerical algorithm requires computing a linear term with the Laplace operator. We propose a novel Laplace operator in triangular meshes to numerically compute the AC equation. Firstly, we introduce the discrete LBO on triangular mesh and propose a novel discrete Laplace operator. Because the plane is a special case of the curved surfaces, let us consider the discretization of LBO ^[19] on curved surfaces. Let $ N $ be the number of node points on the triangular mesh, $ {\bf x}_k = (x_k, y_k) $, $ k = 1, 2, \dots, N $ be node points on the triangular mesh, $ \Delta t $ be temporal step size, and $ u_k^n $ be an approximation of phase-field function $ u({\bf x}_k, n\Delta t) $. In Figure 1, $ {\bf x}_{k} $ denotes $ k $-th node point of a triangular mesh and $ \{ {\bf x}_{k_1}, {\bf x}_{k_2}, {\bf x}_{k_3}, {\bf x}_{k_4}, {\bf x}_{k_5}\} $ are one-ring node points neighboring $ {\bf x}_k $. For simplicity, we assume that there are five node points neighboring $ {\bf x}_k $. The discrete LBO at node point $ {\bf x}_{k} $ is given as

where $ A({\bf x}_{k}) $ is the sum of areas of triangles $ T_m $ sharing the node point $ {\bf x}_{k} $, $ N_k $ is number of node points neighboring $ {\bf x}_k $, $ {\bf x}_{k_0} = {\bf x}_{k_{N_k}} $ and $ {\bf x}_{k_{N_k+1}} = {\bf x}_{k_{1}} $. We also define angles $ \alpha_{k_{m}} = \angle {\bf x}_k{\bf x}_{k_{m-1}}{\bf x}_{k_{m}} $ and $ \beta_{k_{m}} = \angle {\bf x}_k{\bf x}_{k_{m+1}}{\bf x}_{k_{m}} $ for $ m = 1, \cdots, N_k $. For example, when $ m = 5 $, $ \alpha_{k_{5}} = \angle {\bf x}_k{\bf x}_{k_{4}}{\bf x}_{k_{5}} $ and $ \beta_{k_{5}} = \angle {\bf x}_k{\bf x}_{k_{6}}{\bf x}_{k_{5}} $ in triangles $ T_{5} $ and $ T_{6} $, respectively, as shown in Figure 1.

Figure 1.  Schematic illustration of the proposed method.

DownLoad: Full-Size Img PowerPoint

Now, we propose the novel discrete Laplace operator. Let $ N $ be the total number of node points in the triangular mesh and $ N_k $ be the number of one-ring neighboring node points of each node $ {\bf x}_k $, for $ k = 1, 2, \dots, N $ in the triangular mesh. The neighboring node points of each node point $ {\bf x}_k $ are denoted as $ {\bf x}_{k_m} $, $ m = 1, 2, \dots, N_k $. We define the net force $ {{\bf{F}}}_k $ for each node point $ {\bf x}_k $, $ k = 1, 2, \dots, N $, which plays a crucial role in achieving symmetry within the proposed discrete Laplace operator.

The net force at each points in a symmetric object or system is zero because the forces acting on the object are symmetrically distributed. Additionally, the net force of the general five and nine stencil Laplace operators is $ {\bf 0} $. Therefore, we propose the following sub-algorithm using virtual nodes such that the net force $ {{\bf{F}}}_k $ of each node $ {\bf x}_k $, $ k = 1, 2, \dots, N_k $ becomes $ {\bf 0} $. We solve the AC equation on the irregular domain by using a sub-algorithm using virtual nodes to make each node point symmetric for the asymmetric triangular mesh generated on the irregular domain in 2D space.

The sub-algorithm of the proposed method is illustrated in Figure 2. As presented in Figure 2, for each node $ {\bf x}_k $, one-ring neighboring node points are denoted by $ {\bf x}_{k_1} $, $ {\bf x}_{k_2} $, $ \dots $, $ {\bf x}_{k_6} $. First, find all the one-ring neighboring node points for a given point $ {\bf x}_k $. Figure 2(a) presents a given point $ {\bf x}_k $ and one-ring neighboring node points. Second, find the net force $ {{\bf{F}}}_k $ and the triangular domain containing $ {{\bf{F}}}_k $. Here, we use the following method to find the triangle. In the case of Figure 2(b), let the vectors from $ {\bf x}_k $ of the found triangle to the other two node points $ {\bf x}_{k_3} $ and $ {\bf x}_{k_2} $ be $ {{\bf{u}}} $ and $ {{\bf{v}}} $, respectively. Therefore, net force $ {\bf{F}}_k $ is the linear combination of $ {\bf{u}} $ and $ {\bf{v}} $ with constants $ a $ and $ b $ as follows:

where $ {{\bf{u}}} $ and $ {{\bf{v}}} $ are 2D vectors. Let $ {{\bf{u}}} = (u_1, \; u_2)^T $ and $ {{\bf{v}}} = (v_1, \; v_2)^T $. For the third step, we consider the following the matrix form to find values of $ a $ and $ b $.

Figure 2.  Adjusting the net force $ {{\bf{F}}} $ to be $ {\bf 0} $: (a)–(e) represent each step.

DownLoad: Full-Size Img PowerPoint

Then, subtract $ a{{\bf{u}}} $ and $ b{{\bf{v}}} $ from $ {{\bf{u}}} $ and $ {{\bf{v}}} $ to reduce the net force $ {{\bf{F}}}_k $ to $ {\bf 0} $ at $ {\bf x}_k $. Let $ x^*_{k3} $ and $ x^*_{k2} $ be virtual node points corresponding to $ {{\bf{u}}}-a{{\bf{u}}} $ and $ {{\bf{v}}}-b{{\bf{v}}} $, respectively. In the context of the triangular mesh, it is important to note that the physical positions of the node points remain unchanged. Instead, a theoretical construct of imaginary node points is utilized to ensure that the net force acting on the mesh is balanced to zero. By introducing these virtual node points, it is possible to accurately calculate the net force and achieve a state of equilibrium, without the need to manipulate the physical node points. Fourth, we use interpolation for the points where the extension of the modified and original meshes meet on the opposite side of $ {\bf{u}} $ and $ {\bf{v}} $, respectively. In Figure 2(e), the value of $ {\bf x}_{k_2}^{**} $ is obtained by the linear interpolation of the values of $ {\bf x}_{k_5} $ and $ {\bf x}_{k_6} $ as follows:

Here, the virtual node point $ {\bf x}^{**}_{k_2} $ is the intersection point of the line passing through $ {\bf x}_{k} $ and $ {\bf x}_{k_2} $, and the line passing through $ {\bf x}_{k_5} $ and $ {\bf x}_{k_6} $. In addition, $ {\bf x}_{k_3}^{**} $ is obtained analogously. The net force at a point is zero in a symmetric object or system, as the forces acting on it are distributed symmetrically, ensuring that they cancel out each other. Therefore, when virtual node points are used, each node point of the triangular mesh becomes symmetrical with respect to one-ring neighboring points. Finally, the values of modified mesh points are obtained by quadratic interpolation. In Figure 2(f), the value of virtual node $ {\bf x}_{k_3}^{*} $ is obtained by the quadratic interpolation using the values of $ {\bf x}_{k_3}^{**} $, $ {\bf x}_{3} $, and $ {\bf x}_{k_3} $. The value of virtual node $ {\bf x}_{k_2}^{*} $ is similarly obtained. Next, we define the discrete Laplace operator with the values of the neighboring node points set such that the net force $ {{\bf{F}}}_k $ at a given point is zero. Let $ {\bf x}_{k_m} $ for $ m = 1, \dots, N_k $ be the one-ring neighboring node points at a point $ {\bf x}_k $ with net force $ {{\bf{F}}}_k = {\bf 0} $, then we propose the discrete Laplace operator as follows:

where $ S_k $ is a constant given for each $ {\bf x}_k $. The general Laplace operator satisfies $ \Delta f(x, y) = 1 $ if $ f(x, y) = 0.25(x^2+y^2) $. Because the proposed Laplace operator should satisfy the characteristics of the general Laplace operator, we define $ S_k $ as

Here, $ f(x, y) = 0.25(x^2+y^2) $ and $ {\bf d}_{k_m} = {\bf x}_{k_m}-{\bf x}_k $. The reason for using this function is that it is a symmetric function with respect to the $ x $- and $ y $-axes, and unlike a linear function, a non-zero solution exists, and unlike a cubic function, the solution is uniquely determined. Therefore, the proposed discrete Laplace operator satisfies the characteristic of general Laplace operator, $ u(x, y) = 0.25(x^2+y^2) $. Then we discretize Eq (1.1) in temporal and spatial as follows

where $ S_k = \sum_{m = 1}^{N_k}f({\bf d}_{k_m}) $. Therefore, we obtain the following equation by rearranging Eq (2.1) for $ u_k^{n+1} $.

3.   Numerical experiments

3.1.   Discrete Laplace operator

We consider the general five stencil, nine stencil and exact Laplace operators to compare the proposed Laplace operator with the conventional Laplace operators.

which are the well-known standard discrete Laplace operators. Figures 3(a), (b), (c) represent the numerical stencils for Eqs (3.1), (3.2) and (3.3), respectively.

Figure 3.  Stencils for (a) $ \Delta_1 u_{ij} $ using five points, (b) $ \Delta_2 u_{ij} $ using rotated five points, and (c) $ \Delta_3 u_{ij} $ using nine points.

DownLoad: Full-Size Img PowerPoint

3.1.1.   Convergence test on rectangular mesh

In this section, each discrete Laplace operator is compared with the $ C^2 $ function to verify the performance of the proposed discrete Laplace operator. The proposed discrete Laplace operator $ {\tilde \Delta} $ at each node is made by neighboring nodes as given in Figure 3(c). This is one of the benefits of our proposed discrete Laplace operator on a triangular mesh. Therefore, we obtain the discrete Laplace equation as follows using the proposed method.

We consider the following function on the domain $ \Omega = (0, 1)\times (0, 1) $:

with zero Dirichlet boundary condition. Then, the exact solution of $ \Delta u $ is

We use different parameter values $ N_x = 9, 17, 33, 65 $ to consider the convergence rates for different spatial mesh $ h_x = h_y = 1/(N_x-1) $ and discrete domain $ x_i = (i-1)h_x $ for $ i = 1, 2, \dots, N_x $ and $ y_j = (j-1)h_y $ for $ j = 1, 2, \dots, N_y $. Next, to calculate the convergence rate for $ \Delta u_{ij} $, we define the discrete $ l_2 $-error and discrete maximum-error in 2D space, respectively:

where $ {e}_{ij} = {\hat{\Delta}} u_{ij}-{\Delta u_{ij}} $ for $ i = 1, 2, \dots, N_x $ and $ j = 1, 2, \dots, N_y $. Figures 4(a), 4(b) present the $ l_2 $-error and the maximum-error for different Laplace operators, respectively. We make a list of the convergence rates for the maximum-error and $ l_2 $-error of the Laplace operator in Tables 1 and 2, respectively. We observed that the proposed discrete Laplace operator achieved second-order accuracy in spatial directions on a rectangular mesh. The proposed discrete Laplace operator can consider triangular meshes on irregular domains rather than rectangular meshes, while other discrete Laplace operators can only be used on rectangular meshes.

Figure 4.  The value of error for each Laplace operator.

DownLoad: Full-Size Img PowerPoint

Table 1.  Maximum-error of different Laplace operators with respect to different spatial steps.

spatial step	$ h(=1/8) $	rate	$ h/2 $	rate	$ h/4 $	rate	$ h/8 $
$ \Delta_1 u $	$ 2.5237\times10^{-1} $	1.99	$ 6.3336\times10^{-2} $	2.00	$ 1.5849\times10^{-2} $	2.00	$ 3.9633\times10^{-3} $
$ \Delta_2 u $	$ 9.9404\times10^{-1} $	1.98	$ 2.5237\times10^{-1} $	1.99	$ 6.3336\times10^{-2} $	2.00	$ 1.5849\times10^{-2} $
$ \Delta_3 u $	$ 4.9959\times10^{-1} $	1.98	$ 1.2635\times10^{-1} $	2.00	$ 3.1678\times10^{-2} $	2.00	$ 7.9253\times10^{-3} $
$ {\tilde \Delta} u $	$ 7.4682\times10^{-1} $	1.98	$ 1.8936\times10^{-1} $	1.99	$ 4.7507\times10^{-2} $	2.00	$ 1.1887\times10^{-2} $

---

 | Show Table

DownLoad: CSV

Table 2.  $ l_2 $-error of different Laplace operators with respect to different spatial steps.

spatial step	$ h(=1/8) $	rate	$ h/2 $	rate	$ h/4 $	rate	$ h/8 $
$ \Delta_1 u $	$ 1.1216\times10^{-1} $	1.91	$ 2.9805\times10^{-2} $	1.96	$ 7.6845\times10^{-3} $	1.98	$ 1.9511\times10^{-3} $
$ \Delta_2 u $	$ 4.4180\times10^{-1} $	1.90	$ 1.1876\times10^{-1} $	1.95	$ 3.0708\times10^{-2} $	1.98	$ 7.8027\times10^{-3} $
$ \Delta_3 u $	$ 2.2204\times10^{-1} $	1.90	$ 5.9457\times10^{-2} $	1.95	$ 1.5359\times10^{-2} $	1.98	$ 3.9017\times10^{-3} $
$ {\tilde \Delta} u $	$ 3.3192\times10^{-1} $	1.90	$ 8.9110\times10^{-2} $	1.95	$ 2.3034\times10^{-2} $	1.98	$ 5.8522\times10^{-3} $

---

 | Show Table

DownLoad: CSV

3.1.2.   Laplace operator for functions on triangular mesh

The proposed method is for 2D triangular mesh, therefore we need to define some conditions. Let $ N $ be the number of nodes in the triangular mesh and $ {\bf x}_k = (x_k, y_k) $ be $ k $-th node in the triangular mesh for $ k = 1, \dots, N $. Set $ I $ be the set of indices of the interior nodes and $ N_k $ is the number of node points neighboring $ {\bf x}_k $. The average of spatial mesh is

where $ {\bf x}_{k_m} $ is neighboring node points of $ {\bf x}_k $ for $ m = 1 $, $ \dots $, $ N_k $.

Let the initial function $ u(x, y) = ax+by+c $ for $ \; a, b, c \in {\mathbb R} $ be a plane, then $ \Delta u = 0 $. Therefore, we consider the initial condition on the domain $ \Omega = (0, 1)\times(0, 1) $:

Here, we apply Dirichlet boundary conditions. To compute the discrete Laplace operator of $ u $, we use a triangular mesh with $ h_{\rm ave} = 0.1013 $. Figures 5(a), (b), (c) present the triangular mesh, $ u(x, y) = x+2y+3 $, and proposed discrete Laplace operator for $ u $, respectively.

Figure 5.  Discrete Laplace operator for plane equation: (a) triangular mesh, (b) $ u(x, y) = x+2y+3 $, and (c) $ {\tilde \Delta}u(x, y) $.

DownLoad: Full-Size Img PowerPoint

Next, we consider the initial function on the domain $ \Omega = (0, 1)\times(0, 1) $:

Then we have

We apply the zero Dirichlet boundary condition. To compute the discrete Laplace operator of $ u $, we use the parameter $ h_{\rm ave} = 0.1013 $. Figures 6(a), (b), (c) show the triangular mesh, $ u(x, y) = \sin(\pi x)\sin(\pi y) $, and the proposed discrete Laplace operator for $ u $, i.e., $ {\tilde \Delta}u(x, y) $, respectively.

Figure 6.  (a) Triangular mesh, (b) $ u(x, y) = \sin(\pi x)\sin(\pi y) $, and (c) $ {\tilde \Delta}u(x, y) $.

DownLoad: Full-Size Img PowerPoint

3.2.   Diffusion equations on triangular mesh

3.2.1.   Convergence test on square domain

In this section, we compare the numerical and analytical solutions to verify the proposed method. We consider a diffusion equation (1.2) with the following initial condition on the domain $ \Omega = (0, 1)\times(0, 1) $:

with zero Dirichlet boundary condition. Then, the exact solution is

Let $ u_k^n $ be an approximation of $ u(x_k, y_k, n\Delta t) $, where $ \Delta t $ is a time step. Let $ N $ be the number of nodes of the triangular mesh, $ {\bf x}_k = (x_k, y_k) $ be $ k $-th node in triangular mesh for $ k = 1, \dots, N $. For this test, we use different triangular meshes with $ h_{\rm ave} = 0.1254, \; 0.1013, \; 0.0762 $. Using the explicit Euler method and the proposed discrete Laplace operator, Eq (1.2) is as follows:

Let us define the $ l_2 $-error and maximum-error as follows:

where $ {e}^{N_t}_{k} = u_k^{N_t}-u_{exact}(x_k, y_k, T) $, $ k = 1, 2, \dots, N $. Figure 7 presents the triangular mesh and numerical solution using the proposed method for each $ h_{\rm ave} $. We list the maximum-error and $ l_2 $-error of the diffusion equation with $ T = 0.02 $, $ N_t = 2048 $, $ \Delta t = T/N_t $ in Table 3.

Figure 7.  Triangular mesh and numerical solution using the proposed method for each different spatial mesh with $ h_{\rm ave} $. (a) $ h_{\rm ave} = 0.1254 $, (b) $ h_{\rm ave} = 0.1013 $, (c) $ h_{\rm ave} = 0.0762 $.

DownLoad: Full-Size Img PowerPoint

Table 3.  Maximum-error and $ l_2 $-error for proposed discrete Laplace operator with respect to different spatial steps.

Spatial step	$ h_{\rm ave} =0.1254 $	rate	$ h_{\rm ave}=0.1013 $	rate	$ h_{\rm ave} =0.0762 $
Maximum-error	$ 1.8338\times10^{-2} $	1.19	$ 1.4246\times10^{-2} $	1.92	$ 8.2534\times10^{-3} $
$ l_2 $-error	$ 4.2098\times10^{-3} $	1.64	$ 2.9665\times10^{-3} $	1.77	$ 1.7904\times10^{-3} $

---

 | Show Table

DownLoad: CSV

Next, we perform numerical experiments for a larger time evolution with a final time $ T = 0.1 $ using the same triangular meshes as in the previous test. Figure 8 shows the numerical solutions for different spatial steps at the final time $ T = 0.1 $. Table 4 lists the maximum-error and $ l_2 $-error for the proposed discrete Laplace operator. We observed that the $ l_2 $-error and maximum-error of the proposed discrete Laplace operator decreased as the spatial step decreased for a larger final time $ T = 0.1 $.

Figure 8.  Triangular mesh and numerical solution for a larger final time $ T = 0.1 $.

DownLoad: Full-Size Img PowerPoint

Table 4.  Maximum-error and $ l_2 $-error for proposed discrete Laplace operator with respect to different spatial steps.

Spatial step	$ h_{\rm ave} =0.1254 $	rate	$ h_{\rm ave}=0.1013 $	rate	$ h_{\rm ave} =0.0762 $
Maximum-error	$ 4.6855\times10^{-3} $	2.07	$ 3.0121\times10^{-3} $	1.98	$ 1.7142\times10^{-3} $
$ l_2 $-error	$ 2.0184\times10^{-3} $	2.89	$ 1.0894\times10^{-3} $	1.09	$ 7.9989\times10^{-4} $

---

 | Show Table

DownLoad: CSV

We perform a convergence test on the discrete LBO to compare the discrete LBO and the proposed discrete Laplace operator. The initial conditions given are the same as in the previous convergence test with the final time $ T = 0.1 $. Figure 9 shows the numerical solutions obtained using the discrete LBO for different spatial steps. Table 5 lists the maximum-error and $ l_2 $-error for the discrete LBO.

Figure 9.  Triangular mesh and numerical solution using the discrete LBO.

DownLoad: Full-Size Img PowerPoint

Table 5.  Maximum-error and $ l_2 $-error for the discrete LBO with respect to different spatial steps.

Spatial step	$ h_{\rm ave} =0.1254 $	rate	$ h_{\rm ave}=0.1013 $	rate	$ h_{\rm ave} =0.0762 $
Maximum-error	$ 5.3441\times10^{-3} $	2.11	$ 3.4101\times10^{-3} $	–7.49	$ 2.8734\times10^{-2} $
$ l_2 $-error	$ 2.3112\times10^{-3} $	2.01	$ 1.5061\times10^{-3} $	–7.16	$ 1.1569\times10^{-2} $

---

 | Show Table

DownLoad: CSV

We observed that the error of the discrete LBO is larger than that of the proposed discrete Laplace operator and does not converge as the spatial step becomes smaller.

3.2.2.   Diffusion equation on circular domain

In this section, we consider a diffusion equation on a circular domain using the proposed discrete Laplace operator. We use the circular domain of radius $ 1 $ and with center at the origin, i.e., $ \Omega = \{ (x, y) | x^2 + y^2 < 1 \} $. We use the zero Dirichlet boundary condition and $ u(x, y, 0) = 1-x^2-y^2 $ as the initial condition.

Let us consider a benchmark problem using a manufactured solution ^[28], i.e., we assume that the solution of the benchmark problem is

Then, we consider the following benchmark problem:

where $ g(x, y, t) = \partial u(x, y, t) / \partial t - \Delta u(x, y, t) = (3+x^2+y^2)e^{-t} $ is the source term. To construct the numerical solution, we use explicit Euler method as follows:

where $ g_k^{n+\frac{1}{2}} = (3 + x_k^2 + y_k^2)e^{-(n+\frac{1}{2})\Delta t} $.

Let us consider three different spatial meshes with $ h_{\rm ave} = 0.5315, \; 0.2133 $ and $ 0.1020 $. The total time is $ T = 0.1 $ and time step $ \Delta t = 0.001 h_{\rm ave}^2 $ with total iteration $ N_t = T/\Delta t $. Using the proposed discrete Laplace operator, we present the numerical solution at final time $ T $ and compare it with the exact solution, Eq (3.12). Figure 10 presents the triangular mesh, exact solution and numerical results, and the absolute error between the results.

Figure 10.  Triangular mesh, numerical and exact solutions, and absolute error, here top to bottom $ h_{\rm ave} = 0.5315, \; 0.3138, \; 0.2133 $, respectively.

DownLoad: Full-Size Img PowerPoint

3.3.   Two-dimensional Allen-Cahn equation on triangular mesh

In this section, we consider the AC Eq (1.1) on plane with triangular mesh. Let $ \Omega $ be the domain in 2D space and $ N $ be the number of nodes in domain. For $ k = 1, \dots, N $, let $ {\bf x}_k = (x_k, y_k) $ be the $ k $-th node in domain. On the computational triangular mesh, let $ u_k^n $ be an approximation of $ u(x_k, y_k, n\Delta t) $, where $ \Delta t $ is a time step.

3.3.1.   Motion by mean curvature on square domain

In Figure 12, we verify the mean curvature flow of the AC equation using following initial condition on $ \Omega = (-1, 1)\times(-1, 1) $:

where $ R_0 = 0.4 $ is the initial radius with Dirichlet boundary condition $ u(x, y, t) = -1 $ on $ \partial \Omega $. The analytic radius is given as $ R(t) = \sqrt{{R_0}^2-2t} $ for time $ t $ ^[29]. The numerical radius is computed using the average of distances between zero level points of the numerical solution and the center point $ (0, 0) $. As shown in Figure 12, the numerical radii shrink along with the analytic radius.

Next, we consider three triangular meshes generated by different spatial steps, $ h_{\rm ave} = 0.1524 $, $ 0.0997 $, $ 0.0499 $. The initial condition is Eq (3.15) and interface parameter is $ \epsilon = 1.3h_{\rm ave} $. Figures 11(a), (b), (c) show the effect of motion by mean curvature for the AC equation with the different spatial steps. Figure 12 shows temporal evolutions of analytic and numerical radii for different spatial steps. The time steps are used as $ \Delta t = 0.2h_{\rm ave}^2 $ for each $ h_{ave} $. We can observe that the numerical results agree with the analytic exact solution as the spatial step decreases to zero.

Figure 11.  Triangular mesh and zero contour for each different parameter $ h_{\rm ave} $. (a) $ h_{\rm ave} = 0.1524 $, (b) $ h_{\rm ave} = 0.0997 $, (c) $ h_{\rm ave} = 0.0499 $.

DownLoad: Full-Size Img PowerPoint

Figure 12.  Motion by mean curvature for different spatial step $ h_{\rm ave} = 0.1524, \; 0.0997, \; 0.0499 $.

DownLoad: Full-Size Img PowerPoint

3.3.2.   Dynamics of the AC equation on various domains

In this section, we observe the dynamics of the AC equation on various domains. We consider rectangular domain with circular hole $ \Omega_{rectangle} = \big((-1.3, 1.3)\times(-1.3, 1.3)\big)\cap B_2(0.5)^c $ and circular domain with circular hole $ \Omega_{circle} = B_1(1.3)\cap B_2(0.5)^c $ to verify the performance of the proposed method according to the shape of the domain, where $ B_1(r) = \{(x, y)\; |\; x^2+y^2 < r\} $ and $ B_2(r) = \{(x, y)\; |\; x^2+y^2\le r\} $. The initial condition is that all internal node points in the domain are $ 1 $ with Dirichlet boundary condition $ u(x, y, t) = 1 $ in $ \partial\Omega_{inner} $ and $ u(x, y, t) = -1 $ in $ \partial\Omega_{outer} $. Here, $ \partial \Omega_{inner} $ and $ \partial \Omega_{outer} $ are the inner and outer boundaries of the domain, respectively. We use the parameters $ h_{\rm ave} \approx 0.2, \; 0.15, \; 0.1 $, $ \epsilon = 0.15 $, $ \Delta t = 0.001 $, and $ N_t = T/ \Delta t $ to compute Eq (2.2). Here, we use $ T = 5.6 $, which is a sufficient time to reach an numerical equilibrium state for all parameters. The value of $ h_{\rm ave} $ is below each figure. Figures 13 and 14 show the dynamics of the AC equation on $ \Omega_{rectangle} $ and $ \Omega_{circle} $, respectively, for the given initial conditions.

Figure 13.  Dynamics of the AC equation on the rectangle domain with circular hole. (a) Triangular mesh, (b) Initial condition, (c) $ t = 0.4 $, (d) $ t = 5.6 $.

DownLoad: Full-Size Img PowerPoint

Figure 14.  Dynamics of the AC equation on the circular domain with circular hole. (a) Triangular mesh, (b) Initial condition, (c) $ t = 0.4 $, (d) $ t = 5.6 $.

DownLoad: Full-Size Img PowerPoint

4.   Conclusions

In this study, we proposed a novel, effective, simple, and explicit numerical method for solving the diffusion and AC equations on a triangular mesh. In the proposed algorithm, we defined a novel discrete Laplace operator. The proposed discrete Laplace operator has a property that the adjusted net vector of each point using virtual nodes is zero vector on the triangular mesh. This means that the forces acting at each point of the triangular mesh are symmetrically distributed, thus demonstrating the symmetry of the proposed method. To investigate the efficiency and stability of the proposed numerical method, we presented the numerical experiments for solving the diffusion and AC equations in various domains on triangular mesh. We verified the proposed Laplace operator through numerical experiments on general discrete Laplace operators and discrete LBOs. Additionally, we verified the proposed operator division method through numerical experiments involving motion by mean curvature using the AC equation and solving the AC equation in irregular domains. Results indicated that the proposed method can be served in various applications in 2D domain with triangular mesh. In future work, we will solve the AC equation on complex surfaces in 3D space and conduct theoretical analysis on the convergence of the proposed method.

Use of AI tools declaration

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

Acknowledgments

The corresponding author (J. S. Kim) expresses thanks for the support from the BK21 FOUR program. The authors express their sincere gratitude to the reviewers for providing valuable feedback on this revised version. Their constructive comments have significantly contributed to the improvement of the manuscript.

Conflict of interest

The authors declare there is no conflicts of interest.