Semi-Implicit Scheme to Solve Allen-Cahn Equation with Different Boundary Conditions

Abstract

The aim of this paper is to give an appropriate numerical method to solve Allen-Cahn equation, with Dirichlet or Neumann boundary condition. The time discretization involves an explicit scheme for the nonlinear part of the operator and an implicit Euler discretization of the linear part. Finite difference schemes are used for the spatial part. This finally leads to the numerical solution of a sparse linear system that can be solved efficiently.

Share and Cite:

Alqanawi, B. and Aigo, M. (2023) Semi-Implicit Scheme to Solve Allen-Cahn Equation with Different Boundary Conditions. American Journal of Computational Mathematics, 13, 122-135. doi: 10.4236/ajcm.2023.131005.

1. Introduction

Allen-Cahn equation was introduced by Allen and Cahn [1] [2] and can be viewed as a simple model to study the phase separation, such as ferromagnets with positively and negatively magnetized domains and viewed as a reaction diffusion equation in material sciences. The Allen-Cahn equation is a nonlinear kind of the heat equation. Similarly, to the heat equation, it contains a term that tends to make the field constant, given by a Laplace operator accomplishing a local average around any fixed point, and send this average towards zero. Many methods can be used to solve nonlinear PDEs. In [3], the author use Adomian Decomposition Method coupled to the Lesnic’s approach to solve boundary value problems and initial boundary value problems with mixed boundary conditions for linear and nonlinear partial differential equations. Allen-Cahn equation has been applied to many problems, such as image analysis [4] [5] the motion by mean curvature flows [6]. In particular, it has become a essential model equation for the diffuse interface technique proposed to study phase transitions dynamics in materials science [7]. Furthermore, an efficient algorithm of this equation has proposed. Several numerical schemes have been suggested to solve Allen-Cahn equation, and including the finite difference method, the finite element method. In [8], the authors propose some efficient and accurate numerical methods to compute the steady-state of variable coefficients space fractional Cahn-Allen equations. The approach combines an adaptive time stepping semi-implicit gradient flow method to minimize the fractional energy functional and pseudo-spectral approximation schemes. Based on the use of a preconditioned GMRES, the space fractional Cahn-Allen equation is then solved efficiently. In [9], the author proposed and analyzed for the Cahn Allen equation by dividing the equation into linear part and nonlinear part based on the idea of operator splitting. For the linear part, it is discretized by using the Crank-Nicolson scheme and solved by finite element method. The nonlinear part is solved accurately.

In [10], the authors split the Allen Cahn equation into the linear heat and nonlinear equations; and then solve the linear part using the Fourier spectral method and the nonlinear part using an analytic closed-form solution. These steps are unconditionally stable. However, if a large time step is used, then the nonlinear part dominates the evolution and results in a sharp interfacial transition layer.

In this paper, we will look at the two dimensional form of the Allen-Cahn Equation corresponding to Equation (7.9) in [11] with g = 0 and W ( u ) = 1 4 ( 1 u 2 ) 2 , where u = u ( x , y , t ) satisfies

ε u t ( x , y , t ) = ε Δ u ( x , y , t ) ( u 3 u ) / ε + g in Ω × ( 0 , ) , (1)

supplemented with the Dirichlet or Neumann boundary conditions and the initial conditions

u ( x , y , 0 ) = f ( x , y ) in Ω . (2)

Here ( x , y ) Ω = ( 1 , 1 ) 2 . Here, f and g are two given functions, that will be fixed later for numerical study. It is a parabolic, nonlinear and nonhomogeneous PDEs in two dimensions. Moreover, the parameter t > 0 is the time and ε > 0 is a parameter to fix, which can small leading to a stiff problem (and possibly unstable discretization schemes). To get a well-posed problem, we need to add a boundary condition on the boundary of the domain, i.e. Γ : = Ω , that we choose as the non-homogeneous Dirichlet boundary condition u ( x , y , t ) = 1 on Γ . Finally, since we have an initial value problem, we provide an initial data u ( x , y , t = 0 ) : = u 0 ( x , y ) in Ω . Many numerical simulations of Allen-Cahn equation was presented, see [11] [12] for a non-exhaustive list. The paper is organized as follows. In section 2, we recall some properties related to five points approximations of Δ and the tensor product that will be used in our MATLAB code. In section 3, a semi-implicit finite difference will be formulated for homogeneous Allen Cahn equation with Dirichlet boundary condition. In section, 4, a semi-implicit finite scheme for homogeneous Allen Cahn equation with Neumann boundary condition will be presented. In section 5, we will consider the Non-homogeneous Allen Cahn equation with Neumann boundary condition as good application to valid our numerical approach, numerical tests will be presented for different choice of initial data to validate our algorithms and our code of calculation. Then we conclude.

2. Preliminary

We consider the square domain Ω = ( a ; b ) 2 , discretized by N segments in each direction x and y. The number of interior points is m : = N 1 , and the total

number of interior points is M = m 2 . The discretization step is h = b a N . The

boundary Γ is discretized by using 4N points. The boundary is composed of four parts Γ S , W , N , E , for South, West, North, East.

If ones considers a Dirichlet boundary condition, then we have u = g on the boundary Γ for solving the Laplace equation Δ u = f in Ω . We recall that under a suitable mathematical setting, the well-posedness of the problem can be proved. For the interior points, we use a five points stencil approximation of the negative Laplacian Δ

Δ u ( x i , y j ) 1 h 2 ( u i 1 , j + u i , j 1 4 u i , j + u i + 1 , j + u i , j + 1 ) (3)

on the square lattice 2 i , j N with Dirichlet boundary condition u h = 1 on the boundary grid points [13] [14] [15]. In addition we have a right hand side (rhs) f D , h M 1 . The corresponding matrix A D , h of the Laplacian is computed by using (3) for the interior points, resulting in a M × M positive definite matrix. Nevertheless, in the case where Equation (3) holds for a point ( x i , y j ) where one of the points with index ( i 1, j ) , ( i + 1, j ) , ( i , j 1 ) , ( i , j + 1 ) is on the boundary, one must use the fact that we now the boundary value. As an example, if one writes Equation (3) for j = 2 , then the point ( i ,1 ) is on the south boundary and we have u i , 1 = g i , 1 . This directly implies that considering Equation (3) one gets, for 2 i m ,

1 h 2 ( u i 1 , 2 + u i , 1 4 u i , 2 + u i + 1 , 2 + u i , 3 ) = f i , 1

If one assumes that i 2 or m, then we have

1 h 2 ( u i 1 , 2 4 u i , 1 + u i + 1 , 2 + u i , 3 ) = f i , 1 + 1 h 2 g i , 1

For the corner’s points, we must also handle the value of the left ( i = 2 ) and right ( i = m ) endpoints. From above, we see that the matrix is modified as well as the right-hand side. Globally, the approximation of the Laplacian with Dirichlet boundary condition leads to solving a linear system that we latter write

A D , h u h = f D , h .

The resulting matrix A D , h u h is sparse, positive definite and of size M × M . Let us remark that the above system can be solved efficiently by using a LU factorization or any well adapted linear system solver. In practice, the matrix A D , h can be built in MATLAB by a tensor operation, based on the function kron as

(4)

Here, kron is the Kronecker Product function computes the Kronecker tensor product of two matrices, more precisely if we have two matrices A and B the Kronecker tensor product is defined as:

A = [ a 1 a 2 a 3 a 4 a 5 a 6 ] , B = [ x 1 x 2 x 3 x 4 ] ,

then

kron ( A , B ) = [ a 1 x 1 a 1 x 2 a 2 x 1 a 2 x 2 a 3 x 1 a 3 x 2 a 1 x 3 a 1 x 4 a 2 x 3 a 2 x 4 a 3 x 3 a 3 x 4 a 4 x 1 a 4 x 2 a 5 x 1 a 5 x 2 a 6 x 1 a 6 x 2 a 4 x 3 a 4 x 4 a 5 x 3 a 5 x 4 a 6 x 3 a 6 x 4 ] .

3. Discretization for the Dirichlet Boundary Condition

3.1. Semi-Implicit Backward Method: System of Non-Linear Equation

Semi-implicit finite difference schemes are used for the Allen Cahn equation. Namely, the time derivative is approximated by the backward-Euler and a second-order backward differentiation formula scheme is applied for the Laplacian operator (see Equation (3)) whereas the semi-implicit technique is used to deal with nonlinear part. This means we write all the terms implicitly except for the nonlinear term which is treated explicitly. Let h x , h y be the node spacings in the x and y directions and δ t is the step time. If we denote by u i j k u ( x i , y + j , t k ) , here i, j is location node numbers, k is the time time step number, then after discretization get

ε ( u i j k + 1 u i j k ) δ t = ε [ u i 1 j k + 1 2 u i j k + 1 + u i + 1 j k + 1 h x 2 + u i j 1 k + 1 2 u i j k + 1 + u i j + 1 k + 1 h y 2 ] + ( u i j k ) 3 u i j k ε (5)

with initial condition: u i j 0 = f ( x i , y j ) and boundary condition:

u 0 j k = 1 , u N x j k = 1 , u 0 j k = 1 , u i N y k = 1

Rearranging the Allen Cahn equation

( 1 + 2 δ t ( 1 h x 2 + 1 h y 2 ) ) u i j k + 1 δ t h x 2 ( u i 1 j k + 1 + u i + 1 j k + 1 ) δ t h y 2 ( u i j 1 k + 1 + u i j + 1 k + 1 ) = u i j k + δ t ( u i j k ( u i j k ) 3 ) ε 2 (6)

For simplicity, we can assume h x = h y = h we have the final form of discretized Allen Cahn equation

( 1 δ t + 4 h 2 ) u i j k + 1 1 h 2 ( u i 1 j k + 1 + u i + 1 j k + 1 ) 1 h 2 ( u i j 1 k + 1 + u i j + 1 k + 1 ) = u i j k δ t + ( u i j k ( u i j k ) 3 ) ε 2 (7)

In the next subsection, a more compact form of Equation (7) will be presented.

3.2. Semi-Implicit Method: Matrix Formulation

Concerning the discretization, use a uniform time step δ t > 0 and a uniform spatial discretization size h in both the x- and y-directions. We call A D , h the matrix that represents the five points stencil discretization of the 2D negative Laplacian ( Δ ), given by Equation (3), on the square lattice 1 i , j N with Dirichlet boundary condition u h = 1 on the boundary grid points. These values are directly replaced into the system since they are given values. The matrix A D , h is therefore a positive definite matrix of size N 2 × N 2 , where N 2 is the number of interior grid points. Since we use a inhomogeneous Dirichlet boundary condition, and will have a right hand side b N 2 latter related to a source, we denote by b D the modified right hand side based on the inhomogeneous boundary condition and the five points stencil discretization of the Laplacian.

The nonlinearities in the PDE are trivial ( u 3 u ) to deal with if we choose an explicit time integration method for the nonlinear part and backward time integration for the linear part, all the mathematical details for the nonlinear discretization of Allen Cahan equation is presented in subsection 3.1, we call this method a semi-implicit backward scheme.

If we choose an explicit time integration method for Allen-Cahn equation, such as the Forward Euler finite difference, there is a strict restriction on the time step size to get convergence, for further details see [16]. Since the problem is a nonlinear dynamical system, for stability reasons, we propose taken explicitly a semi-implicit backward Euler scheme with the nonlinear part (as presented in subsection 3.1).

Reorganizing the equation gives a PDE for u n + 1 (see Equation (7)), leading then to the following linear system for each time step t n : = n δ t ,

( 1 δ t I + A ) u n + 1 = b n : = 1 δ t u n + 1 ε 2 ( u n ( u n ) 3 ) (8)

where I is the identity operator and A is the discrete two-dimensional matrix for Δ (see Equation (4)), with initial condition u 0 = u 0 and with the inhomogeneous boundary condition u n + 1 = 1 on Γ = Ω × ( 0 , ) , and the boundary condition becomes

u 0 j k = 1 , u N j k = 1 , u i 0 k = 1 , u i N k = 1 , i , j = 0 , , N .

The right-hand side b n must be updated at each time step for evolving the nonlinear system. When passing at the full discretization (both in space and time), then one gets the following linear system at each time step t n to solve

( 1 δ t I h + A D , h ) u D n + 1 : = b D n (9)

for the interior grid points. The matrix I h is the identity matrix of size ( N 1 ) 2 × ( N 1 ) 2 .

3.3. Solvers

When solving (9), we need to solve a sparse positive definite linear system at each time step. For the project, we coded a LU factorization after entering the time loop, to optimize the system solution. In addition, we also wrote Jacobi and Gauss-Seidel iterative methods that work well since the linear system is positive definite [17]. The linear system solution is fast and provides a real-time evolution of the solution. Finally, concerning the post-processing, we plot the solution onto the square at each time step.

3.4. Example

In our numerical example, we take the initial data

u ( x , y , 0 ) = { 1 if x 2 + y 2 < 0.5 1 otherwise (10)

and we work with the Dirichlet boundary condition u = 1 on Ω × ( 0 , ) .

• We will use the Gauss Seidel iterative method to compute the solution of the matrix equation. For that, we fix the tolerance (tol) = 106 and we compute the discrete solution until the stopping criteria is verified. Figure 1 and Figure 2 represent the discrete solution obtained via Gauss Seidel iteration method. The error term in Gauss Seidel iteration is defined as the difference between two consecutive iterate and more precisely our stopping test based on norm ( u n + 1 u n ) < tol . The circle shrinks with time and disappears.

Figure 1. Solution using Gauss Seidel iteration method at different times.

Figure 2. Numerical solution of Allen Cahn equation at different times.

• Like the Gauss Seidel technique, the Jacobi method is considered to solve the matrix equation associated with the Allen Cahn equation. The same phenomena is observed.

• Since, the algorithm proposed is a semi-implicit finite difference, there is no restriction on δ t and h to get convergence. In our numerical implementation, we tested different choice of δ t and h. For practical simulation, we take a small space step size h = 0.01 and a time step δ t = h 2 / 10 . Moreover, the epsilon parameter is chosen very small, ε = 10 h .

• The changes in the circle over time are shown, in Figure 1, Figure 2 with h = 0.01 , δ t = h 2 / 10 , ε = 10 h . The inner circle has a larger curvature, it shrinks faster than the outer circle, and after the inner circle disappears.

4. Discretization for the Neumann Boundary Condition

4.1. Goals

A second variation concerning the problem consists in considering a homogeneous Neumann boundary condition at the interface Γ . We can expect then a very different behavior of the numerical solution because of this fundamentally different boundary condition. Basically, the main point that must be modified is the way the discrete Laplacian matrix is computed. Here we call this matrix A N , h , which is of size N × N . To approximate the normal derivative, we use the approximation

x u n + 1 ( x N + 1 , y j ) u N + 1 , j n + 1 u N , j n + 1 h = 0

for the west side of the square, i.e. u N + 1 , j n + 1 = u N , j n + 1 (the extension to the three other sides is trivial). When considering the five points stencil approximation of the Laplacian at point ( N ,0 ) (except at the upper and lower left corners)

Δ u n + 1 ( x N , y j ) 1 h 2 ( u N 1 , j n + 1 + u N , j 1 n + 1 4 u N , j n + 1 + u N + 1 , j n + 1 + u N , j + 1 n + 1 ) (11)

then, since u N + 1 , j n + 1 = u N , j n + 1 , this simplifies as

Δ u n + 1 ( x N , y j ) 1 h 2 ( u N 1 , j n + 1 + u N , j 1 n + 1 3 u N , j n + 1 + u N , j + 1 n + 1 ) (12)

For the Neumann boundary condition, the right hand side b N is an element of N 2 . The source term is built from the right hand side of Equation (8) at the discrete points. Finally, we have to solve a linear system in the form

( 1 δ t I h + A N , h ) u N n + 1 : = b N n (13)

for the interior grid points. This can clearly be done as for the Dirichlet case, i.e. by a LU factorization (the one which is coded here), a Jacobi or Gauss-Seidel fixed point iteration. We have only retained here the LU factorization from the numerical study performed for the Dirichlet situation. Indeed, the LU factorization is done only once before entering into the time loop, and then a forward-backward substitution is applied at each time step at a cost O ( N 2 ) . The L and U matrices remain with a low bandwidth because of this property on A N , h (as well as A D , h ).

4.2. Example

In our numerical example, we take the initial data (10), namely

u ( x , y , 0 ) = { 1 if x 2 + y 2 < 0.5 1 otherwise (14)

• We set a homogeneous Neumann boundary condition u . n = 0 on Ω , where Ω is the square domain [ 1,1 ] 2 .

• Using Direct method and considering a uniform step-size h x = h y = h = 0.01 and the time step δ t = h 2 / 10 , ε = 10 h .

Figure 3(a) shows the graph of initial condition given by Equation (10).

• The large-time behavior of solutions to Allen Cahn equation is established numerically for initial condition given by Equation (10). It’s observed that the circle shrinks slowly with time and disappears, see Figure 3(b).

Figure 3. Allen Cahn problem with homogeneous Neumann boundary condition and δ t = h 2 / 10 , ε = 10 h .

If we solve the Allen Cahn equation with h = 1 / 200 , ε = 8 h with homogeneous Neumann boundary condition, the solution is unstable It’s necessary that δ t = O ( h 2 ) , for further details see [18]. So taking δ t = h and ε = 8 h , the scheme explodes and is unstable.

5. Non-Homogeneous Allen Cahn Equation with Neumann Boundary Condition

5.1. Formulation and Discretization

Considering the Non-Homogeneous Allen Cahn equation:

ε u t ( x , y , t ) = ε Δ u ( x , y , t ) + ( u u 3 ) / ε + g ( x , y ) in Ω × ( 0 , ) , (15)

where g is a given function and supplemented with the homogeneous Neumann boundary conditions

u ( x , y , t ) n = 0 on Ω × ( 0 , ) , (16)

and the initial conditions

u ( x , y , 0 ) = u 0 ( x , y ) in Ω . (17)

Here ( x , y ) Ω = ( 1 , 1 ) 2 . The same strategy applied for homogeneous Allen equation and more precisely using Equation (7), we deduce that:

( 1 δ t + 4 h 2 ) u i j k + 1 1 h 2 ( u i 1 j k + 1 + u i + 1 j k + 1 ) 1 h 2 ( u i j 1 k + 1 + u i j + 1 k + 1 ) = u i j k δ t + ( u i j k ( u i j k ) 3 ) ε 2 + g i j ε

where g i j = g ( x i , y j ) .

5.2. Numerical Results

In our simulation we consider two different initial conditions: the first one is given by Equation (10) and the second one is given by (18)

Figure 4. Neumann condition: Contour of discrete solution (CDS) for Allen Cahn equation at different times.

Figure 5. Contour of discrete solution (CDS) at different times. The straight line move to the right.

u 0 ( x , y ) = { 1 if x < 0 1 if 1 < x (18)

For numerical implementation, we consider δ t = h 2 / 10 , ε = 10 h . We make different numerical tests to see how the initial data and the function g has an effect on the behavior of the solution.

1) Considering the space domain (‒5, 5)2. If one consider as initial data (10) and g = 2 it’s observed that the circle grow instead of shrink as in the case of g = 0 , see Figure 4.

• Considering the case where g = 1 , and using Equation (18) as initial condition the line curve move to the right direction, see Figure 5.

Figure 6. Contour solution at different times. The straight line move to the right and become curved.

• Note that of g = 1 , one can that the line curve move to the left direction.

• Considering the case where g = 1 y 2 , the straight line move to the right and become curved, see Figure 6.

6. Conclusion

Using finite difference method, we can solve nonlinear PDEs, more precisely we have obtained a numerical solution to the Allen Cahn equation. We have formulated using semi-implicit finite difference method and coded the corresponding algorithm to the Allen-Cahn equation. Three methods were used with the Allen-Cahn equation: the Direct method, the Jacobi and Gauss-Seidel iteration Methods. For small mesh size and small-time steps, it is better to use iterative methods, since the size of the matrix will be very large. Dirichlet boundary condition is considered, the solution is obtained and analyzed sing Direct and iterative technique. We also solved the Allen-Cahn equation using homogenous Neumann boundary conditions, this involved using ghost points next to the boundary. In the case of Neumann boundary condition the circle with g = 2 grow instead of shrink as in the case g = 0 . Moreover, for Neumann boundary condition, if we consider an initial data set u = 1 for x < 0 and u = 1 for x > 1 . If you take g = 0 nothing was happened, but if we take g = 1 the line move in the right direction and if you take g = 1 , the line move in the other direction (left direction). Considering the case where g = 1 y 2 , the straight line moves to the right and become curved. All the numerical tests confirm that the initial data and the function g has an effect on the behavior of the solution. All the methods presented were coded using MATLAB.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1] Allen, S.M. and Cahn, J.W. (1979) A Microscopic Theory for Antiphase Boundary Motion and Its Application to Antiphase Domain Coarsening. Acta Metallurgica, 27, 1085-1095.
https://doi.org/10.1016/0001-6160(79)90196-2
[2] Allen, S.M. and Cahn, J.W. (1972) Ground State Structures in Ordered Binary Alloys with Second Neighbor Interactions. Acta Metallurgica, 20, 423-433.
https://doi.org/10.1016/0001-6160(72)90037-5
[3] Abdullah Alzaid, N. and Omar Bakodah, H. (2018) Numerical Treatment of Initial-Boundary Value Problems with Mixed Boundary Conditions. American Journal of Computational Mathematics, 8, 153-174.
https://doi.org/10.4236/ajcm.2018.82012
[4] Benĕ, M. Chalupecky, V. and Mikula, K. (2004) Geometrical Image Segmentation by the AllenCahn Equation. Applied Numerical Mathematics, 51, 187-205.
https://doi.org/10.1016/j.apnum.2004.05.001
[5] Dobrosotskaya, J.A. and Bertozzi, A.L. (2008) A Wavelet Laplace Variational Technique for Image Deconvolution and Inpainting. The IEEE Transactions on Image Processing, 17, 657-663.
https://doi.org/10.1109/TIP.2008.919367
[6] Feng, X.B. and Prohl, A. (2003) Numerical Analysis of the Allen Cahn Equation and Approximation for Mean Curvature Flows. Numerische Mathematik, 94, 33-65.
https://doi.org/10.1007/s00211-002-0413-1
[7] Chen, L.Q. (2002) Phase Field Models for Microstructure Evolution. Annual Review of Materials Research, 32, 113-140.
https://doi.org/10.1146/annurev.matsci.32.112001.132041
[8] Alzahrani, S.M. and Chokri, C. (2022) Preconditioned Pseudo-Spectral Gradient Flow for Computing the Steady-State of Space Fractional Cahn-Allen Equations with Variable Coefficients. Frontiers in Physics, 10, Article 844294.
https://doi.org/10.3389/fphy.2022.844294
[9] Huang, Y.Q., Yang, W., Wang, H. and Cui, J.T. (2019) Adaptive Operator Splitting Finite Element Method for Allen-Cahn Equation. Numerical Methods for Partial Differential Equations, 35, 1290-1300.
https://doi.org/10.1002/num.22350
[10] Yang, J.X., Li, Y.B., Lee, C., Choi, Y. and Kim, J. (2023) Fast Evolution Numerical Method for the Allen-Cahn Equation. Journal of King Saud University—Science, 35, Article 102430.
https://doi.org/10.1016/j.jksus.2022.102430
[11] Deckelnick, K., Dziuk, G. and Elliott, C.M. (2005) Computation of Geometric Partial Differential Equations and Mean Curvature Flow. Acta Numerica, 14, 139-232.
https://doi.org/10.1017/S0962492904000224
[12] Elliott, C.M. (1997) Approximation of Curvature Dependent Interface Motion. In: Duff, I. and Watson, G.A. Eds., State of the Art in Numerical Analysis, Clarendon Press, Oxford, 407-440.
[13] Causon, D.M. and Mingham, C.G. (2010) Introductory Finite Difference Methods for PDEs. Bookboon, 1-144.
[14] LeVeque, R.J. (2007) Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems. Society for Industrial and Applied Mathematics, Philadelphia, PA.
https://doi.org/10.1137/1.9780898717839
[15] Larsson, S. and Thomee, V. (2009) Finite Difference Methods for Elliptic Equations. In: Bloch, A., Epstein, C.L., Goriely, A. and Greengard, L. Texts in Applied Mathematics, SpringerLink, Heidelberg, 43-49.
[16] Lee, C., Choi, Y. and Kim, J. (2022) An Explicit Stable Finite Difference Method for the Allen-Cahn Equation. Applied Numerical Mathematics, 182, 87-99.
https://doi.org/10.1016/j.apnum.2022.08.006
[17] William, M.K. (1958) Gauss Seidel Methods of Solving Large Systems of Linear Equations. Ph.D. Thesis, University of Toronto, Toronto.
[18] Bao, W. and Dong, X. (2012) Analysis and Comparison of Numerical Methods for the Klein Gordon Equation in the Nonrelativistic Limit Regime. Numerische Mathematik, 120, 189-229.
https://doi.org/10.1007/s00211-011-0411-2

Copyright © 2024 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.