Application of DRBEM for 2D sine-Gordon equation

This research paper introduces an application of dual reciprocity boundary element method (DRBEM) for the solution of sine-Gordon equation (SGE) in two-space dimension. Initially, the time derivatives are expanded using central difference schemes. After inserting the finite difference approximations into the governing equation, the pattern of the modified Helmholtz equation is obtained. The fundamental solution of modified Helmholtz equation is employed in the integral equation formulation. The inhomogeneous terms of the equation cause a domain integral in the integral equation formulation which leads to loss of advantage of the method. The DRBEM provides to transform the domain integral into the boundary integral by approximating the inhomogeneous term of the equation with thin plate spline ( ). First, code validation of the procedure is done using a test problem and then proposed method is applied for several cases involving line and ring solitons to display its capacity to treat the problem. Presented numerical results are observed to be in good agreement with other numerical results available in the literature.


Introduction
The soliton model can be used for describing for many physical phenomena. A soliton is a solitary wave and it behaves like a particle such that whose shape, amplitude and velocity are not changed after interacting or colliding with another soliton [1]. There are many causes of soliton formation. One of them is the balance between the dispersion and nonlinearity [1,2]. One of these models is SGE, which is a nonlinear partial differential equation with a nonlinear sine source term and its solution has the soliton like structure [3]. This equation has great importance in physical phenomena and some important applications include numerical forecasting, optical fibres, nonlinear optics, quantum field theory, etc. [4][5][6][7][8][9][10][11][12]. Particularly, two-dimensional SGE can be used for identifying Josephson-junction model in quantum tunnelling [13].
Time-dependent nonlinear SGE in two space variables is modelled by where z = z(x 1 , x 2 , t) and D = {(x 1 , x 2 ) : |x 1 | ≤ L x 1 , |x 2 | ≤ L x 2 }. Initial conditions associated with Equation (1) will be assumed to be of the form z(x 1 , x 2 , 0) = w 1 (x 1 , x 2 ), (x 1 , x 2 ) ∈ D, z t (x 1 , x 2 , 0) = w 2 (x 1 , x 2 ), (x 1 , and Neumann boundary conditions Here, the function χ(x 1 , x 2 ) named as Josephson current density, the parameter ν is a positive real number and it is known as dissipative term w 1 (x 1 , x 2 ) and w 2 (x 1 , x 2 ) functions are wave modes or kinks and velocity, respectively. is the boundary of the domain D and n is the unit outward normal to the boundary. When ν = 0, Equation (1) is reduced to the undamped sine-Gordon equation, while ν > 0 reduces Equation (1) to the damped sine-Gordon equation.
In the related literature, it is seen that both numerical and exact solutions can be obtained for twodimensional SG equation. Hirota method was used to obtain the exact solution of two-dimensional SGE in [14] and to obtain soliton solutions for more general nonlinear equation in [15]. Another study for the exact solution of two-and three-dimensional time-dependent SGE was presented using Backlund method in [16]. Particular solutions of the SGE were obtained using Lamb's method in [17][18][19]. To date, although many exact solutions of SGE have been found, numerical methods are needed since there are still unsolved soliton solutions. When looking at the studies for the numerical solution, it is seen that many different methods were used. Finite difference method (FDM) is used by Djidjeli et al. in [20], by Duncan in [21], by Wong in [12], by Minzoni et al. in [9], by Bratsos in [22,23], by Cui in [24] and by Liang et al. in [25]. Also, Jiwari [26] and Pekmen [13] used differential quadrature method (DQM) where the intervals are discretized using Chebyshev-Gauss-Lobatto points. Another numerical study with DQM was presented by Shukla et al. in [27] where modified cubic Bspline was used as a basis function. Asgari and Hosseini considered SGE in [28]. In this article, Fourier-spectral method was used for the solution of spatial variables and fourth-order time stepping schemes were used for the time derivatives. Another work with Fourier-spectral method was presented by Jiang et al. in [29]. In this paper, the fourth-order average vector field method was utilized for the discretization in time. Time-spaced pseudo-spectral method was applied for the solution of one-dimensional SGE using Chebyshev-Gauss-Lobatto points in [30]. In [3], the same method was used for the solution of two-dimensional SGE by presenting the stability analysis of the method. Another numerical study for the SGE was done with spectral method using Legendre wavelets as basis in [31]. In addition mesh-free methods were applied for the solution of the SGE presented in [32][33][34][35]. Dehghan and Shokri used radial basis functions for the solution of twodimensional SGEs in [32], Cheng and Liew used meshfree reproducing kernel particle Ritz method in [33] and Li et al. used element-free Galerkin method in [34], LingDe Su used localized method of approximate particular solutions in [35]. Other presented numerical studies for the SGE were given by Argyris et al. in [36] using finite element method, by Sheng et al. in [37] using split cosine scheme, by Bratsos in [38] using method of lines, by Han and Zhang in [39] using operator splitting method, by Suarez in [40] using Chebyshev split-step scheme, by Adak and Natarajan in [41] using virtual element method, by Guo et al. in [42] using local Kriging meshless method and by Jiwari in [43] using barycentric rational interpolation and local radial basis functions based numerical algorithms, etc.
Another preferred numerical solution procedure for the solution of problems modelled with partial differential equations is the boundary element method (BEM) which has found an important place among the other numerical methods. In the BEM solution procedure, the differential equation is transformed into an integral equation that contains only boundary integrals [44]. Comparing to the domain discretization methods, the BEM eliminates the need of the discretization of the whole domain which causes large quantities of data. Thus the dimension of the problem can be reduced and obtained smaller linear systems need less computational effort. On the other hand, the use of BEM as a solution technique of a differential equation depends on the existence of the fundamental solution of that differential equation. However, obtaining the fundamental solution of a differential equation such as nonlinear and time-dependent equations is often not straightforward. Therefore the method should be extended for these cases but this has some difficulties. When this differential equations are transformed into the boundary integral equation, the domain integrals occur and so the method loses its advantage. Thus DRBEM was developed to deal with the difficulty in BEM caused by the domain integral [45].
DRBEM is based on the idea of transforming domain integrals to boundary integrals using fundamental solutions of the Laplace equation. All terms in the original governing equation apart from the Laplace operator are treated as nonhomogeneous terms and expanded with the help of interpolating functions [45]. The DRBEM idea can also be extended to other differential equations for which fundamental solutions exist or can be obtained. In [46], DRBEM was built by the utilization of the fundamental solution of modified Helmholtz equations and used for the solution of many partial differential equations. DRBEM based on the fundamental solution of Laplace equation was used as numerical technique for the solution of SGE by Dehghan and Mirzaei in [47,48]. Two-dimensional SGE was solved by using constant and linear boundary elements in [47] and [48], respectively. In these studies, inhomogeneous terms were approximated by using linear radial basis functions 1 + r.
The main purpose of this study is to show that DRBEM which was built using the fundamental solution of modified Helmholtz equation is suitable as a numerical solution procedure for the solution of the SGE. In this work, DRBEM is used for the solution of SGE after converting it into a modified Helmholtz equation. To obtain the modified Helmholtz equation, the time derivatives are first approximated with central finite difference schemes. Then these approximations are inserted into the governing equation and the governing equation is written in the form of nonhomogeneous modified Helmholtz equation. This enables us to use more information from the original governing equation. Inhomogeneous terms are approximated using thin plate spline r 2 ln r. Constant boundary elements are used for the discretization of the boundary of the domain and thus obtained small size algebraic system can be solved effectively at a reasonable cost. In addition, the requirement for special attention caused by the geometric discontinuity of the boundary nodes that arise in the use of other element types such as linear elements has been eliminated by using constant boundary elements. The code validation of method is carried out by solving a test problem which has an exact solution. Later DRBEM is utilized for both line and ring solitons. Several examples are used to evaluate the capability of the present procedure and the consequences show that the proposed method captures the behaviour of each types of soliton very well.

Mathematical formulation
In order to solve SGE numerically, DRBEM based on fundamental solution of modified Helmholtz equation is adopted. Thus as an initial step, it is necessary to write the original equation in modified Helmholtz equation form. The transformation procedure begins by approximating time derivatives in equation (1). Thus the approximations are where Relaxing the solution z at three time levels using the parameters 0 < σ 1 , σ 2 , σ 3 < 1, σ 1 + σ 2 + σ 3 = 1 as and substituting approximations (4)-(6) into Equation (1), we have where ψ(z (ξ ) ) = χ(x 1 , x 2 ) sin(z (ξ ) ).
Equation (13) can be written in a general form as The solution procedure starts with the solution of Equation (10) by considering it as a modified Helmholtz equation as (15) to obtain the unknown values of z at ξ = 1. The computation of the right-hand side of Equation (15) requires information at −τ . This difficulty is easily overcome by the initial conditions given in (2) as and Then, Equation (11) is solved by considering it as a modified Helmholtz equation as by using the known values of z (0) and z (1) . Then the iterative solution procedure begins with the solution of Equation (14) by taking ξ = 2 and by using the known values of z at ξ = −1, 0, 1, 2. This procedure continues up to the stopping criteria is satisfied.

DRBEM formulation of the problem
In the previous section, the iterative form of the sine-Gordon equation is written as the modified Helmholtz equations given in (15), (18) and (14). Thus the fundamental solution of modified Helmholtz equation z * = 1 2π K 0 (λr) can be used to transform of the domain integrals caused by the inhomogeneous terms to the boundary integrals. Here, K 0 (λr) is the second kind modified Bessel function of order 0. First, for simplicity let us write Equations (15), (18) and (14) in the compact form as where Here By applying the Greens second identity to the resulting weighted residual statements, the left-hand side of the governing equations (26)-(28) is transformed into boundary integral equations as in [44] e j z where e j (j = 1, . . . , N B + N I ) are coefficients and their values change depend on the location of the collocation points (x j , y j ). For constant elements, e j are defined as In order to transform the domain integrals caused by the inhomogeneities to the boundary integrals, the functions b i (x 1 , x 2 ) (i = 1, 2, 3) are expanded as where f s (x, y) = r 2 s ln r s and The values of functions f s can be represented as f js at the collocation point for j = 1, . . . , N B + N I and F can be constructed as a (N B + N I ) × (N B + N I ) square matrix such that F(j, s) = f js . Therefore, α is (i = 1, 2, 3) can be written as The set of functions f s provide the following relationship: whereẑ i,s , (i = 1, 2, 3) are the particular solutions of the corresponding equations. The particular solutions and their normal derivatives can be given aŝ where r x 1 = x 1 − x 1s , r x 2 = x 2 − x 2s and n = ∂x 1 ∂n , ∂x 2 ∂n . γ = 0.57721566490153 is Euler constant and K 1 (x) is the second kind Bessel function of order 1.
Inserting Equations (38)-(40) into the corresponding inhomogeneous terms in Equations (32)-(34), we obtain the modified Helmholtz operators in the domain integrals of Equations (29)- (31). Then using the Green's second identity for the right-hand side of the equations, we have Equations (41)- (43) can be described as the matrixvector form by using the coefficient matrices H i and G i as where their components for i = 1, 2, 3 are such that s is the sth constant boundary element and Also, the dimension of the α i vectors for i = 1, 2, 3 are (N B + N I ) × 1 and they are defined as Z 1 , Z 2 , Z 3 , Q 1 , Q 2 , and Q 3 are all square matrices with the dimensions of (N B + N I ) × (N B + N I ) and they are constructed by using the vectors z 1 , z 2 , z 3 , q 1 , q 2 , and q 3 as column vectors, respectively.

Implementation of the numerical method
In this section, we report the numerical simulations obtained by using the DRBEM for the solution of twodimensional sine-Gordon equation. First, the test problem is solved by using DRBEM for code validation and the results are compared with the exact solution.
Then five problems of the line and ring solitons are considered in order to present numerical solutions of equation (1) with the initial conditions (2) and Neumann boundary conditions (3).

Test problem
The problem defined by sine-Gordon equation [32,43,49] is solved as a test problem. Here, z = z(x 1 , x 2 , t) and D = {(x 1 , x 2 ) : |x 1 | ≤ 7, |x 2 | ≤ 7}, χ = 1 and ν = 0. Initial conditions will be considered as and Neumann boundary conditions are exp(2(x 1 + x 2 )) , exp(2(x 1 + x 2 )) , The exact solution of equations (48)- (50) is N B + N I = 80 + 400 points and τ = 0.1 time step are used in the computation of the test problem. Also, in order to accelerate the convergence the relaxation parameters are chosen as σ 1 = σ 2 = 1/9, and σ 3 = 7/9 which are determined by trial and error. These selected values of the relaxation parameters are used for the rest of this article. In Table 1, the L ∞ and L 2 -errors with the related convergence orders computed at the time t = 1 and with different time step are presented. From the table it can be concluded that the method is of secondorder accuracy. The exact and DRBEM solutions fort = 5 andt = 7 are drawn in the left and middle columns of Figure 1, respectively. Also, the computational results and exact values at the points located x 2 = x 1 are presented at the right column of Figure 1 for the same time levels. From these figures, it is observed that exact  solutions agree very well with the computational results obtained by DRBEM as time advances.

Example 1: Superposition of two line solitons
In the first problem, the numerical solution of Equation (1) is obtained for χ(x 1 , x 2 ) = 1 and ν = 0 such that this case is called superposition of two line solitons. Initial conditions are considered as The numerical results are given at t = 0, 3, 7, 10, 15 and 20 in Figures 2 and 3. The computations are done by using N B + N I = 80 + 400 and τ = 0.05. Initially (t = 0), there are two orthogonal line solitons  occur which are parallel to the diagonal x 2 = −x 1 . As t increases, these two orthogonal line solitons move away from each other in the direction of x 2 = x 1 without losing their form. When t = 7, a deformation is observed on the graphs. All behaviours seen in the graphs are in agreement with those published in [27,32,38,43,47] and [49]. In addition, the effect of dissipative term is evaluated at t = 2 and t = 4 when ν = 0 and ν = 0.5 by giving the contourlines. The results are presented in Figure 4, for ν = 0 and ν = 0.5, by using full and dashed contours, respectively. From the figures, the delaying effect of dissipative term on the propagation of the solitons is seen clearly. As time increases, this delaying effect of dissipative term becomes more evident. This conclusion agrees with [47] and [38].

Example 2: Perturbation of a line soliton
In this example, the proposed method is used for the perturbation of a single soliton which is calculated for χ(x 1 , x 2 ) = 1 and ν = 0. Initial conditions are considered as − 2 sech (x 2 − 7))), The numerical results are obtained by using N B + N I = 100 + 625, τ = 0.5 and are presented in terms of sin(z/2) at t = 0, 2, 4, 7, 10 and 14 in Figure 5. It is seen that two symmetrical dents occur at t = 0 and as time increases, these dents move towards to each other. At t = 7, they reach each other and an interaction occurs between them. After that they are moving away from each other. Looking at the surfaces given fort = 10 and t = 14, it can be concluded that there is no change on the shape of dents after the collision. If the obtained results are compared with the results in [32,38,47] and [49], it can be seen that the agreement of these results is very well.

Example 3: Line soliton in an inhomogeneous medium
In the next problem, the numerical solution of Equation (1) is obtained for χ(x 1 , x 2 ) = 1 + sech 2 x 1 2 + x 2 2 and ν = 0.05 which is called perturbation of a line soliton. Initial conditions are considered as   6 and 7 show the results of the problem at time levels t = 0, 6, 12, 18 in terms of sin(z/2) using N B + N I = 92 + 529 and τ = 0.4 for ν = 0.05. At t = 0, the soliton is in the form of a straight line which is parallel to the x 2 axis. As time increases, it moves slowly in x 1 direction and a curl caused by inhomogeneity of the medium occurs in its straightness. These results show that a good agreement is obtained with those given in [27,32,43,47] and [49]. Also when these results are compared with the results given for ν = 0 in [38], it can be expressed that the propagation of the line soliton is slowing down.

Example 4: Circular ring soliton
The circular ring solitons for the case χ(x 1 , x 2 ) = 1 and ν = 0.05 are analysed with the initial conditions as w 1 (x 1 , x 2 ) = 4 tan −1 exp 3 − x 1 2 + x 2 2 , Numerical computations are carried out by using N B + N I = 92 + 529 and τ = 0.4 and obtained results are evaluated by giving surfaces and corresponding contourlines at t = 0, 2.8, 5.6, 8.4 11.2 and 12.8 in terms of sin(z/2) in Figures 8 and 9, respectively. From the figures, it is observed that at the initial case (t = 0), a single ring soliton occurs and this soliton shrinks until t = 2.8. From t = 5.6 to t = 8.4, an expansion occurs on the ring soliton and this expansion phase starts with the radiation and oscillations. Finally att = 11.2 and t = 12.8, a single ring soliton appears in the cavity which has almost the same form with the initial phase. These results are good agreement with the results given in [27,32,38,43,47] and [49].

Example 5: Elliptic ring soliton
In this example elliptic ring soliton with χ(x 1 , x 2 ) = 1 and ν = 0 is considered with initial conditions The DRBEM results are obtained by using N B + N I = 88 + 484 and τ = 0.4. Figures 10 and 11 present the results for the elliptic ring soliton as surface and corresponding contour lines, respectively, at t = 0, 1.6, 3.2, 4.8, 6.4, 8.0, 9.6 and 11.2. The elliptical circular soliton exhibits shrinkage and expansion behaviour as time progresses. A permanent change in the orientation of the major and minor axes of the ellipse creates wave pulsation. Moreover, it is observed that at t = 8, the elliptical ring soliton begins to take the form of a circular soliton, and when the time is reached t = 11.2, an almost circular ring soliton is formed. The obtained results are in good agreement with those presented in [27,38,43,47] and [50].

Conclusion
In this study, we have employed DRBEM used fundamental solution of modified Helmholtz equation. The proposed method has been applied to the examples of sine-Gordon equation which have both line and ring soliton solutions. The main advantage of the DRBEM is to solve the problem by using small number of points since it does not need to discretize the whole domain. Also, in this study, since the DRBEM has been constructed by using the fundamental solution of modified Helmholtz equation, it has been possible to use more information of the original equation. The obtained results have been compared with the corresponding results for damped and undamped SG equations in [27], [32,38,43,47,49] and [50]. These results reveal that good results are obtained for both line and ring solitons and the behaviour of sine-Gordon equation is very well captured with DRBEM. Therefore, it can be concluded that proposed method is applicable to find numerical solution of sine-Gordon equation.

Disclosure statement
No potential conflict of interest was reported by the author(s).