More accurate results for nonlinear generalized Benjamin-Bona-Mahony-Burgers (GBBMB) problem through spectral meshless radial point interpolation (SMRPI)

Engineering Analysis with Boundary Elements 72 (2016) 42–54

Contents lists available at ScienceDirect

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

More accurate results for nonlinear generalized Benjamin-Bona-Mahony-Burgers (GBBMB) problem through spectral meshless radial point interpolation (SMRPI) Elyas Shivanian n, Ahmad Jafarabadi Department of Mathematics, Imam Khomeini International University, Qazvin 34149-16818, Iran

art ic l e i nf o

a b s t r a c t

Article history: Received 16 February 2016 Received in revised form 3 July 2016 Accepted 8 August 2016

In this article, the spectral meshless radial point interpolation (SMRPI) technique is applied to the nonlinear generalized Benjamin-Bona-Mahony-Burgers (GBBMB) in two-dimension with initial and Dirichlet-type boundary conditions. This method is based on erudite combination of meshless methods and spectral collocation techniques. The point interpolation method with the help of radial basis functions is used to construct shape functions which play as basis functions in the frame of SMRPI. To treat the nonlinearity part, a kind of predictor-corrector scheme combined with Crank-Nicolson technique is adopted. We prove that the time discrete scheme is stable respect to the time variable in H1 and convergent with convergence order 6 (δt ) . To show the high accuracy of SMRPI, a comparison study of the present method and recently applied interpolating element-free Galerkin technique is given through applying on GBBMB equation. The results reveal that the method is more accurate and possesses low complexity. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Spectral meshless radial point interpolation (SMRPI) method Nonlinear generalized Benjamin-BonaMahony-Burgers equation Radial basis function Forward finite difference method

1. Introduction There are no generally applicable methods to solve nonlinear PDEs. Therefore, the numerical methods are important in the field of nonlinear PDEs. Finite difference (FDM), finite volume (FVM), and finite element (FEM) methods have been historically used to model a wide variety of engineering problems in complex geometries that may require extensive meshing which is a major difficulty of these methods. In order to overcome the mentioned shortcoming some methods, the so-called meshless methods, have been proposed [1–3]. Meshless methods are uniquely simple, yet provide solution accuracies for certain classes of equations that rival those of finite elements and boundary elements without requiring the need for mesh connectivity. Ease in programming, no domain or surface discretization, no numerical integration, and similar formulations for 2-D and 3-D make these methods very attractive. There exist several types of meshless methods such as meshless methods based on weak forms such as the element free Galerkin (EFG) method (and its developments like for example improved EFG (IEFG), complex variable EFG (CVEFG) and improved complex variable EFG (ICVEFG)) [4–13], meshless local Petrov-Galerkin

(MLPG) method [14–17] and MLPG based on the particular solutions (MLPG-PS), meshless techniques based on collocation techniques (strong forms) such as the meshless collocation technique based on radial basis functions (RBFs) [18–21] and the method of approximate particular solutions (MAPS) [22], and finally meshless techniques based on the combination of weak forms and collocation technique [23–27]. One of the methods that has been used to solve numerically differential equations is spectral collocation method. The idea is to choose a finite-dimensional space of candidate solutions (usually, polynomials up to a certain degree) and a number of points in the domain (called collocation points), and to select that solution which satisfies the given equation at the collocation points. In spectral meshless radial point interpolation (SMRPI) method [28– 30] the point interpolation method with the help of radial basis functions is proposed to construct shape functions which have Kronecker delta function property and are used as basis function in the frame of SMRPI. In this article, we apply the SMRPI for solving initial boundary value problem of the generalized Benjamin-Bona-Mahony-Burgers (GBBMB) equations in the following form: ut − Δut − β Δu + α∇ . u = γ ∇ . F (u) + f ,

x = (x, y )T ∈ Ω ,

t ∈ [0, T ], (1)

with initial condition n

Corresponding author. E-mail address: [email protected] (A. Jafarabadi).

http://dx.doi.org/10.1016/j.enganabound.2016.08.006 0955-7997/& 2016 Elsevier Ltd. All rights reserved.

u (x, 0) = g (x),

x ∈ Ω,

(2)

E. Shivanian, A. Jafarabadi / Engineering Analysis with Boundary Elements 72 (2016) 42–54

and Dirichlet-type boundary condition

u (x, t ) = h (x, t ),

x ∈ ∂Ω,

and the polynomial moment matrix is

(3)

where F(u) is the nonlinear function, f is a function of space and time variables, Ω ⊆ 2, also Δ and ∇are operators of Laplacian and gradient, respectively. In (1), we assume α and β are positive constants and γ is real number. In Eq. (1), if we put β = 0, α = 1 and γ = − 1 then we obtain the regularized long- wave equation that was studied by Benjamin et al. [32] as an improvement of the Korteweg-de Vries equation (KdV equation) for modeling long surface gravity waves of small amplitude propagating uni-directionally in 1 þ1 dimensions [31–33]. To see the numerical works on GBBMB refer to [34,35]. Recently in [36,37], the GBBMB was solved by interpolating element-free Galerkin technique and the meshless method of radial basis functions. Our aim is to treat the GBBMB via the SMRPI in order to obtain more accurate results in comparison to these two late references. Furthermore, we prove stability and convergence of the present method with respect to the time variable in H1 with convergence order 6 (δt ). Therefore, the framework of this paper arranged as follows: In Section 2, the SMRPI is introduced and the high order operational matrices are obtained. In Section 3, we obtain a time discrete scheme using the finite difference for Eq. (1). In Section 4, it is presented the stability and the convergence theorem. Three numerical examples are solved by the present method with error and comparison analysis in Section 5. Finally, a concluding remarks are given in Section 6.

To prevent the singularity difficulty in the polynomial point interpolation method (PIM), the radial basis function (RBF) is used to develop the radial point interpolation method (RPIM) [38–40] for constructing shape functions as basis functions in SMRPI. Here continuous function u (x, t ) defined in a domain Ω, is approximated as follows:

u (x) =

m

∑ Ri (x) ai + ∑ Pj (x) bj = RT (x) a + PT (x) b, i=1

j=1

(4)

(5)

where the vector of function values Us is

Us = {u1, u2 , u3, …, un }T ,

(6)

the RBFs moment matrix is

⎡ R1 (r1) R2 (r1) ⎢ R (r ) R2 (r2 ) Rn = ⎢ 1 2 ⎢ ⋮ ⋮ ⎢ ⎣ R1 (rn ) R2 (rn )

⋯ Rn (r1) ⎤ ⎥ ⋯ Rn (r2 ) ⎥ , ⋱ ⋮ ⎥ ⎥ ⋯ Rn (rn )⎦n × n

x1 y1 ⋯ Pm (x1) ⎤ ⎥ x2 y2 ⋯ Pm (x2) ⎥ ⋮ ⋮ ⋱ ⋮ ⎥ ⎥ xn yn ⋯ Pm (x n)⎦

, (8)

n×m

Also, the vector of unknown coefcients for RBFs is

aT = {a1, a2, …, an } ,

(9)

and the vector of unknown coefcients for polynomial is

bT = {b1, b2 , …, bm } .

(10)

Note that in Eq. (7), rk in Ri (rk ) is defined as

rk =

(xk − xi )2 + (yk − yi )2 .

(11)

Because in Eq. (5) there are m + n unknown variables, for guarantee uniqueness of approximation, add m condition as follows n

∑ Pj (x i) ai = PTm a = 0,

j = 1, 2, …, m.

i=1

(12)

Now by combining Eqs. (5) and (12) we obtain the following system of equations in the matrix form:

⎡ U ⎤ ⎡ R n Pm ⎤ ⎡ a ⎤ ∼, ⎥ ⎢ ⎥ = Ga U͠ s = ⎢ s ⎥ = ⎢ T ⎢⎣ 0 ⎥⎦ ⎣⎢ Pm 0 ⎥⎦ ⎣⎢ b ⎥⎦

T U͠ s = {u1, u2 , …, un , 0, 0, …, 0} , ∼ aT = {a , a , … , a , b , b , … , b } . 1

2

n

1

2

m

(13)

(7)

(14)

Solving Eq. (13), we obtain

⎡ a⎤ ∼ a = ⎢ ⎥ = G−1U͠ s . ⎣ b⎦

(15)

Now by rewrite of Eq. (4) have

u (x) = {RT (x), PT (x)} G−1U͠ s = Φ͠ T (x) U͠ s ,

where Ri (x) is a radial basis function (RBF), n is the number of RBFs, Pj (x) is monomial in the space coordinate x, and m is the number of polynomial basis functions. When m ¼0, only RBFs are used. Otherwise, the RBF is augmented with m polynomial basis functions. Coefficients ai and bj are unknown which should be determined. To determine the ai and bj in Eq. (4), a support domain is formed for the point of interest at x , and n field nodes are included in the support domain (support domain is usually a disk with radius rs). Coefcients ai and bj in Eq. (4) can be determined by enforcing Eq. (4) to be satisfied at these n nodes surrounding the point of interest x . Therefore we have a linear system algebraic equations, Corresponding each node. these equations in matrix form are as follows

Us = R na + Pm b,

⎡1 ⎢ ⎢1 Pm = ⎢ ⋮ ⎢ ⎣1

where

2. Spectral meshless radial point interpolation (SMRPI) scheme

n

43

(16)

where

Φ͠ T (x) = {RT (x), PT (x)}G−1 = {ϕ1, ϕ2, …, ϕn, ϕn + 1, ϕn + 2, …, ϕn + m } .

(17)

The first n functions of the above vector function are called the RPIM shape functions corresponding to the nodal displacements and we show by the vector ΦT (x), that's mean

ΦT (x) = {ϕ1, ϕ2, …, ϕn } .

(18)

Now Eq. (16) Becomes n

u (x) = ΦT (x) Us =

∑ ϕi (x) ui . i=1

(19)

R−n 1

Because usually exists for arbitrary scattered nodes and therefore the completed matrix G is theoretically non-singular [41,42]. It is significant that the RPIM shape functions have the Kronecker delta function property, that is

⎧ 1 for i = j, j = 1, 2, …, n ϕi (x j ) = ⎨ ⎩ 0 for i ≠ j, i, j = 1, 2, …, n

(20)

This is because the RPIM shape functions are created to pass through nodal values. In the current work, we assume that the number of total nodes covering Ω = (Ω ∪ ∂Ω ) is N. Also, consider the nx , instead n, the

44

E. Shivanian, A. Jafarabadi / Engineering Analysis with Boundary Elements 72 (2016) 42–54

number of nodes included in support domain Ω x corresponding to the point of interest x (for example Ω x can be a disk centered at x with radius rs). By rewriting the Eq. (19), we have

=∂ sϕj (x) /∂y s = 0, s = 1, 2, … , due to Eq. (22). To see a practical example about this feature consider [29].

N

u (x) = ΦT (x) Us =

∑ ϕj (x) uj .

(21)

j=1

Since corresponding to node xj there is a shape function ϕj (x) , j = 1, 2, … , N , we define Ω xc = {xj: xj ∉ Ω x } then obviously from Eq. (20)

∀ x j ∈ Ω xc : ϕj (x) = 0.

(22)

Now the derivatives of u (x) respect to x and y are

∂u (x) = ∂x

N

∑

∂ϕj (x) ∂x

j=1

uj ,

∂u (x) = ∂y

N

∑

∂ϕj (x) ∂y

j=1

uj ,

(23)

and for high derivatives of u (x)

∂ su (x) = ∂x s

N

∑

∂ sϕj (x) ∂x s

j=1

∂ su (x) = ∂y s

uj ,

N

∑ j=1

∂ sϕj (x) ∂y s

∂ϕN (x1) ⎞ ⎟ ∂x ⎟ ⎛ u1 ⎞ ∂ϕN (x2) ⎟ ⎜ ⎟ ⋯ ⎟ u2 ∂x ⎟ ⎜ ⋮ ⎟, ⎜ ⎟ ⋱ ⋮ ⎟ ⎜⎝ uN ⎟⎠ ⎟ ∂ϕN (x N ) ⎟ ⋯ ∂x ⎠

⎛ ∂ϕ1 (x1) ∂ϕ2 (x1) ⎜ ∂y ⎛ u ′y ⎞ ⎜ ∂y 1 ⎜ ⎟ ⎜ ∂ϕ1 (x2) ∂ϕ2 (x2) ⎜ u ′y2 ⎟ ⎜ ∂y ⎜ ⋮ ⎟ = ⎜ ∂y ⎜⎜ ⎟⎟ ⎜ ⋮ ⋮ ⎝ u ′yN ⎠ ⎜ ∂ϕ (x ) ∂ϕ (x ) N 2 N ⎜⎜ 1 ∂y ⎝ ∂y

∂ϕN (x1) ⎞ ⎟ ∂y ⎟ ⎛u ⎞ ∂ϕN (x2) ⎟ ⎜ 1 ⎟ ⋯ ⎟ u2 ∂y ⎟ ⎜ ⎟. ⋮ ⎟ ⎜⎜ ⎟⎟ ⋱ ⋮ ⎟ ⎝ uN ⎠ ∂ϕN (x N ) ⎟ ⋯ ⎟ ∂y ⎠

⋯

D y(sij)

=

∂x s

(25)

=λ

∂ sϕj (x i ) , ∂y s

−

⎞ ⎛ ∂uk (x) ∂uk (x) ⎟ + ⎜ ⎟ ∂y ⎟ ⎝ ∂x ⎠

) + ∂F ( u^ (x)) ⎞⎟⎟ + f k

∂y

⎟ ⎠

k + 1 (x)

+ f k (x), (33)

2 . Consider N nodes with arbitrary distribution on the δt boundary and domain of the problem, i.e. Ω = Ω ∪ ∂Ω . Assuming that u (x i, kδt ), i¼1,2,…, N, are known, our aim is to find unknowns u (x i, (k + 1) δt ), i¼1,2,…,N. For nodes which are located in the interior of the domain, i.e. for x i ∈ Ω , to obtain the discrete equations from (32), substituting approximation formulas (21), (23) and (24) into Eq. (34) yields

(26)

)T

⎛ N λ ⎜ ∑ ϕj (x) ukj + 1 − ⎜ ⎝ j =1 N

(28)

+

∑

∂ 2ϕj (x) ∂y 2

(29)

(31) ∀ xj ∈ Ω xc : ∂ sϕj (x) /∂xs

∂ 2ϕj (x) ∂x 2

j =1

N

∑ j =1

⎛ N ∂ 2ϕ (x) j + β⎜ ∑ ukj + 2 ⎜ ⎝ j = 1 ∂x +

∑

N

ukj + 1 −

∑

∂ 2ϕj (x) ∂y 2

j =1

⎛ N ∂ϕ (x) ⎞ j ukj + 1⎟ + α ⎜ ∑ ukj + 1 + ⎜ ⎟ ⎝ j = 1 ∂x ⎠

⎛ N = λ ⎜ ∑ ϕj (x) ukj − ⎜ ⎝ j =1

N

(30)

N

∑

∂ 2ϕj (x) ukj ∂x 2 N

∑ j =1

N

−

∑ j =1

N

∑ j =1

⎛ N ∂ 2ϕ (x) ⎞ j ukj + 1⎟ − β ⎜ ∑ ukj + 1 2 ⎜ ⎟ ⎝ j = 1 ∂x ⎠

∂ϕj (x) ∂y

⎞ ukj + 1⎟ ⎟ ⎠

∂ 2ϕj (x) ⎞⎟ ukj ⎟ ∂y 2 ⎠

⎛ N ∂ϕ (x) j − α⎜ ∑ ukj ⎜ ⎝ j = 1 ∂x ⎠

∂ 2ϕj (x) ⎞⎟ ukj ⎟ ∂y 2

⎛ ⎛ N ⎞ ⎛ N ⎞⎞ ∂ϕj (x) ⎞⎟ k k ∂ ∂ ⎜ ukj + 2γ ⎜ F ⎜⎜ ∑ ϕj (x) u^ j ⎟⎟ + F ⎜ ∑ ϕj (x) u^ j ⎟⎟ ⎟ ⎟ ⎜ ⎟ ∂y ⎠ ∂y ⎝ j = 1 ⎠⎠ ⎝ ∂x ⎝ j = 1 ⎠

+ f k + 1 (x) + f k (x),

At the end of this section, we mention

)+

βΔuk (x)−α

where λ =

j =1

U = (u1, u2 , …, uN )T .

uk (x)

Δuk (x)

(

(27)

Ux(s) = u y(s1) , u y(s2) , …, u y(sN)

(

uk (x)

⎛ ^k ⎜ ∂F u (x) + 2γ ⎜ ∂x ⎜ ⎝

U y(s) = D(s) yU

,

(32)

tk,

)

j =1

∂ sϕj (x i )

)

f k + 1 (x) + f k (x) , 2

(

where

Dx(sij) =

(

where δt = is the time step size, − = u (x, kδt ) and, Δ and ∇ are operators of Laplacian and gradient, respectively. Also, in the above equation, to eliminate the nonlinearity, an iterative approach in this case a simple predictor-corrector scheme combined with Crank-Nicolson technique is considered for ∇ . F (u) [43]. Eq. (32) can be written as:

t k+1

⋯

(

)

⎛ ⎞ ⎜ ∂uk + 1 (x) ∂uk + 1 (x) ⎟ + λ uk + 1 (x) − Δuk + 1 (x) − βΔuk + 1 (x) + α ⎜ ⎟ ∂x ∂y ⎟ ⎜ ⎝ ⎠

We denote the above coefcients matrices as Dx, Dy and also for high order derivatives we have the following matrix-vector multiplications:

)T

(

(24)

⎛ ∂ϕ1 (x1) ∂ϕ2 (x1) ⎜ ∂x ⎛ u x′ ⎞ ⎜ ∂x ⎜ 1 ⎟ ⎜ ∂ϕ (x2) ∂ϕ (x2) 1 2 ⎜ u x′ 2 ⎟ ⎜ ∂x ⎜ ⋮ ⎟ = ⎜ ∂x ⎜⎜ ⎟⎟ ⎜ ⋮ ⋮ ⎝ u x′ N ⎠ ⎜ ∂ϕ (x N ) ∂ϕ (x N ) 1 2 ⎜ ⎝ ∂x ∂x

(

uk + 1 (x) − uk (x) Δuk + 1 (x) − Δuk (x) Δuk + 1 (x) + Δuk (x) − −β 2 δt δt ⎛ ^k (x) ^k (x) ⎞ F u F u ∂ ∂ ⎜ ⎟ ∇. uk + 1 (x) + ∇. uk (x) +α = γ⎜ + ⎟ 2 x y ∂ ∂ ⎜ ⎟ ⎝ ⎠ +

∂ s (. ) ∂ s (. ) and are s′th derivatives with respect to x and y. s ∂x ∂y s s ∂ (. ) ∂ s (. ) , uy(s) ( . ) = and setting x = x i in Eq. Denoting ux(s) ( . ) = s ∂x ∂y s (23) then the following matrix form is given

Ux(s) = u x(s1) , u x(s2) , …, u x(sN) ,

In this section, a time discretization approach based on the forward finite difference will be applied for approximating the first-order derivative on time variable of the problem (1) at two consecutive time levels k and k + 1 as follows:

uj ,

where

Ux(s) = D(s) xU ,

3. The time discretization approximation and numerical implementation for SMRPI method

(34)

now, setting x = x i, i = 1, 2, … , NΩ (NΩ is the number of nodes in Ω) in Eq. (34) and applying notations (27)–(31) implies

E. Shivanian, A. Jafarabadi / Engineering Analysis with Boundary Elements 72 (2016) 42–54

to remove the non-linear term in Eq. (45), an iterative approach in this case a simple predictor-corrector scheme combined with Crank-Nicolson technique is performed as

N

λuik + 1 −

∑ ( λDx(2ij) + λD y(2ij) + βDx(2ij) + βD y(2ij) − αDx(1ij) − αD y(1ij) ) ukj + 1 j=1 N

= λuik −

∑ ( λDx(2ij) + λD y(2ij) − βDx(2ij) − βD y(2ij) + αDx(1ij) + αD y(1ij) ) ukj

AUlk + 1 = BU k +

j=1

⎛ ∂F u^ (x ) i ⎜ + 2γ ⎜ ∂u ⎝

(

) ∑D N

(1) ^ k x ij uj

+

(

∂F u^ (x i ) ∂u

j=1

+ f k + 1 (x i ) + f k (x i ),

) ∑D N

j=1

⎞

(1) ^ k ⎟ yij uj ⎟

+

⎠

λuik + 1 −

∑ j=1

( (λ + β)D N

= λuik −

∑ j=1

(2) x ij

(35)

(

+ (λ + β ) D y(2ij) − α Dx(1ij) + D y(1ij)

( (λ − β)D

⎛ ∂F u^ (x ) i ⎜ + 2γ ⎜ ∂u ⎝

(

(2) x ij

) ∑D N

(

∂F u^ (x i ) ∂u

j=1

+ f k + 1 (x i ) + f k (x i ),

))u

k+1 j

(

+

) ∑D N

j=1

))u

k j

⎞

i = 1, 2, …, NΩ.

⎠ (36)

i = 1, 2, …, N∂Ω.

(37)

By collecting Eqs. (36) and (37) the following matrix form can be obtained

^ ^ ^ ^ AU k + 1 = BU k + V *Ux + V *Uy + F k + 1,

(38)

T

(

)

U

+

k

(

^ ^ ^ ^ V *Ux + V *Uy

)

⎞ ⎟ ⎟ Ulk−+11 ⎠ (45)

F k + 1.

∥ Ulk + 1 − Ulk−+11 ∥∞ < ϵ,

(46)

ϵ is a

4. The stability and convergence analysis This content has been adapted from Ref. [36]. Consider modification of Eq. (32) with truncation error 6 (δt ):

(1) ^ k ⎟ yij uj ⎟

For nodes which are located on the boundary, using (3) we have the following relation:

uk + 1 (x i ) = h (x i, (k + 1) δt ),

(

where l is the number of iterations in each time level and small positive number.

+ (λ − β ) D y(2ij) + α Dx(1ij) + D y(1ij)

(1) ^ k x ij uj

⎛ 1⎜ ^ ^ ^ ^ V *Ux + V *Uy 2 ⎜⎝

This process is iterated until satisfying the following condition

i = 1, 2, …, NΩ,

or equivalently N

45

)

where U k + 1 = u1k + 1, u2k + 1, … , uNk+ 1 , and N × N matrices A , B and

uk + 1 (x) − uk (x) Δuk + 1 (x) − Δuk (x) Δuk + 1 (x) + Δuk (x) − −β δt δt 2 f k + 1 (x) + f k (x) ∇. uk + 1 (x) + ∇. uk (x) k +α = γ ∇. F (u ) + + 6 (δt ). (47) 2 2 Multiplying δt yields

βδt k + 1 αδt Δu (x) + ∇. uk + 1 (x) 2 2 βδt k αδt = uk (x) − Δuk (x) + Δu (x) − ∇. uk (x) + γδt∇. F (uk ) 2 2 δt + (f k + 1 (x) + f k (x)) + δt 6 (δt ), 2

uk + 1 (x) − Δuk + 1 (x) −

6 (δt 2),

and omitting the small term

^ ^ (1) ^ (1) N-vectors V , Ux , Uy , F k + 1 are defined as follows

(48)

we get the following equation

⎧ ⎛ ⎞ ⎪ λδ ij − ⎜ (λ + β ) D x(2) + (λ + β ) D y(2) − α D x(1) + D y(1) ⎟, i = 1, 2, … , NΩ ij ij ij ij ⎠ ⎝ Aij = ⎨ ⎪ i = NΩ + 1, … , N ⎩ δ ij,

(39)

⎧ ⎛ ⎞ ⎪ λδ ij − ⎜ (λ − β ) D x(2) + (λ − β ) D y(2) + α D x(1) + D y(1) ⎟, i = 1, 2, … , NΩ ij ij ij ij ⎠ ⎝ B ij = ⎨ ⎪ i = NΩ + 1, … , N ⎩ 0,

βδt αδt ΔU k + 1 (x) + ∇. U k + 1 (x) U k + 1 (x) − ΔU k + 1 (x) − 2 2 βδt αδt = U k (x) − ΔU k (x) + ΔU k (x) − ∇. U k (x) + γδt∇. F (U k ) 2 2 δt + (f k + 1 (x) + f k (x)). (49) 2

(40)

In Eqs. (48) and (49) uk (x) and U k (x) are exact and approximate solutions of Eq. (1), respectively. Let L2 (Ω ) be a square integrable function space with inner product and norm as follows:

(

(

)

)

⎧ ∂F u^ (x i ) ^ ⎪ , i = 1, 2, …, NΩ Vi = ⎨ 2γ ∂u ⎪ i = NΩ + 1, …, N , ⎩ 0,

(

〈u, v〉 =

)

⎧ N ⎪∑ D (1) u^j , i = 1, 2, …, NΩ ^ (Ux )i = ⎨ j = 1 xij ⎪ i = NΩ + 1, …, N , ⎩ 0,

⎧ N ⎪∑ D (1) u^j , i = 1, 2, …, NΩ ^ (Uy )i = ⎨ j = 1 yij ⎪ i = NΩ + 1, …, N , ⎩ 0,

Also

(41)

H1 (Ω )

H 1 (Ω)

∥ u ∥L2 (Ω) =

〈u, u〉 .

(50)

denotes the Sobolev-type space defined as

{

}

= u ∈ L2 (Ω): ∂ αu ∈ L2 (Ω) ,

(51)

∂|α|u , |α| ≤ 2, and α is multi-index notation for ∂x α1 ∂y α2 partial derivatives. The inner product, norm and seminorm of H1 (Ω ) are defined by where ∂ αu =

(42)

〈u, v〉1 = 〈u, v〉 + 〈∇. u, ∇. v〉, |u|1 =

(43)

H01 (Ω )

∥ u ∥H1(Ω) =

〈u, u〉1 ,

〈∇. u, ∇. u〉 .

is a subspace of

H01 (Ω) =

⎧ f k + 1 (x ) + f k (x ) i = 1, 2, …, N Ω i i Fik + 1 = ⎨ ⎩ h (x i, (k + 1) δt ) i = NΩ + 1, …, N .

∫Ω u (x) v (x) dΩ,

{ u ∈ H (Ω): u| 1

(52)

H1 (Ω ) ∂Ω

that defined as follows

}

=0 .

(53)

⎪

⎪

(44)

Note that, in Eq. (38), the operator “*” denotes Hadamard product (component by component production). Once again, we note that,

Theorem 4.1. Let U k ∈ H1 (Ω ) and we assume that F ( . ) has the Lipschitz condition as follows

|F (v1) − F (v2 )| ≤ 3|v1 − v2 | ,

∀ v1, v2 ∈ H 1 (Ω)

46

E. Shivanian, A. Jafarabadi / Engineering Analysis with Boundary Elements 72 (2016) 42–54

then for δt ≤

1 the scheme defined by (49) is stable in 3β + sign (γ ) γ 3 + α 2

H1, where sign ( . ) is Sign function.

Also, for right coefficient of (55), clearly

1 αδt γ 3δt 1 3βδt γ 3δt αδt + + ≤ + + sign (γ ) + , 2 4 2 2 4 2 2

^ k + 1 are exact and approximate soluProof. Assuming U k + 1 and U tions of (49) respectively. The round off error is to the following form

and

⎛ βδt ⎞ k + 1 αδt ⎟ Δe ek + 1 − ⎜ 1 + + ∇. e k + 1 ⎝ 2 ⎠ 2 ⎛ ⎛ βδt ⎞ k αδt ^k ⎞ ⎟ Δe − = ek − ⎜ 1 − ∇. e k + γδt∇. ⎜ F (U k ) − F (U ) ⎟, ⎝ ⎝ ⎠ 2 ⎠ 2

therefore

⎛ βδt ⎞ αδt ⎟ 〈Δe k + 1, e k + 1〉 + 〈e k + 1, e k + 1〉 − ⎜ 1 + 〈∇. e k + 1, e k + 1〉 ⎝ 2 ⎠ 2 ⎛ βδt ⎞ αδt ⎟ 〈Δe k , e k + 1〉 − = 〈e k , e k + 1〉 − ⎜ 1 − 〈∇. e k , e k + 1〉 ⎝ 2 ⎠ 2 ⎛ ^k ⎞ + γδt ∇. ⎜ F (U k ) − F (U ) ⎟, e k + 1 . ⎝ ⎠

⎛ βδt ⎞ αδt k + 1 ⎟ 〈∇. e k + 1, ∇. e k + 1〉 − + ⎜1 + 〈e , ∇. e k + 1〉 ⎝ 2 ⎠ 2 ⎛ βδt ⎞ ⎟ 〈∇. e k , ∇. e k + 1〉 = 〈e k , e k + 1〉 + ⎜ 1 − ⎝ 2 ⎠ +

αδt k ^k 〈e , ∇. e k + 1〉 − γδt〈(F (U k ) − F (U )), ∇. e k + 1〉. 2

Now using properties of inner product and Lipschitz condition, we obtain

⎛ βδt ⎞ αδt k + 1 2 ⎟ ‖∇. e k + 1‖22 + ⎜1 + − ‖e ‖ 2 L (Ω ) L (Ω ) ⎝ 2 ⎠ 4 1 1 αδt − ‖∇. e k + 1‖22 ≤ ‖e k + 1‖22 + ‖e k ‖22 L (Ω ) L (Ω ) L (Ω ) 4 2 2 1⎛ βδt ⎞ ⎟ ‖∇. e k ‖22 + ⎜1 − L (Ω ) 2⎝ 2 ⎠

‖e k + 1‖22

L (Ω )

0≤

3β + sign (γ ) γ 3 + α 2

≤⋯ ⎤k + 1 ⎡ 3βδt + sign (γ ) γ 3δt + αδt ⎥ ⎢1 + 2 ⎥ ‖e 0‖2 1 . ≤⎢ H (Ω ) ⎢ 1 − 3βδt − sign (γ ) γ 3δt − αδt ⎥ ⎦ ⎣ 2

1 3βδt γ 3δt αδt 1 3βδt γ 3δt αδt − − sign (γ ) − ≤ + − − . 2 4 2 2 2 4 2 2

(56)

Theorem 4.2. The time discrete solution is H1 - convergent with convergence order 6 (δt ). Proof. If subtract (48) and (49) from each other, obtain

⎛ βδt ⎞ k + 1 αδt ⎟ Δξ + ∇. ξ k + 1 ξk+1 − ⎜ 1 + ⎝ 2 ⎠ 2 ⎛ ⎛ βδt ⎞ k αδt ^k ⎞ ⎟ Δξ − = ξk − ⎜ 1 − ∇. ξ k + γδt∇. ⎜ F (U k ) − F (U ) ⎟ ⎝ ⎠ ⎝ ⎠ 2 2 + δ tRk + 1, ξk

and

0≤

⎡ ⎤2 3βδt + sign (γ ) γ 3δt + αδt ⎥ ⎢1 + 2 ⎥ ‖e k − 1‖2 1 ≤⎢ H (Ω ) ⎢ 1 − 3βδt − sign (γ ) γ 3δt − αδt ⎥ ⎣ ⎦ 2

Let {u (x, tk )}kn= 0 = uk ∈ H1 (Ω ) be the exact solution of (1)–(3) and {U (x, tk )}kn= 0 = U k ∈ H1 (Ω ) be approximate solution with same initial and boundary conditions. Now, the following theorem holds for the convergence of time discrete.

(55)

1 3βδt γ 3δt αδt 1 αδt − − sign (γ ) − ≤ − , 2 4 2 2 2 4

H (Ω )

Therefore, substituting in (56) and the square root of both sides, proof is complete. □

⎛1 ⎛1 3βδt αδt ⎞ k + 1 2 γ 3δt αδt ⎞ ⎜ − ⎟ ‖e ‖ 2 ⎟ ‖∇. e k + 1‖22 +⎜ + − − L (Ω ) L (Ω ) ⎝2 ⎝2 4 ⎠ 4 2 2 ⎠ ⎛1 αδt γ 3δt ⎞ k 2 ⎟ ‖e ‖ 2 ≤⎜ + + L (Ω ) ⎝2 4 2 ⎠

, we conclude

‖e k + 1‖2 1

⎤ ⎡ 3βδt ⎢1 + + sign (γ ) γ 3δt + αδt ⎥ 2 ⎥ ‖e k ‖2 ≤⎢ H1(Ω) 3βδt ⎥ ⎢ 1 3 sign t t − − ( γ ) γ δ − αδ ⎥ ⎢⎣ 2 ⎦

⎛ ⎡ 3β ⎤ ⎞ exp ⎜ ⎢ + sign (γ ) γ 3 + α⎥ T ⎟ ⎛ ⎡ 3β ⎦ ⎠ ⎤ ⎞ ⎝⎣ 2 = exp ⎜ 2 ⎢ + sign (γ ) γ 3 + α⎥ T ⎟. ⎛ ⎡ 3β ⎦ ⎠ ⎤ ⎞ ⎝ ⎣ 2 exp ⎜ −⎢ + sign (γ ) γ 3 + α⎥ T ⎟ ⎦ ⎠ ⎝ ⎣ 2

By simplification

1

⎛1 3βδt γ 3δt αδt ⎞ k + 1 2 ⎜ − ⎟ ‖e ‖ 1 − sign (γ ) − H (Ω ) ⎝2 4 2 2 ⎠ ⎛1 3βδt γ 3δt αδt ⎞ k 2 ⎟ ‖e ‖ 1 . ≤⎜ + + sign (γ ) + H (Ω ) ⎝2 4 2 2 ⎠

is equal to

1⎛ βδt ⎞ ⎜1 − ⎟ ‖∇. e k + 1‖22 L (Ω ) 2⎝ 2 ⎠ αδt k 2 αδt + ‖e ‖ 2 + ‖∇. e k + 1‖22 L (Ω ) L (Ω ) 4 4 γ 3δt k 2 γ 3δt + ‖e ‖ 2 + ‖∇. e k + 1‖22 . L (Ω ) L (Ω ) 2 2

From δt ≤

According to norm of H1 we have

When k + 1 approaches infinity, the coefficient of ‖e0‖2 1 in (56) H (Ω )

+

⎛1 βδt ⎞ ⎟ ‖∇. e k ‖22 . +⎜ − L (Ω ) ⎝2 4 ⎠

{

} }.

By dividing the coefficients and repeat inequality

Because ek + 1, ek ∈ H01 (Ω ), using divergence theorem and a special case of divergence theorem [44], yields

e k + 1〉

⎛1 3βδt γ 3δt αδt ⎞ ⎜ − ⎟ ‖e k + 1‖22 − sign (γ ) − + ‖∇. e k + 1‖22 L (Ω ) L (Ω ) ⎝2 4 2 2 ⎠ ⎛1 ⎞ 3βδt γ 3δt αδt ⎟ ‖e k ‖22 ≤ ⎜ + + sign (γ ) + + ‖∇. e k ‖22 L (Ω ) L (Ω ) ⎝2 4 2 2 ⎠

{

(54)

^ k+1 where ek + 1 = U k + 1 − U . Multiplying Eq. (54) by ek + 1 and integrating on Ω we have

〈e k + 1,

1 βδt 1 3βδt γ 3δt αδt − ≤ + + sign (γ ) + , 2 4 2 4 2 2

(57) uk

Uk

where and Rk + 1 is truncation error such that = − |Rk + 1| ≤ Cδt , for positive constant C. Now multiplying both sides of the above equation by ξk þ 1 and using properties of inner product and also, Lipschitz condition, we obtain

E. Shivanian, A. Jafarabadi / Engineering Analysis with Boundary Elements 72 (2016) 42–54

⎛ βδt ⎞ αδt k + 1 2 ⎟ ‖∇. ξ k + 1‖22 + ⎜1 + − ‖ξ ‖ 2 L (Ω ) L (Ω ) ⎝ 2 ⎠ 4 1 1 αδt − ‖∇. ξ k + 1‖22 ≤ ‖ξ k + 1‖22 + ‖ξ k‖22 L (Ω ) L (Ω ) L (Ω ) 4 2 2 1⎛ βδt ⎞ ⎟ ‖∇. ξ k‖22 + ⎜1 − L (Ω ) 2⎝ 2 ⎠

‖ξ k + 1‖22

L (Ω )

‖ξ k + 1‖2 1 H (Ω )

⎡ ⎢ ⎢ ≤⎢ ⎢ ⎣

1⎛ βδt ⎞ αδt k 2 ⎜1 − ⎟ ‖∇. ξ k + 1‖22 + ‖ξ ‖ 2 L (Ω ) L (Ω ) 2⎝ 2 ⎠ 4 αδt γ 3δt k 2 + ‖∇. ξ k + 1‖22 + ‖ξ ‖ 2 L (Ω ) L (Ω ) 4 2 t γ 3δt δ + ‖∇. ξ k + 1‖22 + ‖ξ k + 1‖22 L (Ω ) L (Ω ) 2 2 δt + ‖Rk + 1‖22 . L (Ω ) 2

⎡ 1 sign (γ ) γ 3δt δt αδt 3βδt + + + + ⎢ 2 2 2 2 ≤⎢ 2 sign (γ ) γ 3δt δt αδt 3βδt ⎢ 1 − − − − ⎢⎣ 2 2 2 2 2 k+1

δt

∑ j =1

If we assume δt ≤

⎤ ⎥ ⎥ 2 ⎥ ‖R k + 1‖L2 (Ω) ⎥ ⎦

≤⋯

By simplification

⎛1 δt αδt ⎞ k + 1 2 ⎜ − ⎟ ‖ξ ‖ 2 − L (Ω ) ⎝2 2 4 ⎠ ⎛1 3βδt γ 3δt αδt ⎞ ⎟ ‖∇. ξ k + 1‖22 +⎜ + − − L (Ω ) ⎝2 4 2 2 ⎠ ⎞ ⎛ ⎛1 1 αδt γ 3δt βδt ⎞ ⎟ ‖ξ k‖22 ⎟ ‖∇. ξ k‖22 ≤⎜ + + +⎜ − L ( Ω ) L (Ω ) ⎝2 ⎠ ⎝ 4 2 2 4 ⎠ δt + ∥ Rk + 1 ∥2 . 2

⎤ ⎥ ⎥ k 2 ⎥ ‖ξ ‖H1(Ω) ⎥ ⎦

⎡ ⎢ 1 + δt + αδt + 3βδt + sign (γ ) γ 3δt ⎢ 2 2 2 2 + δt ⎢ 2 sign (γ ) γ 3δt δt αδt 3βδt 1 − − − ⎢ 2 − 2 2 2 2 ⎣

+

0≤

sign (γ ) γ 3δt δt αδt 1 3βδt + + + + 2 2 2 2 2 sign (γ ) γ 3δt δt αδt 1 3βδt − − − − 2 2 2 2 2

47

⎡ ⎢ ⎢ ⎢ ⎢⎣

⎤k + 1 ⎥ ⎥ ‖ξ 0‖2 H1(Ω) ⎥ ⎥⎦

sign (γ ) γ 3δt δt αδt 1 3βδt + + + + 2 2 2 2 2 sign (γ ) γ 3δt δt αδt 1 3βδt − − − − 2 2 2 2 2

⎤j ⎥ ⎥ ‖R k − j + 2 ‖2 . L2 (Ω) ⎥ ⎥⎦

(59)

Clearly, taking into account ξ 0 = 0 and definition of δt , similar to the Theorem 4.1, we obtain

ξ k +1

2 H1(Ω)

k +1

≤ δt

∑ j =1

⎡ ⎢ ⎢ ⎢ ⎢⎣

sign (γ ) γ 3δt δt αδt 1 3βδt + + + + 2 2 2 2 2 sign (γ ) γ 3δt δt αδt 1 3βδt − − − − 2 2 2 2 2

⎤j ⎥ ⎥ ‖R k − j + 2 ‖2 L2 (Ω) ⎥ ⎥⎦

⎡ 1 sign (γ ) γ 3δt δt αδt 3βδt + + + + ⎢ 2 2 2 2 ≤ δt (k + 1) ⎢ 2 sign (γ ) γ 3δt δt αδt 3βδt ⎢ 1 ⎢⎣ 2 − 2 − 2 − 2 − 2

(58)

1 , it yields 1 + α + 3β + sign (γ ) γ 3

2

2 2T (1 + α + 3β + sign (γ ) γ 3)

≤ TC δt e

⎤k + 1 ⎥ ⎥ C 2δt 2 ⎥ ⎥⎦

.

(60)

1 δt αδt 3βδt sign (γ ) γ 3δt 1 δt αδt − − − − ≤ − − 2 2 2 2 2 2 2 4

By taking the square root of both sides (60) the proof is complete. □

1 δt αδt 3βδt sign (γ ) γ 3δt 1 3βδt γ 3δt αδt − − − − ≤ + − − . 2 2 2 2 2 2 4 2 2

5. Numerical experiments and comparison

and

0≤

Also, for right coefficients of (58) clearly

1 αδt γ 3δt 1 δt αδt 3βδt sign (γ ) γ 3δt + + ≤ + + + + , 2 4 2 2 2 2 2 2 and

1 βδt 1 δt αδt 3βδt sign (γ ) γ 3δt + ≤ + + + + . 2 4 2 2 2 2 2

⎛1 3βδt sign (γ ) γ 3δt ⎞ δt αδt ⎜ − ⎟ − − − ⎝2 ⎠ 2 2 2 2

{

L (Ω )

+ ‖∇. ξ k + 1‖22

L (Ω )

}

⎛1 3βδt sign (γ ) γ 3δt ⎞ δt αδt ⎟ ≤⎜ + + + + ⎝2 ⎠ 2 2 2 2

{ ‖ξ ‖

k 2 L2 (Ω)

L ∞ = ∥ uexact − uapprox ∥∞ = max { |uexact (x i, t ) − uapprox (x i, t )| : i = 1, 2, …, N},

Therefore we have

‖ξ k + 1‖22

In this section, we show the results obtained for three examples (taken from Refs. [36,37]) using the meshless methods described above. To show the accuracy and convergence of the method and also for comparing with the interpolating elementfree Galerkin method [36], two kinds of error measures, maximum absolute error L∞ and RMS error:

+ ‖∇. ξ k‖22

L (Ω )

}

⎛1 3βδt sign (γ ) γ 3δt ⎞ δt αδt ⎟ ‖Rk + 1‖22 . + δt ⎜ + + + + L (Ω ) ⎝2 ⎠ 2 2 2 2 According to norm of H1, we have

⎛1 3βδt sign (γ ) γ 3δt ⎞ k + 1 2 δt αδt ⎜ − ⎟ ‖ξ ‖ 1 − − − H (Ω ) ⎝2 ⎠ 2 2 2 2 ⎛1 ⎞ 3βδt sign (γ ) γ 3δt δt αδt ⎟ ‖ξ k‖2 1 ≤⎜ + + + + H (Ω ) ⎝2 ⎠ 2 2 2 2 ⎛1 3βδt sign (γ ) γ 3δt ⎞ δt αδt ⎟ ‖Rk + 1‖22 . + δt ⎜ + + + + L (Ω ) ⎝2 ⎠ 2 2 2 2 By dividing the coefficients and repeat inequality

RMS =

1 N

N

∑ |uexact (x i, t ) − uapprox (x i, t )| i=1

(61)

are used. In the examples, N is the number of total nodal points covering Ω is regularly distributed. Also in order to implement the SMRPI method, the radius of support domain (that is a disk) to construct basis functions is rs = 4.2 h, where h is the distance between the nodes in x or y direction. This size is significant enough to have sufficient number of nodes (nx ) and gives an appropriate basis functions. In Eq. (4) set m ¼ 15, i.e. the quartic basis functions as follows:

P T (x) = {1, x, y , x2 , xy , y2 , x3, x2y , xy2 , y 3 , x 4 , x3y , x2y2 , xy 3 , y 4 } , also here, we use the thin plate spline (TPS) as radial basis functions in Eq. (4). This RBF is defined as follows:

R (x) = r 2m ln (r )

m = 1, 2, 3, …,

where in examples let m ¼ 2, i.e. R (x) = r 4 ln (r ). We note that all tables are on [0, 1]2. In Fig. 1, four domains which are used in the examples, are shown with N ¼900, N ¼1225, N ¼632 and N ¼804 nodes for Ω1, Ω2, Ω3 and Ω4, respectively.

48

E. Shivanian, A. Jafarabadi / Engineering Analysis with Boundary Elements 72 (2016) 42–54

Ω1

1

Ω2

0.8

0.5

y

y

0.6 0

0.4 −0.5

0.2 0

0

0.2

0.4

0.6

0.8

1

−1

−0.5

0 x

x Ω

0.5

1

Ω

3

4

1

0.6 0.8

0.4

0.6

0

y

y

0.2

0.4

−0.2 −0.4

0.2

−0.6 −0.5

0 x

0.5

0

1

0

0.2

0.4

0.6

0.8

1

x

Fig. 1. Different considered domains for the problem (1)–(3).

Example 1. For the first test problem, consider nonlinear GBBMB equation as the following form ∂u (x, y, t ) ∂ 2u (x, y, t ) ⎞ ∂ 2u (x, y, t ) ∂ 2u (x, y, t ) ∂ ⎛ ∂ 2u (x, y, t ) − + − ⎜ ⎟− ∂t ∂t ⎝ ∂x 2 ∂y 2 ∂x 2 ∂y 2 ⎠ +

∂u (x, y, t ) ∂u (x, y, t ) ∂u (x, y, t ) ∂u (x, y, t ) + = u (x , y , t ) + u (x , y , t ) + f (x , y , t ) , ∂x ∂y ∂x ∂y

with initial condition

u (x, y , 0) = 0,

x, y ∈ Ω = [0, 1]2 ,

and boundary condition

u (x, y , t ) = t sin (x + y),

x, y ∈ ∂Ω, 0 < t ≤ T ,

moreover, set f (x, y, t ) as

f (x, y , t ) = (3 + 2t −

2t 2

cos (x + y)) sin (x + y) + 2t cos (x + y).

The exact solution of this example is

u (x, y , t ) = t sin (x + y). It is evident from Table 1 that errors does not increase much while the time is increasing, in other words the method is rather stable in according to Theorem 4.1. Table 2 shows L∞ and RMS errors for h = 1/8, h = 1/10, and h = 1/20, i.e. N ¼81, N ¼121 and N ¼441, respectively, with two values of δt = 1/100 and δt = 1/1000. The comparison of obtained numerical results between the SMRPI technique and those of [36] is reported in Table 3 at different h and δt . As it is seen, the results are better than those obtained by [36].

Table 1 The L∞ and RMS with the time T ¼1 with fixed N ¼441 and δt = 0.01 for Example 1. t

RMS

L∞

t¼ 0.000000 t¼ 0.100000

0.000000

0.000000

9.618478 × 10−11

9.196118 × 10−9

t¼ 0.200000

1.893497 × 10−10

1.853817 × 10−8

t¼ 0.300000

2.799549 × 10−10

2.801201 × 10−8

t¼ 0.400000

3.683949 × 10−10

3.760483 × 10−8

t¼ 0.500000

4.550066 × 10−10

4.730500 × 10−8

t¼ 0.600000

5.400775 × 10−10

5.710199 × 10−8

t¼ 0.700000

6.238520 × 10−10

6.698624 × 10−8

t¼ 0.800000

7.065381 × 10−10

7.694910 × 10−8

t¼ 0.900000

7.883119 × 10−10

8.698270 × 10−8

t¼ 1.000000

8.693228 × 10−10

9.707991 × 10−8

The computed numerical solutions obtained by the SMRPI method and the absolute errors of these computed solutions with δt = 0.01 and N ¼ 441, i.e. h = 1/20, at some time levels are shown in Fig. 2. Fig. 3 is dedicated to the approximate solution and absolute error using the SMRPI with δt = 1/100, N ¼1225 and T ¼1, on domain Ω2. Also, Fig. 4 shows the approximate solution and absolute error using the SMRPI with δt = 1/100, N ¼632, T ¼1, on domain Ω3. Example 2. Consider nonlinear GBBMB equation as the following form

E. Shivanian, A. Jafarabadi / Engineering Analysis with Boundary Elements 72 (2016) 42–54

49

Table 2 Errors for SMRPI method with T ¼ 1 for Example 1.

δt = 1/100

δt = 1/1000

h

L∞

RMS

L∞

RMS

1/8

7.856140 × 10−6

2.254799 × 10−7

7.903713 × 10−6

2.274534 × 10−7

1/10

2.724002 × 10−6

6.076007 × 10−8

2.957150 × 10−6

5.995971 × 10−8

1/20

9.707991 × 10−8

8.693228 × 10−10

3.952248 × 10−6

2.195382 × 10−8

Table 3 Comparison between errors of present method and method of [36] for Example 1. Method of[36]

SMRPI method

RMS

L∞

RMS

L∞

h = δt = 1/4

1.6729 × 10−3

4.0579 × 10−3

8.5742 × 10−6

1.5251 × 10−4

h = 1/8, δt = 1/64

5.4396 × 10−4

1.1254 × 10−3

2.2550 × 10−7

7.8568 × 10−6

h = 1/16, δt = 1/1024

3.2228 × 10−4

9.1254 × 10−4

2.7754 × 10−8

4.2856 × 10−6

−7

x 10

0.5 0 1

0.5 y

0

0

0.5

0.5 0 1

1

y −7

(u) L

u approx

1

y

0

0

0.5

1

0 1

0

0.5

1

0 1

0

u approx

0

−7

(u) y

0

0.5

x 10

L 0.5

0.5

1

Fig. 2. uapprox (x, y, t ) and L∞ profiles for different time level on [0,

1]2

1

x

0

0.5

1

x

time=4

5

0 1

0.5 y

x

0.5

2

y

2

0

time=3

x

4

x

4

time=4

0 1

0

−7

(u)

u approx

y

0

0.5

x 10

L 0.5

1

2

y

2

0.5

time=2

x

4

0

4

time=3

0 1

0

x 10

2

0.5

0.5

x

time=2

0 1

time=1

1

(u)

1

L

u approx

time=1

0

0

0.5 x

with N ¼ 441 and δt = 0.01 for Example 1 by SMRPI.

1

50

E. Shivanian, A. Jafarabadi / Engineering Analysis with Boundary Elements 72 (2016) 42–54 −3

2.5

2

0.6

1.5

0.5

−0.5

−1

−0.4

−1.5

−2

−0.8

−2.5 −1

0 x

0.4

−1.5

−0.6

−2

0.6

−0.5

−0.2

−1

0.8

0

y

0

0

1

1

0.2

0.5

1.2

1.5

0.4

1

y

2.5

0.8

2

x 10 1.4

0.2

−2.5 −1

1

0 x

0

1

Fig. 3. Graphs of approximate solution and absolute error using the SMRPI with δt = 1/100 , N ¼ 1225, T ¼ 1, on domain Ω2 for Example 1.

8

0.5

6 L (u)

1

0

u

approx

x 10

−6

4

−0.5

2

−1 0.5

0 0.5 1 0 y

1 0

0 −0.5

−1

y

x

0 −0.5

x

−1

Fig. 4. Graphs of approximate solution and absolute error using the SMRPI with δt = 1/100 , N ¼ 632, T ¼1, on domain Ω3 for Example 1.

∂u (x, y, t ) ∂ 2u (x, y, t ) ⎞ ∂ 2u (x, y, t ) ∂ 2u (x, y, t ) ∂ ⎛ ∂ 2u (x, y, t ) ⎜ ⎟− − + − ∂t ∂t ⎝ ∂x 2 ∂y 2 ∂x 2 ∂y 2 ⎠ ∂u (x, y, t ) ∂u (x, y, t ) ∂u (x, y, t ) + + = cos (u (x, y, t )) + cos (u (x, y, t )) ∂x ∂y ∂x ∂u (x, y, t ) + f (x, y, t ), ∂y

with initial condition

u (x, y , 0) = exp (2x + 2y),

x, y ∈ Ω = [0, 1]2 ,

boundary condition

u (x, y , t ) = exp (2x + 2y + 2t ),

x, y ∈ ∂Ω,

0 < t ≤ T,

and source function

Table 4 The L∞ and RMS with the time T ¼1 with fixed N ¼441 and δt = 0.01 for Example 2. t

RMS

L∞

t¼ 0.000000 t¼ 0.100000

0.0000

0.000000

9.425612 × 10−7

4.386349 × 10−5

t¼ 0.200000

2.002017 × 10−6

9.338881 × 10−5

t¼ 0.300000

3.213299 × 10−6

1.502943 × 10−4

t¼ 0.400000

4.616618 × 10−6

2.158037 × 10−4

t¼ 0.500000

6.260532 × 10−6

2.926027 × 10−4

t¼ 0.600000

8.205497 × 10−6

3.849974 × 10−4

t¼ 0.700000

1.052082 × 10−5

4.940170 × 10−4

t¼ 0.800000

1.329085 × 10−5

6.241214 × 10−4

t¼ 0.900000

1.660330 × 10−5

7.804222 × 10−4

t¼ 1.000000

2.063719 × 10−5

9.702541 × 10−4

f (x, y , t ) = − 2 exp (2x + 2y + 2t )(9 + 2 cos (exp (2x + 2y + 2t ))). Then the exact solution is

u (x, y , t ) = exp (2x + 2y + 2t ).

E. Shivanian, A. Jafarabadi / Engineering Analysis with Boundary Elements 72 (2016) 42–54

values of δt = 1/100 and δt = 1/1000. The comparison of obtained numerical results between the SMRPI technique and [37] is reported in Table 6 at different h and δt . From Table 6, it is realized that with the increasing number of nodes, errors is reduced significantly. Therefore compared with method [37] the present method, has the more accuracy. Graphs of approximate solution and absolute error using the present method with h = 1/40 and δt = 0.01 on Ω = [0, 1]2, for two time level, are displayed in Fig. 5. Fig. 6 shows approximate solution and absolute error using the SMRPI with δt = 1/100, N ¼804, T¼ 1, on domain Ω4.

Table 5 Errors L∞ and RMS obtained for SMRPI method with T ¼1 in Example 2.

δt = 1/100 h

δt = 1/1000

L∞

RMS

L∞

RMS

1/8

4.243384 × 10−2

1.006548 × 10−3

4.285819 × 10−2

9.880416 × 10−4

1/10

1.666383 × 10−2

3.036335 × 10−4

1.687995 × 10−2

2.903117 × 10−4

1/20

9.702541 × 10−4

2.063719 × 10−5

7.681819 × 10−4

5.954839 × 10−6

Table 6 Comparison between L∞ error of present method and method of [37] with T ¼1 for Example 2.

L∞ for method [37]

Example 3. For the third (last) test problem consider nonlinear GBBMB equation as the following form

L∞ for SMRPI method

∂u (x, y , t ) ∂ ⎛ ∂ 2u (x, y , t ) ∂ 2u (x, y , t ) ⎞ ∂ 2u (x, y , t ) − ⎜ + ⎟− 2 ∂x2 ∂t ∂t ⎝ ∂x ∂y2 ⎠

δt

h = 1/5 (c = 2.5)

h = 1/20 (c = 0.3)

h = 1/5

1/40

4.4301 × 10−1

5.3785 × 10−1

2.3572 × 10−1

4.5437 × 10−3

1/80

2.0445 ×

10−1

10−1

10−1

10−3

1/160

7.2946 × 10−2

1.4291 × 10−1

2.6186 × 10−1

8.1881 × 10−4

1/320

3.7628 × 10−2

7.2687 × 10−2

2.6237 × 10−1

7.7906 × 10−4

2.8259 ×

2.5846 ×

h = 1/20

1.3221 ×

51

∂ 2u (x, y , t ) ∂u (x, y , t ) ∂u (x, y , t ) + + ∂y2 ∂x ∂y ∂u (x, y , t ) ∂u (x, y , t ) = u (x, y , t ) + u (x, y , t ) + f (x, y , t ), ∂x ∂y −

with initial condition Obviously from Table 4 it is inspected that errors does not increase much while the time is increasing, in other words the method is relatively stable. Table 5 shows L∞ and RMS for h = 1/8, h = 1/10, and h = 1/20, i.e. N ¼81, N ¼ 121 and N ¼441, respectively, with two

u (x, y , 0) = 0,

x, y ∈ Ω = [0, 1]2 ,

and boundary condition

time=1

−3

x 10

500

0.5

L

(u)

u approx

1

0 1 0.5 0 0

y

0.5

0 1

1

0.5 y

x

time=2

0 0

0.5

1

x

time=2 0.1

(u)

4000 2000

0.05

L

u approx

time=1

0 1 0.5 y

0 0

0.5

0 1

1

0.5 y

x

Fig. 5. uapprox (x, y, t ) and L∞ profiles for different time level on [0,

1]2

0 0

0.5 x

with N ¼ 1681 and δt = 0.01 for Example 2 by SMRPI.

1

52

E. Shivanian, A. Jafarabadi / Engineering Analysis with Boundary Elements 72 (2016) 42–54

Fig. 6. Graphs of approximate solution and absolute error using the SMRPI with δt = 1/100 , N¼ 804, T¼ 1, on domain Ω4 for Example 2. Table 7 The L∞ and RMS with the time T ¼ 1 with fixed N ¼ 441, δt = 0.01 and r¼ 3 for Example 3. t

RMS

L∞

t = 0.000000 t = 0.100000

0.000000

0.000000

6.504434 × 10−10

1.096219 × 10−7

t = 0.200000

1.295228 × 10−9

2.189869 × 10−7

t = 0.300000

1.935343 × 10−9

3.281422 × 10−7

t = 0.400000

2.571670 × 10−9

4.371308 × 10−7

t = 0.500000

3.204998 × 10−9

5.459915 × 10−7

t = 0.600000

3.836030 × 10−9

6.547599 × 10−7

t = 0.700000

4.465398 × 10−9

7.634678 × 10−7

t = 0.800000

5.093668 × 10−9

8.721444 × 10−7

t = 0.900000

5.721348 × 10−9

9.808161 × 10−7

t = 1.000000

6.348896 × 10−9

1.089507 × 10−6

accordingly. From Table 7, it is clearly understood that errors does not increase much while the time is increasing, in other words the method is rather stable. Table 8 shows L∞ and RMS errors for h = 1/8, h = 1/10, and h = 1/20, i.e. N ¼ 81, N ¼121 and N ¼441, respectively, with two values of δt = 1/100 and δt = 1/1000. The comparison of obtained numerical results between the SMRPI technique and [36] with r ¼3 is reported in Table 9 at different h and δt . As it is seen, the results are better than those obtained by [36]. The computed numerical solutions obtained by the SMRPI method and the absolute errors of these computed solutions with δt = 0.01 and N ¼441, i.e. h = 1/20, at some time levels are shown in Fig. 7. Also, Fig. 8 is dedicated to the approximate solution and absolute error using the SMRPI with δt = 1/100, N = 632, T = 1, on domain Ω3.

6. Concluding remarks Table 8 Errors for SMRPI method with T ¼ 1 and r ¼3 for Example 3.

δt = 1/100

δt = 1/1000

h

L∞

RMS

1/8

5.263742 × 10−5

1.093370 × 10−6

5.263823 × 10−5

10−5

10−7

10−5

3.373302 × 10−7

1.089181 × 10−6

7.012235 × 10−9

1/10

2.240466 ×

1/20

1.089507 × 10−6

3.374145 ×

RMS

L∞

6.348896 × 10−9

2.240457 ×

u (x, y , t ) = t ( sech2 (x + y − r ) + sech2 (x + y + r ) ),

1.093599 × 10−6

x, y ∈ ∂Ω,

0 < t ≤ T. The exact solution of this example is given

u (x, y , t ) = t ( sech2 (x + y − r ) + sech2 (x + y + ) ), and more, f (x, y, t ) is obtained by setting them into the Eq. (1)

In this article, the spectral meshless radial point interpolation (SMRPI) technique has been applied to the nonlinear generalized Benjamin-Bona-Mahony-Burgers (GBBMB) in two-dimension with initial and Dirichlet-type boundary conditions. We have discretized the main equation using forward finite difference through Crank-Nicolson technique. For the spatial variable, it has been used shape functions which is constructed locally by the help of the combination of thin plate radial basis functions and complete set of monomials via interpolation technique, as a result the shape functions have Kronecker delta function property and are convenient to used for Dirichlet-type boundary conditions. Furthermore, they are free of shape parameter which is tedious and annoying step in other RBFs based technique. It has been proved that the time discrete scheme is relatively stable with respect to the time variable in H1 space and convergent with convergence order 6 (δt ). Also, a comparison study of the present method and recently applied interpolating element-free Galerkin technique has been given through applying on GBBMB

Table 9 Comparison between errors of present method and method of [36] with for Example 3. Method of [36] RMS

h = δt = 1/4

SMRPI method RMS

L∞

L∞

4.0900 × 10−4

1.5252 × 10−3

2.183913 × 10−5

4.177233 × 10−4

10−4

10−4

1.093372 × 10−6

5.263745 × 10−5

7.3766 × 10−4

2.419361 × 10−8

3.021340 × 10−6

h = 1/8, δt = 1/64

2.0866 ×

h = 1/16, δt = 1/1024

1.1328 × 10−4

6.5184 ×

E. Shivanian, A. Jafarabadi / Engineering Analysis with Boundary Elements 72 (2016) 42–54

time=1

53

−6

x 10

time=1

0 1

L

u

(u)

approx

0.5

0.5

0

y

1

0.5

0

2 1 0 1

0.5

0

y

x

0.5

0

1

x

−6

x 10

time=2

0 1

(u) L

0.5

u

approx

1

0.5

0

y

1

0.5

0

time=2

4 2 0 1

0.5

0

y

x

0.5

0

1

x

−6

x 10

time=3

(u)

1 0 1

L

u approx

2

0.5

0

y

1

0.5

0

time=3

4 2 0 1

0.5

0

y

x

0.5

0

1

x

−6

x 10

time=4

(u)

1 0 1

L

u approx

2

0.5

0

y

1

0.5

0

time=4

5

0 1

0.5

0

y

x

0.5

0

x

Fig. 7. uapprox (x, y, t ) and L∞ profiles for different time level on [0, 1]2 with N ¼ 441 and δt = 0.01 for Example 3 by SMRPI.

x 10 4

0.1

3

0.06

L (u)

approx

0.08

u

−6

0.04

2 1

0.02

0 0.5

0 0.5 1 0 y

1 0

0 −0.5

−1

x

y

0 −0.5

−1

x

Fig. 8. Graphs of approximate solution and absolute error using the SMRPI with δt = 1/100 , N ¼632, T ¼ 1, on domain Ω3 for Example 3.

1

54

E. Shivanian, A. Jafarabadi / Engineering Analysis with Boundary Elements 72 (2016) 42–54

equation. Numerical results confirmed error and convergence analysis as well as the efficiency of the SMRPI.

Acknowledgments The authors are very grateful to anonymous reviewers for carefully reading the paper and for their comments and suggestions which have improved the paper very much.

References [1] Liu G-R, Gu Y-T. An introduction to meshfree methods and their programming. Springer Science & Business Media, Berlin; 2005. [2] Atluri SN, Shen S. The meshless local petrov-galerkin (mlpg) method: a simple/& less-costly alternative to the finite element and boundary element methods. CMES: Comput Model Eng Sci 2002;3(1):11–52. [3] Podlubny I. Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications, 198. Academic Press, San Diego, California; 1998. [4] Belytschko T, Lu YY, Gu L. Element-free galerkin methods. Int J Numer Methods Eng 1994;37(2):229–56. [5] Belytschko T, Lu YY, Gu L. Element free Galerkin methods for static and dynamic fracture. Int J Solids Struct 1995;32:2547–70. [6] Duan Y, Tan Y. A meshless galerkin method for dirichlet problems using radial basis functions. J Comput Appl Math 2006;196(2):394–401. [7] Assari P, Adibi H, Dehghan M. A meshless discrete galerkin (MDG) method for the numerical solution of integral equations with logarithmic kernels. J Comput Appl Math 2014;267:160–81. [8] Dehghan RSM. A meshless local petrov-galerkin method for the time-dependent maxwell equations. J Comput Appl Math 2014;268:93–110. [9] Mirzaei D, Dehghan M. Meshless local petrov-galerkin (MLPG) approximation to the two dimensional sine-Gordon equation. J Comput Appl Math 2010;233(10):2737– 54. [10] Fili A, Naji A, Duan Y. Coupling three-field formulation and meshless mixed galerkin methods using radial basis functions. J Comput Appl Math 2010;234(8):2456–68. [11] Dai B, Cheng Y. An improved local boundary integral equation method for twodimensional potential problems. Int J Appl Mech 2010;2(2):421–36. [12] Bai F, Li D, Wang J, Cheng Y. An improved complex variable element-free Galerkin method for two-dimensional elasticity problems. Chin Phys B 2012;21(2):0202041–10. [13] Atluri S, Shen S. The meshless local petrov-galerkin (MLPG) method: a simple and less costly alternative to the finite element and boundary element methods. Comput Model Eng Sci 2002;3:11–51. [14] Belytschko T, Lu Y, Gu L, Tabbara M. Element-free galerkin methods for static and dynamic fracture. Int J Solids Struct 1995;32(17):2547–70. [15] Atluri SN, Shen S. The meshless local petrov-galerkin (mlpg) method: a simple/& less-costly alternative to the finite element and boundary element methods. CMES: Comput Model Eng Sci 2002;3(1):11–52. [16] Shivanian E. Meshless local petrov-galerkin (MLPG)method for three-dimensional nonlinear wave equations via moving least squares approximation. Eng Anal Bound Elem 2015;50:249–57. [17] Shivanian E, Abbasbandy S, Alhuthali MS, Alsulami HH. Local integration of 2-d fractional telegraph equation via moving least squares approximation. Eng Anal Bound Elem 2015;56(0):98–105. [18] Kansa EJ. Multiquadricsa scattered data approximation scheme with applications to computational fluid-dynamicsii solutions to parabolic, hyperbolic and elliptic partial differential equations. Comput Math Appl 1990;19(8):147–61. [19] Lin J, Chen W, Sze K. A new radial basis function for helmholtz problems. Eng Anal Bound Elem 2012;36(12):1923–30.

[20] Abbasbandy S, Ghehsareh HR, Hashim I. A meshfree method for the solution of twodimensional cubic nonlinear schrödinger equation. Eng Anal Bound Elem 2013;37 (6):885–98. [21] Aslefallah M, Shivanian E. Nonlinear fractional integro-differential reaction-diffusion equation via radial basis functions. Eur Phys J 2015;130(47):1–9. [22] Abbasbandy S, Ghehsareh HR, Alhuthali MS, Alsulami HH. Comparison of meshless local weak and strong forms based on particular solutions for a non-classical 2-d diffusion model. Eng Anal Bound Elem 2014;39:121–8. [23] Shirzadi A, Sladek V, Sladek J. A local integral equation formulation to solve coupled nonlinear reaction-diffusion equations by using moving least square approximation. Eng Anal Bound Elem 2013;37(1):8–14. [24] Shivanian E. Analysis of meshless local radial point interpolation (MLRPI) on a nonlinear partial integro-differential equation arising in population dynamics. Eng Anal Bound Elem 2013;37:1693–702. [25] Shivanian E, Khodabandehlo H. Meshless local radial point interpolation (MLRPI) on the telegraph equation with purely integral conditions. Eur Phys J 2014;129:241–51. [26] Hosseini VR, Shivanian E, Chen W. Local integration of 2-D fractional telegraph equation via local radial point interpolant approximation. Eur Phys J 2015;130 (33):1–21. [27] E. Shivanian, On the convergence analysis, stability, and implementation of meshless local radial point interpolation on a class of three‐dimensional wave equations, International Journal for Numerical Methods in Engineering 105 (2) (2016) 83–110. [28] Shivanian E. Analysis of meshless local and spectral meshless radial point interpolation (mlrpi and smrpi) on 3-d nonlinear wave equations. Ocean Eng 2014;89:173–88. [29] Shivanian E. A new spectral meshless radial point interpolation (smrpi) method: a well-behaved alternative to the meshless weak forms. Eng Anal Bound Elem 2015;54:1–12. [30] Fatahi H, Saberi-Nadjafi J, Shivanian E. A new spectral meshless radial point interpolation (smrpi) method for the two-dimensional fredholm integral equations on general domains with error analysis. Journal of Computational and Applied Mathematics. 〈http://dx.doi.org/10.1016/j.cam.2015.08.018〉. [31] Peregrine D. Calculations of the development of an undular bore. J Fluid Mech 1966;25(02):321–30. [32] Benjamin TB, Bona JL, Mahony JJ. Model equations for long waves in nonlinear dispersive systems. Philos Trans R Soc Lond A: Math Phys Eng Sci 1972;272 (1220):47–78. [33] Bona JL, Pritchard W, Scott LR. Solitary-wave interaction. Physics of Fluids (1958– 1988), vol. 23(3), 1980 p. 438–41. [34] Achouri T, Khiari N, Omrani K. On the convergence of difference schemes for the benjamin-bona-mahony (bbm) equation. Appl Math Comput 2006;182(2):999– 1005. [35] Omrani K. Optimal l − ∞ error estimates for finite element galerkin methods for nonlinear evolution equations. J Appl Math Comput 2008;26(1–2):247–62. [36] Dehghan M, Abbaszadeh M, Mohebbi A. The use of interpolating element-free galerkin technique for solving 2d generalized benjamin-bona-mahony-burgers and regularized long-wave equations on non-rectangular domains with error estimate. J Comput Appl Math 2015;286:211–31. [37] Dehghan M, Abbaszadeh M, Mohebbi A. The numerical solution of nonlinear high dimensional generalized benjamin-bona-mahony-burgers equation via the meshless method of radial basis functions. Comput Math Appl 2014;68(3):212–37. [38] Liu G, Gu Y. A local radial point interpolation method (lrpim) for free vibration analyses of 2-d solids. J Sound Vib 2001;246(1):29–46. [39] Wang J, Liu G. A point interpolation meshless method based on radial basis functions. Int J Numer Methods Eng 2002;54(11):1623–48. [40] Wang J, Liu G. On the optimal shape parameters of radial basis functions used for 2-d meshless methods. Comput Methods Appl Mech Eng 2002;191(23):2611–30. [41] Powell MJ. The theory of radial basis function approximation in 1990. Univ Camb Dep Appl Math Theor Phys 1990. [42] Wendland H. Error estimates for interpolation by compactly supported radial basis functions of minimal degree. J Approx Theory 1998;93(2):258–72. [43] Dehghan M, Ghesmati A. Numerical simulation of two-dimensional sine-gordon solitons via a local weak meshless technique based on the radial point interpolation method (rpim). Comput Phys Commun 2010;181(4):772–86. [44] Kreyszig E. Advanced engineering mathematics. John Wiley & Sons, New York; 1988.