Analysis of meshless local and spectral meshless radial point interpolation (MLRPI and SMRPI) on 3-D nonlinear wave equations

Ocean Engineering 89 (2014) 173–188

Contents lists available at ScienceDirect

Ocean Engineering journal homepage: www.elsevier.com/locate/oceaneng

Analysis of meshless local and spectral meshless radial point interpolation (MLRPI and SMRPI) on 3-D nonlinear wave equations Elyas Shivanian n 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 6 November 2013 Accepted 8 August 2014 Available online 28 August 2014

In this paper, the meshless local radial point interpolation (MLRPI) method is applied to simulate three space dimensional nonlinear wave equation of the form utt þ αut þ βu ¼ uxx þ uyy þ uzz þ δgðuÞut þ f ðx; y; z; tÞ, 0 o x; y; z o 1; t 40 subject to given appropriate initial and Dirichlet boundary conditions and then a new spectral meshless radial point interpolation (SMRPI) method is proposed to solve the mentioned problem. The main drawback of methods in fully 3-D problems is the large computational costs. In the MLRPI method, all integrations are carried out locally over small quadrature domains of regular shapes such as spheres or cubes. In innovative SMRPI method, it is not needed to any integration. The radial point interpolation method is proposed to construct shape functions for MLRPI and basis functions for SMRPI. A weak formulation with a Heaviside step function transforms the set of governing equations into local integral equations on local subdomains in MLRPI whereas the operational matrices converts easily the governing equations (even high order) into linear system of equations in SMRPI. A two-step time discretization method is employed to approximate the time derivatives. To treat the nonlinearity part, a kind of predictor–corrector scheme combined with one-step time discretization and Crank–Nicolson technique is adopted. A comparison study of the efficiency and accuracy of the MLRPI and SMRPI method is given by applying on mentioned problem. Convergence studies in the numerical examples show that the SMRPI method possesses excellent rates of convergence. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Spectral meshless radial point interpolation (SMRPI) Local weak formulation MLRPI Radial basis function 3-D wave equation

1. Introduction The wave equations model vibrations of structures such as buildings, beams and machines, and more are the basis for fundamental equations of atomic physics. The wave equation usually describes water waves, the vibrations of a string, the propagation of sound waves, and the transmission of electric signals in a cable, etc. The typical model of wave equation is an initial boundary value problem being valid in a bounded domain or an initial value problem being valid in an unbounded domain. It is interesting to mention that two initial conditions should be prescribed, namely the initial displacement uðx; 0Þ ¼ φðxÞ and the initial velocity ut ðx; 0Þ ¼ ψ ðxÞ which describe the initial displacement and the initial velocity at the starting time t¼0, respectively. Unlike the heat equation, the wave equation includes the term utt which represents the vertical acceleration at point x. The wave equation plays a significant role in various physical problems and is required in diverse areas of science and engineering (Selvadurai, 2000; Nettel, 2003; Hirose, 1985; Wazwaz, 2002; Dai et al., 2006; Grilli et al., 2001; Yan and Konotop, 2009; Senturk, 2011).

n

Tel./fax: þ98 9126825371. E-mail address: [email protected]

http://dx.doi.org/10.1016/j.oceaneng.2014.08.007 0029-8018/& 2014 Elsevier Ltd. All rights reserved.

The numerical solution of nonlinear wave equation with high order of accuracy in space and time is still a challenging task, especially if 3-dimensional computational domains are involved. The popular technique to solve wave equations is to discretize time and spatial derivatives by second-order finite difference schemes (Dehghan, 2006a, 2005, 2006b). Kamruzzaman has applied a collocation method to hyperbolic equations (Kamruzzaman, 2008), and Xiong has used the finite element method (FEM) for semilinear hyperbolic equations (Xiong, 2008). Dehghan et al. have applied the combination of collocation, finite difference, and multigrid methods as well as radial basis functions to some wave (Dehghan and Shokri, 2009; Dehghan and Mohebbi, 2008). The numerical methods for hyperbolic problems are mainly based on difference methods, or the FEM (Avilez-Valente and Seabra-Santos, 2004; Gamallo and Astley, 2006). However, the meshless methods have been seldom applied in this field especially in three dimensional waves (Carlos et al., 2005; Lu et al., 1995; Bouillard et al., 2004; Dehghan and Salehi, 2012a, 2012b; Dehghan and Ghesmati, 2010a). Consider the following three-space dimensional nonlinear hyperbolic partial differential equation: ∂2 u ∂u þ α þ βu ¼ Δu þ δgðuÞut þ f ðx; tÞ; ∂t ∂t 2

x A Ω; t A ½0; T;

ð1Þ

174

E. Shivanian / Ocean Engineering 89 (2014) 173–188

subject to the initial conditions uðx; 0Þ ¼ φðxÞ;

∂u ðx; 0Þ ¼ ψ ðxÞ; x A Ω; ∂t

ð2Þ

and the boundary conditions uðx; tÞ ¼ a0 ðy; z; tÞ

for x ¼ 0;

uðx; tÞ ¼ a1 ðy; z; tÞ uðx; tÞ ¼ b0 ðx; z; tÞ

for x ¼ 1; x A Ω; t 4 0;

ð3Þ

for y ¼ 0;

uðx; tÞ ¼ b1 ðx; z; tÞ

for y ¼ 1; x A Ω; t 4 0;

uðx; tÞ ¼ c0 ðx; y; tÞ

for z ¼ 0;

uðx; tÞ ¼ c1 ðx; y; tÞ

for z ¼ 1; x A Ω; t 4 0;

ð4Þ

ð5Þ

where xT ¼ ½x; y; z is spatial variable, Ω ¼ ½0; 13 ¼ fðx; y; zÞ : 0 r x; y; z r1g and α; β A R are constant. Mohanty and Gopal (2013) have 2 2 2 4 applied a numerical method of order Oðk þ k h þ h Þ based on off-step discretization for the solution of the above equation. Liew and Cheng (2013) have considered the special case of the problem (1)–(5) when it is linear and applied a technique based on the mesh-free kp-Ritz method. Also, some other special cases of the mentioned problem have been considered by some researchers (Wen, 2010; Zhang et al., 2012; Han and Atluri, 2004). The main shortcoming of mesh-based methods such as the finite element method (FEM), the finite volume method (FVM) and the boundary element method (BEM) is that these numerical methods rely on meshes or elements. In the two last decades, in order to overcome the mentioned difficulties some techniques so-called meshless methods have been proposed (Liu and Gu, 2005). This method is used to establish system of algebraic equations for the whole domain of the problem without the use of predefined mesh for the domain discretization so that it is used a set of nodes scattered within the domain of the problem as well as sets of nodes on the boundaries of the domain to represent (but not to discretize) the domain of the problem and its boundaries. These sets of scattered nodes are usually called field nodes. There are three types of meshless methods: Meshless methods based on weak forms such as the element free Galerkin (EFG) method (Belytschko et al., 1994, 1995), meshless methods based on collocation techniques (strong forms) such as the meshless collocation method based on radial basis functions (RBFs) (Kansa, 1990; Dehghan and Shokri, 2008) and meshless methods based on the combination of weak forms and collocation technique. Due to the ill-conditioning of the resultant linear systems in the RBFcollocation method, various approaches are proposed to circumvent this problem, Ling et al. (2006), Ling and Schaback (2008), Libre et al. (2008), Libre et al. (2009) being among them. The weak forms are used to derive a set of algebraic equations through a numerical integration process using a set of quadrature domain that may be constructed globally or locally in the domain of the problem. In the global weak form methods, global background cells are needed for numerical integration in computing the algebraic equations. To avoid the use of global background cells, a so-called local weak form is used to develop the meshless local Petrov–Galerkin (MLPG) method (Atluri and Zhu, 1998a, 1998b, 2000a, 2000b; Dehghan and Mirzaei, 2008, 2009; Gu and Liu, 2001; Abbasbandy and Shirzadi, 2010, 2011; Shirzadi et al., 2012, 2013). When a local weak form is used for a field node, the numerical integrations are carried out over a local quadrature domain defined for the node, which also be the local domain where the test (weight) function is defined. The local domain usually has a regular and simple shape for an internal node such as sphere and cube and the integration is done numerically within the local domain. Hence the domain and boundary integrals in the weak form methods can easily be evaluated over the regularly shaped sub-domains and their boundaries.

In the literature, several meshless weak form methods have been reported such as: diffuse element method (DEM) (Nayroles et al., 1992), smooth particle hydrodynamic (SPH) (Bratsos, 2008; Clear, 1998), reproducing kernel particle method (RKPM) (Liu et al., 1995), boundary node method (BNM) (Mukherjee and Mukherjee, 1997), partition of unity finite element method (PUFEM) (Melenk and Babuska, 1996), finite sphere method (FSM) (De and Bathe, 2000), boundary point interpolation method (BPIM) (Gu and Liu, 2002) and boundary radial point interpolation method (BRPIM) (Gu and Liu, 2003). G.R. Liu applied the concept of MLPG and developed meshless local radial point interpolation (MLRPI) method (Liu et al., 2002; Liu and Gu, 2001; Dehghan and Ghesmati, 2010b). In this paper, we concentrate on the numerical solution of Eqs. (1)–(5) using the meshless local radial point interpolation (MLRPI) method, and then when observing the results are not satisfactory, we propose an innovative spectral meshless radial point interpolation (SMRPI) method. In the spectral meshless radial point interpolation (SMRPI) technique, the point interpolation method with the help of radial basis functions is proposed to construct basis functions which have Kronecker delta function property. Based on spectral methods, evaluation of high order derivatives of given differential equation is not difficult by constructing and using operational matrices. In the SMRPI method, it does not require any kind of integration locally over small quadrature domains. Therefore, computational costs of the SMRPI method are less expensive. Furthermore, four illustrative examples namely 3-dimensional telegraphic equation, Van der Pol type nonlinear wave equation, dissipative nonlinear wave equation and an example with the exponential growth respect to the time are presented to compare MLRPI and SMRPI. The results show that SMRPI is superior to MLRPI.

2. The modified radial point interpolation scheme In the conventional point interpolation method (PIM) there is a main difficulty that inverse of the polynomial moment matrix does not often exist. This condition could always be possible depending on the locations of the nodes in the support domain and the terms of monomials used in the basis. If an inappropriate polynomial basis is chosen for a given set of nodes, it may yield in a badly conditioned or even singular moment matrix (Liu and Gu, 2005). In order to avoid the singularity problem in the polynomial point interpolation method (PIM), radial basis function (RBF) is used to develop the radial point interpolation method (RPIM) for meshless weak form techniques (Liu and Gu, 2001; Wang and Liu, 2002a, 2002b). The combination of RPIM and polynomial PIM is described as following: Consider a continuous function uðxÞ defined in a domain Ω, which is represented by a set of field nodes. The uðxÞ at a point of interest x is approximated in the form of n

m

i¼1

j¼1

uðxÞ ¼ ∑ Ri ðxÞai þ ∑ pj ðxÞbj ¼ RT ðxÞa þ PT ðxÞb;

ð6Þ

where Ri ðxÞ is a radial basis function (RBF), n is the number of RBFs, pj ðxÞ is monomial in the space coordinate xT ¼ ½x; y; z, and m is the number of polynomial basis functions. The pj ðxÞ in Eq. (6) is built using Pascal's triangle and a complete basis is usually preferred. The linear basis functions are given by PT ðxÞ ¼ f1; x; y; zg;

m ¼ 4;

ð7Þ

the quadratic basis functions are PT ðxÞ ¼ f1; x; y; z; x2 ; y2 ; z2 ; xy; xz; yzg;

m ¼ 10;

ð8Þ

E. Shivanian / Ocean Engineering 89 (2014) 173–188

and the cubic basis functions are presented as T

2

2

2

3

3

3

2

P ðxÞ ¼ f1; x; y; z; x ; y ; z ; xy; xz; yz; x ; y ; z ; x y; x2 z; xy2 ; y2 z; xz2 ; yz2 ; xyzg; m ¼ 20:

ð9Þ

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. There are a number of types of RBFs, and the characteristics of RBFs have been widely investigated (Kansa, 1990; Franke and Schaback, 1997; Sharan et al., 1997). In the current work, we have chosen the thin plate spline (TPS) as radial basis functions in Eq. (6). This RBF is defined as follows: RðxÞ ¼ r 2m lnðrÞ;

m ¼ 1; 2; 3; …

ð10Þ

Since RðxÞ in Eq. (10) belongs to C 2m  1 (all continuous function to the order 2m  1), so higher-order thin plate splines must be used for higher-order partial differential operators. For the secondorder partial differential equation (1), m ¼ 2 is used for thin plate splines (i.e. second-order thin plate splines). In the radial basis function Ri ðxÞ, the variable is only the distance between the point of interest x and a node at qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi xi ¼ ðxi ; yi ; zi Þ, i.e. r ¼ ðx  xi Þ2 þðy  yi Þ2 þ ðz zi Þ2 . In order to determine ai and bj in Eq. (6), a support domain is formed for the point of interest at x, and n field nodes are included in the support domain (see Fig. 1). Coefficients ai and bj in Eq. (6) can be determined by enforcing Eq. (6) to be satisfied at these n nodes surrounding the point of interest x. This leads to a linear system of n equations, one for each node. The matrix form of this system of linear equations can be expressed as Us ¼ R n a þ Pm b

ð11Þ

where the vector of function values Us is Us ¼ fu1 u2 u3 … un gT ;

ð12Þ

the RBFs moment matrix is 0 1 R1 ðr 1 Þ R2 ðr 1 Þ … Rn ðr 1 Þ B C B R1 ðr 2 Þ R2 ðr 2 Þ … Rn ðr 2 Þ C C Rn ¼ B B ⋮ ⋮ ⋱ ⋮ C @ A R1 ðr n Þ R2 ðr n Þ … Rn ðr n Þ

175

and the polynomial moment matrix is 0 1 1 x1 y1 … pm ðx1 Þ B C B 1 x2 y2 … pm ðx2 Þ C C Pm ¼ B : B⋮ ⋮ ⋮ ⋱ ⋮ C @ A 1 xn yn … pm ðxn Þ

ð14Þ

nm

Also, the vector of unknown coefficients for RBFs is aT ¼ fa1 a2 … an g

ð15Þ

and the vector of unknown coefficients for polynomial is T

b ¼ fb1 b2 …bm g:

ð16Þ

We notify that in Eq. (13), rk in Ri ðr k Þ is defined as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r k ¼ ðxk  xi Þ2 þ ðyk  yi Þ2 þðzk  zi Þ2 :

ð17Þ

we mention that there are m þ n variables in Eq. (11). The additional m equations can be added using the following m constraint conditions: n

∑ pj ðxi Þai ¼ PTm a ¼ 0;

i¼1

j ¼ 1; 2; …; m:

ð18Þ

Combining Eqs. (11) and (18) yields the following system of equations in the matrix form !   R n Pm  a  ~ s ¼ Us ¼ ~ U ¼ Ga; ð19Þ PTm 0 0 b where ~ s ¼ fu1 u2 … un 0 0 … 0g; U

T a~ ¼ fa1 a2 … an b1 b2 … bm g:

ð20Þ Because the matrix Rn is symmetric, the matrix G will also be symmetric. Solving Eq. (19), we obtain   a ~ s: ð21Þ a~ ¼ ¼ G  1U b Eq. (6) can be re-written as uðxÞ ¼ RT ðxÞa þ PT ðxÞb ¼ fR T ðxÞ PT ðxÞg

; nn

ð13Þ

  a : b

ð22Þ

Now using Eq. (21) we obtain: ~ T ðxÞU ~ s; ~s¼Φ uðxÞ ¼ fR T ðxÞ PT ðxÞgG  1 U

Fig. 1. Local support and local quadrature domains for an arbitrary nodal point xi ¼ ðx; y; zÞ.

ð23Þ

176

E. Shivanian / Ocean Engineering 89 (2014) 173–188

where u~ is the latest available approximation of u. Using the above discussion, Eq. (1) can be written as

~ T ðxÞ can be rewritten as where Φ ~ ðxÞ ¼ fR T ðxÞ PT ðxÞgG  1 Φ ¼ fϕ1 ðxÞ ϕ2 ðxÞ … ϕn ðxÞ ϕn þ 1 ðxÞ … ϕn þ m ðxÞg: T

ð24Þ

1 ðΔtÞ2

ðuk þ 1 ðxÞ  2uk ðxÞ þ uk  1 ðxÞÞ

α

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 Φ ðxÞ so that it is

þ ðuk þ 1 ðxÞ þ uk ðxÞ þ uk  1 ðxÞÞ 3 1 ¼ ðΔuk þ 1 ðxÞ þ Δuk ðxÞ þ Δuk  1 ðxÞÞ 3 " # δ u~ k þ 1  u~ k u~ k þ 1  u~ k  1 u~ k  u~ k  1 þ gðu~ k Þ þ gðu~ k  1 Þ þ gðu~ k þ 1 Þ 3 Δt 2 Δt Δt

ð25Þ

then Eq. (23) is converted to the following one n

uðxÞ ¼ Φ ðxÞUs ¼ ∑ ϕi ðxÞui : T

ð26Þ

i¼1

þ

The derivatives of uðxÞ are easily obtained as

 1 f ðx; ðk þ 1ÞΔtÞ þ f ðx; kΔtÞ þ f ðx; ðk  1ÞΔtÞ : 3

ð34Þ

Suppose that λ ¼ 3=ðΔtÞ2 , η ¼ 3α=2Δt and Fðx; kÞ ¼ f ðx; ðk þ1Þ ΔtÞ þ f ðx; kΔtÞ þ f ðx; ðk  1ÞΔtÞ, then we obtain

n ∂ϕ ðxÞ ∂uðxÞ i ¼ ∑ ui ; ∂x ∂x i¼1

ðλ þ η þ β Þuk þ 1  Δuk þ 1

n ∂ϕ ðxÞ ∂uðxÞ i ¼ ∑ u; ∂y ∂y i i¼1 n ∂ϕ ðxÞ ∂uðxÞ i ¼ ∑ ui : ∂z ∂z i¼1

ðuk þ 1 ðxÞ uk  1 ðxÞÞ

β

T

ΦT ðxÞ ¼ fϕ1 ðxÞ ϕ2 ðxÞ … ϕn ðxÞg;

2 Δt

¼ ð2λ  βÞuk þ Δuk þðη  λ  βÞuk  1 þ Δuk  1

Note that R n 1 usually exists for arbitrary scattered nodes (Powell, 1992; Wendland, 1998). In addition, the order of polynomial used in Eq. (6) is relatively low (m ¼4 or 10). It is remarkable that the RPIM shape functions have the Kronecker delta function property, that is ( 1; i ¼ j; j ¼ 1; 2; …; n; ϕi ðxj Þ ¼ ð28Þ 0; i aj; i; j ¼ 1; 2; …; n: This is because the RPIM shape functions are created to pass thorough nodal values.

3. The time discretization approximation In the current work, we employ a time-stepping scheme to overcome the time derivatives. For this purpose, the following finite difference approximations for the time derivative operators are used: ∂2 uðx; tÞ 1 ffi ðuk þ 1 ðxÞ  2uk ðxÞ þ uk  1 ðxÞÞ; ∂t 2 ðΔtÞ2

ð29Þ

∂uðx; tÞ 1 ffi ðuk þ 1 ðxÞ uk  1 ðxÞÞ: ∂t 2 Δt

ð30Þ

Also we employ the following approximations by using the Crank– Nicolson technique: uðx; tÞ ffi 13 ðuk þ 1 ðxÞ þ uk ðxÞ þ uk  1 ðxÞÞ;

ð31Þ

Δuðx; tÞ ffi 13 ðΔuk þ 1 ðxÞ þ Δuk ðxÞ þ Δuk  1 ðxÞÞ;

ð32Þ

where uk ðxÞ ¼ uðx; kΔtÞ and x ¼ ðx; y; zÞ. To eliminate the nonlinearity (the term gðuÞut ), an iterative approach in this case a kind of predictor–corrector scheme combined with one-step time discretization and Crank–Nicolson technique (Dehghan and Ghesmati, 2010b) is performed as follows: " # u~ k þ 1  u~ k u~ k þ 1  u~ k  1 u~ k  u~ k  1 1 gðuÞut ¼ gðu~ k þ 1 Þ þ gðu~ k Þ þ gðu~ k  1 Þ ; 3 Δt 2 Δt Δt ð33Þ

δ

½2gðu~ k þ 1 Þðu~ k þ 1  u~ k Þ þ gðu~ k Þðu~ k þ 1  u~ k  1 Þ 2 Δt þ2gðu~ k  1 Þðu~ k  u~ k  1 Þ þ Fðx; kÞ þ

ð27Þ

ð35Þ

Remark. We notice that derivative of unknown function uðx; tÞ at t¼0 is known in the problem (1)–(5). Therefore, we are encouraged to use the central-difference formulas (29) and (30) for time derivatives which have a leading error of order h2. Hence, having order of convergence two in time dimension we are led to have a two-step time discretization approximation method. 4. The meshless local weak form formulation for MLRPI Let us now before continuing recall a certain Green's formula in 3-dimension which will be of fundamental importance in what follows. Let us start from the divergence theorem (in three dimensions): Z Z div A dx ¼ A  n ds; ð36Þ Ω

Γ

where A ¼ ðA1 ; A2 ; A3 Þ is a vector-valued function defined on x ¼ ðx; y; zÞ, div A ¼

Ω,

∂A1 ∂A2 ∂A3 þ þ ∂x ∂y ∂z

and n ¼ ðn1 ; n2 ; n3 Þ is the outward unit normal to Γ . Here dx denotes the element of volume in R3 and ds the element of area on Γ . If we apply the divergence theorem to A ¼ ðvw; 0; 0Þ, A ¼ ð0; vw; 0Þ and A ¼ ð0; 0; vwÞ we find that Z Z Z ∂v ∂w w dx þ v dx ¼ vwn1 ds; ð37Þ Ω ∂x Ω ∂x Γ Z Z ∂v ∂w w dx þ v dx ¼ vwn2 ds; Ω ∂y Ω ∂y Γ

ð38Þ

Z Z ∂v ∂w w dx þ v dx ¼ vwn3 ds; Ω ∂z Ω ∂z Γ

ð39Þ

Z

Z

Denoting by ∇v the gradient of v, i.e. ∇v ¼ ð∂v=∂x; ∂v=∂y; ∂v=∂zÞ we get from (37) to (39) the following 3-dimensional Green's formula:  Z  Z ∂v ∂w ∂v ∂w ∂v ∂w þ þ dx ∇v  ∇w dx ¼ Ω Ω ∂x ∂x ∂y ∂y ∂z ∂z

E. Shivanian / Ocean Engineering 89 (2014) 173–188

Z 

 ∂w ∂w ∂w n1 þ v n2 þ v n3 ds ∂x ∂y ∂z Γ  Z  2 ∂ w ∂2 w ∂2 w þ þ dx  v ∂x2 ∂y2 ∂z2 Z Z Ω ∂w ¼ v ds  vΔw dx; Γ ∂n Ω ¼

v

ð40Þ

where ∂w ∂w ∂w ∂w ¼ n1 þ n2 þ n3 ∂n ∂x ∂y ∂z is the normal derivative, i.e. the derivative in the outward normal direction to the boundary Γ . Instead of giving the global weak form, the MLRPI method constructs the weak form over local quadrature cell such as Ωq, which is a small region taken for each node in the global domain Ω. The local quadrature cells over lap each other and cover the whole global domain Ω. The local quadrature cells could be of any geometric shape and size. In this paper they are taken to be of spherical shape (see Fig. 1). The local weak form of Eq. (35) for i xi ¼ ðxi ; yi ; zi Þ A Ωq can be written as Z ððλ þ η þ βÞuk þ 1  Δuk þ 1 ÞvðxÞ dΩ Ωiq

Z

¼

Ωiq

ðð2λ  β Þuk þ Δuk ÞvðxÞ dΩ

Z

þ Z þ Z þ

Ωiq Ω

i q

ððη  λ  βÞu

k1

þ Δu

k1

δ

ÞvðxÞ dΩ

½2gðu~ k þ 1 Þðu~ k þ 1  u~ k Þ þgðu~ k Þðu~ k þ 1  u~ k  1 Þ

þ 2gðu~ k  1 Þðu~ k  u~ k  1 ÞvðxÞ dΩ

ð41Þ

where Ω is the local quadrature domain associated with the point i, i.e. it is a sphere centered at xi of radius rq and, vðxÞ is the Heaviside step function (Hu et al., 2006; Liu et al., 2006), ( 1; x A Ωq ; ð42Þ vðxÞ ¼ 0; x 2 = Ωq ; i q

as the test function in each local quadrature domain. Using the 3-dimensional Green's formula (40), Eq. (41) yields the following expression: Z Z ðλ þ η þ β Þ uðk þ 1Þ vðxÞ dΩ þ ∇uðk þ 1Þ ∇v dΩ Z  Z 

Ωiq

Ωiq

ðk þ 1Þ

∂u v i ∂n ∂Ωq Ω

i q

dΓ ¼ ð2λ  β Þ

þ ðη  λ  β Þ Z

Ωiq

Ωiq

uðkÞ vðxÞ dΩ

i q

v

ðk  1Þ

i

ð43Þ

where ∂Ωq is the boundary of Ωiq, n ¼ ðn1 ; n2 ; n3 Þ is the outward i unit normal to the boundary ∂Ωq , and i

∂u ∂u ∂u ∂u ¼ n1 þ n2 þ n3 ∂n ∂x ∂y ∂z

In this section, we consider Eq. (44) to see how to obtain discrete equations. Consider N regularly located points on the boundary and domain of the problem so that the distance between consecutive nodes in each direction is constant and equal to h (see Fig. 1). Assuming that uðxi ; ðk  1ÞΔtÞ and uðxi ; kΔtÞ; i ¼ 1; 2; …; N are known, our aim is to compute uðxi ; ðk þ 1ÞΔtÞ; i ¼ 1; 2; …; N. So, we have N unknowns and to compute these unknowns we need N equations. As it will be described, corresponding to each node we obtain one equation. For nodes which are located in the interior of the domain, i.e. for xi A interior Ω, to obtain the discrete equations from the locally weak forms (44), substituting approximation formulas (26) and (27) into local integral equations (44) yields ! ! Z Z N N ∂ϕj ðk þ 1Þ þ 1Þ d ðλ þ η þ β Þ ∑ ϕ d Ω  ∑ Γ u uðk j j j i Ωiq ∂Ωq ∂n j¼1 j¼1 ! ! Z Z N N ∂ ϕj ðkÞ dΓ uðkÞ ¼ ð2λ  βÞ ∑ ϕj dΩ uj þ ∑ j i Ωiq ∂Ωq ∂n j¼1 j¼1 ! ! Z Z N N ∂ϕj ðk  1Þ  1Þ d þðη  λ  βÞ ∑ ϕ d Ω þ ∑ Γ u uðk j j j i Ωiq ∂Ωq ∂n j¼1 j¼1

Ωiq

Fðx; kÞ dΩ þ γ ½2gðu~ ki þ 1 Þðu~ ki þ 1  u~ ki Þ

or equivalently " Z N ðλ þ η þ β Þ ∑

∂Ωq

Ωiq

ð44Þ

where γ ¼ 2πδr 3q =3Δt. Applying the radial point interpolation (RPI) for the unknown functions, the local integral equation (44) is transformed into a system of algebraic equations with used unknown quantities, as described in the next section.

þgðu~ ki Þðu~ ki þ 1  u~ ki  1 Þ þ 2gðu~ ki  1 Þðu~ ki  u~ ki  1 Þ;

Ωiq

þ gðu~ k Þðu~ k þ 1  u~ k  1 Þ þ 2gðu~ k  1 Þðu~ k  u~ k  1 ÞvðxÞ dΩ Z þ Fðx; kÞvðxÞ dΩ;

þgðu~ ki Þðu~ ki þ 1  u~ ki  1 Þ þ 2gðu~ ki  1 Þðu~ ki  u~ ki  1 Þ;

þ

∂u dΓ v ∂n Z δ þ ½2gðu~ k þ 1 Þðu~ k þ 1  u~ k Þ 2Δt Ωiq þ

Ωiq

Z

∂uðkÞ dΓ ∂n ∂Ω Z uðk  1Þ vðxÞ dΩ  ∇uðk  1Þ ∇v dΩ

∇uðkÞ ∇v dΩ þ Z

Z

Z

is the normal derivative, i.e. the derivative in the outward normal i direction to the boundary ∂Ωq . Because the derivative of the Heaviside step function vðxÞ is equal to zero, then the local weak form equation (43) is transformed into the following simple integral equation: Z Z ∂uðk þ 1Þ ðλ þ η þ β Þ dΓ uðk þ 1Þ dΩ  i i ∂n Ωq ∂Ωq Z Z ∂uðkÞ dΓ ¼ ð2λ  βÞ uðkÞ dΩ þ i Ωiq ∂Ωq ∂n Z Z ∂uðk  1Þ dΓ þðη  λ  βÞ uðk  1Þ dΩ þ i ∂n Ωiq ∂Ωq Z þ Fðx; kÞ dΩ þ γ ½2gðu~ ki þ 1 Þðu~ ki þ 1  u~ ki Þ

5. Discretization and numerical implementation for the MLRPI method

Fðx; kÞvðxÞ dΩ

Ωiq 2Δt

177

!

ð45Þ

!# ∂ ϕj þ 1Þ dΓ uðk j i Ωq ∂Ωq ∂n j¼1 j¼1 " ! !# Z Z N N ∂ ϕj d ¼ ð2λ  βÞ ∑ ϕ d Ω Γ þ ∑ uðkÞ j j i Ωiq ∂Ωq ∂n j¼1 j¼1 " ! ! # Z Z N N ∂ ϕj  1Þ d þðη  λ  βÞ ∑ ϕ d Ω Γ þ ∑ uðk j j i Ωiq ∂Ωq ∂n j¼1 j¼1 Z þ Fðx; kÞ dΩ þ γ ½2gðuki þ 1 Þðuki þ 1  uki Þ N

ϕj dΩ  ∑ i

Z

Ωiq

þgðuki Þðuki þ 1  uki  1 Þ þ 2gðuki  1 Þðuki  uki  1 Þ;

ð46Þ

178

E. Shivanian / Ocean Engineering 89 (2014) 173–188

because u~ is the latest available approximation of u. For nodes which are located on the boundary, we have the following cases:

(b1) For nodes on the boundary 0 r y; z r 1 and x¼ 0, using Eq. (3) we have uk þ 1 ðxi Þ ¼ a0 ðxi ; ðk þ 1ÞΔtÞ:

ð47Þ

(b2) For nodes on the boundary 0 r y; z r1 and x ¼1, using Eq. (3) we have uk þ 1 ðxi Þ ¼ a1 ðxi ; ðk þ 1ÞΔtÞ:

ð48Þ

(b3) For nodes on the boundary 0 r x; z r 1 and y¼ 0, using Eq. (4) we have uk þ 1 ðxi Þ ¼ b0 ðxi ; ðk þ 1ÞΔtÞ:

ð49Þ

(b4) For nodes on the boundary 0 r x; z r 1 and y¼1, using Eq. (4) we have uk þ 1 ðxi Þ ¼ b1 ðxi ; ðk þ 1ÞΔtÞ:

ð50Þ

(b5) For nodes on the boundary 0 r x; y r 1 and z ¼0, using Eq. (5) we have uk þ 1 ðxi Þ ¼ c0 ðxi ; ðk þ 1ÞΔtÞ:

ð51Þ

(b6) For nodes on the boundary 0 r x; y r 1 and z¼ 1, using Eq. (5) we have uk þ 1 ðxi Þ ¼ c1 ðxi ; ðk þ 1ÞΔtÞ:

ð52Þ

G½U ðk þ 1Þ i ¼ G½U ðkÞ i ¼ G½U ðk  1Þ i ¼ 0; ( 1; j ¼ i 8 j : Bij ¼ Cij ¼ 0; Aij ¼ 0; ja i and for the vector Fk , we assign (1) If xi is belong to the boundary then Fki ¼ a0 ðxi ; ðk þ 1ÞΔtÞ. (2) If xi is belong to the boundary then Fki ¼ a1 ðxi ; ðk þ 1ÞΔtÞ. (3) If xi is belong to the boundary then Fki ¼ b0 ðxi ; ðk þ 1ÞΔtÞ. (4) If xi is belong to the boundary then Fki ¼ b1 ðxi ; ðk þ 1ÞΔtÞ. (5) If xi is belong to the boundary then Fki ¼ c0 ðxi ; ðk þ 1ÞΔtÞ. (6) If xi is belong to the boundary then Fki ¼ c1 ðxi ; ðk þ 1ÞΔtÞ.

(b1) described in Eqs. (47)–(52) (b2) described in Eqs. (47)–(52) (b3) described in Eqs. (47)–(52) (b4) described in Eqs. (47)–(52) (b5) described in Eqs. (47)–(52) (b6) described in Eqs. (47)–(52)

At first step, when k ¼0, according to the initial conditions that were introduced in Eq. (2), we apply the following assumptions: U ð0Þ ¼ Φ ¼ ½φðx1 Þ; φðx2 Þ; …; φðxN ÞT

ð57Þ

U ð  1Þ ¼ U ð1Þ  2Δt Ψ ¼ U ð1Þ  2Δt½ψ ðx1 Þ; ψ ðx2 Þ; …; ψ ðxN ÞT

ð58Þ

Therefore Eq. (55) is converted to the following AU ð1Þ ¼ 2γ G½U ð1Þ nðU ð1Þ  U ð0Þ Þ þ γ G½U ð0Þ nU ð1Þ þ BU ð0Þ þ CðU ð1Þ 2Δt Ψ Þ þ F0 þ 2γ G½U ð1Þ  2Δt Ψ 

nðU

 U ð1Þ þ2Δt Ψ Þ  γ G½U ð0Þ nðU ð1Þ  2Δt Ψ Þ

ð0Þ

¼ 2γ G½U ð1Þ nðU ð1Þ  U ð0Þ Þ þ BU ð0Þ þCðU ð1Þ  2Δt Ψ Þ þ F0 þ 2γ G½U ð1Þ  2Δt Ψ nðU ð0Þ  U ð1Þ þ 2Δt Ψ Þ þ 2γΔt Ψ nG½U ð0Þ ;

The matrix forms of Eqs. (46)–(52) for all N nodal points in domain and on boundary of the problem are given below: " # N

N

j¼1

j¼1

þ 1Þ ðλ þ η þ β Þ ∑ Aij  ∑ Bij uðk j

¼ γ ½2gðuki þ 1 Þðuki þ 1  uki Þ þ gðuki Þuki þ 1  "

N

N

ð2λ  β Þ ∑ Aij þ ∑ j¼1

"

j¼1

#

N

j¼1

j¼1

þ ðη  λ  β Þ ∑ Aij þ ∑

 1Þ Bij uðk j

Ω

i q

Z Bij ¼

∂ ϕj dΓ ; i ∂n ∂Ωq

F ðkÞ i ¼

Z Ω

i q

Fðx; kÞ dΩ:

ð60Þ

ð53Þ

It can be shown that if h, distance between to consecutive nodes in each direction is sufficiently small, then rq in γ ¼ 2π r 3q =3Δt would be small enough, consequently the iterates U ð1Þ ½j þ 1 will converge to U ð1Þ as j-1. In general, we have the following iteration process for each time level:

ð54Þ

þ 1Þ ðk þ 1Þ þ 1Þ AU ðk nðU ðk  U ðkÞ Þ ½j þ 1 ¼ 2γ G½U ½j ½j

Assuming Aij ¼ ðλ þ η þ β ÞAij  Bij , Bij ¼ ð2λ  β ÞAij þ Bij , Cij ¼ ðη  λ  ðkÞ ðkÞ T βÞAij þ Bij , Fk ¼ ½F ðkÞ 1 ; F 2 ; …; F N  , U ¼ ðui ÞN1 and G½U ¼ ½gðui ÞN1 ,

Eq. (53) yields AU ðk þ 1Þ ¼ 2γ G½U ðk þ 1Þ nðU ðk þ 1Þ  U ðkÞ Þ þ γ G½U ðkÞ nU ðk þ 1Þ þ BU ðkÞ þ CU ðk  1Þ þ Fk þ 2γ G½U ðk  1Þ nðU ðkÞ  U ðk  1Þ Þ  γ G½U ðkÞ nU ðk  1Þ ;

we need to solve the above nonlinear system of To obtain U equations. Traditional rootfinding methods such as Newton's method can be applied to Eq. (59) to find U ð1Þ ; but often that is very time-consuming process. Instead, Eq. (59) is usually solved by ð1Þ a simple iteration technique. Given an initial guess U ½0 ¼ U ð0Þ , ð1Þ define U ð1Þ ; U , etc., by ½1 ½2   ð1Þ ð1Þ ð1Þ ð0Þ AU ½j þ 1 ¼ 2γ G½U ð1Þ  2Δt Ψ Þ þF0 þ BU ð0Þ þ CðU ½j ½j n U ½j  U

j ¼ 0; 1; 2; …:

#

k1 þ F ðkÞ Þðuki  uki  1 Þ  gðuki Þuki  1 ; i þ γ ½2gðui

where Z ϕj dΩ; Aij ¼

ð59Þ

ð1Þ

ð1Þ ð1Þ þ 2γ G½U ½j  2Δt Ψ nðU ð0Þ  U ½j þ 2Δt Ψ Þ þ2γΔt Ψ nG½U ð0Þ ;

Bij uðkÞ j

N

ð56Þ

ð55Þ

where the operator “n” denotes Hadamard product (component by component production). Furthermore, to satisfy Eqs. (47)–(52), for all nodes belong to the boundary, i.e. xi A ∂Ω, we set for each step

þ 1Þ þBU ðkÞ þ CU ðk  1Þ þ γ G½U ðkÞ nU ðk ½j

þ Fk þ 2γ G½U ðk  1Þ nðU ðkÞ  U ðk  1Þ Þ  γ G½U ðkÞ nU ðk  1Þ ;

ð61Þ

where j denotes the number of iterations in each time level. This process is iterated until the unknown quantity converges to within a prescribed number of tolerance. In other words in each time level, these iterations are continued until satisfying the following condition:

þ 1Þ ðk þ 1Þ ‖U ðk ‖1 r ϵ; ½j þ 1  U ½j

ð62Þ

E. Shivanian / Ocean Engineering 89 (2014) 173–188

where ϵ is a fixed number, in this case is equal to 10  10 and ‖  ‖1 is infinity norm.

N ∂s ϕ ðxÞ ∂s uðxÞ j ¼ ∑ uj ; s ∂z ∂zs j¼1

Suppose that the number of total nodes covering Ω ¼ ½0; 13 is N as previous section. On the other hand, we know that n is dependent on point of interest x (so, after that we call it nx ) in Eq. (26) which is the number of nodes included in support domain Ωx corresponding to the point of interest x (for example Ωx can be a sphere centered at x with radius rs (see Fig. 1)). Therefore, we have nx r N and Eq. (26) can be modified as

∂x

N

uðxÞ ¼ Φ ðxÞUs ¼ ∑ ϕj ðxÞuj :

ð63Þ

0

j¼1

c

ð64Þ

0

The derivatives of uðxÞ are easily obtained as

ð65Þ

N ∂s ϕ ðxÞ ∂s uðxÞ j ¼ ∑ uj ; ∂ys ∂ys j¼1

1

1

∂ϕ1 ðxN Þ ∂z

∂ϕ2 ðxN Þ ∂z

⋯ ⋯

⋮

⋱ ⋯

ðsÞ U ðsÞ y ¼ D yU;

0 1 u1 C ∂ϕN ðx2 Þ CB u C CB 2 C ∂y C; CB @ ⋮ A ⋮ C A ∂ϕN ðxN Þ uN

z

0.5

0.5 x

1

0

∂y

0.12

0.8

ðsÞ U ðsÞ z ¼ D zU;

z

0.08

0.4

0.06

0.3

0.6

0.2

0.3

0.4

0.15

0.2

0.2

0.2

0 1

0.1

0.4

0.5 y

0

0

0.5 x

1

0 1

0

0.1 0.05 0.5 y

0

0.06

1

0.8

0.05

0.8

0.6

0.04

0.6

0.4

0.03 0.02

0.04

0.2

0 1

0.02

0 1

1

0

0.01 0.5 y

0

0

0

0.5 x

1

0

t=2

1

0.2

0.5 x

ð70Þ

0.25

z

0.1

0.6

ð69Þ

∂z

t=1.6 0.14

1

0 1 u1 C CB u2 C C CB C; CB @ ⋮ A ⋮ C A ∂ϕN ðxN Þ uN ∂ϕN ðx1 Þ ∂z ∂ϕN ðx2 Þ ∂z

0.8

t=1.2 1

ð68Þ

0.5

0.4

0.2

1

1

0.6

0.6

0.4

0

∂ϕ2 ðx1 Þ ∂z ∂ϕ2 ðx2 Þ ∂z

⋱ ⋯

∂ϕN ðx1 Þ ∂y

t=0.8

z

0.6

0

⋯

⋮

0.7

0.8

0.5 y

∂ϕ2 ðx2 Þ ∂y

∂ϕ1 ðx1 Þ ∂z ∂ϕ1 ðx2 Þ ∂z

⋮

⋯

z

0.8

0

∂ϕ2 ðx1 Þ ∂y

∂ϕ2 ðxN Þ ∂y

0

C B C B C¼B C B A B @

1

1.5

0

∂y

t=0.4

1

0.5 y

1 ðx1 Þ

∂y

ðsÞ U ðsÞ x ¼ D xU;

t=0

0 1

0 ∂ϕ

∂x

we denote above coefficients matrices as Dx, Dy and Dz, respectively. Also, for high order derivatives we have the following matrix-vector multiplications:

and also high derivatives of uðxÞ are easily given as N ∂s ϕ ðxÞ ∂s uðxÞ j ¼ ∑ uj ; ∂xs ∂xs j¼1

u0z1

B u0 B z2 B B ⋮ @ u0zN

N ∂ϕ ðxÞ ∂uðxÞ j ¼ ∑ u; ∂y ∂y j j¼1

N ∂ϕ ðxÞ ∂uðxÞ j ¼ ∑ uj ; ∂z ∂z j¼1

1

∂x

x

0.5

1

0

0.025 0.02 0.015

z

N ∂ϕ ðxÞ ∂uðxÞ j ¼ ∑ uj ; ∂x ∂x j¼1

u0y1

B u0 C B ∂ϕ1 ðx2 Þ B y2 C B B C B ∂y B ⋮ C¼B B @ A @ ⋮ ∂ϕ1 ðxN Þ u0yN

In fact, corresponding to node x j there is a shape function ϕj ðxÞ, c j ¼ 1; 2; 3; …; N, we define Ωx ¼ xj : xj 2 = Ωx then it is clear from the previous section that 8 xj A Ωx : ϕj ðxÞ ¼ 0:

ð66Þ

where ∂s ðÞ=∂xs , ∂s ðÞ=∂ys and ∂s ðÞ=∂zs are s’th derivative with s s respect to x, y and z, respectively. Denoting uðsÞ x ðÞ ¼ ∂ ðÞ=∂x and ðsÞ s s uy ðÞ ¼ ∂ ðÞ=∂y , and setting x ¼ xi in Eq. (65), then we can represent the discrete differentiations process as a matrix-vector multiplications 1 0 0 1 0 ∂ϕ1 ðx1 Þ ∂ϕ2 ðx1 Þ ∂ϕN ðx1 Þ 0 1 ux1 ⋯ u1 ∂x ∂x ∂x B C B u0 C B ∂ϕ1 ðx2 Þ ∂ϕ2 ðx2 Þ ∂ ϕ ðx Þ C B C 2 N B x2 C B ⋯ CB u2 C ∂x ∂x ∂x B C ð67Þ CB C; B ⋮ C¼B @ ⋮ A ⋮ ⋮ ⋱ ⋮ C @ A B @ A ∂ϕ1 ðxN Þ ∂ϕ2 ðxN Þ u0xN uN ⋯ ∂ϕN ðxN Þ

6. High order derivatives by operational matrices for SMRPI

T

179

0.4

0.01

0.2 0 1

0.005 0.5 y

0

0

0.5 x

1

0

Fig. 2. Initial condition and numerical contours of Example 1 by MLRPI, at times t ¼ 0; 0:4; 0:8; 1:2; 1:6 and 2 with Δt ¼ 0:1 and h ¼ 0.1 (N ¼1331).

180

E. Shivanian / Ocean Engineering 89 (2014) 173–188

where U xðsÞ

ðsÞ ðsÞ T ðuðsÞ x1 ; ux2 ; …; uxN Þ ;

ð71Þ

ðsÞ ðsÞ T U yðsÞ ¼ ðuðsÞ y1 ; uy2 ; …; uyN Þ ;

ð72Þ

ðsÞ ðsÞ T U zðsÞ ¼ ðuðsÞ z1 ; uz2 ; …; uzN Þ ;

ð73Þ

DðsÞ xij ¼

∂s ϕj ðxi Þ ; ∂xs

ð74Þ

DðsÞ yij ¼

∂ ϕj ðxi Þ ; ∂ys

¼

s

ð75Þ

∂s ϕj ðxi Þ ; DðsÞ zij ¼ ∂zs

ð76Þ

U ¼ ðu1 ; u2 ; …; uN ÞT :

ð77Þ

We conjecture here that there is a relation between the matrix Dx and DðsÞ x (also between Dy and DðsÞ y and, between Dz and DðsÞ z) or how to choose nodes distribution so that there exists a relation between them, then it would be a great advantage for SMRPI, because we can easily compute the DðsÞ x from Dx (or DðsÞ y from Dy, or DðsÞ z from Dz) and apply to high order partial differential equation without difficulty. Overall, these conjectures would be challenge opportunity for the researchers in near future. c It is definitely right if we say 8 xj A Ωx : ∂s ϕj ðxÞ=∂xs ¼ ∂s ϕj ðxÞ=∂ys ¼ ∂s ϕj ðxÞ=∂zs ¼ 0, s ¼ 1; 2; 3; … due to Eq. (64), as the result of this fact all DðsÞ y, DðsÞ x and DðsÞ z are banded matrices.

7. Discretization and numerical implementation for SMRPI method In this section, we consider Eq. (35) and explain how to implement spectral meshless radial point interpolation (SMRPI) method to obtain discrete equations. Eq. (35) can be rewritten as  2 kþ1  ∂ u ðxÞ ∂2 uk þ 1 ðxÞ ∂2 uk þ 1 ðxÞ þ þ ðλ þ η þ β Þuk þ 1  ∂x2 ∂y2 ∂z2  2 k  2 k 2 k ∂ u ðxÞ ∂ u ðxÞ ∂ u ðxÞ ¼ ð2λ  βÞuk þ þ þ ∂x2 ∂y2 ∂z2  2 k1  2 k1 ∂ u ðxÞ ∂ u ðxÞ ∂2 uk  1 ðxÞ þ ðη  λ  β Þuk  1 þ þ þ ∂x2 ∂y2 ∂z2 þ

δ ½2gðu~ k þ 1 Þðu~ k þ 1  u~ k Þ þ gðu~ k Þðu~ k þ 1  u~ k  1 Þ 2Δt

þ 2gðu~ k  1 Þðu~ k  u~ k  1 Þ þ Fðx; kÞ:

ð78Þ

0.8 0.7

u(x,y,0.5)

0.6 0.5 0.4 0.3 0.2 0.1 0 0

0.2

0.4

0.6

x

0.5 0.8

1

0

1

y

Fig. 3. Numerical behavior of uðx; y; 0:5Þ by passing the time with h ¼ 0.1 and time step Δt ¼ 0:1 for Example 1.

Consider N nodes with arbitrary distribution on the boundary and domain of the problem (see Fig. 1). Assuming that uðxi ; kΔtÞ; i ¼ 1; 2; …; N are known, our aim is to compute uðxi ; ðk þ1Þ ΔtÞ; i ¼ 1; 2; …; N. So, we have N unknowns and to compute these unknowns we need N equations. As it will be described, corresponding to each node we obtain one equation. For nodes which are located in the interior of the domain, to obtain the discrete equations from (78), substituting approximation formulas (63) and (66) into Eq. (78) yields ðλ þ η þ β Þuik þ 1 

N

∂2 ϕj ðxÞ

j¼1

∂x2

∑

¼ ð2λ  βÞuki þ

N

∂2 ϕj ðxÞ

j¼1

∂x2

∑

þðη  λ  βÞuki  1 þ N

∂2 ϕj ðxÞ

j¼1

∂z2

þ ∑

N

∂2 ϕj ðxÞ

j¼1

∂y2

ukj þ 1 þ ∑

N

∂2 ϕj ðxÞ

j¼1

∂y2

ukj þ ∑

N

∂2 ϕj ðxÞ

j¼1

∂x2

∑

N

∂2 ϕj ðxÞ

j¼1

∂z2

ukj þ 1 þ ∑

N

∂2 ϕj ðxÞ

j¼1

∂z2

ukj þ ∑

N

∂2 ϕj ðxÞ

j¼1

∂y2

ukj  1 þ ∑

! ukj þ 1

! ukj

ukj  1

! ukj  1

Table 1 The ARE error for different time levels and different number of nodal points for Example 1 with Δt ¼ 0:1 by MLRPI. t

N ¼ 125 ðh ¼ 0:25Þ

N ¼ 216 ðh ¼ 0:2Þ

N ¼ 1331 ðh ¼ 0:1Þ

0.000000 0.200000 0.400000 0.600000 0.800000 1.000000 1.200000 1.400000 1.600000 1.800000 2.000000 2.200000 2.400000 2.600000 2.800000 3.000000

0.000000eþ 000 3.232892e  003 1.075105e 002 1.370814e  002 7.598482e 003 7.597848e  003 1.388343e  002 1.246365e  002 2.924713e  002 3.718199e  002 1.699281e  002 5.330926e  002 5.789303e  002 4.844782e  002 1.263714e  001 1.177617e  001

0.000000eþ000 5.227482e 003 1.652042e  002 2.067552e  002 1.134996e 002 1.172441e  002 2.107833e  002 1.892462e  002 4.406778e  002 5.608845e  002 2.866204e 002 8.042301e  002 9.481466e  002 8.873924e 002 1.920666e  001 2.635218e 001

0.000000e þ000 1.746704e  002 5.028264e  002 6.027377e  002 3.204873e  002 3.595841e  002 6.088789e  002 6.001171e  002 1.276347e  001 1.619998e  001 9.539429e  002 2.222164e  001 2.681669e  001 2.244535e  001 5.314648e  001 5.643713e 001

Table 2 The ARE error for different time levels and different number of nodal points for Example 1 with Δt ¼ 0:01 by SMRPI. t

N ¼ 125 ðh ¼ 0:25Þ

N ¼ 216 ðh ¼ 0:2Þ

N ¼ 1331 ðh ¼ 0:1Þ

0.000000 0.200000 0.400000 0.600000 0.800000 1.000000 1.200000 1.400000 1.600000 1.800000 2.000000 2.200000 2.400000 2.600000 2.800000 3.000000 3.200000 3.400000 3.600000 3.800000 4.000000

0.000000eþ 000 9.209334e  005 4.440123e  005 6.282209e  005 9.649442e  005 8.402105e  005 1.593789e  004 1.192284e 004 1.644851e  004 2.210102e  004 2.859203e  004 3.748244e  004 3.078280e  004 3.333790e  004 4.332384e  004 8.005493e  004 1.122433e 003 1.102172e  003 1.073087e  003 1.251636e  003 1.616210e  003

0.000000eþ000 5.249165e  005 3.624698e  005 2.920665e 005 6.949653e  005 5.320248e  005 7.872472e  005 8.399698e  005 1.205528e  004 1.340773e  004 1.493039e  004 1.589036e  004 2.910061e 004 3.078231e  004 3.399931e  004 3.191275e  004 6.812933e  004 6.750575e 004 7.011417e  004 8.047384e  004 1.384964e 003

0.000000e þ000 1.711119e 005 1.791014e  005 1.311040e  005 2.145122e  005 3.007383e  005 1.798631e  005 5.193769e 005 4.566172e  005 5.006974e  005 8.442648e  005 8.317509e  005 7.304358e  005 1.773245e  004 1.255156e  004 1.890204e  004 2.909437e  004 2.207709e  004 3.540545e  004 5.518299e  004 2.995422e  004

E. Shivanian / Ocean Engineering 89 (2014) 173–188

δ

½2gðu~ k þ 1 Þðu~ k þ 1  u~ k Þ þgðu~ k Þðu~ k þ 1  u~ k  1 Þ 2 Δt þ 2gðu~ k  1 Þðu~ k  u~ k  1 Þ þ Fðx; kÞ;

ð79Þ

now, setting x ¼ xi , i ¼ 1; 2; 3; …; NΩ (N Ω is the number of nodes in interior of Ω) in the above equation and applying notations (71)– (77) implies ! ðλ þ η þ βÞuki þ 1 

N

N

because u~ is the latest available approximation of u. For nodes which are located on the boundary, we have the same Eqs. (47)– (52). The matrix forms of Eqs. (46)–(52) for all N nodal points in domain and on boundary of the problem are given below: " #

N

j¼1

N

¼ ð2λ  βÞuki þ

j¼1

N

!

N

∑ Dð2Þ xij ukj þ ∑ Dð2Þ yij ukj þ ∑ Dð2Þ zij ukj

j¼1

þ ðη  λ  β Þuki  1 þ

j¼1

j¼1

1

!

N

N

N

j¼1

j¼1

j¼1

∑ Dð2Þ xij ukj  1 þ ∑ Dð2Þ yij ukj  1 þ ∑ Dð2Þ zij ukj  1

0.9 0.8

þ γ ½2gðu~ ki þ 1 Þðu~ ki þ 1  u~ ki Þ þgðu~ ki Þðu~ ki þ 1  u~ ki  1 Þ

0.7

where γ ¼ δ=2Δt, or equivalently ðλ þ η þ β

N

∑ D

ð2Þ

j¼1

xij ukj þ 1 þ

N

¼ ð2λ  βÞuki þ

N

∑ D

j¼1

ð2Þ

yij ukj þ 1 þ

N

u(x,0.5,0.5)

þ 2gðu~ ki  1 Þðu~ ki  u~ ki  1 Þ þ Fðxi ; kÞ;

Þuki þ 1 

!

N

∑ D

ð2Þ

j¼1

zij ukj þ 1

!

N

þ ðη  λ  β Þuki  1 þ !

N

j¼1

j¼1

N

N

j¼1

j¼1

0.6 0.5 0.4 0.3

∑ Dð2Þ xij ukj þ ∑ Dð2Þ yij ukj þ ∑ Dð2Þ zij ukj

j¼1

N

j¼1

¼ γ ½2gðuki þ 1 Þðuki þ 1 uki Þ þgðuki Þuki þ 1 

∑ Dð2Þ xij ukj þ 1 þ ∑ Dð2Þ yij ukj þ 1 þ ∑ Dð2Þ zij ukj þ 1

j¼1

N

j¼1

þ 1Þ ðλ þ η þ βÞ ∑ δij  ∑ Bij uðk j

0.2

∑ Dð2Þ xij ukj  1 þ ∑ Dð2Þ yij ukj  1

0.1 0

þ ∑ Dð2Þ zij ukj  1 þ γ ½2gðuki þ 1 Þðuki þ 1  uki Þ þ gðuki Þðuki þ 1 uki  1 Þ

0

0.2

0.4

0.6

j¼1

0.8

0.4

0.2

y

0 1 0.5

0.2 0

0

x

0.5

1

time=1.2s

1 0.8

0.4

0.6

0.4

0.3

0.4

0.2

0.2

0.2

0.2

0.15

0.25

0.8

0.2

0.6

0.4

0.15

0.4

0.2

0.1

0.2

y

0 1 0.5

0.05 0

0

x

0.5

1

0

0 1 0.5 0 y

0 1 0.5 y 0

0.1 0.5 x

0

1

0

time=1.6s

0.35 0.3 0.25

0.1

0.05

0.5 x

0.05 0

1

0

0.5 x

1

time=2.0s 0.12

0.8

0.15

0

0.1

1

z

z

0.6

0.6

1

0.3

0.4

0.8

0 1 0.5 0 y

0

time=0.8s

0.5

z

z

0.6

0.4

1

0.6

0.8

0.8

0.6

time=0.4s

1

1

z

time=0.0s

1

Fig. 5. Numerical behavior of uðx; 0:5; 0:5Þ by passing the time with h ¼ 0.1 and time step Δt ¼ 0:1 for Example 2.

þ 2gðuki  1 Þðuki  uki  1 Þþ Fðxi ; kÞ;

1

0.8

x

0.1

0.6

0.08

0.4

0.06

0.2

0.04

z

þ

181

0 1 0.5 y 0

0.02 0

0.5 x

1

Fig. 4. Initial condition and numerical contours of Example 2 by MLRPI, at times t¼ 0, 0.4, 0.8, 1.2, 1.6 and 2 with Δt ¼ 0:1 and h ¼ 0.1 (N ¼ 1331).

182

E. Shivanian / Ocean Engineering 89 (2014) 173–188

"

N

N

#

ð2λ  β Þ ∑ δij þ ∑ Bij uðkÞ j j¼1

"

j¼1

N

N

j¼1

j¼1

#

 1Þ þ ðη  λ  β Þ ∑ δij þ ∑ Bij uðk j

þ F ðkÞ i þ where (

δij ¼

γ

½2gðuki  1 Þðuki  uki  1 Þ  gðuki Þuki  1 ;

1;

i ¼ j;

0;

i a j;

Bij ¼ Dð2Þ xij þ Dð2Þ yij þ Dð2Þ zij ;

F ðkÞ i ¼

Z Ωiq

ð80Þ

Fðx; kÞ dΩ:

ð81Þ

Suppose that Aij ¼ ðλ þ η þ β Þδij  Bij , Bij ¼ ð2λ  βÞδij þ Bij , Cij ¼ ðη  λ  β Þδij þ Bij , then all equations (55)–(62) hold for the SMRPI method as well.

8. Numerical experiments In this section, we show the results obtained for four examples using the meshless methods described above. In these examples,

the domain integrals are evaluated with 8 points Gaussian quadrature rule while the boundary integrals are evaluated with 9 points Gaussian quadrature rule in the MLRPI method. To show the behavior of the solution and the efficiency of the MLRPI and proposed SMRPI methods, the following average relative error (ARE) and relative error with infinity norm (REI) are applied to make comparison vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N u∑ ðU exact ðxi Þ  U approx ðxi ÞÞ2 ARE ¼ t i ¼ 1 2 ∑N i ¼ 1 ðU exact ðx i Þ  U exact ðx i ÞÞ max1 r i r N jU exact ðxi Þ  U approx ðxi Þj REI ¼ max1 r i r N jU exact ðxi Þ  U exact ðxi Þj where U exact ðxi Þ and U approx ðxi Þ are achieved by exact and approximate solutions on xi points, N is number of nodal points and the over bar denotes the average. In the following examples (except the last one) the regular node distribution is used. Also in order to implement the meshless local weak form (MLRPI), the radius of the local quadrature domain (sphere) r q ¼ 0:7h is selected, where h is the distance between the nodes in x, y or z direction. The radius of support domain (sphere) to local radial point interpolation

Table 3 The ARE error for different time levels and different number of nodal points for Example 2 with Δt ¼ 0:1 by MLRPI. t

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0

N ¼ 216 ðh ¼ 0:2Þ

N ¼ 1331 ðh ¼ 0:1Þ

κ¼1

κ¼2

κ¼3

κ¼1

κ¼2

κ¼3

0.000000eþ 000 4.009916e  002 8.162033e  002 4.220947e 002 3.060119e 002 5.608115e  002 1.033977e 001 1.205526e 002 6.866501e  002 1.162832e 001 1.319023e  001 7.852188e  002 8.550022e  002 1.648103e 001 1.742357e  001 1.609675e  001

0.000000eþ 000 4.442574e 002 9.298299e  002 4.444208e  002 6.222293e  002 8.699010e  002 1.808660e 001 8.176070e 002 2.018298e  001 2.943552e  001 3.453811e  001 4.079489e  001 4.455365e  001 6.917634e  001 7.139160e 001 1.104397e þ000

0.000000eþ000 4.958507e  002 1.031252e  001 4.581301e  002 1.291325e  001 1.525989e  001 3.462113e  001 3.061031e 001 5.671946e  001 8.555966e  001 1.123776eþ 000 1.689115e þ000 1.980682eþ 000 3.369709eþ000 3.659398e þ000 6.969116eþ 000

0.000000e þ000 5.475741e  002 1.021499e  001 2.608464e  002 2.940503e  002 8.433469e  002 1.317274e  001 2.294165e 002 5.546278e 002 1.459881e  001 1.554682e  001 1.057924e  001 8.138474e  002 2.463378e  001 1.646810e  001 2.342819e  001

0.000000e þ000 6.155718e  002 1.167794e 001 1.377169e  002 7.096877e  002 1.297843e 001 2.181628e 001 1.282529e  001 2.016637e 001 3.676951e 001 4.099992e  001 5.578408e  001 4.975496e 001 1.071231e þ000 6.809457e  001 1.721901eþ 000

0.000000eþ 000 7.266897e  002 1.402382e  001 2.451617e 002 1.490538e  001 2.321455e 001 3.967552e 001 4.055053e  001 6.358388e 001 1.051686e þ000 1.379142e þ000 2.510127e þ000 2.660900eþ 000 5.417346e þ000 4.620693eþ 000 1.135441eþ 001

Table 4 The ARE error for different time levels and different number of nodal points for Example 2 with Δt ¼ 0:1 by SMRPI. t

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0

N ¼ 216 ðh ¼ 0:2Þ

N ¼ 1331 ðh ¼ 0:1Þ

κ¼1

κ¼2

κ¼3

κ¼1

κ¼2

κ¼3

0.000000eþ 000 7.934007e  003 1.399029e 002 1.418313e 002 1.012631e  002 1.622013e 002 7.559992e 003 7.625957e  003 1.935730e  002 2.628794e 002 2.047448e  002 6.658618e  003 3.769284e 002 7.696695e 003 3.277411e  002 4.912468e  002

0.000000eþ 000 9.405547e  003 1.286754e  002 1.553698e  002 1.806678e 002 2.771652e  002 1.141390e  002 3.111662e  002 3.728080e  002 4.918937e 002 9.124677e  002 2.424747e  002 1.469487e 001 1.075185e  001 1.540705e  001 2.525204e 001

0.000000eþ000 1.022862e  002 1.587702e  002 2.320886e  002 1.749547e  002 2.565103e 002 3.979277e  002 2.237179e  002 7.031231e 002 1.282080e  001 4.197244e 002 1.858540e  001 2.641417e  001 7.214897e  002 4.442677e  001 7.098267e  001

0.000000e þ000 9.906430e 004 1.513334e  003 8.082008e 004 4.622793e 004 1.795780e  003 6.502347e  004 1.203809e  003 2.512350e 003 1.748893e  003 1.666536e 003 2.291035e  003 3.038723e 003 2.187851e 003 3.857794e  003 4.762359e  003

0.000000e þ000 5.897250e  003 9.924226e 003 5.139951e  003 8.032542e  003 1.703924e 002 7.974590e  003 1.736484e  002 3.522403e 002 1.848321e  002 3.149295e  002 1.426062e  002 3.917556e  002 5.930850e  002 1.260325e  001 8.868599e  002

0.000000eþ 000 1.017353e 003 1.223196e 003 1.257071e  003 1.649337e  003 2.041814e  003 3.039767e  003 3.345790e  003 5.628526e  003 8.895086e  003 6.913095e  003 1.303440e  002 2.141918e  002 1.396238e  002 3.441102e 002 5.423630e  002

E. Shivanian / Ocean Engineering 89 (2014) 173–188

method is r s ¼ 4r q and r s ¼ 3:2h for MLRPI and SMRPI, respectively. This size is significant enough to have sufficient number of nodes (nx ) and gives an appropriate approximation. Also, in Eq. (6), we set m ¼ 10 and m ¼20, i.e. the quadratic and cubic basis functions ((8) and (10)) are used for MLRPI and SMRPI, respectively.

Example 2 (Mohanty and Gopal, 2013, Van der Pol type nonlinear wave equation). In this Example, we take α ¼ δ ¼ κ , β ¼ 0 and gðuÞ ¼ u2 . We assume that uðx; y; z; tÞ ¼ sin ðπ xÞ sin ðπ yÞ sin ðπ zÞ expð  κ tÞ is the exact solution of the problem (1)–(5) where φðx; y; zÞ, ψ ðx; y; zÞ, a0 ðy; z; tÞ, a1 ðy; z; tÞ, b0 ðx; z; tÞ, b1 ðx; z; tÞ, c0 ðx; y; tÞ, c1 ðx; y; tÞ and f ðx; y; z; tÞ are defined accordingly.

Example 1 (Mohanty and Gopal, 2013, Three space dimensional linear telegraphic hyperbolic equation). We set α ¼ β ¼ 2 and δ ¼ 0. The exact solution of the problem (1)–(5) is taken as

As previous example, the numerical contour solutions of this problem, with time step Δt ¼ 0:1, at t¼0, 0.4, 0.8, 1.2, 1.6 and 2 are shown in Fig. 4. It is found that the solution in the whole of cube tends to zero which is in full agreement with exact solution. In Fig. 5, the numerical simulation of uðx; y; z; tÞ when y ¼ z ¼ 12, i.e. uðx; 0:5; 0:5Þ has been plotted for a large range of times (t ¼ 0:2i; i ¼ 0; 1; 2; …; 15). As it is clear from this figure, uðx; 0:5; 0:5Þ goes to zero slowly by passing time. In Tables 3 and 4, as previous example, ARE error of numerical solutions for different time levels, different number of nodal points and different values of κ have been reported with MLRPI and SMRPI, respectively. It is observed again that ARE does not converge with respect to N for MLRPI (Table 3) whereas convergence does occur for SMRPI (Table 4). Moreover, the accuracy of these two methods shows that SMRPI is obviously superior to MLRPI.

uðx; y; z; tÞ ¼ sinhðxÞsinhðyÞsinhðzÞexpð  2tÞ; where φðx; y; zÞ, ψ ðx; y; zÞ, a0 ðy; z; tÞ, a1 ðy; z; tÞ, b0 ðx; z; tÞ, b1 ðx; z; tÞ, c0 ðx; y; tÞ, c1 ðx; y; tÞ and f ðx; y; z; tÞ are defined accordingly. The numerical contour solutions of this problem, with time step Δt ¼ 0:1, at t ¼ 0; 0:4; 0:8; 1:2; 1:6 and 2 are shown in Fig. 2. The numerical simulation of uðx; y; z; tÞ when z ¼ 12, i.e. uðx; y; 0:5Þ has been shown in Fig. 3 for a large range of times (t ¼ 0:2i; i ¼ 0; 1; 2; …; 15). As it is clear from this figure, uðx; y; 0:5Þ goes to zero slowly by passing time, this is in full agreement with the exact solution. In Tables 1 and 2, ARE error of numerical solutions for different time levels and different number of nodal points have been reported with MLRPI and SMRPI, respectively. It can be seen that ARE does not converge with respect to N for MLRPI (Table 1) whereas convergence does occur for SMRPI (Table 2). Moreover, the accuracy of these two methods shows that SMRPI is obviously superior to MLRPI.

time=0.1s

1

0.1

183

Example 3 (Mohanty and Gopal, 2013, Dissipative nonlinear wave equation). We set α ¼ β ¼ 0, δ ¼ 2 and gðuÞ ¼ u. The exact

time=0.4s

1

time=0.8s

1

0.6

0

0

0.5

x

1

0 1 y 0.5

0

0.2

0 0

x

0.5

1

1

z 0

0.8 0.6

0.6

x

0.5

1

time=2.8s

1

0.6

0.2 0

0

time=2.0s

0.5

0 1 y 0.5 0

time=2.4s

1

x

0.5

0.5

0.4

0.4

0.2 1

0

1

0.4

0.5

0

0.8

0.6

x

0

time=1.6s

0.5

0 0

0.2 0 1 0.5 y

1

0.8

0 1 0.5 y

0.4

0.5

0.1

time=1.2s

1

z

0.5

z

0 1 y 0.5

0.05

z

0.5

z

z

0.3

0 1 y 0.5

0

1

0.3

0.2 0

0

0.5

1

0

x

time=3.1s

0.2

0.5

z

0.5

z

z

0.03 0.4

0

0

x

0.5

1

0.02

0.1

0.2 0 1 y 0.5

0.5

0 1 0.5 y

0

0

0.5

x

1

0 1 0.5 y 0

0.01

0

x

0.5

1

0

Fig. 6. Diagram of numerical contours of Example 3 by MLRPI, at times t ¼0.1, 0.4, 0.8, 1.2, 1.6, 2, 2.4, 2.8 and 3.1 with Δt ¼ 0:1 and h ¼0.1 (N ¼ 1331).

184

E. Shivanian / Ocean Engineering 89 (2014) 173–188

Example 4 (Exponential growth with respect to the time). In this Example, we take α ¼ 20, β ¼ 200 and δ ¼ 0. We assume that

solution of problem (1)–(5) is taken as uðx; y; z; tÞ ¼ sin ðπ xÞ sin ðπ yÞ sin ðπ zÞ sin ðtÞ;

uðx; y; z; tÞ ¼ sin ðπ xÞ sin ðπ yÞ sin ðπ zÞ expð2tÞ

where φðx; y; zÞ, ψ ðx; y; zÞ, a0 ðy; z; tÞ, a1 ðy; z; tÞ, b0 ðx; z; tÞ, b1 ðx; z; tÞ, c0 ðx; y; tÞ, c1 ðx; y; tÞ and f ðx; y; z; tÞ are defined accordingly. The numerical contour solutions of this problem, with time step Δt ¼ 0:1, at t¼ 0.1, 0.4, 0.8, 1.2, 1.6, 2, 2.4, 2.8 and 3.1 are shown in Fig. 6. The numerical simulation of uðx; y; z; tÞ when y ¼ z ¼ 12, i.e. uðx; 0:5; 0:5Þ has been shown in Figs. 7 and 8 for a large range of times: at t ¼ 0:1i; i ¼ 0; 1; 2; …; 15 in Fig. 7 and at t ¼ 0:1i; i ¼ 16; 17; 18; …; 47 in Fig. 8. As it is clear from these figures, uðx; 0:5; 0:5Þ oscillates between  1 and 1 by passing time which is in full agreement with the exact result. In Tables 5 and 6, as previous examples, ARE error of numerical solutions for different time levels and different number of nodal points have been reported with MLRPI and SMRPI, respectively. An inspection of Table 5 shows that ARE does not converge with respect to N for MLRPI whereas convergence happens for SMRPI (Table 6). Moreover, the comparison of these two methods shows that SMRPI is obviously superior to MLRPI and it seems to be unconditionally stable. 1 0.9 0.8

u(x,0.5,0.5)

f ðx; y; z; tÞ ¼ ð244 þ 3π 2 Þ sin ðπ xÞ sin ðπ yÞ sin ðπ zÞexpð2tÞ: We consider two types of covering domain for this example first by regularly distributed nodes as previous examples and second by nonregularly distributed nodes as follows: (i) Regularly distributed nodes: As previous examples, the numerical contour solutions of this problem, with time step Δt ¼ 0:1, at t ¼0, 0.4, 0.8, 1.2, 1.6 and 2 are shown in Fig. 9. It is found that the solution in the whole of cube tends to infinity which is in full agreement with exact solution. In Fig. 10, the numerical simulation of uðx; y; z; tÞ when y ¼ z ¼ 12, i.e. uðx; 0:5; 0:5Þ has been plotted for a large range of times (t ¼ 0:1i; i ¼ 0; 1; 2; …; 20). As it is clear from this figure, uðx; 0:5; 0:5Þ goes to infinity exponentially by passing time. In Tables 7 and 8, REI error of numerical solutions for different time levels and different number of nodal points have been reported with MLRPI and SMRPI, respectively. An inspection of Table 7 shows that REI does not converge with respect to N for MLRPI whereas convergence happens for SMRPI (Table 8). Table 5 The ARE error for different time levels and different number of nodal points for Example 3 with Δt ¼ 0:1 by MLRPI.

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

is the exact solution of the problem (1)–(5) where φðx; y; zÞ, ψ ðx; y; zÞ, a0 ðy; z; tÞ, a1 ðy; z; tÞ, b0 ðx; z; tÞ, b1 ðx; z; tÞ, c0 ðx; y; tÞ and c1 ðx; y; tÞ are defined accordingly and f ðx; y; z; tÞ is as follows:

0

0.2

0.4

0.6

0.8

1

x Fig. 7. Numerical behavior of uðx; 0:5; 0:5Þ by passing the time (t ¼ 0; o:1; …; 1:5) with h ¼0.1 and time step Δt ¼ 0:1 for Example 3.

t

N ¼ 125 ðh ¼ 0:25Þ

N ¼ 216 ðh ¼ 0:2Þ

N ¼ 1331 ðh ¼ 0:1Þ

0.000000 0.200000 0.400000 0.600000 0.800000 1.000000 1.200000 1.400000 1.600000 1.800000 2.000000 2.200000 2.400000 2.600000 2.800000 3.000000

0.000000eþ 000 5.416840e  003 2.449643e  002 4.105962e 002 4.951100e  002 5.050371e  002 4.956679e  002 5.194146e  002 5.897987e  002 6.858276e  002 7.718151e 002 8.228449e  002 8.436246e  002 8.631984e  002 9.113406e 002 1.039664e  001

0.000000eþ000 4.929363e  003 2.121046e  002 3.348623e  002 3.921548e  002 4.115983e  002 4.280265e  002 4.640587e  002 5.292397e  002 6.194720e  002 7.075515e  002 7.676188e  002 8.002492e  002 8.294663e  002 8.787686e 002 9.803009e  002

0.000000e þ000 1.105437e 002 3.652967e  002 4.431229e  002 3.956246e  002 4.018692e  002 4.700923e  002 5.170772e  002 5.461558e  002 6.212952e  002 7.193451e  002 7.848981e  002 8.185022e  002 8.625674e  002 9.376719e  002 1.005838e 001

1

Table 6 The ARE error for different time levels and different number of nodal points for Example 3 with Δt ¼ 0:1 by SMRPI.

0.8

u(x,0.5,0.5)

0.6 0.4

t

N ¼ 125 ðh ¼ 0:25Þ

N ¼ 216 ðh ¼ 0:2Þ

N ¼ 1331 ðh ¼ 0:1Þ

0.2

0.000000 0.200000 0.400000 0.600000 0.800000 1.000000 1.200000 1.400000 1.600000 1.800000 2.000000 2.200000 2.400000 2.600000 2.800000 3.000000

0.000000eþ 000 3.692590e  003 1.421025e  002 2.239623e  002 2.629020e  002 2.519844e  002 2.209853e  002 1.997922e  002 2.075454e  002 2.352843e  002 2.549821e  002 2.538855e  002 2.361819e  002 2.149194e  002 2.077941e  002 2.611472e  002

0.000000eþ000 2.357024e  003 6.491011e 003 8.176596e  003 8.460333e  003 8.343195e  003 7.896826e  003 6.956851e  003 6.921330e  003 7.922469e  003 8.370498e  003 7.999323e  003 7.698898e  003 7.537323e  003 6.649593e 003 7.188876e  003

0.000000e þ000 1.846129e 003 1.057290e  003 5.269951e 004 5.460865e  004 4.976281e 004 5.447249e  004 6.719116e  004 6.951108e  004 6.485115e  004 5.765768e  004 5.823181e 004 8.009381e 004 1.071039e  003 1.168595e  003 8.560389e  004

0 −0.2 −0.4 −0.6 −0.8 −1

0

0.2

0.4

0.6

0.8

1

x Fig. 8. Numerical behavior of uðx; 0:5; 0:5Þ by passing the time (t ¼ 1:6; 1:7; …; 4:7) with h ¼0.1 and time step Δt ¼ 0:1 for Example 3.

E. Shivanian / Ocean Engineering 89 (2014) 173–188

time=0.4s

1

1

time=0.8s 2.5

1 0.8

0.8 0.6

0.8

1.5

z

z 0.4

0.4

0

0

1

0.5 x

1

0 1

0

0.5 y

0

0

0.5 x

1

0.8

10

0.8

0.6

8

0.6

0.4

6

0.4

0.2

4

0.2

0 1

2

0 1

1

0

1 0.5 y

1

0.5 x

0

25 20 15 10

y

5 0.5

0

0

0.5 x

0

0

0.5 x

1

0

time=2.0s

z

z

12

0

0 1

time=1.6s

1

0.5

2

0.2

0.5

time=1.2s

y

3

1

0

1

60

0.8

50

0.6

40

0.4

30

0.2

20

0 1

10

z

0.5 y

4

0.6 0.4

0.2 0.2

5

0.8

0.4

0.2 0 1

1

2

0.6

0.6

6

z

time=0.0s

185

0.5 y

0

0

0.5

1

x

Fig. 9. Initial condition and numerical contours of Example 4 by MLRPI, at times t¼ 0, 0.4, 0.8, 1.2, 1.6 and 2 with Δt ¼ 0:1 and h ¼ 0.1 (N ¼ 1331). 60

Table 7 The REI error for different time levels and different number of nodal points for Example 4 with Δt ¼ 0:001 by MLRPI.

50

t

N ¼ 125 ðh ¼ 0:25Þ

N ¼ 216 ðh ¼ 0:2Þ

N ¼ 1331 ðh ¼ 0:1Þ

0.000000 0.100000 0.200000 0.300000 0.400000 0.500000 0.600000 0.700000 0.800000 0.900000 1.000000 1.100000 1.200000 1.300000 1.400000 1.500000 1.600000 1.700000 1.800000 1.900000 2.000000

0.000000eþ 000 1.387420e  002 1.999864e  002 1.899463e  002 1.843618e  002 1.849459e  002 1.854540e  002 1.854328e  002 1.853856e  002 1.853851e  002 1.853895e  002 1.853897e  002 1.853893e  002 1.853893e  002 1.853893e  002 1.853893e  002 1.853893e  002 1.853893e  002 1.853893e  002 1.853893e  002 1.853893e  002

0.000000eþ 000 5.703190e  003 8.990607e  003 8.981378e  003 8.692393e  003 8.652422e 003 8.671244e  003 8.676840e 003 8.675997e  003 8.675487e  003 8.675505e 003 8.675542e 003 8.675543e 003 8.675540e 003 8.675540e 003 8.675540e 003 8.675540e 003 8.675540e 003 8.675540e 003 8.675540e 003 8.675540e 003

0.000000e þ000 6.506800e  003 1.012395e  002 9.986951e  003 9.689304e  003 9.659557e  003 9.674535e  003 9.680928e  003 9.680863e  003 9.680039e  003 9.679954e  003 9.680049e  003 9.680048e  003 9.680047e  003 9.680044e  003 9.680045e  003 9.680045e  003 9.680045e  003 9.680045e  003 9.680045e  003 9.680045e  003

u(x,0.5,0.5)

40 30 20 10 0 −10

0

0.2

0.4

0.6

0.8

1

x Fig. 10. Numerical behavior of uðx; 0:5; 0:5Þ by passing the time with h ¼ 0.1 and time step Δt ¼ 0:001 for Example 4.

(i) Non-regularly distributed nodes: Let us choose the set of nodes for interior domain to be Halton points in ½0; 13 . Furthermore, it is used regularly distributed nodes for the boundary of the domain. We use the following commands in MATLAB to generate such sets:

%to generate boundary nodes (152 nodes) [x,y,z] ¼meshgrid(0:0.2:1,0:0.2:1,0:0.2:1); xx ¼x(:);yy ¼y(:);zz¼z(:);

186

E. Shivanian / Ocean Engineering 89 (2014) 173–188

Table 8 The REI error for different time levels and different number of nodal points for Example 4 with Δt ¼ 0:001 by SMRPI. t

N ¼ 125 ðh ¼ 0:25Þ

N ¼ 216 ðh ¼ 0:2Þ

N ¼ 1331 ðh ¼ 0:1Þ

0.000000 0.200000 0.400000 0.600000 0.800000 1.000000 1.200000 1.400000 1.600000 1.800000 2.000000 2.200000 2.400000 2.600000 2.800000 3.000000 3.200000 3.400000 3.600000 3.800000 4.000000

0.000000eþ 000 6.823826e  003 6.264095e  003 6.300513e  003 6.299577e  003 6.299413e  003 6.299443e  003 6.299440e  003 6.299441e  003 6.299441e  003 6.299441e  003 6.299441e  003 6.299441e  003 6.299441e  003 6.299441e  003 6.299441e  003 6.299441e  003 6.299441e  003 6.299441e  003 6.299441e  003 6.299441e  003

0.000000eþ000 3.920191e  003 3.811504e  003 3.800478e  003 3.802276e  003 3.802157e  003 3.802159e  003 3.802159e  003 3.802159e  003 3.802159e  003 3.802159e  003 3.802159e  003 3.802159e  003 3.802159e  003 3.802159e  003 3.802159e  003 3.802159e  003 3.802159e  003 3.802159e  003 3.802159e  003 3.802159e  003

0.000000eþ 000 2.632109e  003 2.685450e  003 2.657973e  003 2.660674e  003 2.660585e  003 2.660574e  003 2.660576e  003 2.660576e  003 2.660576e  003 2.660576e  003 2.660576e  003 2.660576e  003 2.660576e  003 2.660576e  003 2.660576e  003 2.660576e  003 2.660576e  003 2.660576e  003 2.660576e  003 2.660576e  003

Table 9 The ARE and REI errors for different time levels and different number of nodal points for Example 4 with Δt ¼ 0:001 by SMRPI method. t

ARE Error

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0

REI Error

N ¼252

N ¼352

N¼ 452

N ¼ 252

N ¼ 352

N ¼ 452

0.000000eþ 000 3.794128e 003 3.594094e 003 3.591029e 003 3.591293e  003 3.591339e  003 3.591334e  003 3.591334e  003 3.591334e  003 3.591334e  003 3.591334e  003 3.591334e  003 3.591334e  003 3.591334e  003 3.591334e  003 3.591334e  003

0.000000eþ 000 2.938477e  003 2.882871e  003 2.869468e  003 2.871015e 003 2.870960e  003 2.870954e  003 2.870955e  003 2.870955e  003 2.870955e  003 2.870955e  003 2.870955e  003 2.870955e  003 2.870955e  003 2.870955e  003 2.870955e  003

0.000000eþ 000 2.615877e  003 2.631230e 003 2.609978e  003 2.612181e  003 2.612118e  003 2.612108e 003 2.612109e 003 2.612109e 003 2.612109e 003 2.612109e 003 2.612109e 003 2.612109e 003 2.612109e 003 2.612109e 003 2.612109e 003

0.000000eþ000 4.410879e  003 4.188327e  003 4.194650e  003 4.195165e  003 4.195079e  003 4.195079e  003 4.195080e  003 4.195080e  003 4.195080e  003 4.195080e  003 4.195080e  003 4.195080e  003 4.195080e  003 4.195080e  003 4.195080e  003

0.000000e þ000 3.707623e  003 3.644006e 003 3.631960e  003 3.633117e  003 3.633136e 003 3.633123e  003 3.633122e  003 3.633122e  003 3.633122e  003 3.633122e  003 3.633122e  003 3.633122e  003 3.633122e  003 3.633122e  003 3.633122e  003

0.000000e þ000 2.645540e  003 2.708867e 003 2.691070e 003 2.692250e  003 2.692165e  003 2.692167e 003 2.692170e  003 2.692169e  003 2.692169e  003 2.692169e  003 2.692169e  003 2.692169e  003 2.692169e  003 2.692169e  003 2.692169e  003

N=252

N=452

N=352

1

1

1

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0 1

1

0.5 0 0

0.5

0 1 0.5 0 0

0.5

0 1

1

0.5 3

Fig. 11. Non-regularly distributed nodes covering the domain ½0; 1 for Example 4.

0 0

0.5

1

E. Shivanian / Ocean Engineering 89 (2014) 173–188

a ¼ findð xx ¼ ¼ 1jyy ¼ ¼ 1jzz ¼ ¼ 1j xx ¼ ¼ 0jyy ¼ ¼ 0jzz ¼ ¼ 0Þ ; %to generate interior nodes p ¼haltonset(3,’Skip’,1e3,’Leap’,1e2); XX ¼net(p,100);YY ¼net(p,200);ZZ¼ net(p,300); %to show N ¼100 þ152 nodes subplot(1,3,1), scatter3([XX(:,1);xx(a)],[XX(:,2);yy(a)],[XX (:,3);zz(a)],’ro’,’filled’), axis tight, title(’N ¼252’) set(gca,’box’,’on’,’ticklength’,[0.01 0.01]) %to show N ¼200 þ152 nodes subplot(1,3,2), scatter3([YY(:,1);xx(a)],[YY(:,2);yy(a)],[YY (:,3);zz(a)],’ro’,’filled’), axis tight, title(’N ¼352’) set(gca,’box’,’on’,’ticklength’,[0.01 0.01]) %to show N ¼300 þ152 nodes subplot(1,3,3), scatter3([ZZ(:,1);xx(a)],[ZZ(:,2);yy(a)],[ZZ (:,3);zz(a)],’ro’,’filled’), axis tight, title(’N ¼452’) set(gca,’box’,’on’,’ticklength’,[0.01 0.01])

The results are depicted in Fig. 11. In this case, we have reported ARE and REI errors of numerical solutions obtained by the SMRPI method for different time levels and different number of nodal points in Table 9.

9. Conclusions In this paper meshless local radial point interpolation (MLRPI) method has been applied to solve a class of three space dimensional nonlinear wave equations. The main drawback of almost all methods in fully 3-dimensional problems is the large computational cost. The main challenges of current work are the improvement of meshless techniques, which relies on the robust and efficient implementation of effective integration rules, to the high-dimensional nonlinear waves so that the time derivative is a part of nonlinear term. However, it has been shown that traditional MLRPI method is not an effective technique to solve these kinds of the problems so that the results respect to the nodal points are divergence. Then, a reliable and innovative method namely spectral meshless radial point interpolation (SMRPI) has been proposed. The present method does not need any element or mesh as well as any integration. In this method, the shape functions have been constructed by the radial point interpolation method to construct basis functions in SMRPI. Some time stepping schemes were employed to approximate the time derivatives. In order to eliminate the nonlinearity, a simple predictor–corrector scheme combined with one-step time discretization and Crank–Nicolson technique has been adopted.

Acknowledgments The author is very grateful to two anonymous reviewers for carefully reading the paper and for their comments and suggestions which have improved the paper very much. References Abbasbandy, S., Shirzadi, A., 2010. A meshless method for two-dimensional diffusion equation with an integral condition. Eng. Anal. Boundary Elem. 34 (12), 1031–1037.

187

Abbasbandy, S., Shirzadi, A., 2011. MLPG method for two-dimensional diffusion equation with Neumann's and non-classical boundary conditions. Appl. Numer. Math. 61, 170–180. Atluri, S., Zhu, T., 1998a. A new meshless local Petrov–Galerkin (MLPG) approach in computational mechanics. Comput. Mech. 22, 117–127. Atluri, S., Zhu, T., 1998b. A new meshless local Petrov–Galerkin (MLPG) approach to nonlinear problems in computer modeling and simulation. Comput. Model. Simul. Eng. 3 (3), 187–196. Atluri, S., Zhu, T., 2000a. New concepts in meshless methods. Int. J. Numer. Meth. Eng. 13, 537–556. Atluri, S., Zhu, T., 2000b. The meshless local Petrov–Galerkin (MLPG) approach for solving problems in elasto-statics. Comput. Mech. 25, 169–179. Avilez-Valente, P., Seabra-Santos, F., 2004. A Petrov–Galerkin finite element scheme for the regularized long wave equation. Comput. Mech. 34, 256–270. Belytschko, T., Lu, Y.Y., Gu, L., 1994. Element-free Galerkin methods. Int.J. Numer. Methods Eng. 37 (2), 229–256. Belytschko, T., Lu, Y.Y., Gu, L., 1995. Element free Galerkin methods for static and dynamic fracture. Int. J. Solids Struct. 32, 2547–2570. Bouillard, P., Lacroix, V., DeBel, E., 2004. A wave-oriented meshless formulation for acoustical and vibro-acoustical applications. Wave Motion 39, 295–305. Bratsos, A., 2008. An improved numerical scheme for the sine-Gordon equation in 2þ1 dimensions. Int. J. Numer. Meth. Eng. 75, 787–799. Carlos, J., Alves, S., Valtchev, S., 2005. Numerical comparison of two meshfree methods for acoustic wave scattering. Eng. Anal. Boundary Elem. 29, 371–382. Clear, P., 1998. Modeling conned multi-material heat and mass flows using SPH. Appl. Math. Model. 22, 981–993. Dai, W., Song, H., Su, S., Nassar, R., 2006. A stable finite difference scheme for solving a hyperbolic two-step model in a 3D micro sphere exposed to ultrashort-pulsed lasers. Int. J. Numer. Meth. H 16, 693–717. De, S., Bathe, K., 2000. The method of finite spheres. Comput. Mech. 25, 329–345. Dehghan, M., 2005. On the solution of an initial-boundary value problem that combines Neumann and integral condition for the wave equation. Numer. Meth. Partial Diff. Eq. 21, 24–40. Dehghan, M., 2006a. Finite difference procedures for solving a problem arising in modeling and design of certain optoelectronic devices. Math. Comput. Simul. 71, 16–30. Dehghan, M., 2006b. A computational study of the one-dimensional parabolic equation subject to nonclassical boundary specifications. Numer. Meth. Partial Diff. Eq. 22, 220–257. Dehghan, M., Ghesmati, A., 2010a. Combination of meshless local weak and strong (MLWS) forms to solve the two dimensional hyperbolic telegraph equation. Eng. Anal. Boundary Elem. 34, 324–336. Dehghan, M., Ghesmati, A., 2010b. Numerical simulation of two-dimensional sineGordon solitons via a local weak meshless technique based on the radial point interpolation method (RPIM). Comput. Phys. Commun. 181, 772–786. Dehghan, M., Mirzaei, D., 2008. The meshless local Petrov–Galerkin (MLPG) method for the generalized two-dimensional non-linear Schrödinger equation. Eng. Anal. Boundary Elem. 32, 747–756. Dehghan, M., Mirzaei, D., 2009. Meshless local Petrov–Galerkin (MLPG) method for the unsteady magnetohydrodynamic (MHD) flow through pipe with arbitrary wall conductivity. Appl. Numer. Math. 59, 1043–1058. Dehghan, M., Mohebbi, A., 2008. The combination of collocation, finite difference, and multigrid methods for solution of the two-dimensional wave equation. Numer. Methods Partial Diff. Eq. 24, 897–910. Dehghan, M., Salehi, R., 2012a. A method based on meshless approach for the numerical solution of the two-space dimensional hyperbolic telegraph equation. Math. Method Appl. Sci. 35, 1120–1233. Dehghan, M., Salehi, R., 2012b. A meshless based numerical technique for traveling solitary wave solution of Boussinesq equation. Appl. Math. Model. 36, 1939–1956. Dehghan, M., Shokri, A., 2008. A numerical method for solution of the two dimensional sine-Gordon equation using the radial basis functions. Math. Comput. Simul. 79, 700–715. Dehghan, M., Shokri, A., 2009. Numerical solution of the nonlinear Klein–Gordon equation using radial basis functions. J. Comput. Appl. Math 230, 400–410. Franke, C., Schaback, R., 1997. Solving partial differential equations by collocation using radial basis functions. Appl. Math. Comput. 93, 73–82. Gamallo, P., Astley, R., 2006. The partition of unity finite element method for short wave acoustic propagation on non-uniform potential flows. Int. J. Numer. Meth. Eng. 65, 425–444. Grilli, S.T., Guyenneb, P., Diasc, F., 2001. A fully non-linear model for threedimensional overturning waves over an arbitrary bottom. Int. J. Numer. Meth. Fluids 35, 829–867. Gu, Y., Liu, G., 2001. A meshless local Petrov–Galerkin (MLPG) method for free and forced vibration analyses for solids. Comput. Mech. 27, 188–198. Gu, Y., Liu, G., 2002. A boundary point interpolation method for stress analysis of solids. Comput. Mech. 28, 47–54. Gu, Y., Liu, G., 2003. A boundary radial point interpolation method (BRPIM) for 2-d structural analyses. Struct. Eng. Mech. 15, 535–550. Han, Z., Atluri, S., 2004. A meshless local Petrov–Galerkin (MLPG) approach for 3-dimensional elasto-dynamics. CMC: Comput. Mater. Continua 1 (2), 129–140. Hirose, A., 1985. Introduction to Wave Phenomena. Wile, New York. Hu, D., Long, S., Liu, K., Li, G., 2006. A modified meshless local Petrov–Galerkin method to elasticity problems in computer modeling and simulation. Eng. Anal. Boundary Elem. 30, 399–404.

188

E. Shivanian / Ocean Engineering 89 (2014) 173–188

Kamruzzaman, M., 2008. A new meshless collocation method for partial differential equations. Commun. Numer. Meth. Eng. 24, 1617–1639. Kansa, E., 1990. Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics. I. Surface approximations and partial derivative estimates. Comput. Math. Appl. 19 (8–9), 127–145. Libre, N., Emdadi, A., Kansa, E., Shekarchi, M., Rahimian, M., 2008. A fast adaptive wavelet scheme in RBF collocation for nearly singular potential PDEs. CMESComput. Modeling Eng. Sci. 38 (3), 263–284. Libre, N., Emdadi, A., Kansa, E., Shekarchi, M., Rahimian, M., 2009. A multiresolution prewavelet-based adaptive refinement scheme for RBF approximations of nearly singular problems. Eng. Anal. Boundary Elem. 33 (7), 901–914. Liew, K., Cheng, R., 2013. Numerical study of the three-dimensional wave equation using the mesh-free kp-Ritz method. Eng. Anal. Boundary. Elem. 37, 977–989. Ling, L., Schaback, R., 2008. Stable and convergent unsymmetric meshless collocation methods. SIAM. J. Numer. Anal. 46 (3), 1097–1115. Ling, L., Opfer, R., Schaback, R., 2006. Results on meshless collocation techniques. Eng. Anal. Boundary Elem. 30 (4), 247–253. Liu, G., Gu, Y., 2001. A local radial point interpolation method (LR-PIM) for free vibration analyses of 2-D solids. J. Sound Vib. 246 (1), 29–46. Liu, G., Gu, Y., 2005. An Introduction to Meshfree Methods and Their Programming. Springer, Berlin. Liu, W., Jun, S., Zhang, Y., 1995. Reproducing kernel particle methods. Int. J. Numer. Meth. Eng. 20, 1081–1106. Liu, G., Yan, L., Wang, J., Gu, Y., 2002. Point interpolation method based on local residual formulation using radial basis functions. Struct. Eng. Mech. 14, 713–732. Liu, K., Long, S., Li, G., 2006. A simple and less-costly meshless local Petrov–Galerkin (MLPG) method for the dynamic fracture problem. Eng. Anal. Boundary Elem. 30, 72–76. Lu, Y., Belytschko, T., Tabbara, M., 1995. Element-free Galerkin method for wave propagation and dynamic fracture. Comput. Meth. Appl. Mech. Eng. 26, 131–153. Melenk, J., Babuska, I., 1996. The partition of unity finite element method: basic theory and applications. Comput. Meth. Appl. Mech. Eng. 139, 289–314. Mohanty, R., Gopal, V., 2013. A new off-step high order approximation for the solution of three-space dimensional nonlinear wave equations. Appl. Math. Model. 37, 2802–2815.

Mukherjee, Y., Mukherjee, S., 1997. Boundary node method for potential problems. Int. J. Numer. Meth. Eng. 40, 797–815. Nayroles, B., Touzot, G., Villon, P., 1992. Generalizing the finite element method: diffuse approximation and diffuse elements. Comput. Mech. 10, 307–318. Nettel, S., 2003. Wave Physics: Oscillations-Solitons-Chaos. Springer, Berlin. Powell, M., 1992. Theory of radial basis function approximation in 1990. In: Light, F.W. (Ed.), Advanced Numerical Analysis, pp. 303–322. Selvadurai, A., 2000. Partial Differential Equations in Mechanics. Springer, Berlin. Senturk, U., 2011. Modeling nonlinear waves in a numerical wave tank with localized meshless RBF method. Comput. Fluids 44, 221–228. Sharan, M., Kansa, E., Gupta, S., 1997. Application of the multiquadric method for numerical solution of elliptic partial differential equations. Appl. Math. Comput. 84, 275–302. Shirzadi, A., Ling, L., Abbasbandy, S., 2012. Meshless simulations of the twodimensional fractional-time convection-diffusion-reaction equations. Eng. Anal. Boundary Elem. 36, 1522–1527. Shirzadi, A., Sladek, V., Sladek, J., 2013. A local integral equation formulation to solve coupled nonlinear reaction-diffusion equations by using moving least square approximation. Eng. Anal. Boundary Elem. 37, 8–14. Wang, J., Liu, G., 2002a. A point interpolation meshless method based on radial basis functions. Int. J. Numer. Meth. Eng. 54, 1623–1648. Wang, J., Liu, G., 2002b. On the optimal shape parameters of radial basis functions used for 2-D meshless methods. Comput. Meth. Appl. Math. 191, 2611–2630. Wazwaz, A., 2002. Partial Differential Equations: Methods and Applications. Swets & Zeitlinger, B.V., Lisse, Netherlands. Wen, P., 2010. Meshless local Petrov–Galerkin (MLPG) method for wave propagation in 3D poroelastic solids. Eng. Anal. Boundary Elem. 34, 315–323. Wendland, H., 1998. Error estimates for interpolation by compactly supported radial basis functions of minimal degree. J. Approx. Theory 93, 258–396. Xiong, Z., 2008. Convergence of FEM with interpolated coefficients for semilinear hyperbolic equation. J. Comput. Appl. Math. 24, 313–317. Yan, Z., Konotop, V.V., 2009. Exact solutions to three-dimensional generalized nonlinear schrdinger equations with varying potential and nonlinearities. Phys. Rev. E 80, 036607:1–036607:9. Zhang, Z., Li, D., Cheng, Y., Liew, K., 2012. The improved element-free Galerkin method for three-dimensional wave equation. Acta Mech. Sin. 28 (3), 808–818.