Meshless numerical analysis of a class of nonlinear generalized Klein–Gordon equations with a well-posed moving least squares approximation

Accepted Manuscript

Meshless numerical analysis of a class of nonlinear generalized Klein-Gordon equations with a well-posed moving least squares approximation Xiaolin Li PII: DOI: Reference:

S0307-904X(17)30258-5 10.1016/j.apm.2017.03.063 APM 11704

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

29 July 2015 19 February 2017 28 March 2017

Please cite this article as: Xiaolin Li, Meshless numerical analysis of a class of nonlinear generalized Klein-Gordon equations with a well-posed moving least squares approximation, Applied Mathematical Modelling (2017), doi: 10.1016/j.apm.2017.03.063

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • A meshless method is developed for some generalized Klein-Gordon equations. • A convergent iterative scheme is performed to tackle the nonlinearity. • Examples involving five kinds of nonlinear equations and solitons are given.

AC

CE

PT

ED

M

AN US

CR IP T

• Numerical results indicate that the method is about second-order convergence in both time and space.

1

ACCEPTED MANUSCRIPT

Meshless numerical analysis of a class of nonlinear generalized Klein-Gordon equations with a well-posed moving least squares approximation Xiaolin Li∗

CR IP T

College of Mathematics Science, Chongqing Normal University, Chongqing 400047, P R China

Abstract

This paper presents a meshless method for the numerical solution of a class of nonlinear generalized Klein-Gordon equations. In this method, a time discrete technique is first adopted to discretize the time derivatives, and then a well-posed moving least squares (WP-MLS) approximation using shifted and scaled

AN US

orthogonal basis functions is developed to approximate the spatial derivatives. To deal with the nonlinearity, an iterative scheme is presented and the corresponding convergence is discussed theoretically. Numerical examples involving Klein-Gordon, Dodd-Bullough-Mikhailov, sine-Gordon, double sine-Gordon and sinhGordon equations, and line and ring solitons are provided to illustrate the performance and efficiency of the method.

M

Keywords: Meshless; Well-posed moving least squares approximation; Nonlinear generalized

ED

Klein-Gordon equation; Sine-Gordon; Convergence; Soliton

1. Introduction

PT

The Klein-Gordon equation arises in a wide variety of modeling situations in optics, quantum mechanics, solid mechanics, kink and fluid dynamics, nonlinear differential geometry, and plays an important role in the modeling of a number of physical problems such as the fluxon propagation between two superconductors,

CE

the behavior of elementary particles, and the dislocation propagation in crystals [1–4]. In this paper, we consider the following n-dimensional nonlinear generalized Klein-Gordon equation

AC

∂ 2 u (x, t) ∂u (x, t) = α∆u (x, t) − g (u (x, t)) + ω (x, t) , +β ∂t2 ∂t

x ∈ Ω,

t > 0,

(1)

with initial conditions u (x, 0) = φ1 (x) , x ∈ Ω, ∂u (x, t) = φ2 (x) , x ∈ Ω, ∂t t=0

(2) (3)

∗ Corresponding

author. Email address: [email protected] (Xiaolin Li)

Preprint submitted to Applied Mathematical Modelling

April 4, 2017

ACCEPTED MANUSCRIPT

and appropriate boundary conditions x ∈ ΓD ,

u (x, t) = u ¯ (x, t) , q (x, t) :=

∂u (x, t) = q¯ (x, t) , ∂n

t ≥ 0,

x ∈ ΓN ,

(4) t ≥ 0,

(5) T

where Ω is a domain of Rn (n = 1, 2, · · · ) with boundary Γ = ΓD + ΓN , x = (x1 , x2 , · · · , xn ) is the space

CR IP T

variable, t is the time variable, u (x, t) is the unknown function, g (u) is a nonlinear function of u, ω (x, t) is a given function, φ1 (x) and φ2 (x) are prescribed initial functions, u ¯ (x, t) and q¯ (x, t) are prescribed boundary T

functions, n = (n1 , n2 , · · · , nn ) is the unit outward normal on Γ, and α and β are known constants. In the physical applications, the nonlinear force g (u) takes many forms [5–7] as shown in Table 1. Table 1: Forms of the nonlinear force g (u).

Name of Eq. (1)

sin u

Sine-Gordon equation

AN US

g (u)

sinh u

Sinh-Gordon equation

sin u + sin 2u

Double sine-Gordon equation

sinh u + sinh 2u

Double sinh-Gordon equation

exp u

Liouville equation

Dodd-Bullough-Mikhailov equation Tzitzeica-Dodd-Bullough equation

ED

exp (−u) + exp (−2u)

M

exp u + exp (−2u)

In addition, the nonlinear force g (u) can be chosen as the following polynomial (6)

PT

g (u) = γ1 u + γ2 um ,

where γ1 and γ2 are given constants. In this situation, Eq. (1) is called the Klein-Gordon equation with

CE

quadratic nonlinearity if m = 2, and with cubic nonlinearity if m = 3 [6–8]. The solutions of these nonlinear equations are important in mathematical physics. Some analytical methods, such as B¨acklund transformation [9], Lamb’s method [10], Painlev`e transcendents [11], Jacobi

AC

elliptic-function method [12], the sine-cosine method [6], the tanh method [13, 14], the auxiliary equation method [15] and Hirota’s method [16] have been proposed to obtain analytical solutions of these equations. It should be stressed that analytical solutions are difficult to obtain in the general cases. For the sine-Gordon equation, the study of analytical solutions has been focused on the special undamped case (β = 0). In the past three decades, considerable attention has been paid to the numerical solutions of these

nonlinear generalized Klein-Gordon equations with the finite difference method (FDM) [17–21], the finite element method (FEM) [22], the boundary element method (BEM) [23–25], the spline difference method [26], the differential quadrature method (DQM) [27], and so on. In the FDM, the FEM and the BEM, the 3

ACCEPTED MANUSCRIPT

accuracy of the numerical solutions relies on the quality of meshes, which can be arduous and computationally expensive. Meshless (or meshfree) methods [28, 29] alleviate the difficulty of meshing the computational domain by using a scattered set of distributed nodes. Table 2 gives some previously reported applications of meshless methods for generalized Klein-Gordon equations. In these works, each individual equation has been studied by a special meshless method. Besides, compared with the FDM and the FEM, not many meshless methods

meshless methods for the numerical solution of these equations.

CR IP T

have been applied to nonlinear generalized Klein-Gordon equations. Therefore, it is important to develop

Table 2: Some previously reported applications of meshless methods for generalized Klein-Gordon equations.

Meshless method

1D Klein-Gordon equation

Radial basis function (RBF) method [7, 30], Chebyshev tau meshless method [8]

1D sine-Gordon equation

AN US

Equation

Chebyshev tau meshless method [8], multiquadric quasi-interpolation method [31], meshless Ritz method [32]

2D sine-Gordon equation

RBF method [30, 33], meshfree reproducing kernel particle Ritz (MRKPR) method [34], radial point interpolation method (RPIM)

M

[35], meshless local Petrov-Galerkin (MLPG) method [36], moving least square radial reproducing polynomial (MLSRRP) meshless method [37], element-free Galerkin (EFG) method [38] RBF method [39], RBF-pseudospectral method [39], MLS method [39],

ED

2D sinh-Gordon equation

element-free Galerkin method [38] RBF method [40]

PT

Time fractional sine-Gordon

CE

and Klein-Gordon equations

The moving least squares (MLS) approximation [41] is a frequently used meshless methods. As in the traditional least squares method, the moment matrix in the MLS maybe ill-conditioned or singular. Then,

AC

the inversion of the moment matrix leads to the decrease in stability and the increase in computational cost and error. To overcome this drawback, Refs. [42–45] introduce weighted orthogonal basis functions into the MLS approximation and produce the improved moving least squares (IMLS) approximation. The moment matrix in the IMLS is diagonal, and therefore the final system of algebra equations in the IMLS is solved without the inversion of the moment matrix. Hence, the IMLS improves the computational efficiency. In this paper, a meshless method is developed for solving nonlinear generalized Klein-Gordon equations. In this method, the time derivatives are discretized by a time discrete technique. Then, a new implementation of the IMLS approximation is developed to establish the full discretization system. In our implementation, 4

ACCEPTED MANUSCRIPT

here called the well-posed moving least squares (WP-MLS) approximation, shifted and scaled orthogonal basis functions are used to obtain more accurate and stable results. Besides, an iterative scheme is performed to tackle the nonlinearity. The convergence of the iterative process is proved theoretically. The present method is truly meshless, i.e. absolutely no elements or cells are required in the problem domain. As a result, the present method is easy to apply for multidimensional problems and is expected to have higher computational speed and efficiency. Numerical examples involving Klein-Gordon, Dodd-Bullough-Mikhailov,

CR IP T

sine-Gordon, double sine-Gordon and sinh-Gordon equations, and line and ring solitons are given to show the capability of the method.

The rest of this paper is organized as follows. Section 2 develops a time discrete technique to discretize the time derivatives. Section 3 describes and analyzes the present WP-MLS approximation. Then, detailed numerical implementation of our meshless method for nonlinear generalized Klein-Gordon equations is given

AN US

in Section 4. Section 5 presents some numerical examples. Finally, Section 6 contains conclusions. 2. Discretization of time

According to the Taylor series, the time derivatives can be obtained by the difference quotient as [46]
u(k+1) (x) − u(k−1) (x) ∂u (x, kτ ) = + O τ2 , ∂t 2τ

M


∂2u u(k+1) (x) − 2u(k) (x) + u(k−1) (x) (x, kτ ) = + O τ2 , 2 2 ∂t τ

(7) (8)

ED

where τ > 0 is the time-step size and u(k) (x) = u (x, kτ ).

(9)

PT

In addition, according to the Crank-Nicolson scheme, we have i
1 h (k+1) ∆u (x, kτ ) = ∆u (x) + ∆u(k) (x) + ∆u(k−1) (x) + O τ 2 . 3

AC

and

CE

Furthermore, we can expand g (u (x, t)) in t to obtain     
∂u dg  (k) g u(k+1) (x) = g u(k) (x) + τ (x, kτ ) u (x) + O τ 2 , ∂t du

Consequently,

    
∂u dg  (k) g u(k−1) (x) = g u(k) (x) − τ (x, kτ ) u (x) + O τ 2 . ∂t du

g (u (x, kτ )) =

    i
1 h  (k+1) g u (x) + g u(k) (x) + g u(k−1) (x) + O τ 2 . 3

(10)

By considering Eq. (1) in point (x, kτ ) and applying Eqs. (7)-(10), we have

u(k+1) − u(k−1) u(k+1) − 2u(k) + u(k−1) + β τ2 2τ i 1h      i
α h (k+1) (k) (k−1) = ∆u + ∆u + ∆u − g u(k+1) + g u(k) + g u(k−1) + ω (k) + O τ 2 , 3 3 5

(11)

ACCEPTED MANUSCRIPT

where ω (k) (x) = ω (x, kτ ) is known at the kth time level. Then, omitting the small term in Eq. (11) yields

where

  −2ατ 2 ∆u(k+1) + 3 (τ β + 2) u(k+1) + 2τ 2 g u(k+1) = f (k) ,

k = 0, 1, 2, . . .

(12)


Clearly, the truncation error of Eq. (12) is O τ 2 .

CR IP T

h    i f (k) = 2ατ 2 ∆u(k) + 12u(k) + 2ατ 2 ∆u(k−1) + 3 (τ β − 2) u(k−1) − 2τ 2 g u(k) + g u(k−1) + 6τ 2 ω (k) . Eq. (12) requires u(0) and u(−1) when k = 0. From Eq. (2) it follows that u(0) (x) = φ1 (x) ,

x ∈ Ω.

Additionally, from Eqs. (3) and (7) we obtain

AN US


u(−1) (x) = u(1) (x) − 2τ φ2 (x) + O τ 3 ,

x ∈ Ω.

(14)

Letting k = 0 in Eq. (11), invoking Eqs. (13) and (14), and omitting the small term, we have     −2ατ 2 ∆u(1) + 6u(1) + τ 2 g u(1) + τ 2 g u(1) − 2τ φ2 = f (0) , where

(13)

(15)

M

f (0) = ατ 2 ∆φ1 + 6φ1 − 2ατ 3 ∆φ2 − 3τ (τ β − 2) φ2 − τ 2 g (φ1 ) + 3τ 2 ω (0) . Finally, according to Eqs. (12) and (15), the time-dependent Klein-Gordon problem (1)-(5) is transformed

ED

into the following time-independent equations     −2ατ 2 ∆u(1) + 6u(1) + τ 2 g u(1) + τ 2 g u(1) − 2τ φ2 = f (0) , 



PT

−2ατ 2 ∆u(k+1) + 3 (τ β + 2) u(k+1) + 2τ 2 g u(k+1) = f (k) ,

CE

with boundary conditions

∆

u(k+1) (x) = u ¯(k+1) (x) = u ¯ (x, (k + 1) τ ) ,

x ∈ ΓN ,

k = 1, 2, . . .

k = 0, 1, 2, . . . k = 0, 1, 2, . . .

(16) (17)

(18) (19)

AC

∂u(k+1) (x) ∆ = q¯(k+1) (x) = q¯ (x, (k + 1) τ ) , ∂n

x ∈ ΓD ,

x ∈ Ω,

x ∈ Ω,

3.

A well-posed moving least squares (WP-MLS) approximation using shifted and scaled

orthogonal basis functions N

To approximate the unknown function u(k+1) (x) in Eqs. (16)-(19), let {xi }i=1 be the set of all nodes in

¯ = Ω ∪ Γ. Then, we use Ω

h = max

min

1≤i≤N 1≤j≤N,j6=i

to denote the nodal spacing. 6

|xi − xj |

ACCEPTED MANUSCRIPT

3.1. Formulation of the improved moving least squares (IMLS) approximation As in the MLS [41], the local approximation of u(k+1) (x) at x in the IMLS can be defined as u

(x) ≈

(k+1) uh

m X

(x, x ¯) =

(k+1)

qj (¯ x) aj

(x),

j=1

(k+1)

where qj (¯ x) are constructed basis functions, aj

x ∈ Γ,

(20)

(x) are unknown coefficients, m is the number of un-

knowns, and x ¯ can either be x or xi .

CR IP T

(k+1)

The basis functions qj (¯ x) are the following weighted orthogonal functions [42–45], q1 (¯ x) = 1,

qj (¯ x) = pj (¯ x) −

j−1 X (pj , qk )

x

k=1

(qk , qk )x

qk (¯ x),

j = 2, 3, · · · m,

(21)

where the inner product (·, ·)x on nodes {xi }i∈∧(x) is defined for functions v1 and v2 as X

wi (x) v1 (xi ) v2 (xi ),

AN US

(v1 , v2 )x =

i∈∧(x)

where xi , i ∈ ∧ (x), are nodes whose influence domains cover x, and ∧ (x) = {I1 , I2 , · · · , I` } represents the global sequence numbers of these nodes.

A number of weight functions, such as Gaussian, exponential, cubic and quartic splines have been extensively investigated for meshless methods [29]. In this study, the weight function is chosen as the cubic

M

spline [29], i.e.,

ED

   2/3 − 4d2 + 4d3 , d ≤ 1/2,     wi (x) = 4/3 − 4d + 4d2 − 4d3 3, 1/2 < d ≤ 1,      0, d > 1,

PT

 ¯ d¯ is the radius of the influence domain of xi . In all examples presented in Section 5, where d = |x − xi | d,

d¯ is taken to be 2.5h, with h as the nodal spacing.

The functions pj (¯ x) in Eq. (21) can be taken as the polynomial bases used in the MLS [45]. Let T

AC

CE

p (¯ x) = [p1 (¯ x) , p2 (¯ x) , · · · , pm (¯ x)], then a linear basis is

and a quadratic basis is

pT (¯ x) = [1, x ¯1 , x ¯2 , · · · , x ¯n ] ,   pT (¯ x) = 1, x ¯1 , x ¯2 , · · · , x ¯n , x ¯21 , x ¯1 x ¯2 , · · · , x ¯2n .

Besides, we can define pj (¯ x) as [42–44] p1 (¯ x) = 1,

pj (¯ x) = rj−1 (¯ x) ,

where r (¯ x) = x ¯1 for one-dimensional problems, and r (¯ x) = problems. 7

j = 2, 3, · · · , m,

Pn

i=1

x ¯i or r (¯ x) =

pPn

i=1

x ¯2i for n-dimensional

ACCEPTED MANUSCRIPT

(k+1)

The unknowns aj

J=

X

wi (x)

i∈∧(x)

as

(x) in Eq. (20) are solved minimizing the L2 norm h

(k+1) uh

(x, xi ) −

(k+1) ui

i2

X

=

i∈∧(x)

 2 m X (k+1) (k+1)  wi (x)  qj (xi ) aj (x) − ui j=1

h iT (k+1) (k+1) a(k+1) (x) = a1 (x) , a2 (x) , · · · , a(k+1) (x) = A−1 (x) B (x) u(k+1) , m

(22)

B (x) are defined as

i, j = 1, 2, · · · , m,

Aij (x) = (qi , qj )x ,

j = 1, 2, · · · , m;

Bji (x) = wIi (x) qj (xIi ) ,

CR IP T

h iT (k+1) (k+1) (k+1) (k+1) where u(k+1) = uI1 , uI2 , · · · , uI` with ui are the nodal values, and the matrices A (x) and (23)

i = 1, 2, · · · , `.

Aij (x) = Consequently, (x) =

P`

  (qj , qj ) , x 

0,

(k+1)

i=1

Bji (x) uIi (qj , qj )x

=

i = j,

i, j = 1, 2, · · · , m.

i 6= j,

P

(k+1)

i∈∧(x)

wi (x) qj (xi ) ui

j = 1, 2, · · · , m.

,

(qj , qj )x

M

(k+1) aj

AN US

The basis functions given in Eq. (21) are weighted orthogonal, i.e., (qi , qj )x = 0 for i 6= j. Thus,

(24)

By inserting Eq. (24) into Eq. (20), we have the following local approximation (x, x ¯) =

m X

qj (¯ x)

j=1

P

(k+1)

wi (x) qj (xi ) ui

ED

(k+1) uh

i∈∧(x)

(qj , qj )x

=

X

wi (x)

m X qj (xi ) qj (¯ x) j=1

i∈∧(x)

(qj , qj )x

(k+1)

ui

.

PT

Finally, the global approximation of u(k+1) (x) is (k+1)

(k+1)

(x) = uh

(x, x ¯) |x¯=x =

CE

u(k+1) (x) ≈ uh

X

(k+1)

Φi (x) ui

i∈∧(x)

=

N X

(k+1)

Φi (x) ui

,

(25)

i=1

AC

where the IMLS shape function is  m X qj (xi ) qj (x)   ,  wi (x) (qj , qj )x j=1 Φi (x) =    0,

i ∈ ∧ (x) ,

i = 1, 2, · · · , N.

i∈ / ∧ (x) ,

In Refs. [47–51], Li et al. analyzed theoretically the properties and errors of the IMLS and the MLS, and then obtained error estimates and convergence results of the associated element-free Galerkin method and Galerkin boundary node method in Sobolev spaces.

8

ACCEPTED MANUSCRIPT

3.2. Analysis of the moment matrix In Eq. (21), the functions pj (¯ x) are used to generate the orthogonal basis functions qj (¯ x). Thus, the inversion of the moment matrix A (x) is not required in the IMLS. If pj (¯ x) are directly used as basis functions in the MLS, then it is well known that the associated moment ˜ (x) in the MLS can be written as matrix A i, j = 1, 2, · · · , m.

CR IP T

˜ ij (x) = (pi , pj ) 6= 0, A x

˜ (x) is not diagonal, the inversion of A ˜ (x) in the MLS leads to the decrease in stability and the Since A increase in computational cost. In Refs. [47, 51, 52], it is proved theoretically that

AN US

  P ˆ ˜ (x) = Cd (x, n, m) h2 m j=1 j , det A

where Cd (x, n, m) is a bounded and computable number, and ˆj expresses the largest degree of pj (¯ x). Besides, there exists a lower unitriangular matrix Q (x) such that [47] ˜ (x) QT (x) . A (x) = Q (x) A

M

Clearly, det (Q (x)) = 1, then we have

(26)

  P ˆ ˜ (x) = Cd (x, n, m) h2 m j=1 j . det (A (x)) = det A

ED

˜ (x) in the MLS approximation, there is a bounded and computable Moreover, for the moment matrix A

PT

˜ (x) is [52] number C˜c (x, n, m) independent of h such that the 2-norm condition number of A   ˆ ˜ (x) = C˜c (x, n, m) h−2m cond A .

Because the lower unitriangular matrix Q (x) satisfies cond (Q (x)) = 1, from Eq. (26) we can verify that

AC

CE

the 2-norm condition number of the moment matrix A (x) in the IMLS approximation can be written as ˆ cond (A (x)) = Cc (x, n, m) h−2m ,

where Cc (x, n, m) is a bounded and computable number. Thus, for h small enough, the IMLS may form an ill-conditioned moment matrix A (x). And then, the (k+1)

accuracy of computing the unknowns aj

(x) decreases.

3.3. Shifted and scaled orthogonal basis functions In Ref. [53], a shifted and scaled polynomial basis function is used in the MLS to improve the condition number of the associated moment matrix. This technique is verified to be very efficient [52]. 9

ACCEPTED MANUSCRIPT

In this study, the following shifted and scaled functions are first defined using pj (¯ x) as   x ¯ − xe psj (¯ x ) = pj , h where xe is a fixed point located on the influence domain of the evaluation point x. In particular, we can psj (¯ x) as q1s

(¯ x) = 1,

qjs

psj

(¯ x) =

(¯ x) −

j−1 X psj , qks

k=1




x s qk

(qks , qks )x

(¯ x),

CR IP T

take xe = x. Then, as in Eq. (21), shifted and scaled orthogonal basis functions qjs (¯ x) are constructed using

j = 2, 3, · · · m.

Finally, substituting qjs (x) for qj (x) in the IMLS approximation, the WP-MLS approximation can be obtained directly. The corresponding moment matrix is denoted as As (x) and can be written as



0,

Besides, the WP-MLS shape function is

i = j, i 6= j,

i, j = 1, 2, · · · , m.

AN US

Asij (x) =


 qjs , qjs , x

 m X qjs (xi ) qjs (x)  
,  wi (x) qjs , qjs x j=1 Φi (x) =    0,

i ∈ ∧ (x) ,

i = 1, 2, · · · , N.

(27)

M

i∈ / ∧ (x) ,

For the moment matrix As (x) in the WP-MLS approximation, it can be proved as in the MLS approxdet (As (x)) = Cds (x, n, m) ,

(28)

cond (As (x)) = Ccs (x, n, m) ,

(29)

ED

imation [52] that

PT

where Cds (x, n, m) and Ccs (x, n, m) are bounded and computable numbers. From Eqs. (28) and (29) we can observe that both the determinant and the condition number of the

CE

moment matrix As (x) in the WP-MLS approximation are independent of h. Consequently, the present WP-MLS always form a well-conditioned moment matrix As (x) even if h decreases. And then, the present

AC

WP-MLS yields more accurate and stable results.

4. The WP-MLS approximation for generalized Klein-Gordon equations 4.1. Discrete formulas Using the WP-MLS shape function Φi (x) given in Eq. (27), the unknown function u(k+1) (x) can be approximated as u(k+1) (x) =

N X i=1

10

(k+1)

ui

Φi (x).

(30)

ACCEPTED MANUSCRIPT

Similarly, we can approximate u(k) (x) and u(k−1) (x) respectively as u(k) (x) =

N X

ui Φi (x),

N X

ui

(k)

(31)

i=1

u(k−1) (x) =

(k−1)

Φi (x).

(32)

Differentiating Eq. (30) yields N X

∆u(k+1) (x) =

(k+1)

ui

CR IP T

i=1

x ∈ Ω,

∆Φi (x),

i=1

q (k+1) (x) =

N X

(k+1) ∂Φi

(x) , ∂n

ui

i=1

(k+1)

(34)

and to simplify the representation, we further assume that

AN US

To obtain the unknown coefficient ui

x ∈ Γ.

(33)

N

N +N

I I D the first NI nodes {xj }j=1 belong to Ω, the next ND nodes {xj }j=N belong to ΓD and the remaining I +1

N

N − NI − ND nodes {xj }j=NI +ND +1 belong to ΓN .

The boundary value problem given by Eqs. (17)-(19) can be solved iteratively. The details follow. Firstly,

M

using Eqs. (30)-(33) and collocating the governing equation (17) for all internal nodes xj ∈ Ω yield ! N N X X  (k+1)  (k+1) 2 2 ui 3 (τ β + 2) Φi (xj ) − 2ατ ∆Φi (xj ) + 2τ g ui Φi (xj ) i=1

(k)

ui

i=1



 2ατ 2 ∆Φi (xj ) + 12Φi (xj ) + (k) ui Φi

!

(xj )

i=1

− 2τ 2 g

N X

N X

i=1

(k−1)

ui

i=1

(k−1) ui Φi

  2ατ 2 ∆Φi (xj ) + 3 (τ β − 2) Φi (xj )

(xj )

i=1

PT

− 2τ 2 g

N X

ED

=

N X

!

+ 6τ 2 ω (k) (xj ) ,

j = 1, 2, · · · , NI .

(35)

Then, using Eq. (30) and collocating the Dirichlet boundary condition (18) for all boundary nodes xj ∈ ΓD N X

ui

N X

ui

CE

yield

(k+1)

Φi (xj ) = u ¯(k+1) (xj ) ,

i=1

j = NI + 1, · · · , NI + ND .

(36)

Finally, using Eq. (34) and collocating the Neumann boundary condition (19) for all boundary nodes

AC

xj ∈ ΓN yield

i=1

(k+1) ∂Φi

(xj ) = q¯(k+1) (xj ) , ∂n

j = NI + ND + 1, · · · , N.

(37)

Eqs. (35)-(37) contain a set of coupled N algebraic equations and can be represented in matrix form as       Au(k+1) + f u(k+1) = Bu(k) + Cu(k−1) − f u(k) − f u(k−1) + b(k) ,

where

h iT (k) (k) (k) u(k) = u1 , u2 , · · · , uN , 11

k = 0, 1, 2, · · · ,

k = 1, 2, · · · ,

(38)

ACCEPTED MANUSCRIPT

AN US

CR IP T

   3 (τ β + 2) Φi (xj ) − 2ατ 2 ∆Φi (xj ) , j = 1, 2, · · · , NI ,    [A]ij = Φi (xj ) , i = 1, 2, · · · , N, j = NI + 1, · · · , NI + ND ,      ∂Φ (x )/∂n, j = NI + ND + 1, · · · , N, i j   2ατ 2 ∆Φi (xj ) + 12Φi (xj ) , j = 1, 2, · · · , NI , [B]ij = i = 1, 2, · · · , N,  0, j = NI + 1, · · · , N,   2ατ 2 ∆Φi (xj ) + 3 (τ β − 2) Φi (xj ) , j = 1, 2, · · · , NI , [C]ij = i = 1, 2, · · · , N,  0, j = NI + 1, · · · , N,      2τ 2 g Φ (xj ) u(k) , j = 1, 2, · · · , NI , i h  = k = 0, 1, 2, · · · , f u(k)  j  0, j = NI + 1, · · · , N,    6τ 2 ω (k) (xj ) , j = 1, 2, · · · , NI ,    h i (k) k = 1, 2, · · · . b = u ¯(k+1) (xj ) , j = NI + 1, · · · , NI + ND ,  j     q¯(k+1) (x ) , j = NI + ND + 1, · · · , N, j

(39)

(40)

(41)

(42)

Similarly, at the initial stage k = 0, the boundary value problem given by Eqs. (16), (18) and (19) can be discretized as

M

   6Φi (xj ) − 2ατ 2 ∆Φi (xj ) , j = 1, 2, · · · , NI ,      ¯ i = 1, 2, · · · , N, A = Φi (xj ) , j = NI + 1, · · · , NI + ND , ij      ∂Φ (x )/∂n, j = NI + ND + 1, · · · , N, i j        τ 2 g Φ (xj ) u(1) + τ 2 g Φ (xj ) u(1) − 2τ φ2 (xj ) , j = 1, 2, · · · , NI , i h  ¯ f u(1) =  j  0, j = NI + 1, · · · , N,    f (0) (xj ) , j = 1, 2, · · · , NI ,    h i ¯ (0) = u b ¯(1) (xj ) , j = NI + 1, · · · , NI + ND ,  j     q¯(1) (x ) , j = N + N + 1, · · · , N, j I D

(43)

(44)

(45)

(46)

AC

CE

PT

ED

where

  ¯ (0) , ¯ (1) + ¯ Au f u(1) = b

with

f (0) (xj ) = ατ 2 ∆φ1 (xj ) + 6φ1 (xj ) − 2ατ 3 ∆φ2 (xj ) − 3τ (τ β − 2) φ2 (xj ) − τ 2 g (φ1 (xj )) + 3τ 2 ω (0) (xj ) .

12

ACCEPTED MANUSCRIPT

4.2. Algorithms
Because g (u) is a nonlinear function of u, from Eqs. (42) and (45) we conclude that f u(k+1) and
¯ f u(1) are nonlinear of u(k+1) and u(1) , respectively. Thus, the algebraic systems given in Eqs. (38) and

(43) are nonlinear. In this research, an iterative algorithm is used in each time level to solve these systems. The following iterative algorithm is expound on the nonlinear algebraic system (38). Obviously, this

CR IP T

algorithm can be used directly for system (43). Algorithm 4.1. Select a tolerance ε > 0.

1. Let l = 0 and u(k+1),l = u(k) .

2. Let f u(k+1) = f u(k+1),l , then Eq. (38) is transformed into the following linear system

      Au(k+1),l+1 = −f u(k+1),l + Bu(k) + Cu(k−1) − f u(k) − f u(k−1) + b(k) ,

AN US

from which we can solve u(k+1),l+1 with ease.

(47)



3. Stop the algorithm and return u(k+1) = u(k+1),l+1 if u(k+1),l+1 − u(k+1),l u(k+1),l+1 ≤ ε. Otherwise, update l to l + 1, and go to 2.

Theorem 4.1. Assume that the nonlinear term g (u) satisfies the Lipschitz condition with respect to u,

M

|g (u1 ) − g (u2 )| ≤ L |u1 − u2 | ,

∀u1 , u2 ∈ L∞ (Ω) ,

i = 1, 2, · · · , N.

(49)

PT

ED

where L is a Lipschitz constant. Let the matrix A0 be defined as    6Φi (xj ) , j = 1, 2, · · · , NI ,    [A0 ]ij = Φi (xj ) , j = NI + 1, · · · , NI + ND ,      ∂Φ (x )/∂n, j = N + N + 1, · · · , N, i j I D

(48)

CE

If A0 is invertible, then there exists a positive constant δ such that if τ < δ, Algorithm 4.1 generates a  ∞ convergent sequence u(k+1),l l=0 , i.e., lim u(k+1),l = u(k+1) .

l→∞

AC

Proof. From Eq. (47) it follows that       Au(k+1),l+2 = −f u(k+1),l+1 + Bu(k) + Cu(k−1) − f u(k) − f u(k−1) + b(k) .

Subtracting Eq. (47) from Eq. (50) yields     Au(k+1),l+2 − Au(k+1),l+1 = f u(k+1),l − f u(k+1),l+1 ,

thus

h    i u(k+1),l+2 − u(k+1),l+1 = A−1 f u(k+1),l − f u(k+1),l+1 , 13

(50)

ACCEPTED MANUSCRIPT

and



(k+1),l+2 − u(k+1),l+1

u

∞

   



≤ A−1 ∞ f u(k+1),l+1 − f u(k+1),l

∞

.

(51)

According to Eqs. (42) and (48), we have



≤ 2τ 2 LCΦ u(k+1),l+1 − u(k+1),l

   

f u(k+1),l+1 − f u(k+1),l

where CΦ = max1≤j≤NI kΦ (xj )k∞ .

CR IP T

∞

∞

(52)

From Eqs. (39) and (49) it follows that limτ →0 A = A0 . Besides, the invertibility of the matrix A0



. Thus, A−1 exists as τ is small enough, and then we have a yields the existence of the norm A−1 0 ∞ ∞ positive constant CA such that

−1

A ≤ CA . ∞

Inserting Eqs. (52) and (53) into Eq. (51) arrives at

∞



≤ 2τ 2 LCΦ CA u(k+1),l+1 − u(k+1),l

AN US



(k+1),l+2

− u(k+1),l+1

u

∞

(53)

.

√ By choosing τ such that 2τ 2 LCΦ CA < 1, namely τ < δ with δ = 1 2LCΦ CA > 0, we finally obtain

(k+1),l+2

− u(k+1),l+1

u

∞



< u(k+1),l+1 − u(k+1),l

∞

.

M

From this error reduction property we can see that liml→∞ u(k+1),l+1 − u(k+1),l ∞ = 0. Therefore, we can  ∞ conclude that the sequence u(k+1),l l=0 converges to u(k+1) . This completes the proof.

ED

Finally, the detailed iterative algorithm of the WP-MLS approximation for the nonlinear generalized

Klein-Gordon problem (1)-(5) is concluded as follows. N

PT

Algorithm 4.2. Select the time-step size τ and N nodes {xi }i=1 in Ω ∪ Γ, and set k = 0.

CE

1. Using Eq. (27), calculate the WP-MLS shape functions for the N nodes.
¯ (0) . ¯ and vectors ¯ 2. Using Eqs. (44)-(46), calculate the matrix A f u(1) and b 3. Solving system (43), obtain u(1) .

AC

4. Using Eqs. (39)-(41), calculate the matrices A, B and C. 5. Update k to k + 1 in system (38). Solving the system, obtain u(k+1) . Repeating this step process, until (k + 1) τ reaches the desired time.

6. Using Eq. (30), calculate the numerical solution of the nonlinear Klein-Gordon problem. Since the coefficient matrices A, B and C are invariable with respect to the time level k, these matrices

can be calculated once. In addition, an efficient technique for solving system (38) is to perform the LU factorization for the matrix A at the initial level. 14

ACCEPTED MANUSCRIPT

5. Numerical experiments To demonstrate the validity of the proposed method, some numerical experiments are presented. In all computations, unless indicated otherwise, regular node distributions are used. 5.1. Validations

CR IP T

Five nonlinear Klein-Gordon, Dodd-Bullough-Mikhailov, sine-Gordon, double sine-Gordon and sinhGordon equations with known analytical solutions are solved in this section to show the capability and convergence of the WP-MLS approximation.

For the error estimation and convergence analysis, the

following L∞ and root-mean-square (RMS) errors are used to evaluate the performance of the numerical method L∞ = max |u (yi ) − uh (yi )| ,

AN US

1≤i≤Ns

v u Ns u 1 X 2 RMS = t [u (yi ) − uh (yi )] , Ns i=1

N

s is the set of test points in Ω, and u (yi ) and uh (yi ) are the exact and approximate values of where {yi }i=1

the function u at yi , respectively. 5.1.1. Klein-Gordon equation

Consider the Klein-Gordon equation with cubic nonlinearity,

M


∂2u ∂2u = α − γ1 u + γ2 u3 , 2 2 ∂t ∂x

t > 0,

ED

subject to the initial conditions

x ∈ (0, 1) ,

PT

u (x, 0) = B tan (Kx) , x ∈ (0, 1) , ∂u (x, t) = BcK sec2 (Kx) , x ∈ (0, 1) , ∂t t=0

CE

and with Dirichlet boundary conditions

u (0, t) = B tan (cKt) ,

t > 0,

u (1, t) = B tan [K (1 + ct)] ,

t > 0,

AC

q p  where B = γ1 /γ2 , K = −0.5γ1 (α + c2 ) and c is a constant. The analytical solution of this problem is [7]

u (x, t) = B tan [K (x + ct)] . In this study, we use α = −2.5, γ1 = 1, γ2 = 1.5 and c = 0.5. Fig. 1 gives the space-time graph of the numerical solution and the graph of the associated error. The results are obtained by using τ = 0.01 and h = 0.01. We can find that the WP-MLS approximation yields very accurate numerical results. 15

ACCEPTED MANUSCRIPT

b

CR IP T

a

Fig. 1. Graphs of (a) the numerical solution and (b) the associated error up to t = 4 with τ = 0.01 and h = 0.01 for

AN US

the Klein-Gordon equation.

Fig. 2 gives the L∞ - and RMS-errors at t = 1, 2, 3 and 4 for the WP-MLS approximation with τ = 0.01 and h = 0.01. For comparison, Fig. 2 also gives the results of the RBF method [7] obtained using τ = 0.001 and h = 0.01. The value of τ used in the WP-MLS approximation is tenfold that used in the RBF method. However, the proposed WP-MLS approximation is more accurate than the RBF method. b

CE

PT

ED

M

a

Fig. 2. Results of (a) L∞ -error and (b) RMS-error obtained using the WP-MLS approximation with τ = 0.01 and

AC

h = 0.01 and the RBF method with τ = 0.001 and h = 0.01 for the Klein-Gordon equation.

To investigate the convergence in time, the L∞ - and RMS-errors are displayed with respect to τ in Fig.

3. The errors are also tabulated in Table 3. These results are obtained by the WP-MLS approximation using h = 0.01. Clearly, both errors decrease as τ decreases, and the experimental convergence form in time
approximates to O τ 2 .

For investigating the convergence in space, the errors are displayed with respect to h in Fig. 4 and Table 16

ACCEPTED MANUSCRIPT

b

AN US

CR IP T

a

PT

ED

M

Fig. 3. Graphs of (a) L∞ -error and (b) RMS-error with respect to τ for the Klein-Gordon equation.

τ

CE

Table 3: Errors of the WP-MLS approximation using h = 0.01 for the Klein-Gordon equation. L∞ −error

RMS-error

t=2

t=3

t=4

t=1

t=2

t=3

t=4

1

6.7199×10−4

1.0236×10−3

4.4296×10−3

1.2042×10−1

4.7679×10−4

7.4219×10−4

3.1683×10−3

8.5058×10−2

1/2

1.4344×10−4

3.0533×10−4

1.7543×10−3

1.9184×10−2

1.0519×10−4

2.1739×10−4

1.2493×10−3

1.2407×10−2

1/4

−5

−5

−4

−3

−5

−5

−4

4.5665×10−3

−5

1.3559×10−3

−6

3.8661×10−4

AC

t=1

1/8

1/16

8.8721×10

−5

1.8214×10

−6

2.5562×10

5.8414×10

−5

2.4813×10

−6

1.2493×10

4.4599×10

−4

1.0151×10

−5

1.3065×10

7.4165×10

−3

2.3084×10

−4

6.9557×10

17

6.0462×10

−5

1.4183×10

−6

1.5498×10

3.9792×10

−5

1.6709×10

−7

8.3358×10

3.1295×10

7.0482×10 7.4653×10

ACCEPTED MANUSCRIPT

4. These errors are obtained by the WP-MLS approximation using τ = 0.01. We can find that the errors
decrease as h decreases, and the experimental convergence form in space approximates to O h2 . b

AN US

CR IP T

a

Fig. 4. Graphs of (a) L∞ -error and (b) RMS-error with respect to h for the Klein-Gordon equation.

Table 4: Errors of the WP-MLS approximation using τ = 0.01 for the Klein-Gordon equation. L∞ −error

h t=2

t=3

1/10

3.8626×10−4

1.0959×10−3

2.8303×10−3

1/20

9.1078×10−5

2.7091×10−4

1/40

−5

−5

1/160

−6

5.0816×10

−6

1.1360×10

6.5086×10

−5

1.5730×10

t=1

t=2

t=3

t=4

1.6956×10−2

2.6262×10−4

8.0855×10−4

2.0372×10−3

5.7566×10−3

6.5676×10−4

3.0947×10−3

6.3985×10−5

1.9812×10−4

4.8490×10−4

8.4410×10−4

−4

−4

−5

−5

−4

1.7915×10−4

−5

7.4076×10−5

−6

2.9885×10−5

1.5120×10

−5

3.5331×10

ED

1/80

2.1484×10

−6

3.7038×10

RMS-error

t=4

M

t=1

−6

7.9686×10

5.8936×10

−4

1.7371×10

−5

6.6497×10

1.5363×10

−6

3.6774×10

−7

8.2860×10

4.7657×10

−5

1.1492×10

−6

2.6940×10

1.1156×10 2.5892×10 5.8033×10

PT

Fig. 5 gives the errors at times t = 1 and 3 for the IMLS approximation and the WP-MLS approximation with τ = 0.01. For the WP-MLS approximation, we can find that the errors in all cases decrease

CE

monotonously as h decreases, and there is no instability. However, when h is sufficiently small, i.e., a large number of nodes are used, the errors of the IMLS approximation do not decrease but increase along with the decrease of h. Consequently, the WP-MLS approximation obtains much more stable and convergent

AC

numerical results than the original IMLS approximation. 5.1.2. Dodd-Bullough-Mikhailov equation In this example, we consider the following nonlinear Dodd-Bullough-Mikhailov equation ∂2u ∂2u − + eu + e−2u = 0, ∂t2 ∂x2

18

x ∈ (0, 1) ,

t > 0.

ACCEPTED MANUSCRIPT

b

CR IP T

a

Fig. 5. Comparison of (a) L∞ -error and (b) RMS-error of the IMLS approximation and the WP-MLS

AN US

approximation for the Klein-Gordon equation.

The initial functions φ1 (x) and φ2 (x) and Dirichlet boundary function u ¯ (x, t) are extracted from the following analytical solution [14]

#!) " r 3 1 (x + ct) , u (x, t) = ln 0.5 1 + 3 tan2 2 1 − c2 (

M

where c is a constant satisfying c2 < 1. In this study, we use c = 0.5.

Fig. 6 gives the space-time graph of the numerical solution and the graph of the associated error. The results are obtained by using τ = 0.1 and h = 0.01. It can be found from Fig. 6 that the WP-MLS

ED

approximation yields very accurate numerical results. b

AC

CE

PT

a

Fig. 6. Graphs of (a) the numerical solution and (b) the associated error up to t = 6 with τ = 0.1 and h = 0.01 for the Dodd-Bullough-Mikhailov equation.

The log-log plot of the L∞ - and RMS-errors at t = 2, 4 and 6 for the WP-MLS approximation with 19

ACCEPTED MANUSCRIPT

respect to τ is shown in Fig. 7. These results are obtained by the WP-MLS approximation using h = 0.0025.
Both errors decrease as τ decreases, and the experimental convergence form in time approximates to O τ 2 . b

AN US

CR IP T

a

Fig. 7. Graphs of (a) L∞ -error and (b) RMS-error with respect to τ for the Dodd-Bullough-Mikhailov equation.

The plot of the errors with respect to h is shown in Fig. 8. These errors are obtained by the WP-MLS

M

approximation using τ = 0.1. We can find from Fig. 8 that the errors decrease as h decreases, and the
experimental convergence form in space approximates to O h2 . b

AC

CE

PT

ED

a

Fig. 8. Graphs of (a) L∞ -error and (b) RMS-error with respect to h for the Dodd-Bullough-Mikhailov equation.

5.1.3. Sine-Gordon equation Consider the nonlinear sine-Gordon equation, ∂2u ∂u ∂2u ∂2u + β = + − sin (u) , ∂t2 ∂t ∂x21 ∂x22 20

x1 , x2 ∈ (−7, 7) ,

t > 0,

ACCEPTED MANUSCRIPT

subject to the initial conditions u (x1 , x2 , 0) = 4 tan−1 exp (x1 + x2 ) , x1 , x2 ∈ (−7, 7) , 4 exp (x1 + x2 ) ∂u (x1 , x2 , t) =− , x1 , x2 ∈ (−7, 7) , ∂t 1 + exp (2x1 + 2x2 ) t=0

and with Neumann boundary conditions

x1 = ±7,

∂u (x1 , x2 , t) 4 exp (x1 + x2 + t) = , ∂x2 exp (2t) + exp (2x1 + 2x2 )

x2 = ±7,

The analytical solution for this problem, in which β = 0, is [33]

x2 ∈ (−7, 7) ,

t > 0,

x1 ∈ (−7, 7) ,

t > 0.

CR IP T

∂u (x1 , x2 , t) 4 exp (x1 + x2 + t) = , ∂x1 exp (2t) + exp (2x1 + 2x2 )

AN US

u (x1 , x2 , t) = 4 tan−1 exp (x1 + x2 − t) .

The WP-MLS approximation is used to solve this problem with τ = 0.1 and h = 0.25. Fig. 9 gives the graphs of the numerical solution and the associated error. b

PT

ED

M

a

CE

Fig. 9. Graphs of (a) the numerical solution and (b) the associated error at t = 7 with τ = 0.1 and h = 0.25 for the sine-Gordon equation.

AC

Fig. 10 gives the L∞ - and RMS-errors at t = 1, 3, 5 and 7 for the WP-MLS approximation with τ = 0.1

and h = 0.25. For comparison, Fig. 10 also gives the results of the RBF method [33] obtained using τ = 0.001 and h = 0.25. The value of τ used in the WP-MLS approximation is 100 times that used in the

RBF method. Nevertheless, Fig. 10 indicates that the WP-MLS gives better accuracy than the RBF. Fig. 11 gives the L∞ - and RMS-errors for the WP-MLS approximation against the time step τ . The results of errors are also tabulated in Table 5. These results are obtained using h = 0.25. Again, both errors
decrease as τ decreases, and the experimental convergence form in time approximates to O τ 2 . 21

ACCEPTED MANUSCRIPT

b

CR IP T

a

Fig. 10. Results of (a) L∞ -error and (b) RMS-error obtained by the WP-MLS approximation with τ = 0.1 and

AN US

h = 0.25 and the RBF method with τ = 0.001 and h = 0.25 for the sine-Gordon equation.

b

CE

PT

ED

M

a

AC

Fig. 11. Graphs of (a) L∞ -error and (b) RMS-error with respect to τ for the sine-Gordon equation.

Table 5: Errors of the WP-MLS approximation using h = 0.25 for the sine-Gordon equation. L∞ −error

τ

RMS-error

t=1

t=3

t=5

t=7

t=1

t=3

t=5

t=7

1

5.3032×10−1

4.7404×10−1

9.6356×10−1

1.7368

8.6102×10−2

1.1642×10−1

2.2588×10−1

3.6475×10−1

1/2

1.2627×10−1

1.8035×10−1

3.1151×10−1

5.4000×10−1

2.0915×10−2

3.2476×10−2

5.9863×10−2

9.5137×10−2

1/4

−2

−2

−2

−2

−3

−2

−2

1.8058×10−2

3.9275×10

8.2267×10

9.6909×10

9.5480×10

22

5.6845×10

1.0545×10

1.5005×10

ACCEPTED MANUSCRIPT

Fig. 12 gives the errors against the nodal spacing h. The results of errors are also tabulated in Table 6. These results are obtained using τ = 0.1. We can find that the experimental convergence form in space
approximates to O h2 . b

AN US

CR IP T

a

Fig. 12. Graphs of (a) L∞ -error and (b) RMS-error with respect to h for the sine-Gordon equation.

L∞ −error

h

M

Table 6: Errors of the WP-MLS approximation using τ = 0.1 for the sine-Gordon equation. RMS-error

t=1

t=3

t=5

t=7

t=1

t=3

t=5

t=7

2

1.9273

2.2295

2.6798

3.1986

3.3487×10−1

4.1167×10−1

5.8630×10−1

6.2923×10−1

1

3.7421×10−1

7.3596×10−1

4.2929×10−2

9.7580×10−2

1.4199×10−1

1.8378×10−1

−1

−1

−3

−2

−2

4.6515×10−2

−3

1.0667×10−2

−2

2.8764×10

2.1361×10

−2

5.6524×10

−1

3.2791×10

−2

8.2637×10

1.4041 −1

4.4798×10

−1

1.0564×10

8.4887×10

−3

2.2170×10

2.4547×10

−3

6.4064×10

3.5767×10 8.7801×10

PT

1/4

1.0861×10

ED

1/2

1.0735

CE

5.1.4. Double sine-Gordon equation

In this example, we consider the following nonlinear double sine-Gordon equation x ∈ (−10, 10) ,

t > 0.

AC

 2 ∂2u ∂2u u 2η sin u − sin = ω (x, t) , − + ∂t2 ∂x2 1 + 4 |η| 2

where η =

1 4

sinh2 R, and R is a positive constant. The initial functions φ1 (x) and φ2 (x), Neumann

boundary function q¯ (x, t) and the right-hand function ω (x, t) are extracted from the following analytical

solution [5]

−1

u (x, t) = 4πm ± 4 tan

 √  sinh (x − ct) 1 − c2 , cosh R

where m is an integer and c is a constant. The signs “+” and “-” correspond to the kink and anti-kink wave solutions, respectively. We can interpret the function u as a superposition of two sine-Gordon solitons, 23

ACCEPTED MANUSCRIPT

separated by the distance 2R [5]. In this study, we use m = 0, c = 0.8 and R = 4, and consider the case of anti-kink for numerical test. Fig. 13 gives the space-time graph of the numerical solution and the graph of the associated error. The results are obtained by the WP-MLS approximation using τ = 0.1 and h = 0.05. We can find that the WP-MLS approximation yields very accurate numerical results. Besides, the profile of the wave consists two kinks.

CR IP T

b

AN US

a

Fig. 13. Graphs of (a) the numerical solution and (b) the associated error up to t = 10 with τ = 0.1 and h = 0.05

M

for the double sine-Gordon equation.

ED

Fig. 14 gives the plot of errors obtained using h = 0.05 with respect to τ , while Fig. 15 gives the plot of errors obtained using τ = 0.1 with respect to h. Again, high rates of convergence are achievable. b

AC

CE

PT

a

Fig. 14. Graphs of (a) L∞ -error and (b) RMS-error with respect to τ for the double sine-Gordon equation.

24

ACCEPTED MANUSCRIPT

b

CR IP T

a

AN US

Fig. 15. Graphs of (a) L∞ -error and (b) RMS-error with respect to h for the double sine-Gordon equation.

5.1.5. Sinh-Gordon equation

Consider the nonlinear sinh-Gordon equation,

∂2u ∂2u ∂2u = + − sinh (u) + ω (x1 , x2 , t) , 2 ∂t ∂x21 ∂x22

T

(x1 , x2 ) ∈ Ω,

t > 0,

(54)

The initial functions φ1 (x1 , x2 ) and φ2 (x1 , x2 ), Dirichlet boundary function u ¯ (x1 , x2 , t) and the right-hand

ED

M

function ω (x1 , x2 , t) are extracted from the following analytical solution [39]   q 1 u (x1 , x2 , t) = 4 tan−1 exp t + 36 + 30x21 + 12x1 x2 + 30x22 . 6

2

First, the sinh-Gordon equation (54) is solved in the rectangular domain Ω = (−6, 6) . For different values of τ , Fig. 16 presents the L∞ -error at time t = 1 obtained by the WP-MLS

PT

approximation with h = 0.2. In Ref. [39], this problem is solved by the MLS, RBF and RBF-pseudospectral (PS) meshless methods. The errors of the three methods with h = 0.2 are also given in Fig. 16. It can be

CE

found that the accuracy and the convergence of the present WP-MLS approximation are much better than the methods in Ref. [39].

For different values of h, Fig. 17 gives the L∞ -error and the condition number of the system matrix

AC

at time t = 1 obtained by the WP-MLS approximation with τ = 0.01. The results of the MLS, RBF and RBF-PS methods with τ = 0.01 are also given for comparison. From Fig. 17 (a), it can be seen again that the error of the present WP-MLS is the best. From Fig. 17 (b) we can observe that the coefficient matrix in the WP-MLS is well-conditioned, and the WP-MLS has a significantly smaller condition number than the MLS and RBF methods. Therefore, the WP-MLS has higher computational precision and better numerical stability. Using the same step lengths h = 0.2 and τ = 0.01 as in Ref. [39], Fig. 18 presents the graphs of the WP-MLS solution and the associated error at t = 1, 5 and 10. Compared with the results in Ref. [39], the 25

AN US

CR IP T

ACCEPTED MANUSCRIPT

M

Fig. 16. Graph of L∞ -error at t = 1 with respect to τ for the sinh-Gordon equation.

b

AC

CE

PT

ED

a

Fig. 17. Graphs of (a) L∞ -error and (b) condition number at t = 1 with respect to h for the sinh-Gordon equation.

26

ACCEPTED MANUSCRIPT

error obtained by the WP-MLS approximation is about 100 times smaller than the errors obtained by the MLS and RBF methods. Thus, the accuracy of the WP-MLS approximation is much better than the two

AC

CE

PT

ED

M

AN US

CR IP T

methods.

Fig. 18. The WP-MLS solutions (left panel) and associated errors (right panel) for the sinh-Gordon equation.

Next, to show the advantages of the meshless nature of the present method, a more complicated geometry is also tested. The boundary Γ of the domain Ω is defined by the parametric equation [54] n o T Γ = (x1 , x2 ) ∈ R2 : x1 (θ) = r (θ) cos θ, x2 (θ) = r (θ) sin θ, θ ∈ [0, 2π] , 27

ACCEPTED MANUSCRIPT

where

p r (θ) = 3 cos2 (5θ) + cos2 (2θ) + cos2 θ.

Tests with both uniformly distributed and Halton-Smith points [53, 54] are performed on the domain. Fig. 19 shows the domain and the distribution of collocation points. Halton-Smith points sets are nested, and one can add new points without generating the whole set from the beginning [54]. Generally, it is

¯ collocation points in Ω.

b

M

AN US

a

CR IP T

expensive to evaluate the value of h exactly for Halton-Smith points. As in Ref. [53], one may use the

¯ N , where mes Ω ¯ denotes the area of Ω ¯ = Ω ∪ Γ and N is the number of approximation h ≈ mes Ω

Fig. 19. Distribution of collocation points in the domain and on its boundary. (a) Uniformly distributed points and

ED

(b) Halton-Smith points.

b

AC

CE

PT

a

Fig. 20. Error of the WP-MLS approximation using (a) uniformly distributed points and (b) Halton-Smith points for the sinh-Gordon equation.

Fig. 20 gives error |u − uh | at time t = 1 obtained using the present method with τ = 0.01. In this 28

ACCEPTED MANUSCRIPT

analysis, the two kinds of collocation points shown in Fig. 19 are used, and the number of collocation points is N = 296. From Fig. 20 we can see that the present method obtains very accurate results in both cases. In addition, Table 7 presents the L∞ - and RMS-errors at time t = 1 obtained by the present method using τ = 0.01 with different number of collocation points. The RMS-error obtained by uniformly distributed points is slightly better than that obtained by Halton-Smith points, but the L∞ -error obtained by uniformly

CR IP T

distributed points is close to that obtained by Halton-Smith points. Table 7: Errors of the present method using uniformly distributed and Halton-Smith points for the sinh-Gordon equation.

Uniformly distributed points

N

L∞ −error

RMS-error

Halton-Smith points L∞ −error

RMS-error

7.3942×10−2

1.4900×10−2

7.2825×10−2

1.8456×10−2

104

1.7768×10−2

2.5988×10−3

2.2981×10−2

3.9494×10−3

296

7.1477×10−3

7.0671×10−4

6.5459×10−3

1.1502×10−3

AN US

42

5.2. Applications

Case studies of more complex nonlinear equations are made to further illustrate the performance and

M

efficiency of the developed WP-MLS approximation.

In this section, we consider the numerical simulation of the solution to the following nonlinear sine-Gordon equation

ED

∂u (x, t) ∂ 2 u (x, t) = ∆u (x, t) − ψ (x) sin u (x, t) , +β 2 ∂t ∂t

x ∈ Ω,

t > 0,

(55)

PT

subject to the initial conditions

u (x, 0) = φ1 (x) , x ∈ Ω, ∂u (x, t) = φ2 (x) , x ∈ Ω, ∂t t=0

∂u (x, t) = 0, ∂n

x ∈ Γ = ∂Ω,

AC

CE

and with Neumann boundary conditions

t ≥ 0.

We can interpret the function ψ (x) in Eq. (55) as a Josephson current density. Besides, β ≥ 0 denotes the dissipative term. Eq. (55) is the damped equation if β > 0, and is the undamped equation if β = 0. In the following numerical experiments, cases including line and ring solitons for the solution of Eq. (55)

are considered.

29

ACCEPTED MANUSCRIPT

5.2.1. Superposition of two orthogonal line solitons This type of soliton can be simulated by choosing ψ (x1 , x2 ) = 1, β = 0 and [17, 19–24, 27, 30, 33–38]
φ1 (x1 , x2 ) = 4 tan−1 exp (x1 ) + tan−1 exp (x2 ) ,

φ2 (x1 , x2 ) = 0,

x1 , x2 ∈ (−6, 6) .

In this example, we use h = 0.5 and τ = 0.1 in computation. The numerical results of the displacement u at times t = 1, 3, 5 and 7 are plotted in Fig. 21, from which

CR IP T

we can see that the two solitons are separated paralleling to the line x1 + x2 = 0. The results agree well with those obtained by the FDM [17, 19–21], the FEM [22], the BEM [23, 24], the DQM [27], the RBF [30, 33],

CE

PT

ED

M

AN US

the MRKPR [34], the RPIM [35], the MLPG [36], the MLSRRP [37] and the EFG [38].

Fig. 21. Numerical results of u for the superposition of two orthogonal line solitons when β = 0.

AC

To further show the separation, Fig. 22 depicts the comparison between the separations at t = 1, 3, 5

and 7. We can observe that the separation appears with no deformation until t = 5, and a deformation is visible at t = 7. The results are in good agreement with those obtained by the FDM [17, 20], the BEM [23, 24], the DQM [27], the MRKPR [34], the RPIM [35], the MLPG [36] and the MLSRRP [37].  Figs. 23 and 24 give the numerical results of the velocity ∂u/∂t and the acceleration ∂ 2 u ∂t2 , respectively. The two figures show the prominent performance of the WP-MLS approximation. Clearly, no reflection

appears on the boundaries of the velocity field, and no perturbation appears on the boundaries of the acceleration field. The results are in good agreement with the results given in Ref. [22]. 30

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Fig. 22. Comparison of u for the superposition of two orthogonal line solitons when β = 0.

Fig. 23. Numerical results of the velocity ∂u/∂t for the superposition of two orthogonal line solitons when β = 0.

31

AN US

CR IP T

ACCEPTED MANUSCRIPT

M

 Fig. 24. Numerical results of the acceleration ∂ 2 u ∂t2 for the superposition of two orthogonal line solitons when

AC

CE

PT

ED

β = 0.

Fig. 25. Comparison of u obtained using β = 0 and 0.5 for the superposition of two orthogonal line solitons.

32

ACCEPTED MANUSCRIPT

To study the influence of the dissipative term, the numerical results of u obtained at t = 3 and 7 with β = 0 and 0.5 are displayed in Fig. 25. It is true that the movement of the solitons is delayed by the dissipative term. At t = 3, the delay is visible. While at t = 7, the delay is very obvious. This conclusion agrees with that given in Refs. [17, 20, 23, 24, 33, 34]. For the sine-Gordon equation given in Eq. (55), it is well known that the energy given by [17, 20, 21,

E (t) = E (0) =

1 2

Z h Ω

i 2 u2t + |∇u| + 2 (1 − cos u) dΩ

CR IP T

23, 24, 34–36, 38]

is conserved when β = 0. The integration in the E (t) can be carried out by the composite trapezoidal rule. To study the conservation of the energy, the values of E (t) obtained using β = 0 and 0.5 are displayed in Fig. 26. We can see that the energy is conserved when β = 0. While for β = 0.5, the energy E (t) decreases

ED

M

AN US

gradually as t increases.

PT

Fig. 26. Conservation of the energy E (t).

CE

5.2.2. Line soliton in an inhomogeneous medium This type of soliton can be simulated by choosing ψ (x1 , x2 ) = 1 + sech2

AC

24, 27, 33]

φ1 (x1 , x2 ) = 4 tan−1 exp



x1 − 3.5 0.954



,

φ2 (x1 , x2 ) = 0.629sech



x1 − 3.5 0.954

p

x21 + x22 , β = 0 and [19–



x1 , x2 ∈ (−7, 7) .

,

Fig. 27 displays the numerical results in term of sin (u/2). The results are obtained at t = 0, 6, 12

and 18 using τ = 0.25 and h = 2/3. We can see that the soliton moves straightly at the beginning, and a deformation appears at t = 6. Then, the movement is impeded when t tends to 12. Finally, the soliton regains the straightness at t = 18.

The graphical results agree well with those obtained by the FDM

[19–21], the FEM [22], the BEM [23, 24], the DQM [27] and the RBF [33]. 33

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 27. Surface plots of sin (u/2) (left panel) and the associated contours (right panel) for the line soliton in an inhomogeneous medium.

34

ACCEPTED MANUSCRIPT

5.2.3. Elliptical breather soliton This type of soliton can be simulated by choosing ψ (x1 , x2 ) = 1, β = 0 and [21, 22, 24, 27, 35, 36] s    2 2 (x − x ) (x + x ) 1 2 1 2  , φ2 (x1 , x2 ) = 0, x1 , x2 ∈ (−7, 7) . φ1 (x1 , x2 ) = 4 tan−1 2sech 0.86 + 3 2

Fig. 28 depicts the numerical results for τ = 0.05 and h = 2/7 at times t = 0, 1.6, 8.0, 11.2, 12.8 and

CR IP T

15.2. We can observe from this figure that the major axes of the breather from the initial direction x2 = x1 are turned and shrank until t = 1.6, and a reflection phase is visible at t = 8.0. Then, at t = 12.8, the major axes nearly regains the initial direction. Finally, an expanding phase is obvious at t = 12.8 and 15.2. The graphical results agree well with those obtained by the FDM [21], the FEM [22], the BEM [24], the DQM

AC

CE

PT

ED

M

AN US

[27], the RPIM [35] and the MLPG [36].

Fig. 28. Surface plots of u (left panel) and the associated contours (right panel) for the elliptical breather soliton.

35

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 28. (continued)

36

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 29. Surface plots of sin (u/2) (left panel) and the associated contours (right panel) for the collision of circular solitons.

37

Fig. 29. (continued)

5.2.4. Collision of four circular solitons

CR IP T

ACCEPTED MANUSCRIPT

φ1 (x1 , x2 ) = 4 tan−1 exp φ2 (x1 , x2 ) = cosh





AN US

This type of soliton can be simulated by choosing ψ (x1 , x2 ) = 1, β = 0.05 and [17, 20–22, 24, 33, 36, 38]

4−

4−

q

q

2

(x1 + 3) + (x2 + 3) 4.13 2

(x1 + 3) + (x2 + 3)

2

2





 0.436 ,

0.436

,

M

in the domain −30 ≤ x1 , x2 ≤ 10. As in Refs. [17, 20–22, 24, 33, 36, 38], the approximate solution is computed by the WP-MLS approximation on one quarter of the domain. Then, it is extended across the

ED

lines x1 = −10 and x2 = −10 using symmetry relations. Fig. 29 depicts the numerical results obtained at t = 0, 1.5, 3, 6 and 9 using h = 2/3 and τ = 0.1. We can observe that the four circular solitons are shrank until t = 3. Then, an expanding phase is obvious and

PT

a bigger ring soliton appears at t = 6 and 9. In this example, the process of shrinking and expanding is interactive and very complex. The graphical results agree well with those obtained by the FDM [17, 20, 21],

CE

the FEM [22], the BEM [24], the RBF [33], the MLPG [36] and the EFG [38]. For a large β, the dissipative term is expected to delay the process of shrinking and expanding. Fig. 30 depicts the numerical results obtained at t = 6 for β = 0.05, 0.4, 2 and 5.

From Figs. 29 and 30, we can

AC

see that with β = 5 the solitons at t = 6 are similar to those shown at t = 1.5 for β = 0.05.

38

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 30. Surface plots of sin (u/2) (left panel) and the associated contours (right panel) for the collision of circular solitons obtained at t = 6 for β = 0.05, 0.4, 2 and 5.

39

ACCEPTED MANUSCRIPT

6. Conclusions A meshless method based on the WP-MLS approximation is developed in this paper for the numerical analysis of a class of nonlinear generalized Klein-Gordon equations. In this method, no mesh is required to discretize the problem domain, and the approximate solution is generated entirely based on a cluster of scattered nodes.

CR IP T

To demonstrate the accuracy, convergence and efficiency of the present method, numerical examples involving Klein-Gordon, Dodd-Bullough-Mikhailov, sine-Gordon, double sine-Gordon and sinh-Gordon equations with known analytical solutions are given. For all examples, by decreasing the time and the space steps, we find that the errors decrease and thus convergent solutions are obtained. Besides, numerical re

sults indicate that the method has a experimental convergence order of about O τ 2 in time and O h2

in space. Therefore, our numerical scheme is of second order in both time and space variables. Moreover,

AN US

the errors and condition numbers are compared with those obtained using the RBF and MLS methods. These comparisons indicate that the present method has higher precision and convergence order, and lower condition numbers.

To further illustrate the performance and efficiency of the present method, numerical experiments involving more complex line and ring solitons are also considered. The numerical results agree well with previously

M

reported numerical results.

The results obtained in this study show that the present method is well-suited for solving nonlinear generalized Klein-Gordon equations. This meshless method is extensible to solve more complicated time

PT

Acknowledgements

ED

dependent nonlinear problems in mathematical physics. Nevertheless, more research work is required.

This research was supported by the National Natural Science Foundation of China under the Grant No.

CE

11471063, the Chongqing Research Program of Basic Research and Frontier Technology under the Grant No. cstc2015jcyjBX0083 and the Scientific and Technological Research Program of Chongqing Municipal

AC

Education Commission under the Grant No. KJ1600330.

References

[1] R.K. Dodd, J.C. Eilbeck, J.D. Gibbon, H.C. Morris, Solitons in Nonlinear Wave Equations, Academic Press, New York, 1982.

[2] P.J. Drazin, R.S. Johnson, Soliton: An Introduction, Cambridge University Press, Cambridge, UK, 1989. [3] W. Greiner, Relativistic Quantum Mechanics-Wave Equations, third ed., Springer, Berlin, 2000. [4] A.M. Wazwaz, New travelling wave solutions to the Boussinesq and the Klein-Gordon equations, Commun. Nonlinear. Sci. 13 (2008) 889-901.

40

ACCEPTED MANUSCRIPT

[5] H. Han, Z. Zhang, Split local absorbing conditions for one-dimensional nonlinear Klein-Gordon equation on unbounded domain, J. Comput. Phys. 227 (2008) 8992-9004. [6] A.M. Wazwaz, The tanh and the sine-cosine methods for compact and noncompact solutions of the nonlinear Klein-Gordon equation, Appl. Math. Comput. 167 (2005) 1179-1195. [7] M. Dehghan, A. Shokri, Numerical solution of the nonlinear Klein-Gordon equation using radial basis functions, J. Comput. Appl. Math. 230 (2009) 400-410.

tau meshless method, Comput. Phys. Comm. 185 (2014) 1399-1409.

CR IP T

[8] W. Shao, X. Wu, The numerical solution of the nonlinear Klein-Gordon and Sine-Gordon equations using the Chebyshev

[9] G. Leibbrandt, New exact solutions of the classical sine-Gordon equation in 2+1 and 3+1 dimensions, Phys. Rev. Lett. 41 (1978) 435-438.

[10] J. Zagrodzinsky, Particular solutions of the sine-Gordon equation in 2+1 dimensions, Phys. Lett. A 72 (1979) 284-286. [11] P. Kaliappan, M. Lakshmanan, Kadomtsev-Petviashvili and two-dimensional sine-Gordon equations: Reduction to Painlev` e transcendents, J. Phys. A: Math. Gen. 12 (1979) 249-252.

[12] E.J. Parkes, B.R. Duffy, P.C. Abbott, The Jacobi elliptic-function method for finding periodic-wave solutions to nonlinear

AN US

evolution equations, Phys. Lett. A 295 (2002) 280-286.

[13] A.M. Wazwaz, Exact solutions for the generalized sine-Gordon and the generalized sinh-Gordon equations, Chaos Solitons Fractals 28 (2006) 127-135.

[14] A.M. Wazwaz, The tanh method: solitons and periodic solutions for the Dodd-Bullough-Mikhailov and the TzitzeicaDodd-Bullough equations, Chaos Solitons Fractals 25 (2005) 55-63.

[15] Sirendaoreji, Auxiliary equation method and new solutions of Klein-Gordon equations, Chaos Solitons Fractals 31 (2007) 943-950.

Math. Lett. 25 (2012) 2354-2358.

M

[16] A.M. Wazwaz, One and two soliton solutions for the sinh-Gordon equation in (1+1), (2+1) and (3+1) dimensions, Appl.

[17] Z. Liang, Y. Yan, G. Cai, A Dufort-Frankel difference scheme for two-dimensional sine-Gordon equation, Discrete Dyn.

ED

Nat. Soc. 2014 (2014) 784387.

[18] M. Lakestani, M. Dehghan, Collocation and finite difference-collocation methods for the solution of nonlinear Klein-Gordon equation, Comput. Phys. Comm. 181 (2010) 1392-1401. [19] P. Christiansen, P. Lomdahl, Numerical solution of 2+1 dimensional sine-Gordon solitons, Physica D 2 (1981) 482-494.

PT

[20] K. Djidjeli, W.G. Price, E.H. Twizell, Numerical solutions of a damped sine-Gordon equation in two space variables, J. Eng. Math. 29 (1995) 347-369.

[21] A.G. Bratsos, The solution of the two-dimensional sine-Gordon equation using the method of lines, J. Comput. Appl.

CE

Math. 206 (2007) 251-277.

[22] J. Argyris, M. Haase, J.C. Heinrich, Finite element approximation to two- dimensional sine-Gordon solitons, Comput. Methods Appl. Mech. Eng. 86 (1991) 1-26.

AC

[23] M. Dehghan, D. Mirzaei, The dual reciprocity boundary element method (DRBEM) for two-dimensional sine-Gordon equation, Comput. Methods Appl. Mech. Eng. 197 (2008) 476-486.

[24] D. Mirzaei, M. Dehghan, Boundary element solution of the two-dimensional sine-Gordon equation using continuous linear elements, Eng. Anal. Bound. Elem. 33 (2009) 12-24.

[25] M. Dehghan, A. Ghesmati, Application of the dual reciprocity boundary integral equation technique to solve the nonlinear Klein-Gordon equation, Comput. Phys. Comm. 181 (2010) 1410-1418. [26] J. Rashidinia, R. Mohammadi, Tension spline approach for the numerical solution of nonlinear Klein-Gordon equation, Comput. Phys. Comm. 181 (2010) 78-91. [27] R. Jiwari, S. Pandit, R. Mittal, Numerical simulation of two-dimensional sine-Gordon solitons by differential quadrature

41

ACCEPTED MANUSCRIPT

method, Comput. Phys. Comm. 183 (2012) 600-616. [28] S. Li, W.K. Liu, Meshfree Particle Methods, Springer, Berlin, 2004. [29] G.R. Liu, Meshfree Methods: Moving beyond the Finite Element Method, second ed., CRC Press, Boca Raton, 2009. [30] A. Hussain, S. Haq, M. Uddin, Numerical solution of Klein-Gordon and sine-Gordon equations by meshless method of lines, Eng. Anal. Bound. Elem. 37 (2013) 1351-1366. [31] Z.W. Jiang, R.H. Wang, Numerical solution of one-dimensional Sine-Gordon equation using high accuracy multiquadric quasi-interpolation, Appl. Math. Comput. 218 (2012) 7711-7716.

Probl. Eng. 2014 (2014) 383219.

CR IP T

[32] Q. Wei, R.J. Cheng, The improved moving least-square Ritz method for the one-dimensional sine-Gordon equation, Math.

[33] M. Dehghan, A. Shokri, A numerical method for solution of the two-dimensional sine-Gordon equation using the radial basis functions, Math. Comput. Simulat. 79 (2008) 700-715.

[34] R.J. Cheng, K.M. Liew, Analyzing two-dimensional sine-Gordon equation with the mesh-free reproducing kernel particle Ritz method, Comput. Methods Appl. Mech. Eng. 245-246 (2012) 132-143.

[35] M. Dehghan, A. Ghesmati, Numerical simulation of two-dimensional sine-Gordon solitons via a local weak meshless

AN US

technique based on the radial point interpolation method (RPIM), Comput. Phys. Comm. 181 (2010) 772-786.

[36] D. Mirzaei, M. Dehghan, Meshless local Petrov-Galerkin (MLPG) approximation to the two dimensional sine-Gordon equation, J. Comput. Appl. Math. 233 (2010) 2737-2754.

[37] R. Salehi, M. Dehghan, A moving least square reproducing polynomial meshless method, Appl. Numer. Math. 69 (2013) 34-58.

[38] X.L. Li, S.G. Zhang, Y. Wang, H. Chen, Analysis and application of the element-free Galerkin method for nonlinear sine-Gordon and generalized sinh-Gordon equations, Comput. Math. Appl. 71 (2016) 1655-1678.

M

[39] M. Dehghan, M. Abbaszadeh, A. Mohebbi, The numerical solution of the two-dimensional sinh-Gordon equation via three meshless methods, Eng. Anal. Bound. Elem. 51 (2015) 220-235. [40] M. Dehghan, M. Abbaszadeh, A. Mohebbi, An implicit RBF meshless approach for solving the time fractional nonlinear

ED

sine-Gordon and Klein-Gordon equations, Eng. Anal. Bound. Elem. 50 (2015) 412-434. [41] P. Lancaster, K. Salkauskas, Surface generated by moving least squares methods, Math. Comput. 37 (1981) 141-158. [42] K.M. Liew, Y.M. Cheng, S. Kitipornchai, Boundary element-free method (BEFM) and its application to two-dimensional elasticity problems, Int. J. Numer. Methods Eng. 65 (2006) 1310-1332.

PT

[43] M.J. Peng, Y.M. Cheng, A boundary element-free method (BEFM) for two-dimensional potential problems, Eng. Anal. Bound. Elem. 33 (2009) 77-82.

[44] Y.M. Cheng, M.J. Peng, Boundary element-free method for elastodynamics, Science in China, Series G 48 (2005) 641-657.

CE

[45] Y.Y. Lu, T. Belytschko, L. Gu, A new implementation of the element free Galerkin method, Comput. Methods Appl. Mech. Eng. 113 (1994) 397-414. [46] R.L. Burden, J.D. Faires, Numerical Analysis, ninth ed., Cengage Learning, Boston, 2011.

AC

[47] X.L. Li, H. Chen, Y. Wang, Error analysis in Sobolev spaces for the improved moving least-square approximation and the improved element-free Galerkin method, Appl. Math. Comput. 262 (2015) 56-78.

[48] X.L. Li, S.G. Zhang, Meshless analysis and applications of a symmetric improved Galerkin boundary node method using the improved moving least-square approximation, Appl. Math. Model. 40 (2016) 2875-2896.

[49] X.L. Li, J.L. Zhu, A Galerkin boundary node method and its convergence analysis, J. Comput. Appl. Math. 230 (2009) 314-328. [50] X.L. Li, Meshless Galerkin algorithms for boundary integral equations with moving least square approximations, Appl. Numer. Math. 61 (2011) 1237-1256. [51] X.L. Li, Error estimates for the moving least-square approximation and the element-free Galerkin method in n-dimensional

42

ACCEPTED MANUSCRIPT

spaces, Appl. Numer. Math. 99 (2016) 77-97. [52] X.L. Li, S.L. Li, On the stability of the moving least squares approximation and the element-free Galerkin method, Comput. Math. Appl. 72 (2016) 1515-1531. [53] D. Mirzaei, R. Schaback, M. Dehghan, On generalized moving least squares and diffuse derivatives, IMA J. Numer. Anal. 32 (2012) 983-1000. [54] A. Bouhamidi, M. Hachedb, K. Jbilou, A meshless RBF method for computing a numerical solution of unsteady Burgers’-

AC

CE

PT

ED

M

AN US

CR IP T

type equations, Comput. Math. Appl. 68 (2014) 238-256.

43