Boundary element solution of the two-dimensional sine-Gordon equation using continuous linear elements

ARTICLE IN PRESS Engineering Analysis with Boundary Elements 33 (2009) 12– 24

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Boundary element solution of the two-dimensional sine-Gordon equation using continuous linear elements Davoud Mirzaei, Mehdi Dehghan  Department of Applied Mathematics, Faculty of Mathematics and Computer Science, Amirkabir University of Technology, No. 424, Hafez Ave., 15914 Tehran, Iran

a r t i c l e in f o

a b s t r a c t

Article history: Received 8 November 2007 Accepted 12 March 2008 Available online 20 June 2008

This article studies the boundary element solution of two-dimensional sine-Gordon (SG) equation using continuous linear elements approximation. Non-linear and in-homogenous terms are converted to the boundary by the dual reciprocity method and a predictor–corrector scheme is employed to eliminate the non-linearity. The procedure developed in this paper, is applied to various problems involving line and ring solitons where considered in references [Argyris J, Haase M, Heinrich JC. Finite element approximation to two-dimensional sine-Gordon solitons. Comput Methods Appl Mech Eng 1991;86:1–26; Bratsos AG. An explicit numerical scheme for the sine-Gordon equation in 2+1 dimensions. Appl Numer Anal Comput Math 2005;2(2):189–211, Bratsos AG. A modified predictor–corrector scheme for the two-dimensional sine-Gordon equation. Numer Algorithms 2006;43:295–308; Bratsos AG. The solution of the two-dimensional sine-Gordon equation using the method of lines. J Comput Appl Math 2007;206:251–77; Bratsos AG. A third order numerical scheme for the twodimensional sine-Gordon equation. Math Comput Simul 2007;76:271–8; Christiansen PL, Lomdahl PS. Numerical solutions of 2+1 dimensional sine-Gordon solitons. Physica D: Nonlinear Phenom 1981;2(3):482–94; Djidjeli K, Price WG, Twizell EH. Numerical solutions of a damped sine-Gordon equation in two space variables. J Eng Math 1995;29:347–69; Dehghan M, Mirzaei D. The dual reciprocity boundary element method (DRBEM) for two-dimensional sine-Gordon equation. Comput Methods Appl Mech Eng 2008;197:476–86]. Using continuous linear elements approximation produces more accurate results than constant ones. By using this approach all cases associated to SG equation, which exist in literature, are investigated. & 2008 Elsevier Ltd. All rights reserved.

Keywords: Sine-Gordon equation Soliton Boundary element method Dual reciprocity method Continuous linear elements

1. Introduction The boundary element method (BEM) is a modern numerical technique, which has enjoyed increasing popularity over the last three decades, and it has become one of the most general and effective numerical tools for solving different engineering problems [1–7]. The main idea in this method is to convert the original partial differential equation to an equivalent boundary integral equation by using Green’s theorem and a fundamental solution. Consequently, the main advantage in this method over the classical domain methods such as finite-element method and finite-difference method, is that only boundary discretization is required due to dimension reduction. But the method becomes less attractive when dealing with non-homogeneous and nonlinear problems and more general linear PDEs that do not have an explicit fundamental solution. The dual reciprocity method

 Corresponding author.

E-mail addresses: [email protected] (D. Mirzaei), [email protected] (M. Dehghan). 0955-7997/$ - see front matter & 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2008.03.011

(DRM), is a good tool for dealing with an equation in BEM for which a complete fundamental solution is either unavailable or inconvenient, leading to one or more terms remaining as domain integrals after the usual operations on the boundary integral equation. The main idea behind this approach is to expand the inhomogeneous and non-linear terms in terms of their values at the nodes, which lie in domain and boundary. The inhomogeneous term is approximated by interpolation in terms of some well-known functions called radial basis functions, where r is the distance between a source point and the field point. Some criteria for selecting this approximation functions can be found in [8]. The method was first introduced by Nardini and Brebbia [9], Brebbia and Nardini [10] and Partridge and Brebbia [11]. By applying the DRM, the problem will be reduced to a boundary only formulation, thus we do not have any domain integration in the boundary integral equation. The DRBEM was originally proposed for the numerical solution of dynamic problems in solid mechanics, but the method has now been successfully extended to a wide range of problems in engineering, such as those involving diffusion processes, inhomogeneous media and non-linearity. For some examples of these problems, see Refs. [12–18] and other

ARTICLE IN PRESS D. Mirzaei, M. Dehghan / Engineering Analysis with Boundary Elements 33 (2009) 12–24

references therein. In this report, we adapt this approach for twodimensional sine-Gordon (SG) equation. We employ the continuous linear elements to approximate the unknown function and its derivative over each element. In this method, to impose the essential (Dirichlet) boundary condition (BC) at corners we have a major problem due to the changing of gradient direction [5]. In this equation, the BC is natural (Neumenn); therefore, this problem does not arise in numerical implementation. Over the last decade, some attentions have been paid to models, which possess soliton-like structures in higher dimensions. Soliton solutions have been found in a variety of non-linear differential equations such as the Korteweg and de Vries equation [41], the Schrodinger equation [42], the SG equation, etc. Physical applications of solitons have been found among others in shallowwater waves, optical fibers, Josephson-junction oscillators, etc. In particular, the Josephson junction model that consists of two layers of superconducting material separated by an isolating barrier. This model is found to have many applications in electronics and can be described by the two-dimensional undamped SG equation. Moreover, it is found to possess solitonlike solutions [19–25]. Consider the two-dimensional SG equation: q2 u qu q2 u q2 u ¼ 2 þ 2  fðx; yÞ sinðuÞ; þb 2 qt qx qy qt

ðx; yÞ 2 O;

tX0,

(1.1)

where O ¼ {(x,y)|apxpb,cpypd}. The BCs associated with (1.1) impose a zero gradient along the boundary G of O: qu ðx; y; tÞ ¼ 0; qt

ðx; yÞ 2 G;

tX0,

(1.2)

while the initial conditions are given by uðx; y; 0Þ ¼ f ðx; yÞ;

qu ðx; y; 0Þ ¼ gðx; yÞ; qt

ðx; yÞ 2 O.

(1.3)

The function f(x, y) can be interpreted as a Josephson current density, and f(x, y) and g(x, y) are wave modes or kinks and velocity, respectively. The parameter b is the so-called dissipative term, which is assumed to be a real number with bX0. When b ¼ 0, Eq. (1.1) reduces to the undamped SG equation in two space variables, while when b 40, to the damped one. For the undamped SG equation in higher dimensions, the exact soliton solutions have obtained in [26] using Hirota’s method, in [27,28] using Lamb’s method, in [29] by Painleve transcendents and in [30] by Blacklund transformation, etc. Numerical solutions for the SG equation have given by Guo et al. [31] using two difference schemes, Christiansen and Lomdahl [24] using a generalized leapfrog method, Argyris et al. [19] by the finite-element approach, Sheng et al. [32] by a split cosine scheme, Djidjeli et al. [25] using a two-step one-parameter leapfrog scheme, Bratsos [20] using a three time-level fourth-order explicit finite-difference scheme, Bratsos [21] using a modified predictor–corrector scheme, Bratsos [22] using the method of lines, Bratsos [23] by a third-order numerical scheme and Dehghan and Mirzaei [33] by the constant boundary element approximation. The organization of the rest of this article is as follows: In Section 2, the boundary element formulation of SG equation is presented, a time stepping method and a predictor–corrector scheme are described for the time derivatives and non-linearity, respectively. In Section 3, numerical results for some problems, involving line and ring solitons, are investigated and the reformations of the solitons from initial condition to final time are given. Also to test the accuracy of the proposed method the value of energy at several time steps is computed. Finally in Section 4, the report ends with a brief conclusion.

13

2. Formulations 2.1. Boundary element formulation We present a description of the dual reciprocity BEM formulation for SG equation using continuous linear element approximation. In elements of this kind, it is assumed that the variables of the integral equation have a linear evaluation along the elements. These will be segments of a straight line, the linear evolution being expressed through the values of the functions at the end of the elements. We start with the following weighted residual statement for the approximate solution of Eq. (1.1): # Z " 2 q u qu 2  r u þ fðx; yÞ sinðuÞ un dO ¼ 0. þb (2.1) 2 qt O qt Applying Green’s theorem, the integral formulation of Eq. (2.1) is concluded by  Z " 2 Z  qun q u qu n qu dG ¼ u u þb ci ui þ 2 qn qt qn qt G O þ fðx; yÞ sinðuÞun dO,

(2.2)

where ci are coefficients depending on the location of the point i, and weighting function u called fundamental solution is given by un ¼

1 ln r. 2p

(2.3)

The above integral equation contains a domain integral corresponding to the time derivatives and non-linear terms. To transfer this domain integral to the boundary, we introduce the producttype approximation of the non-homogeneous terms by a set of coordinate functions mk given by NþL X q2 u qu þ fðx; yÞ sinðuÞ ¼ bðx; y; tÞ ¼ þb mk ðx; yÞak ðtÞ, (2.4) 2 qt qt k¼1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where mk(x, y) ¼ 1+rk and r k ¼ ðx  xk Þ2 þ ðy  yk Þ2 for k ¼ 1, 2,y, N+L. ak(t) are unknown functions of time and N and L represent the number of functions (collocation points) on the boundary and domain, respectively. In fact, we determine N boundary elements on G and take the first N collocation points to be the end of these boundary elements and the rest of L collocation points are chosen in the domain of the problem. The collocation points are denoted by (xj, yj) for j ¼ 1, 2,y,N+L. Functions mk are so chosen that we can easily find u^ k satisfying

r 2 u^ k ¼ mk .

(2.5)

With this approximation, the domain integral in Eq. (2.2) is converted to an equivalent boundary integral by applying the same idea as for the main equation, that is Z

bun dO ¼

O

NþL X k¼1

¼

Z ak ðtÞ

ðr 2 u^ k Þun dO

O

   Z  qun qu^  un k dG , ak ðtÞ ci u^ ik þ u^ k qn qn G k¼1

NþL X

(2.6)

where u^ ik is the known value of u^ k at the source point (xi, yi). Combining Eqs. (2.2) and (2.6), yields  Z  qun qu ci ui þ  un k dG u qn qn G    Z  N þL X qun qu^ ¼  un k dG . (2.7) ak ðtÞ ci u^ ik þ u^ k qn qn G k¼1 If we represent the value of the function mk at collocation point i by mik for i ¼ 1, 2,y,N+L, and set F as a (N+L)  (N+L) matrix that

ARTICLE IN PRESS 14

D. Mirzaei, M. Dehghan / Engineering Analysis with Boundary Elements 33 (2009) 12–24

F(i, k) ¼ mik and E ¼ F1, then we have [6]

ak ðtÞ ¼

NþL X

Eq. (2.10) can be written in the compact form

Ekj bj ðtÞ,

(2.8)

N h X

ci ui ðtÞ þ

where bj(t) ¼ b(xj, yj, t). The discretization of the integral equation (2.7) has the form N X j¼1

Z

uj ðtÞqn dG 

¼ Gj

N Z X j¼1

8 < N þL N X X ak ðtÞ ci u^ ik þ ¼ :

Gm ij ¼ (2.9)

^ It must be observed where q ¼ qu/qn, q ¼ qu/qn and q^ ¼ qu=qn. that the first integrals appearing in left- and right-hand sides of Eq. (2.9) have the symbol of finite part of the integral. The meaning of this symbol being only formal in any case when the collocation point i belongs to the element j (i ¼ j or j+1). This sense of finite part, instead of the sense of Cauchy principle value (CPV), is required due to the positioning of the collocation point at the extreme of the elements. In these circumstances, the integrals along each of the two elements to which one collocation point belongs must be considered as finite part, the set of the two integrations maintaining the sense of CPV that the integral along the whole actual boundary has [5]. In accordance with the linear evaluation of the variables to be taken along the elements, they can be expressed along the generic element j by introducing natural coordinates. Owing to the fact that the variables involving in the formulation can be related to nodes and elements, variable u (or q) will be represented in the form uelement ðor qelement Þ. This notation node node makes the possibility to distinguish between the values of a function at a node, which now belongs to two adjacent elements. But since the function u is a scaler function and not dependent on the direction, its values at source point j are equal for elements j and j+1, but they may be different for q. Therefore, from Eq. (2.9) we obtain N1 ðxÞ

j¼1 G j N Z  X j¼1

"

N1 ðxÞ

Gj

" #  uj ðtÞ Lj N 2 ðxÞ qn dx 2 ujþ1 ðtÞ 2 3 j 6 qj ðtÞ 7 n Lj N 2 ðxÞ 4 dx 5u 2 qjjþ1 ðtÞ

8 N þL N < X X ¼ ak ðtÞ ci u^ ik þ :



# 

ujþ1 ðtÞ

N h X

G1ij

G2ij

i

j¼1

2 6 4

qjj ðtÞ qjjþ1 ðtÞ

3 7 5

Sik ak ðtÞ,

(2.12)

N 1 ðxÞ

N2 ðxÞ



"

Gj

u^ j

# qn

Lj dx 2

(2.10)

where x is the natural coordinate along the elements taking 71 at the end edges and N2 ðxÞ ¼ 12ðx þ 1Þ,

N m ðxÞqn

Lj dx; 2

m ¼ 1; 2,

(2.13)

Nm ðxÞun

Lj dx; 2

m ¼ 1; 2

(2.14)

Gj

Z Gj

and

Sik ¼ ci u^ ik þ

N h X

H1ij

H2ij

i

"

u^ j

# 

u^ jþ1

j¼1

N h X

G1ij

G2ij

j¼1

2 3 i q^ jj 4 5. j q^ jþ1

(2.15)

From (2.8), the right-hand side of Eq. (2.12) is NþL X

Sik ak ðtÞ ¼

k¼1

N þL X k¼1

Sik

NþL X

Ekj bj ðtÞ ¼

j¼1

N þL X

M ij bj ðtÞ,

(2.16)

j¼1

where Mij ¼

N þL X

Sik Ekj .

(2.17)

k¼1

Combining Eqs. (2.12) and (2.16), leads to

ci ui ðtÞ þ

N h X

H1ij

j¼1

¼

N þL X

H2ij

i

"

uj ðtÞ ujþ1 ðtÞ

M ij bj ðtÞ.

# 

N h X j¼1

G1ij

G2ij

i

2 6 4

qjj ðtÞ qjjþ1 ðtÞ

3 7 5 (2.18)

j¼1

This equation can be written for every node i of the N+L nodes on the boundary and domain. The above equations are the discretized version of the SG equation using the continuous linear elements approximation and DRM. Details on the computation of the integrals over Gj in Eqs. (2.13) and (2.14) may be found in Refs. [4,5]. 2.2. Application of the BC

u^ jþ1 j¼1 k¼1 9 2 3 j > ^ = N Z  X  qj 7 n Lj N1 ðxÞ N 2 ðxÞ 6 dx ,  4 j 5u > 2 ; q^ jþ1 j¼1 Gj

N1 ðxÞ ¼ 12ðx  1Þ;

uj ðtÞ

k¼1

Hm ij ¼

Gj

j¼1



i

where

9 Z N Z = X n n ¼ u q^ k dG , u^ k q  ; Gj Gj

N I  X

¼

N þL X

un qðtÞ dG

j¼1

k¼1

ci ui ðtÞ þ

H2ij

j¼1

j¼1

ci ui þ

H1ij

(2.11)

are linear shape functions. In Eq. (2.10), dG is expressed in the natural coordinate as Lj/2 dx. The integrations along all the elements are extended due to the use of natural coordinate x to the normalized interval (1, +1).

It must be pointed out that the number of variables associated to a node j is initially 3, because although only one value of the potential uj is associated to node j, but two fluxes qjj1 and qjj corresponding to those existing at the ends of elements j1 and j, which can be different, are associated to it. It is a problem for applying the BC, when the boundary is not smooth at node j and the BC before node j (at element j1) and after it (at element j) are Dirichlet. Because in this case there are two unknown variables associated to node j. For dealing with this problem, an extra equation is required or else a direct relation between two variables, or between them and third variable must be introduced. Also, the discontinuous linear element approximation can be used. The above problem does not arise in this equation due to the Neumann type BC (1.2). Using the BC the second summation in left-hand side of Eq. (2.18), equals to zero and we have AuðtÞ ¼ MbðtÞ,

(2.19)

ARTICLE IN PRESS D. Mirzaei, M. Dehghan / Engineering Analysis with Boundary Elements 33 (2009) 12–24

where u(t) ¼ [u1(t), u2(t),y,uN+L(t)]T, b(t) ¼ [b1(t), b2(t),y,bN+L(t)]T and Aij ¼ H1ij þ H2ij1 þ ci dij ,

(2.20)

where d is kronecker symbol. Now, (2.19) constitutes a system of (N+L) equations in (N+L) unknown functions of t. The unknown variables are given by uj(t), for j ¼ 1, 2,y,N+L. This system is solved approximately using the iterative scheme, which is described in the following. 2.3. Time difference formulation and a predictor– corrector scheme

15

distribution of nodes is regular in all examples. The routine has written in MATLAB 7.1 and is run on a Pentium 4 PC, with 1 GB RAM. It is known (see for example [20–23,25,31]) that, when b ¼ 0, for the SG equation the energy given by the following expression: ZZ 1 E ¼ EðtÞ ¼ ½u2x þ u2y þ u2t þ 2ð1  cos uÞ dx dy, (3.1) 2 is conserved. We also investigate this property for SG equation. The evaluation of E is performed using the composite trapezoidal rule for integration. 3.1. Line solitons

We make the approximations q2 u uðn1Þ  2uðnÞ þ uðnþ1Þ ’ , qt 2 ðDtÞ2 qu uðn1Þ  uðnþ1Þ ’ , qt 2Dt

(2.21)

uðtÞ ’ Z1 uðn1Þ þ Z2 uðnÞ þ Z3 uðnþ1Þ , where t ¼ nDt, Z1+Z2+Z3 ¼ 1 and Z1, Z2, Z3X0. If x ¼ [x1, x2,y,xn+L] ~ ¼ fðx; yÞ  sinðuÞ, ~ where and y ¼ [y1, y2,y,yN+L]T, by assuming cðuÞ ~ is given by the known approximation of u(x, y, t), as described u earlier, and l ¼ 1/(Dt)2 and g ¼ b/(2Dt), then from (2.19) and (2.21) we have ½Z3 A  ðl þ gÞMuðnþ1Þ ¼ ½2lM  Z2 AuðnÞ ~ þ ½ðl þ gÞM  Z1 Auðn1Þ þ McðuÞ.

(2.22)

(0)

At the first time level, when n ¼ 0, u ¼ f(x, y) and u(1) ¼ u(1)2Dtg(x, y) are determined by the initial conditions (1.3). In each time level at first we put ~ ¼ cðuðnÞ Þ. cðuÞ Having this, Eq. (2.22) is solved as a system of linear algebraic equations for unknowns u(n+1), assuming that u(n1) and u(n) are known, from the two previous time levels. Recompute ~ ¼ cðZ1 uðn1Þ þ Z2 uðnÞ þ Z3 uðnþ1Þ Þ, cðuÞ where u(n+1) is obtained recently from solving Eq. (2.22). We ~ and solving the approximately iterate between calculating cðuÞ values of the unknowns, until all the unknowns quantities converge to within a prescribed number of significant figures, i.e. a predictor–corrector approach is adopted in each time level. Once the prescribed convergence is achieved, we can move on to the following time level. The resulting solutions from the current time level provide the known values for the next iteration. This process is iterated, until reaching to the desirable time t.

3. Numerical results and discussions The above boundary element procedure is applied to various cases involving two dimensional line and ring solitons. The examples are chosen for comparison with the results of [19–24,33] for undamped and damped equations. We need to ~ to recompute c and iterate between finding an estimation for u solving the system (2.22) for a new u that described in the previous section. In all the cases considered, the iteration was stopped when the absolute values of all the unknowns in (2.22) from two consecutive iterations differ by less than 105. Convergence was achieved after less than five iterations in all cases tested. Also in all problems Dt ¼ 0.2 and Z1 ¼ Z2 ¼ Z3 ¼ 1/3 are chosen. The numbers of boundary elements and interior collocation points have stated in each cases. Note that, the

3.1.1. Superposition of two orthogonal line solitons The superposition of two orthogonal line solitons will be obtained for the case of f(x, y) ¼ 1 and b ¼ 0, with initial conditions f ðx; yÞ ¼ 4 tan1 expðxÞ þ 4 tan1 expðyÞ,

(3.2)

gðx; yÞ ¼ 0,

(3.3)

over the region 6px, yp6. The results are presented in Fig. 1, where the initial condition at time t ¼ 0 and the numerical solutions at t ¼ 4 are shown. The results in Fig. 1 show the break up of two orthogonal line solitons which are parallel to the diagonal y ¼ x and are moving away from each other in the direction of y ¼ x, undisturbed. The comparison between the separations at t ¼ 3, 4 and 7 are also depicted in Fig. 2. It is maintained from the relevant contours that until t ¼ 4 this separation occurs without any deformation, while at t ¼ 7 a deformation has already appeared. The present solutions are in good agreement with the corresponding results of [19–24,33]. Also in order to examine the behavior of the dissipative term, the solution is evaluated at t ¼ 3 when b ¼ 0, 0.5, 1.5. The resulting solitons are shown in Fig. 3 when full contours show the solitons when b ¼ 0, while the dashed ones when b ¼ 0.5 and 1.5. As can be seen from this figure, the presence of the dissipative term delays the propagation of the solitons. This delay becomes obvious when the larger value of b (b ¼ 1.5) is used. This conclusion agrees with [22]. Let us check the conservation of energy in this case. Table 1 presents the values of E(t) when b ¼ 0 and 0.5 at some selected times t up to t ¼ 7. It is worth to point out that these solutions are calculated using 192 boundary elements and 2209 interior collocation points. We note that, using the constant boundary element [33] did not produce accurate results for t45. 3.1.2. Perturbation of a line soliton A static line soliton is perturbed to produce two symmetric dents moving towards each other with a constant unit velocity. According to [19,22,24], the dents collide and continue with the same velocity and no shift occurs. Perturbation of a single soliton has been calculated for b ¼ 0 and f(x, y) ¼ 1, in terms of sin(u/2) at t ¼ 2, 7, 11, with the initial conditions f ðx; yÞ ¼ 4 tan1 expðx þ 1  2 sechðy þ 7Þ  2 sechðy  7ÞÞ,

(3.4)

gðx; yÞ ¼ 0,

(3.5)

over the region 7px, yp7. The results in Fig. 4 show, two symmetric dents moving toward each other, collapsing at t ¼ 7 and continuing to move away from each other thereafter. It can be deduced that after the collision the dents retain their shape, which verifies the conclusions of [19,22,24,25,33]. Table 2 presents the values of E(t) at some selected times t, that shows the energy remains constant as time increases.

ARTICLE IN PRESS 16

D. Mirzaei, M. Dehghan / Engineering Analysis with Boundary Elements 33 (2009) 12–24

Fig. 1. Superposition of two orthogonal line solitons: solution at t ¼ 0, 4.

Fig. 2. Comparison of superposition of two orthogonal line solitons at t ¼ 3, 4, 7.

Fig. 3. Comparison of superposition of two orthogonal line solitons at t ¼ 3 for b ¼ 0, 0.5, 1.5.

ARTICLE IN PRESS D. Mirzaei, M. Dehghan / Engineering Analysis with Boundary Elements 33 (2009) 12–24

Note that these solutions were calculated using 224 boundary elements and 3025 interior collocation points.

3.1.3. Line soliton in an inhomogeneous medium A model for an inhomogeneity on large-area Josephson junction is given by the Josephson current density fðx; yÞ ¼ 1 þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2 sech x2 þ y2 and initial conditions

x  3:5 , f ðx; yÞ ¼ 4 tan1 exp 0:954

(3.6)



x  3:5 , gðx; yÞ ¼ 0:629 sech 0:954

(3.7)

over the region 7px, yp7 for b ¼ 0. Numerical results are presented in Fig. 5, in term of sin(u/2), at t ¼ 0, 6, 12 and 18. The results show that the line soliton is moving in direction x as a straight line during the transmission through inhomogeneity. As can be seen from Fig. 5 when t tends to 12 a deformation in its straightness appears. After t ¼ 12 until t ¼ 18 this movement seems to be prevented, while when t ¼ 18 the soliton recovers its straightness. Christiansen and Lomdahl [24], claimed that this was due to the BC used. This phenomenon was observed in [22,24,25,33] and was not observed in [19,20]. The results were calculated using 224 boundary elements and 3025 interior collocation points.

Table 1 The energy of the superposition of two orthogonal line solitons when tA[0,7]

3.2. Ring solitons 3.2.1. Circular ring soliton The behavior of a circular ring quasi-soliton arising from the SG equation was investigated numerically in [34,35], who named these waves pulsons because of their pulsating behavior. Then some studies concerning the behavior of these quasi-soliton waves have been published in [36–40] and, etc. [21,22]. In Fig. 6, circular ring solitons for the case of f(x, y) ¼ 1, b ¼ 0 and initial conditions qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

f ðx; yÞ ¼ 4 tan1 exp 3  x2 þ y2 , (3.8) gðx; yÞ ¼ 0

(3.9)

over the region 7px, yp7 at time t ¼ 0 (initial condition) and the numerical solutions at t ¼ 2.8, 5.6, 8.4, 11.2 and 12.6 are shown in terms of sin(u/2) for corresponding contours. The soliton from its initial position, shrinks until t ¼ 2.8 appearing as a singlering soliton. From t ¼ 5.6, which could be considered as the beginning of the expansion phase, a radiation appears, which is followed by oscillations at the boundaries. This expansion is continued until t ¼ 11.2, where the soliton is almost reformed. Finally, since t ¼ 12.6 it appears to be again in its shrinking phase. During all the above transformations no displacement of the center of the soliton occurred. These behaviors have also appeared in [22,23]. Let e(t) ¼ |E(t)E(0)|. Table 3 presents the values of e at some selected times t, which shows the conservation of energy when dissipative term is zero. Table 2 The energy of the perturbation of a soliton when tA[0,7]

b

Initial

E(1.0)

E(2.0)

E(3.0)

E(4.0)

E(7.0)

0.0 0.5

175.5745 175.5745

175.2215 172.7906

175.1703 165.9170

175.1538 158.1395

175.1369 150.2747

175.2949 123.9263

Time t

Initial

t¼2

t¼5

t¼7

t¼9

t ¼ 11

E(t)

122.6648

121.9386

122.1123

122.1691

122.0431

121.9471

t=2

1

sin (u/2)

sin (u/2)

t = 0 (initialcondition)

0.5 0

1 0.5 0

5 y

5

5

0 -5

0 -5

y

x

5

0

0 -5

sin (u/2)

sin (u/2)

1 0.5 0 5 5

0 -5

-5

-5

x

t = 11

t=7

y

17

0 x

1 0.5 0 5 y

5

0 -5

Fig. 4. Perturbation of a line soliton: solution at t ¼ 0, 2, 7, 11.

-5

0 x

ARTICLE IN PRESS 18

D. Mirzaei, M. Dehghan / Engineering Analysis with Boundary Elements 33 (2009) 12–24

Fig. 5. Line soliton in an inhomogeneous medium: solution at t ¼ 0, 6, 12, 1.

The present solutions are in good agreement with the corresponding results of [19–23]. The computations were performed using 224 boundary elements and 3025 interior collocation points. 3.2.2. Elliptic ring soliton Elliptical ring solitons are obtained for f(x, y) ¼ 1 and b ¼ 0 with initial conditions qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

f ðx; yÞ ¼ 4 tan1 exp 3  ðx  yÞ2 =3 þ ðx þ yÞ2 =2 , (3.10) gðx; yÞ ¼ 0,

(3.11)

over the region 7px, yp7. The contour plots at times t ¼ 0, 4.8, 6.4, 8.0, 9.6 and 11.2 are shown in Fig. 7, in terms of sin(u/2). The temporal behavior of the soliton consists of a shrinking and an expanding phase. The associated pulsation of the wave is accompanied by a permanent change in the orientation of the major and minor axes of the ellipse. At t ¼ 11.2, an essentially circular ring soliton is seen [19,22]. These results are in good agreement with those given in [19,22,33]. The results were performed using 224 boundary elements and 3025 interior collocation points. 3.2.3. Elliptical breather An elliptical breather is obtained for f(x, y) ¼ 1 with initial conditions f ðx; yÞ

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ¼ 4 tan1 exp 2:0 sech 0:866 ðx  yÞ2 =3 þ ðx þ yÞ2 =2 , (3.12)

gðx; yÞ ¼ 0,

(3.13)

over the region 7px, yp7. In Fig. 8, the contours of the elliptical breather at t ¼ 0, 1.6, 8.0, 9.6, 11.2, 12.8, 14.4 and 15.2 are given when b ¼ 0. The major axes of the breather from its initial direction y ¼ x seems to be turned clockwise shrinking until t ¼ 1.6, while at t ¼ 8, 9.6 a reflexion phase is observed. At t ¼ 11.2, the major axes has almost recovered its initial direction y ¼ x, but strong oscillations are observed [22,23]. After t ¼ 12.8, an expansion phase is observed. The graphs are in agreement with those given in [19,20,22,23]. As previous, Table 4 shows that the energy remains constant as time increases. The results were calculated using 224 boundary elements and 3025 interior collocation points.

3.2.4. Collision of two circular solitons The collision of two expanding circular ring solitons is considered with f(x, y) ¼ 1 and initial conditions  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

1 f ðx; yÞ ¼ 4 tan1 exp 4  ðx þ 3Þ2 þ ðy þ 7Þ2 , 0:436

(3.14)

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

1 gðx; yÞ ¼ 4:13 sinh 4  ðx þ 3Þ2 þ ðy þ 7Þ2 , 0:436

(3.15)

over the region 30pxp10 and 21pyp7 using 170 boundary elements and 1666 interior collocation points only in region 10pxp10 and 7pyp7 and then expanding across x ¼ 10 and y ¼ 7 by the symmetry relations. Fig. 9 depicts the surface plots of the collision of two expanding circular ring solitons at t ¼ 0, 4, 8 and 11 in terms of sin(u/2) when b ¼ 0. The results are in good agreement with those given in [19,22].

ARTICLE IN PRESS D. Mirzaei, M. Dehghan / Engineering Analysis with Boundary Elements 33 (2009) 12–24

t = 2.8

6

6

4

4

2

2

0

0

y

y

t = 0 (initial condition)

-2

-2

-4

-4

-6

-6 -6

-4

-2

0 x

2

4

6

-6

-4

-2

6

4

4

2

2

0

0

y

y

6

-2

-2

-4

-4

-6

-6 -4

-2

0 x

2

4

6

-6

-4

-2

6

6

4

4

2

2

0

0

-2

-2

-4

-4

-6

-6 -4

-2

0 x

2

4

6

0 x

2

4

6

2

4

6

t = 12.6

y

y

t = 11.2

-6

0 x t = 8.4

t = 5.6

-6

19

2

4

6

-6

-4

-2

0 x

Fig. 6. Circular ring solitons: solution at t ¼ 0, 2.8, 5.6, 8.4, 11.2, 12.6.

3.2.5. Collision of four circular solitons The collision of four expanding circular ring solitons is considered for f(x, y) ¼ 1 and initial conditions

Table 3 The error e(t) of the circular ring soliton when E(0) ¼ 150.4599 Time t

t ¼ 2.8

t ¼ 5.6

t ¼ 8.4

t ¼ 11.2

t ¼ 12.6

Error e(t)

0.0466

0.1598

0.4421

0.3015

0.2985

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

1 4  ðx þ 3Þ2 þ ðy þ 3Þ2 , f ðx; yÞ ¼ 4 tan1 exp 0:436

(3.16)

ARTICLE IN PRESS 20

D. Mirzaei, M. Dehghan / Engineering Analysis with Boundary Elements 33 (2009) 12–24

t = 4.8

6

6

4

4

2

2

0

0

y

y

t = 0 (initial condition)

-2

-2

-4

-4

-6

-6 -6

-4

-2

0 x

2

4

6

-6

-4

-2

6

6

4

4

2

2

0

0

-2

-2

-4

-4

-6

-6 -6

-4

-2

0 x

2

4

6

-6

-4

-2

6

4

4

2

2

0

0

y

y

6

-2

-2

-4

-4

-6

-6 -4

-2

0 x

4

6

0 x

2

4

6

2

4

6

t = 11.2

t = 9.6

-6

2

t = 8.0

y

y

t = 6.4

0 x

2

4

6

-6

-4

-2

0 x

Fig. 7. Elliptic ring solitons: solution at t ¼ 0, 4.8, 6.4, 8.0, 9.6, 11.2.

gðx; yÞ ¼ 4:13= cosh



 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

1 4  ðx þ 3Þ2 þ ðy þ 3Þ2 0:436

(3.17)

over the region 30px, yp10. For numerical solution we chose 200 boundary elements and 2401 interior collocation points over one quarter of the domain. The solution is expanded

using the symmetry relations. The calculation is carried out to t ¼ 11 and the corresponding surface plots are depicted in Fig. 10 in terms of sin(u/2). The results show an extremely complex interaction with rapidly varying values of u in the center and are in good agreement with corresponding surfaces given in [19].

ARTICLE IN PRESS D. Mirzaei, M. Dehghan / Engineering Analysis with Boundary Elements 33 (2009) 12–24

t = 1.6

6

6

4

4

2

2

0

0

y

y

t = 0 (initial condition)

-2

-2

-4

-4

-6

-6 -6

-4

-2

0 x

2

4

6

-6

-4

-2

6

6

4

4

2

2

0

0

-2

-2

-4

-4

-6

-6 -4

-2

0 x

2

4

6

-6

-4

-2

6

6

4

4

2

2

0

0

-2

-2

-4

-4

-6

-6 -4

-2

0 x

2

4

6

-6

-4

-2

6

6

4

4

2

2

0

0

-2

-2

-4

-4

-6

-6 -4

-2

0 x

4

6

0 x

2

4

6

0 x

2

4

6

2

4

6

t = 15.2

y

y

t = 14.4

-6

2

t = 12.8

y

y

t = 11.2

-6

0 x t = 9.6

y

y

t = 8.0

-6

21

2

4

6

-6

-4

-2

0 x

Fig. 8. Elliptical breather: solution at t ¼ 0, 1.6, 8.0, 9.6, 11.2, 12.8, 14.4, 15.2.

ARTICLE IN PRESS 22

D. Mirzaei, M. Dehghan / Engineering Analysis with Boundary Elements 33 (2009) 12–24

4. Concluding remarks In this paper, the BEM formulated and successfully implemented for the non-linear two-dimensional SG equation. To approximate the unknown function and its derivative, the continuous linear elements employed and the DRM applied to transfer the

Table 4 The error e(t) of the elliptical breather when E(0) ¼ 89.6761 Time t

t ¼ 1.6

t ¼ 8.0

t ¼ 9.6

t ¼ 11.2

t ¼ 12.8

t ¼ 14.4

t ¼ 15.2

Error e(t)

0.0387

0.1319

0.2109

0.8745

0.1006

0.0007

0.1993

volume integrals onto the boundary. A time stepping method employed to deal with the time derivatives and a predictor–corrector scheme presented to eliminate the non-linearity. Some numerical experiments involving line and ring solitons depicted. To demonstrate the accuracy of this method the values of energy were given at several time steps. The number of boundary elements and interior points in each example is reported. The distribution of points was regular in all cases. Using linear elements approximation leads to more accurate numerical results than the constant ones. Also by this method more numerical experiments such as elliptical breather, collision of two and four circular ring solitons investigated and the corresponding table and figures presented. In fact by using this approach, all cases associated to SG equation, which exist in literature, were

Fig. 9. Collision of two line solitons: solution at t ¼ 0, 4, 8, 11.

ARTICLE IN PRESS D. Mirzaei, M. Dehghan / Engineering Analysis with Boundary Elements 33 (2009) 12–24

23

Fig. 10. Collision of four circular solitons: solution at t ¼ 0, 3, 5, 7, 9, 11.

considered. The higher order elements such as quadratic and cubic could be employed also. Finally the authors would like to point out the possible extension of the new approach to solve the two-dimensional version of the time-dependent partial differential equations studied in [43–45].

Acknowledgments The authors thank both reviewers for carefully reading the paper and for their comments and suggestions. References [1] Aliabadi MH. The boundary element method, applications in solids and structures, vol. 2. New York: Wiley; 2002. [2] Ang WT. A beginner’s course in boundary element method. Universal Publisher; 2007. [3] Brebbia CA, Walker S. Boundary element technique in engineering. London: Newnes-Butterworths; 1980. [4] Brebbia CA, Telles JCF, Wrobel LC. Boundary element techniques. Berlin, Heidelberg, New York, Tokyo: Springer; 1984. [5] Paris F, Canas J. Boundary element method: fundamental and applications. Oxford: Oxford University Press; 1997. [6] Ramachandran PA. Boundary element methods in transport phenomena. Southampton, Boston: Computational Mechanics Publication; 1994. [7] Wrobel LC. The boundary element method, applications in thermo-fluids and acoustics, vol. 1. New York: Wiley; 2001.

[8] Partridge PW. Towards criteria for selecting approximation functions in the dual reciprocity method. Eng Anal Bound Elem 2000;24:519–29. [9] Nardini D, Brebbia CA. A new approach for free vibration analysis using boundary elements. In: Brebbia CA, editor. Boundary element methods in engineering. Berlin: Springer; 1992. p. 312–26. [10] Brebbia CA, Nardini D. Dynamic analysis in solid mechanics by an alternative boundary element procedure. Int J Soil Dyn Earthquake Eng 1983;2:228–33. [11] Partridge PW, Brebbia CA. The dual reciprocity boundary element method for the Helmholtz equation. In: Brebbia CA, Choudouet-Miranda A, editors. Proceedings of the international boundary elements symposium. Berlin: Computational Mechanics Publications/Springer; 1990. p. 543–55. [12] Ang WT. The two-dimensional reaction–diffusion Brusselator system: a dualreciprocity boundary element solution. Eng Anal Bound Elem 2003;27: 897–903. [13] Ang WT, Ang KC. A dual reciprocity boundary element solution of a generalized nonlinear Schro¨dinger equation. Numer Meth Part Diff Eq 2004;20:843–54. [14] Krishna KM, Tanaka M. Dual reciprocity boundary element analysis of nonlinear diffusion: temporal discretization. Eng Anal Bound Elem 1999;23: 419–33. [15] Krishna KM, Tanaka M. Dual reciprocity boundary element analysis of inverse heat conduction problems. Comput Meth Appl Mech Eng 2001;190:5283–95. [16] Partridge PW, Sensale B. The method of fundamental solutions with dual reciprocity for diffusion and diffusion–convection using subdomains. Eng Anal Bound Elem 2000;24:633–41. [17] Tanaka M, Matsumoto T, Takakuwa S. Dual reciprocity BEM for time-stepping approach to the transient heat conduction problem in nonlinear materials. Comput Meth Appl Mech Eng 2006;195:4953–61. [18] Wrobel LC, Brebbia CA. The dual reciprocity boundary element formulation for nonlinear diffusion problems. Comput Meth Appl Mech Eng 1987;65: 147–64.

ARTICLE IN PRESS 24

D. Mirzaei, M. Dehghan / Engineering Analysis with Boundary Elements 33 (2009) 12–24

[19] Argyris J, Haase M, Heinrich JC. Finite element approximation to twodimensional sine-Gordon solitons. Comput Meth Appl Mech Eng 1991;86: 1–26. [20] Bratsos AG. An explicit numerical scheme for the sine-Gordon equation in 2+1 dimensions. Appl Numer Anal Comput Math 2005;2(2):189–211. [21] Bratsos AG. A modified predictor–corrector scheme for the two-dimensional sine-Gordon equation. Numer Algorithms 2006;43:295–308. [22] Bratsos AG. The solution of the two-dimensional sine-Gordon equation using the method of lines. J Comput Appl Math 2007;206:251–77. [23] Bratsos AG. A third order numerical scheme for the two-dimensional sineGordon equation. Math Comput Simulat 2007;76:271–8. [24] Christiansen PL, Lomdahl PS. Numerical solutions of 2+1 dimensional sineGordon solitons. Physica D: Nonlinear Phenom 1981;2(3):482–94. [25] Djidjeli K, Price WG, Twizell EH. Numerical solutions of a damped sineGordon equation in two space variables. J Eng Math 1995;29:347–69. [26] Hirota R. Exact three-soliton solution of the two-dimensional sine-Gordon equation. J Phys Soc Jpn 1973;35:1566. [27] Christiansen PL, Olsen OH. Return effect for rotationally symmetric solitary wave solutions to the sine-Gordon equation. Phys Lett A 1978;68(2):185–8. [28] Zagrodzinsky J. Particular solutions of the sine-Gordon equation in 2+1 dimensions. Phys Lett A 1979;72:284–6. [29] Kaliappan P, Lakshmanan M. Kadomtsev–Petviashvili and two-dimensional sine-Gordon equations: reduction to Painleve` transcendents. J Phys A: Math Gen 1979;12:L249–52. [30] Leibbrandt G. New exact solutions of the classical sine-Gordon equation in 2+1 and 3+1 dimensions. Phys Rev Lett 1978;41:435–8. [31] Guo BY, Pascual PJ, Rodriguez MJ, Va´zquez L. Numerical solution of the sineGordon equation. Appl Math Comput 1986;18:1–14. [32] Sheng Q, Q .M. Khaliq A, Voss DA. Numerical simulation of two-dimensional sine-Gordon solitons via a split cosine scheme. Math Comput Simulat 2005;68:355–73.

[33] Dehghan M, Mirzaei D. The dual reciprocity boundary element method (DRBEM) for two-dimensional sine-Gordon equation. Comput Meth Appl Mech Eng 2008;197:476–86. [34] Bogolyubskil˘ IL, Makhankov VG. Lifetime of pulsating solitons in certain classical models. JETP Lett 1976;24(1):12–4. [35] Bogolyubskil˘ IL. Oscillating particle-like solutions of the nonlinear Klein–Gordon equation. JETP Lett 1976;24(10):535–8. [36] Christiansen PL, Grønbech-Jensen N, Lomdahl PS, Malomed BA. Oscillations of eccentric Pulsons. Phys Scripta 1997;55:131–4. [37] Malomed BA. Decay of shrinking solitons in multidimensional sine-Gordon equation. Physica D 1987;24:155–71. [38] Malomed BA. Dynamic of quasi-one-dimensional kinks in the two-dimensional sine-Gordon model. Physica D 1991;52:157–70. [39] Maslov EM. Dynamics of rotationally symmetric solitons in near-SG field model with applications to large-area Josephson junctions and ferromagnets. Physica D 1985;15:433–43. [40] Maslov EM. Rotationally symmetric SG oscillator with tunable frequency. Phys Lett A 1988;131(6):364–7. [41] Dehghan M, Shokri A. A numerical method for KdV equation using collocation and radial basis functions. Nonlinear Dyn 2007;50:111–20. [42] Dehghan M. Finite difference procedures for solving a problem arising in modeling and design of certain optoelectronic devices. Math Comput Simulat 2006;71:16–30. [43] Dehghan M. A computational study of the one-dimensional parabolic equation subject to nonclassical boundary specifications. Numer Meth Partial Differ Equations 2006;22:220–57. [44] Dehghan M. The one-dimensional heat equation subject to a boundary integral specification. Chaos, Soliton Fract 2007;32:661–75. [45] Shakeri F, Dehghan M. Numerical solution of the Klein–Gordon equation via He’s variational iteration method. Nonlinear Dyn 2008;51:89–97.