Numerical Solution of Sine-Gordon Equation with the Local Kriging Meshless Method

'is paper develops a local Kriging meshless solution to the nonlinear 2 + 1-dimensional sine-Gordon equation. 'e meshless shape function is constructed by Kriging interpolationmethod to have Kronecker delta function property for the two-dimensional field function, which leads to convenient implementation of imposing essential boundary conditions. Based on the local Petrov–Galerkin formulation and the center difference method for time discretization, a system of nonlinear discrete equations is obtained.'e numerical examples are presented and the numerical solutions are found to be in good agreement with the results in the literature to validate the ability of the present meshless method to handle the 2 + 1-dimensional sine-Gordon equation related problems.


Introduction
e nonlinear sine-Gordon equation (SGE), a type of hyperbolic partial differential equation, is often used to describe and simulate the physical phenomena in a variety of fields of engineering and science, such as nonlinear waves, propagation of fluxons and dislocation of metals [1][2][3][4]. Because sine-Gordon equation leads to different types of soliton solutions, SGE has been receiving an enormous amount of attention. Soliton solution travels without experiencing any deformation through the medium, even when it collides with another soliton. e solitons, identified in many wave and particle systems, are of importance in the theory of nonlinear differential equations.
As one of the crucial equations in nonlinear science, the sine-Gordon equation has been constantly investigated and solved numerically and analytically in recent years [4][5][6][7][8]. For the one-dimensional sine-Gordon equation, Wazwaz [9] adopted the tanh method to find the traveling wave solutions. Johnson et al. [10] provided three exact solutions to deal with the two-dimensional SGE. Zhong and Belić analytically established special two-soliton solution to solve the generalized SGE based on the Hirota bilinear method and self-similar method [11]. For more related research about exactly solving the sine-Gordon equations, the interested readers may refer to [12][13][14][15]. Despite the exact solutions or analytical analyses of the SGEs that can offer deeper understanding and insights of the associated physical phenomena, numerical methods are still indispensable in pursuit those unknown soliton solutions to the SGEs. Dehghan et al. used the collocation method based on the radial basis function to obtain the numerical solutions to the one-dimensional nonlinear sine-Gordon equation [16]. Besides, the difference approach [17], the boundary integral equation scheme [18], the homotopy-perturbation method [19], the modified Adomain decomposition method [20], and the multilevel augmentation method [21] have also been applied to solve the one-dimensional SGEs. For dealing with the higher dimensional problems, Argyris et al. [22] proposed a finite element algorithm to obtain soliton-type solutions. Djidjeli et al. [23] developed an explicit numerical approach to solve a damped sine-Gordon equation in two space variables. Sheng et al. presented a split cosine scheme [24], and Bratsos used the method of lines [25] for seeking the solitary solutions of the two-dimensional sine-Gordon equation.
During recent decades, meshless methods have also been developed for solving the nonlinear sine-Gordon equations, which include the mesh-free kp-Ritz method [4], the radius basis function (RBF) method [7,16], the radial point interpolation method (RPIM) [8], the meshless local Petrov-Galerkin (MLPG) method [26], and the well-posed moving least-squares (WP-MLS) approximation [27]. Meshless methods are a type of numerical schemes by which the problem domain is represented by a set of scattered nodes, and the mesh or discretization can thus be avoided. ey have been usually used to effectively handle those physical and engineering problems or to solve the partial differential equations that might pose challenges of the traditional numerical approaches, such as the boundary element method [28] and finite element method [29]. A variety of meshless methods have been developed and then applied in scientific and engineering modelling based on different combinations of certain formulation schemes, such as the element-free Galerkin formulation [30], and the approximation/interpolation techniques, such as moving least-squares (MLS) [31,32]. Based on MLS, the element-free Galerkin method has been viewed as one of the predominant meshless methods. Cheng et al. proposed the improved moving least-squares (IMLS) approximation to develop the improved element-free Galerkin (IEFG) method [33][34][35] and the interpolating element-free Galerkin method [36][37][38]. By introducing the dimension splitting method into the reproducing kernel particle method, a new meshless method was developed to solve three-dimensional potential problems [39][40][41][42][43].
In this study, we present the local Kriging meshless method [44][45][46][47] to seek the numerical solutions of the 2 + 1dimensional nonlinear sine-Gordon equation. e present meshless approach is based on the Kriging interpolation technique and the local Petrov-Galerkin formulation, which leads to the shape functions having the Kronecker delta property and thus to a convenient implementation required for imposing the essential boundary conditions, like in the traditional finite element analysis. Combining with the Petrov-Galerkin formulation over the local subdomains, the meshless method is featured by avoiding background cells for the numerical integration of the global weak form. Based on the center difference scheme for making time discretization, we have a global system of nonlinear algebraic equations. Numerical tests are performed to validate the convergence and accuracy of the present meshless method for solving the 2 + 1-dimensional sine-Gordon equation. e numerical results are found to agree well with the existing results reported in the published literature.

Governing Equation.
In the present study, we seek to acquire a meshless approximation to the damped 2 + 1-dimensional sine-Gordon equation [4], within the problem domain of Ω � (x, y) | a ≤ x ≤ b, c ≤ y ≤ d} for t > 0. In equation (1), the parameter β ≥ 0 is the damping factor for the dissipative term. When β � 0, this equation will reduce to an undamped sine-Gordon equation. e functions ϕ(x, y) represent Josephson current density. e boundary conditions linked to the sine-Gordon equation (equation (1)) enforce a zero normal gradient on the boundary Γ of Ω as follows: At time t � 0, we assume the initial conditions to be in the form of where the functions, φ 1 (x, y) and φ 2 (x, y), are specified wave modes and velocity.

Kriging
Interpolation. e problem domain Ω with the boundary Γ can be represented by a set of scattered nodes, the I th of which is assumed to be located at x I (I � 1, 2, . . . , N), where N is the total number of nodes. e field function, such as u(x), is approximated/interpolated via the nodal values u(x I ). rough a weighted linear combination, the field value u(x 0 ) at point x 0 can be evaluated approximately as where u(x I ) stands for the nodal value at x I (I � 1, 2, . . . , n) and Ω s (Ω s ⊂ Ω) represents the support domain associated with the point at x 0 , which contains n nodes that exert influence on the point at x 0 . e weight coefficient λ I will be determined according to the Kriging interpolation method and the Kriging-based shape function can thus be written by e matrices A(x I ) and B(x I ) are expressed, respectively, as where I is the unit matrix of order n, and p(x) is a column vector that contains m monomial basis functions, which are given, as an example, in case of two-dimensional linear and quadratic basis functions, as follows: x y x 2 xy y 2 , (m � 6).

(7)
In terms of the basis function at nodes in the support domain, the matrix P is then written as e matrix R and vector r(x) are, respectively, expressed as each entry of which, c(x I , x J ), is evaluated based on the semivariogram model that determines the relation of any pair of nodes at x I and x J as follows: where E[·] stands for an expected value of a random function. e semivariogram model in equation (11) is independent of the locations of the two nodes but only depends on their lag vector h.
In this numerical investigation, we adopt the Gauss model to construct the shape functions via the Kriging interpolation scheme [48], which is given as In equation (12), h is the Euclidean distance of any point at x and the node at x I , or the distance of any pair of nodes at x I and x J . us, h � ‖x − x I ‖ for equation (10) and h � ‖x I − x J ‖ for equation (9). e two factors a 0 and c 0 are the sill and range that are empirically determined. In the present study, we set c 0 to 1 because its value has little influence on the shape functions, whereas a 0 is evaluated according to the influence domain size (d s ) linked to nodes in which α 0 is a specified factor [46]. In the local Kriging meshless method, the interpolation domain for any point at x is constructed in terms of the influence domain size of d I that is written as where β 0 is a dimensional coefficient and d ave represents the average spacing among nodes, which is determined by considering the balance of computational costs and accuracy [49,50].
In the problem domain Ω, we have the derivatives of the constructed shape functions with respect to x 1 and x 2 based on equation (5), and their first-and second-order partial derivatives are given as follows: in which ( ) ,i and ( ) ,ij represent z( )/zx i and z 2 ( )/zx i zx j , respectively.

Local Weak Form of the Governing Equation.
In the present study, we apply the local Petrov-Galerkin formulation to construct the weak form of the governing equation over the pre-established local subdomains (Ω q ) associated with the nodes in the global problem domain Ω (see Figure 1). e positive definite test functions W I that are assigned to each node need to be predefined within a test function region (Ω w ). Over one local subdomain Ω qI centered at node I, the local weighted residual weak form of equation (1) is then expressed as where we specify a cubic spine function [44] as the test function W I in this study, which is defined as By decomposing the subdomain boundary Γ q into three sub-boundaries and applying the Gauss divergence theorem, equation (16) is re-expressed as in which n x and n y stand for the outward normal direction cosines of the boundary (Γ qI ) of the local subdomain Ω qI Mathematical Problems in Engineering and the boundary contains three disjoint parts According to the Kriging interpolation method and the constructed shape function (equation (5)) presented in the previous section, the field function u(x) to be determined can be approximated in terms of nodal values at the n nodes that are included in the support domain associated with the point at x as follows: e time and space derivatives of the filed function are expressed in the following form of where Substituting the approximate field function u h (x) in equation (19) for u(x) in the local weak form of the governing equation (equation (18)), we obtain the nodal system equation for the Ith node in the matrix form: in which In order to construct the global system equation, the assembling process similar to the finite difference method has been used for all N nodes scattered in the problem domain Ω and it leads to Apparently, equation (26) is a system of nonlinear algebraic equations, and an iterative technique is then applied to cope with the nonlinearity. In equation (25), the first term on the right-hand side needs to be evaluated according to the latest available approximation of the field function u. Figure 1: Local quadrature, support, and test function domains.

e Time Discretization and Iterative
Procedure. In the present study, the time derivatives of the global system equations (equation (25)) is handled by the center difference technique for making time discretization, and the global system equation (equation (25)) is transformed into e above equation can be rewritten as When t � 0, U (0) is determined according to the specified initial condition (equation (3)), that is, U (0) � φ 1 (x, y), and for the next time level, U (1) can be obtained as U (1) � U (0) + 2Δtφ 2 (x, y). For the purpose of handling the nonlinearity of equation (28) using a predictor-corrector technique [4,26], at a time level of t � k + 1, we need the field function value U � U (k) at the last time level (t � k) to calculate the term 2Δt 2 F in the equation, and the system of equations is then solved linearly with U (k+2) � U (k+2), 0 . Afterwards, we recalculate U � (1/3)[U (k+2),0 + U (k+1) + U (k) ]. Equation (28) can be solved with the latest available U to obtain the field function values U (k+2), 1 at the time level of (k + 2). Performing repeatedly this iterative procedure and setting U � (1/3)[U (k+2),l + U (k+1) + U (k) ], we then achieve the nonlinear solutions to equation (28) when t � k + 2 if the solutions can converge to a prescribed number ε. In the present study, the condition for terminating the iteration is given as follows: We can set U (k+2) � U (k+2),l and proceed to the next time level, in case the above condition is met.
By conducting the above iterative procedure to solve the equation (equation (28)) till the desired time level, we achieve the numerical solutions to the 2 + 1-dimensional nonlinear sine-Gordon equation.

Numerical Examples
In this section, we apply the local Kriging meshless method to several numerical examples regarding two-dimensional line and ring solitons. [4,26,51]. In the case of an undamped 2D sine-Gordon equation (β � 0) with ϕ(x, y) � 1, we obtain the superposition of two orthogonal line solitons under the given initial conditions,

Example 1: Superposition of Two Orthogonal Line Solitons (Undamped and Damped SG Equation)
and over a square domain of Ω � (x, y) | − 6 ≤ x ≤ 6, −6 ≤ y ≤ 6}. e numerical results of u(x, y, t) are shown in Figure 2 for t � 1.0, 2.0, 3.0, and 4.0, respectively. For better depicting the breakup of two orthogonal line solitons, we illustrate the contours of the superposition of the two orthogonal line solitons in Figure 3 that move separately from each other along the direction y � −x over time. All numerical solutions are in good agreement with those available results in [4,26,51]. In order to investigate the effects of the dissipative term on the behavior, the same numerical example is used with β � 0.5 and 1.5, respectively. Figures 3(d), 4(a), and 4(b) illustrate the contours of the line solitons with varying damping factors when t � 4, in which exist apparent propagation delay of the line solitons due to the presence of the dissipative term in the SG equation.
In order to quantitatively compare the numerical results with the literature, we adopt the energy for an undamped sine-Gordon equation (β � 0), which is conserved and defined as the following form [4]: e integration of the above equation over the problem domain Ω is calculated according to the trapezoidal rule. Table 1 exhibits the conservation of the energy E(t) given by the present method, which agrees well with the literature [4,26]. e convergence analysis of the present method is conducted based on the energy, E(t), shown in Figure 5. [4,22]. For the numerical case, in which β � 0.05 and ϕ(x, y) � 1, we numerically obtain the circular ring solitons with the specified initial conditions,

Concluding Remarks
In this study, the local Kriging meshless method has been extended to the 2 + 1-dimensional nonlinear sine-Gordon equation. e Kriging interpolation technique is used to approximate the two-dimensional field function, which leads to convenient implementation of imposing essential boundary conditions. A system of nonlinear discrete equations can be established based on the adoption of the local Petrov-Galerkin formulation and the center difference method for time discretization. e nonlinear algebraic equations are solved by applying the iterative technique and a predictor-corrector scheme. e numerical results are in good agreement with those available in the literature. e present meshless method is thus a potential alternative to other numerical methods, such as the finite element method and finite difference method, for dealing with a 2 + 1-dimensional nonlinear sine-Gordon equation.

Data Availability
All results have been obtained by conducting the numerical procedure and the ideas can be shared for the researchers.

Conflicts of Interest
e authors declare that they have no conflicts of interest.