Abstract
This paper develops a local Kriging meshless solution to the nonlinear 2 + 1-dimensional sine-Gordon equation. The meshless shape function is constructed by Kriging interpolation method 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. The 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.
1. Introduction
The 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–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. The 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–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–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. They 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–35] and the interpolating element-free Galerkin method [36–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–43].
In this study, we present the local Kriging meshless method [44–47] to seek the numerical solutions of the 2 + 1-dimensional nonlinear sine-Gordon equation. The 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. The numerical results are found to agree well with the existing results reported in the published literature.
2. The Local Kriging Meshless Formulation for the 2D Sine-Gordon Equation
2.1. 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 for . In equation (1), the parameter is the damping factor for the dissipative term. When , this equation will reduce to an undamped sine-Gordon equation. The functions represent Josephson current density.
The boundary conditions linked to the sine-Gordon equation (equation (1)) enforce a zero normal gradient on the boundary of as follows:
At time , we assume the initial conditions to be in the form ofwhere the functions, and , are specified wave modes and velocity.
2.2. Kriging Interpolation
The problem domain with the boundary can be represented by a set of scattered nodes, the of which is assumed to be located at , where is the total number of nodes. The field function, such as , is approximated/interpolated via the nodal values . Through a weighted linear combination, the field value at point can be evaluated approximately aswhere stands for the nodal value at and represents the support domain associated with the point at , which contains nodes that exert influence on the point at . The weight coefficient will be determined according to the Kriging interpolation method and the Kriging-based shape function can thus be written by
The matrices and are expressed, respectively, aswhere is the unit matrix of order , and is a column vector that contains monomial basis functions, which are given, as an example, in case of two-dimensional linear and quadratic basis functions, as follows:
In terms of the basis function at nodes in the support domain, the matrix is then written as
The matrix and vector are, respectively, expressed aseach entry of which, , is evaluated based on the semivariogram model that determines the relation of any pair of nodes at and as follows:where stands for an expected value of a random function. The semivariogram model in equation (11) is independent of the locations of the two nodes but only depends on their lag vector .
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), is the Euclidean distance of any point at and the node at , or the distance of any pair of nodes at and . Thus, for equation (10) and for equation (9). The two factors and are the sill and range that are empirically determined. In the present study, we set to 1 because its value has little influence on the shape functions, whereas is evaluated according to the influence domain size linked to nodesin which is a specified factor [46]. In the local Kriging meshless method, the interpolation domain for any point at is constructed in terms of the influence domain size of that is written aswhere is a dimensional coefficient and 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 and based on equation (5), and their first- and second-order partial derivatives are given as follows:in which and represent and , respectively.
2.3. 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 associated with the nodes in the global problem domain (see Figure 1). The positive definite test functions that are assigned to each node need to be predefined within a test function region . Over one local subdomain centered at node , the local weighted residual weak form of equation (1) is then expressed aswhere we specify a cubic spine function [44] as the test function in this study, which is defined asin which .
By decomposing the subdomain boundary into three sub-boundaries and applying the Gauss divergence theorem, equation (16) is re-expressed asin which and stand for the outward normal direction cosines of the boundary of the local subdomain 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 to be determined can be approximated in terms of nodal values at the nodes that are included in the support domain associated with the point at as follows:
The time and space derivatives of the filed function are expressed in the following form ofwhere
Substituting the approximate field function in equation (19) for in the local weak form of the governing equation (equation (18)), we obtain the nodal system equation for the 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 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 .
2.4. The 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
The above equation can be rewritten as
When , is determined according to the specified initial condition (equation (3)), that is, , and for the next time level, can be obtained as . For the purpose of handling the nonlinearity of equation (28) using a predictor-corrector technique [4, 26], at a time level of , we need the field function value at the last time level to calculate the term in the equation, and the system of equations is then solved linearly with . Afterwards, we recalculate . Equation (28) can be solved with the latest available to obtain the field function values at the time level of . Performing repeatedly this iterative procedure and setting , we then achieve the nonlinear solutions to equation (28) when 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 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.
3. Numerical Examples
In this section, we apply the local Kriging meshless method to several numerical examples regarding two-dimensional line and ring solitons.
3.1. Example 1: Superposition of Two Orthogonal Line Solitons (Undamped and Damped SG Equation) [4, 26, 51]
In the case of an undamped 2D sine-Gordon equation with , we obtain the superposition of two orthogonal line solitons under the given initial conditions,and over a square domain of . The numerical results of are shown in Figure 2 for , and , 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 over time. All numerical solutions are in good agreement with those available results in [4, 26, 51].
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
In order to investigate the effects of the dissipative term on the behavior, the same numerical example is used with and , respectively. Figures 3(d), 4(a), and 4(b) illustrate the contours of the line solitons with varying damping factors when , in which exist apparent propagation delay of the line solitons due to the presence of the dissipative term in the SG equation.
(a)
(b)
In order to quantitatively compare the numerical results with the literature, we adopt the energy for an undamped sine-Gordon equation , which is conserved and defined as the following form [4]:
The integration of the above equation over the problem domain is calculated according to the trapezoidal rule. Table 1 exhibits the conservation of the energy given by the present method, which agrees well with the literature [4, 26].
The convergence analysis of the present method is conducted based on the energy, , shown in Figure 5.
3.2. Example 2: Circular Ring Soliton (Damped) [4, 22]
For the numerical case, in which and , we numerically obtain the circular ring solitons with the specified initial conditions,and in the problem domain of . Figures 6 and 7 show the numerical solutions and their contours in term of for and , respectively. Comparing the surfaces in Figure 6, a complete agreement can be observed with those published in [4] for the damped SG equation.
(a)
(b)
(a)
(b)
3.3. Perturbation of a Static Line Soliton [4, 52]
The perturbation of a line soliton is studied in the numerical example with and and the specified initial conditions,over the problem domain . In Figure 8 which illustrates the symmetric perturbation of single line soliton in terms of for , two symmetric dents firstly move towards each other, then collide at certain time, and proceed to move away from each other.
(a)
(b)
(c)
(d)
(e)
(f)
4. Concluding Remarks
In this study, the local Kriging meshless method has been extended to the 2 + 1-dimensional nonlinear sine-Gordon equation. The 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. The nonlinear algebraic equations are solved by applying the iterative technique and a predictor-corrector scheme. The numerical results are in good agreement with those available in the literature. The 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
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by Hainan Provincial Natural Science Foundation of China (Grant no. 119MS039).