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

Authors: Pengfei Guo; Ariunkhishig Boldbaatar; Dutao Yi; Pengxiang Dai
Journal: Mathematical Problems in Engineering
 (2020)

---

## 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.


---

## Body

## 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-freekp-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],(1)∂2ux,t∂t2+β∂ux,t∂t=∂2ux,t∂x2+∂2ux,t∂y2−ϕx,ysinux,t,within the problem domain of Ω=x,ya≤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. The functions ϕx,y 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:(2)∂ux,y,t∂n=0,x,y∈Γ.At timet=0, we assume the initial conditions to be in the form of(3)ux,y,0=φ1x,y,∂ux,y,0∂t=φ2x,y,where the functions, φ1x,y and φ2x,y, 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 Ith of which is assumed to be located at xII=1,2,…,N, where N is the total number of nodes. The field function, such as ux, is approximated/interpolated via the nodal values uxI. Through a weighted linear combination, the field value ux0 at point x0 can be evaluated approximately as(4)u^x0=∑I=1nλIuxI,xI∈Ωs,where uxI stands for the nodal value at xII=1,2,…,n and ΩsΩs⊂Ω represents the support domain associated with the point at x0, which contains n nodes that exert influence on the point at x0. The weight coefficient λI will be determined according to the Kriging interpolation method and the Kriging-based shape function can thus be written by(5)u^x,xI=Ψx,xIuI=pxTAxI+rxTBxIuI.The matricesAxI and BxI are expressed, respectively, as(6)AxI=PTR−1P−1PTR−1,BxI=R−1I−PA,where I is the unit matrix of order n, and px 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:(7)pTx=1xy,m=3,pTx=1xyx2xyy2,m=6.In terms of the basis function at nodes in the support domain, the matrixP is then written as(8)P=p1x1⋯pmx1⋮⋱⋮p1xn⋯pmxn.The matrixR and vector rx are, respectively, expressed as(9)R=γx1,x1⋯γx1,xn⋮⋱⋮γxn,x1⋯γxn,xn,(10)rTx=γx,x1γx,x2γx,xn,each entry of which, γxI,xJ, is evaluated based on the semivariogram model that determines the relation of any pair of nodes at xI and xJ as follows:(11)γxI,xJ=12EZbxI−ZbxJ2=γh,where E⋅ 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 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(12)γh=c01−exp−3h2a02h≤a0.In equation (12), h is the Euclidean distance of any point at x and the node at xI, or the distance of any pair of nodes at xI and xJ. Thus, h=x−xI for equation (10) and h=xI−xJ for equation (9). The two factors a0 and c0 are the sill and range that are empirically determined. In the present study, we set c0 to 1 because its value has little influence on the shape functions, whereas a0 is evaluated according to the influence domain size ds linked to nodes(13)a0=α0⋅ds,α0≥1,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 dI that is written as(14)dI=β0⋅dave,where β0 is a dimensional coefficient and dave 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 x1 and x2 based on equation (5), and their first- and second-order partial derivatives are given as follows:(15)Ψ,ix,xI=p,ixTAxI+r,ixTBxI,Ψ,ijx,xI=p,ijxTAxI+r,ijxTBxI,in which ,i and ,ij represent ∂/∂xi and ∂2/∂xi∂xj, 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Ωq associated with the nodes in the global problem domain Ω (see Figure 1). The 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(16)∫ΩqIW^I∂2u∂t2+β∂u∂t−∂2u∂x2−∂2u∂y2+ϕx,ysinudΩ=0,where we specify a cubic spine function [44] as the test function W^I in this study, which is defined as(17)W^I=W^x,xI=23−4r¯I2+4r¯I3,r¯I≤0.5,43−4r¯I+4r¯I2−4r¯I33,0.5<r¯I≤1,0,r¯I>1,in which rI¯=x−xI/rw.Figure 1
Local quadrature, support, and test function domains.By decomposing the subdomain boundaryΓq into three sub-boundaries and applying the Gauss divergence theorem, equation (16) is re-expressed as(18)∫ΩqIW^I∂2u∂t2+β∂u∂tdΩ−∫ΓqIi+ΓqIuW^Inx∂u∂x+ny∂u∂ydΓ+∫ΩqIW^I,x∂u∂x+W^I,y∂u∂ydΩ=∫ΓqItW^Inx∂u∂x+ny∂u∂ydΓ−∫ΩqIW^Iϕx,ysinudΩ,in which nx and ny stand for the outward normal direction cosines of the boundary ΓqI of the local subdomain ΩqI and the boundary contains three disjoint parts ΓqI=ΓqIi∪ΓqIu∪ΓqIt.According to the Kriging interpolation method and the constructed shape function (equation (5)) presented in the previous section, the field function ux 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:(19)uhx,t=∑J=1nψJxuxJ,t=∑J=1nψJuJ.The time and space derivatives of the filed function are expressed in the following form of(20)∂ux,t∂t=∑J=1nψJx∂uxJ,t∂t=ΨxU˙t∂2ux,t∂t2=∑J=1nψJx∂u2xJ,t∂t2=ΨxU¨t,∂ux,t∂xi=∑J=1n∂ψJx∂xiuxJ,t=ΨixUt,where(21)Ψx=ψ1xψ2x…ψnx,Ut=u1tu2t…untT.Substituting the approximate field functionuhx in equation (19) for ux 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:(22)CIU¨+βCIU˙+KIU=FIU,in which(23)CI=∫W^IΨxdΩ,(24)KI=∫ΩqIW^I,xΨ,xx+W^I,yΨ,yxdΩ−∫ΓqIi+ΓqIuW^InxΨ,xx+nyΨ,yxdΓ,(25)FIU=−∫ΩqIW^Iϕx,ysinudΩ+∫ΓqItW^Inx∂u∂x+ny∂u∂ydΓ.In order to construct the global system equation, the assembling process similar to the finite difference method has been used for allN nodes scattered in the problem domain Ωand it leads to(26)CU¨+βCU˙+KU=FU.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.


### 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(27)CUk+2−2Uk+1+UkΔt2+βCUk+1−UkΔt+KUk+1+Uk2=F.The above equation can be rewritten as(28)2CUk+2=4C−2βCΔt−KΔt2Uk+1+−2C+2βCΔt−KΔt2Uk+2Δt2F.Whent=0, U0 is determined according to the specified initial condition (equation (3)), that is, U0=φ1x,y, and for the next time level, U1 can be obtained as U1=U0+2Δtφ2x,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˜=Uk at the last time level t=k to calculate the term 2Δt2F in the equation, and the system of equations is then solved linearly with Uk+2=Uk+2,0. Afterwards, we recalculate U˜=1/3Uk+2,0+Uk+1+Uk. Equation (28) can be solved with the latest available U˜ to obtain the field function values Uk+2,1 at the time level of k+2. Performing repeatedly this iterative procedure and setting U˜=1/3Uk+2,l+Uk+1+Uk, 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:(29)Uk+2,l−Uk+2,l−1≤ε.We can setUk+2=Uk+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.


## 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],(1)∂2ux,t∂t2+β∂ux,t∂t=∂2ux,t∂x2+∂2ux,t∂y2−ϕx,ysinux,t,within the problem domain of Ω=x,ya≤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. The functions ϕx,y 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:(2)∂ux,y,t∂n=0,x,y∈Γ.At timet=0, we assume the initial conditions to be in the form of(3)ux,y,0=φ1x,y,∂ux,y,0∂t=φ2x,y,where the functions, φ1x,y and φ2x,y, 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 Ith of which is assumed to be located at xII=1,2,…,N, where N is the total number of nodes. The field function, such as ux, is approximated/interpolated via the nodal values uxI. Through a weighted linear combination, the field value ux0 at point x0 can be evaluated approximately as(4)u^x0=∑I=1nλIuxI,xI∈Ωs,where uxI stands for the nodal value at xII=1,2,…,n and ΩsΩs⊂Ω represents the support domain associated with the point at x0, which contains n nodes that exert influence on the point at x0. The weight coefficient λI will be determined according to the Kriging interpolation method and the Kriging-based shape function can thus be written by(5)u^x,xI=Ψx,xIuI=pxTAxI+rxTBxIuI.The matricesAxI and BxI are expressed, respectively, as(6)AxI=PTR−1P−1PTR−1,BxI=R−1I−PA,where I is the unit matrix of order n, and px 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:(7)pTx=1xy,m=3,pTx=1xyx2xyy2,m=6.In terms of the basis function at nodes in the support domain, the matrixP is then written as(8)P=p1x1⋯pmx1⋮⋱⋮p1xn⋯pmxn.The matrixR and vector rx are, respectively, expressed as(9)R=γx1,x1⋯γx1,xn⋮⋱⋮γxn,x1⋯γxn,xn,(10)rTx=γx,x1γx,x2γx,xn,each entry of which, γxI,xJ, is evaluated based on the semivariogram model that determines the relation of any pair of nodes at xI and xJ as follows:(11)γxI,xJ=12EZbxI−ZbxJ2=γh,where E⋅ 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 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(12)γh=c01−exp−3h2a02h≤a0.In equation (12), h is the Euclidean distance of any point at x and the node at xI, or the distance of any pair of nodes at xI and xJ. Thus, h=x−xI for equation (10) and h=xI−xJ for equation (9). The two factors a0 and c0 are the sill and range that are empirically determined. In the present study, we set c0 to 1 because its value has little influence on the shape functions, whereas a0 is evaluated according to the influence domain size ds linked to nodes(13)a0=α0⋅ds,α0≥1,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 dI that is written as(14)dI=β0⋅dave,where β0 is a dimensional coefficient and dave 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 x1 and x2 based on equation (5), and their first- and second-order partial derivatives are given as follows:(15)Ψ,ix,xI=p,ixTAxI+r,ixTBxI,Ψ,ijx,xI=p,ijxTAxI+r,ijxTBxI,in which ,i and ,ij represent ∂/∂xi and ∂2/∂xi∂xj, 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Ωq associated with the nodes in the global problem domain Ω (see Figure 1). The 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(16)∫ΩqIW^I∂2u∂t2+β∂u∂t−∂2u∂x2−∂2u∂y2+ϕx,ysinudΩ=0,where we specify a cubic spine function [44] as the test function W^I in this study, which is defined as(17)W^I=W^x,xI=23−4r¯I2+4r¯I3,r¯I≤0.5,43−4r¯I+4r¯I2−4r¯I33,0.5<r¯I≤1,0,r¯I>1,in which rI¯=x−xI/rw.Figure 1
Local quadrature, support, and test function domains.By decomposing the subdomain boundaryΓq into three sub-boundaries and applying the Gauss divergence theorem, equation (16) is re-expressed as(18)∫ΩqIW^I∂2u∂t2+β∂u∂tdΩ−∫ΓqIi+ΓqIuW^Inx∂u∂x+ny∂u∂ydΓ+∫ΩqIW^I,x∂u∂x+W^I,y∂u∂ydΩ=∫ΓqItW^Inx∂u∂x+ny∂u∂ydΓ−∫ΩqIW^Iϕx,ysinudΩ,in which nx and ny stand for the outward normal direction cosines of the boundary ΓqI of the local subdomain ΩqI and the boundary contains three disjoint parts ΓqI=ΓqIi∪ΓqIu∪ΓqIt.According to the Kriging interpolation method and the constructed shape function (equation (5)) presented in the previous section, the field function ux 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:(19)uhx,t=∑J=1nψJxuxJ,t=∑J=1nψJuJ.The time and space derivatives of the filed function are expressed in the following form of(20)∂ux,t∂t=∑J=1nψJx∂uxJ,t∂t=ΨxU˙t∂2ux,t∂t2=∑J=1nψJx∂u2xJ,t∂t2=ΨxU¨t,∂ux,t∂xi=∑J=1n∂ψJx∂xiuxJ,t=ΨixUt,where(21)Ψx=ψ1xψ2x…ψnx,Ut=u1tu2t…untT.Substituting the approximate field functionuhx in equation (19) for ux 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:(22)CIU¨+βCIU˙+KIU=FIU,in which(23)CI=∫W^IΨxdΩ,(24)KI=∫ΩqIW^I,xΨ,xx+W^I,yΨ,yxdΩ−∫ΓqIi+ΓqIuW^InxΨ,xx+nyΨ,yxdΓ,(25)FIU=−∫ΩqIW^Iϕx,ysinudΩ+∫ΓqItW^Inx∂u∂x+ny∂u∂ydΓ.In order to construct the global system equation, the assembling process similar to the finite difference method has been used for allN nodes scattered in the problem domain Ωand it leads to(26)CU¨+βCU˙+KU=FU.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.


## 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(27)CUk+2−2Uk+1+UkΔt2+βCUk+1−UkΔt+KUk+1+Uk2=F.The above equation can be rewritten as(28)2CUk+2=4C−2βCΔt−KΔt2Uk+1+−2C+2βCΔt−KΔt2Uk+2Δt2F.Whent=0, U0 is determined according to the specified initial condition (equation (3)), that is, U0=φ1x,y, and for the next time level, U1 can be obtained as U1=U0+2Δtφ2x,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˜=Uk at the last time level t=k to calculate the term 2Δt2F in the equation, and the system of equations is then solved linearly with Uk+2=Uk+2,0. Afterwards, we recalculate U˜=1/3Uk+2,0+Uk+1+Uk. Equation (28) can be solved with the latest available U˜ to obtain the field function values Uk+2,1 at the time level of k+2. Performing repeatedly this iterative procedure and setting U˜=1/3Uk+2,l+Uk+1+Uk, 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:(29)Uk+2,l−Uk+2,l−1≤ε.We can setUk+2=Uk+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.


## 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β=0 with ϕx,y=1, we obtain the superposition of two orthogonal line solitons under the given initial conditions,(30)ux,y,0=4tan−1expx+4tan−1expy,∂ux,y,0∂t=0,and over a square domain of Ω=x,y−6≤x≤6,−6≤y≤6. The numerical results of ux,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].Figure 2
Superposition of two orthogonal line solitons (β=0 and ϕx,y=1). (a) t=1.0, (b) t=2.0, (c) t=3.0, and (d) t=4.0.
(a)(b)(c)(d)Figure 3
Contour for superposition of two orthogonal line solitons (β=0 and ϕx,y=1). (a) t=1.0, (b) t=2.0, (c) t=3.0, and (d) t=4.0.
(a)(b)(c)(d)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.Figure 4
Contour for superposition of two orthogonal line solitons (t=4.0, ϕx,y=1). (a) β=0.5 and (b) β=1.5.
(a)(b)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]:(31)Et=12∫Ωux2+uy2+ut2+21−cosudΩ.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 Et given by the present method, which agrees well with the literature [4, 26].Table 1
The energyEt of superposition of two orthogonal line solitons (β=0 and ϕx,y=1) for t∈0,7.
tt=0t=1t=2t=3t=4t=5t=6t=7Et [4]–175.3027175.3107175.3289175.3743175.3921175.4131175.7865Et [26]175.5745175.2927175.2761175.3188175.3637175.3899175.4045175.6831Et (present)175.7503175.7515175.8155175.9262176.0511176.2110176.3888176.5849The convergence analysis of the present method is conducted based on the energy,Et, shown in Figure 5.Figure 5
Relationship between the energyEt and nodal distribution in Example 1.


### 3.2. Example 2: Circular Ring Soliton (Damped) [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,(32)ux,y,0=4tan−1exp3−x2+y2,∂ux,y,0∂t=0,and in the problem domain of Ω=x,y−7≤x≤7,−7≤y≤7. Figures 6 and 7 show the numerical solutions and their contours in term of sinux,y,t/2 for t=5.6 and 8.4, respectively. Comparing the surfaces in Figure 6, a complete agreement can be observed with those published in [4] for the damped SG equation.Figure 6
Circular ring solitons in terms ofsinux,y,t/2 (β=0.05 and ϕx,y=1). (a) t=2.8 and (b) t=5.6.
(a)(b)Figure 7
Contours of circular ring solitons in terms ofsinu/2 (β=0.05 and ϕx,y=1). (a) t=2.8 and (b) t=5.6.
(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β=0.05 and ϕx,y=1 and the specified initial conditions,(33)ux,y,0=4tan−1expx+1−2sechy+7−2sechy−7,∂ux,y,0∂t=0,over the problem domain Ω=x,y−7≤x≤7,−7≤y≤7. In Figure 8 which illustrates the symmetric perturbation of single line soliton in terms of sinux,y,t/2 for t=0,2,5,7,9,and11, two symmetric dents firstly move towards each other, then collide at certain time, and proceed to move away from each other.Figure 8
Symmetric perturbation of static line soliton in terms ofsinux,y,t/2 at different times (β=0.05 and ϕx,y=1). (a) t=0, (b) t=2.0, (c) t=5.0, (d) t=7.0, (e) t=9.0, and (f) t=11.0.
(a)(b)(c)(d)(e)(f)


## 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β=0 with ϕx,y=1, we obtain the superposition of two orthogonal line solitons under the given initial conditions,(30)ux,y,0=4tan−1expx+4tan−1expy,∂ux,y,0∂t=0,and over a square domain of Ω=x,y−6≤x≤6,−6≤y≤6. The numerical results of ux,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].Figure 2
Superposition of two orthogonal line solitons (β=0 and ϕx,y=1). (a) t=1.0, (b) t=2.0, (c) t=3.0, and (d) t=4.0.
(a)(b)(c)(d)Figure 3
Contour for superposition of two orthogonal line solitons (β=0 and ϕx,y=1). (a) t=1.0, (b) t=2.0, (c) t=3.0, and (d) t=4.0.
(a)(b)(c)(d)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.Figure 4
Contour for superposition of two orthogonal line solitons (t=4.0, ϕx,y=1). (a) β=0.5 and (b) β=1.5.
(a)(b)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]:(31)Et=12∫Ωux2+uy2+ut2+21−cosudΩ.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 Et given by the present method, which agrees well with the literature [4, 26].Table 1
The energyEt of superposition of two orthogonal line solitons (β=0 and ϕx,y=1) for t∈0,7.
tt=0t=1t=2t=3t=4t=5t=6t=7Et [4]–175.3027175.3107175.3289175.3743175.3921175.4131175.7865Et [26]175.5745175.2927175.2761175.3188175.3637175.3899175.4045175.6831Et (present)175.7503175.7515175.8155175.9262176.0511176.2110176.3888176.5849The convergence analysis of the present method is conducted based on the energy,Et, shown in Figure 5.Figure 5
Relationship between the energyEt and nodal distribution in Example 1.


## 3.2. Example 2: Circular Ring Soliton (Damped) [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,(32)ux,y,0=4tan−1exp3−x2+y2,∂ux,y,0∂t=0,and in the problem domain of Ω=x,y−7≤x≤7,−7≤y≤7. Figures 6 and 7 show the numerical solutions and their contours in term of sinux,y,t/2 for t=5.6 and 8.4, respectively. Comparing the surfaces in Figure 6, a complete agreement can be observed with those published in [4] for the damped SG equation.Figure 6
Circular ring solitons in terms ofsinux,y,t/2 (β=0.05 and ϕx,y=1). (a) t=2.8 and (b) t=5.6.
(a)(b)Figure 7
Contours of circular ring solitons in terms ofsinu/2 (β=0.05 and ϕx,y=1). (a) t=2.8 and (b) t=5.6.
(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β=0.05 and ϕx,y=1 and the specified initial conditions,(33)ux,y,0=4tan−1expx+1−2sechy+7−2sechy−7,∂ux,y,0∂t=0,over the problem domain Ω=x,y−7≤x≤7,−7≤y≤7. In Figure 8 which illustrates the symmetric perturbation of single line soliton in terms of sinux,y,t/2 for t=0,2,5,7,9,and11, two symmetric dents firstly move towards each other, then collide at certain time, and proceed to move away from each other.Figure 8
Symmetric perturbation of static line soliton in terms ofsinux,y,t/2 at different times (β=0.05 and ϕx,y=1). (a) t=0, (b) t=2.0, (c) t=5.0, (d) t=7.0, (e) t=9.0, and (f) t=11.0.
(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.



---

*Source: 9057387-2020-09-08.xml*