An efficient numerical technique for solution of two-dimensional cubic nonlinear Schrödinger equation with error analysis

Engineering Analysis with Boundary Elements 83 (2017) 74–86

Contents lists available at ScienceDirect

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

An efficient numerical technique for solution of two-dimensional cubic nonlinear Schrödinger equation with error analysis Elyas Shivanian∗, Ahmad Jafarabadi Department of Applied Mathematics, Imam Khomeini International University, Qazvin, 34149-16818, Iran

a r t i c l e

i n f o

Keywords: Spectral meshless radial point interpolation (SMRPI) method Radial basis function Cubic nonlinear Schrödinger equation

a b s t r a c t In this paper, the spectral meshless radial point interpolation (SMRPI) technique is applied to the solution of twodimensional cubic nonlinear Schrödinger equations. Firstly, we obtain a time discrete scheme by approximating time derivative via a finite difference formula, then we use the SMRPI approach to approximate the spatial derivatives. This method is based on a 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 act as basis functions in frame of SMRPI. In the current work, the thin plate splines (TPS) are used as the basis functions and in order to eliminate the nonlinearity, a simple predictor-corrector (P-C) scheme is performed. We prove that the time discrete scheme is unconditionally stable and convergent in time variable using the energy method. We show that convergence order of the time discrete scheme is (𝛿𝑡). The aim of this paper is to show that the SMRPI method is suitable for the treatment of the nonlinear Schrödinger equations. Also, the SMRPI has less computational complexity than the other methods that have already solved this problem. The results of numerical experiments are compared with analytical solution to confirm the accuracy and efficiency of the presented scheme. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction The present paper considers the following generalized cubic nonlinear two-dimensional Schrödinger equation with wave function 𝑢(x, 𝑡): 𝜕𝑢(x, 𝑡) + 𝛼Δ𝑢(x, 𝑡) + 𝛽|𝑢(x, 𝑡)|2 𝑢(x, 𝑡) + 𝑤(x)𝑢(x, 𝑡) = 0, 𝜕𝑡 x ∈ Ω, 𝑡 ∈ (0, 𝑇 ], 𝑖

(1)

with initial condition 𝑢(x, 0) = 𝑢0 (x),

x ∈ Ω,

(2)

and Dirichlet boundary condition 𝑢(x, 𝑡) = ℎ(x, 𝑡),

x ∈ 𝜕Ω, 𝑡 > 0,

(3) √ 2 where 𝑢(x, 𝑡) is an unknown complex function, 𝑖 = −1, Ω ⊂ ℝ , 𝜕Ω is the boundary of Ω, 𝑢0 (x) and ℎ(x, 𝑡) are given sufficiently smooth functions and 𝑤(x) is an arbitrary potential function that is real value and bounded on Ω. The cubic nonlinear Schrödinger equation occurs in a variety of areas, including, quantum mechanics [1], nonlinear optics [2,3], electromagnetic wave propagation [4]. In optics, the nonlinear Schrödinger equation occurs in the Manakov system, a model of wave propagation in fiber optics. The function u ∗

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

http://dx.doi.org/10.1016/j.enganabound.2017.07.012 Received 27 October 2016; Received in revised form 15 June 2017; Accepted 17 July 2017 0955-7997/© 2017 Elsevier Ltd. All rights reserved.

represents a wave and the nonlinear Schrödinger equation describes the propagation of the wave through a nonlinear medium. The second-order derivative represents the dispersion, while the 𝛽 term represents the nonlinearity. The equation models many nonlinearity effects in a fiber, including but not limited to self-phase modulation, four-wave mixing, second harmonic generation, stimulated Raman scattering, etc [5–9]. Since the solutions of the Schrödinger equation is very important for describing several phenomena in physics and engineering, therefore solving this equation is necessary. As a short list that numerically investigated the solutions of various aspects of the Schrödinger equation, we refer to some of them. In [10–13], they have applied various types of the finite difference schemes, in [14–20] they have employed some meshfree methods and in [21–23] they have used the alternating direction implicit (ADI) schemes for solving several models of linear and nonlinear Schrödinger equations. As references that have investigated analytical solutions, can be cited to [24,25]. In [24], the author has used an effective method namely the improved tan (Φ(𝜉)/2)-expansion method for constructing a range of exact solutions for the following partial differential equation 𝑖𝑢𝑡 + 𝑢𝑥𝑥 ± 2𝛾|𝑢|2 𝑢 = 0,

(4)

where 𝛾 is a non-zero real constants and 𝑢 = 𝑢(𝑥, 𝑡) is a complex-valued function of two real variables x, t. The authors of [25] have used the

E. Shivanian, A. Jafarabadi

Engineering Analysis with Boundary Elements 83 (2017) 74–86

𝑖𝜆𝑢𝑘+1 (x) + 𝛼Δ𝑢𝑘+1 (x) + 𝑤(x)𝑢𝑘+1 (x)

𝐺′ -expansion 𝐺

method for finding the exact solutions of the five complex nonlinear Schrödinger equations via the nonlinear Schrödinger equations, Eq. (4), the unstable Schrödinger equation 𝑖𝑢𝑡 + 𝑢𝑥𝑥 + 2𝜆|𝑢|2 𝑢 − 2𝛾𝑢 = 0,

= 𝑖𝜆𝑢𝑘 (x) − 𝛼Δ𝑢𝑘 (x) − 2𝛽𝑓 (𝑢𝑘 ) − 𝑤(x)𝑢𝑘 (x)−2𝑖𝑅𝑘+1 , or

(5)

𝑢𝑘+1 (x) − 𝑖𝜆−1 𝛼Δ𝑢𝑘+1 (x) − 𝑖𝜆−1 𝑤(x)𝑢𝑘+1 (x)

and two-dimensional 𝑖𝑢𝑡 + 𝑢𝑥𝑥 + 𝑢𝑦𝑦 + 𝛾|𝑢|2 𝑢 = 0,

= 𝑢𝑘 (x) + 𝑖𝜆−1 𝛼Δ𝑢𝑘 (x) + 2𝑖𝜆−1 𝛽𝑓 (𝑢𝑘 ) + 𝑖𝜆−1 𝑤(x)𝑢𝑘 (x) + 𝑅𝑘∗+1 ,

(6)

(7)

It is well known some traditional methods such as the finite element method (FEM), the finite volume method (FVM) and the boundary element method (BEM) are based on meshes or elements, and this is the main deficiency of these numerical techniques. In the last two decades, in order to overcome the mentioned shortcoming some methods so-called meshless methods have been proposed. There are various types of meshless techniques, for example, meshless techniques based on weak forms such as the element-free Galerkin (EFG) method [26–34], diffuse element method [35], meshless local radial point interpolation method [36–38], meshless local Petrov–Galerkin method [39–42] and including their developments; meshless techniques based on collocation techniques (strong forms) such as the meshless collocation technique based on radial basis functions (RBFs) and finally meshless techniques based on the combination of weak forms and collocation technique [43–48]. Shivanian [49,50] proposed SMRPI method which is based on meshless methods and benefits from spectral collocation ideas. In SMRPI technique, the point interpolation method with the help of those radial basis functions, which were free of shape parameter, has been proposed to construct shape 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.

𝑈 𝑘+1 (x) − 𝑖𝜆−1 𝛼Δ𝑈 𝑘+1 (x) − 𝑖𝜆−1 𝑤(x)𝑈 𝑘+1 (x) = 𝑈 𝑘 (x) + 𝑖𝜆−1 𝛼Δ𝑈 𝑘 (x) + 2𝑖𝜆−1 𝛽𝑓 (𝑈 𝑘 ) + 𝑖𝜆−1 𝑤(x)𝑈 𝑘 (x).

2.1. The stability and convergence analysis Here, we examine the stability and convergence for Eqs. (1)–(3) using the proposed time discrete scheme. To introduce the variational formulation of the Eq. (10), we define some functional spaces endowed with standard norms and inner products that will be used hereafter. Let Ω denote a bounded and open domain in ℝ𝑑 , for 𝑑 = 2 or 3 and let dx be the Lebesgue measure on ℝ𝑑 . For p < ∞, we denote by Lp (Ω) the space of the measurable functions 𝑢 ∶ Ω → ℂ such that ∫ Ω |u(x)|p dx < ∞ that is a Banach space with the standard norm ( )1∕𝑝 ‖𝑢‖𝐿𝑝 (Ω) =

∫Ω

|𝑢(𝑥)|𝑝 𝑑𝑥

.

In particle case, let L2 (Ω) be the set of all measurable functions 𝑢 ∶ Ω → ℂ such that ∫ Ω |u(x)|2 dx < ∞. From the inequality 𝑎𝑏 ≤ 12 (𝑎2 + 𝑏2 ), valid

for all a, b ≥ 0, we see that if u, v ∈ L2 (Ω) then |𝑢𝑣̄ | ≤ 12 (|𝑢|2 + |𝑣|2 ), so that 𝑢𝑣̄ ∈ 𝐿1 (Ω). It follows easily that the formula [51]

Our aim of this paper is the development of spectral meshless radial point interpolation to obtain the solution of two-dimensional cubic nonlinear Schrödinger equation. We will show that the SMRPI method is suitable for the treatment of the nonlinear Schrödinger equations. Also, the SMRPI has less computational complexity than the other methods that have already solved this problem, e.g. [18]. The outline of this paper is as follows: In Section 2, we obtain a time discrete scheme to the mentioned equation. In this section, we prove the unconditional stability and convergence of the time discrete scheme using the energy method. In Section 3, we introduce the spectral meshless radial point interpolation scheme briefly so that the high order operational matrices are obtained. The numerical implementation of SMRPI is given in Section 4. In Section 5, we report the numerical experiments of solving Eq. (1) for three test problems. Finally a conclusion is given in Section 6.

⟨𝑢, 𝑣⟩ =

2. Time discretization

𝑒𝑘+1 (x) − 𝑖𝜆−1 𝛼Δ𝑒𝑘+1 (x) − 𝑖𝜆−1 𝑤(x)𝑒𝑘+1 (x)

∫Ω

𝑢(𝑥)𝑣(𝑥)𝑑𝑥

defines an inner product on L2 (Ω). Now, above inner product induces the norm ( )1∕2 ‖𝑢‖𝐿2 (Ω) =

∫Ω

|𝑢(𝑥)|2 𝑑𝑥

.

Theorem 1. Let Uk ∈ L2 (Ω) and we assume that the solution of Eq. (1) is ̄ × [0, 𝑇 ]. Also, suppose that f(·) has the Lipschitz condition, analytical over Ω i.e. |𝑓 (𝑣1 ) − 𝑓 (𝑣2 )| ≤ |𝑣1 − 𝑣2 |,

∀𝑣1 , 𝑣2 ∈ 𝐿2 (Ω).

Then the scheme defined by (11) is unconditionally stable in L2 (Ω). Proof. The roundoff error is to the following form

= 𝑒𝑘 (x) + 𝑖𝜆−1 𝛼Δ𝑒𝑘 (x) + 2𝑖𝜆−1 𝛽(𝑓 (𝑈 𝑘 ) − 𝑓 (𝑈̂ 𝑘 )) + 𝑖𝜆−1 𝑤(x)𝑒𝑘 (x), (12)

In this section, we discretize the time variable using forward finite difference relation for approximating the first-order derivative on time variable. If we consider Eq. (1) in point (x, 𝑡𝑘+1 ), then we have ) 𝑢𝑘+1 (x) − 𝑢𝑘 (x) 𝛼 ( 𝑘+1 + Δ𝑢 (x) + Δ𝑢𝑘 (x) + 𝛽|𝑢𝑘 (x)|𝑢𝑘 (x) 𝛿𝑡 2 ) 𝑤(x) ( 𝑘+1 + 𝑢 (x) + 𝑢𝑘 (x) + 𝑖𝑅𝑘+1 = 0, 2

(11)

In the above relations 𝑢𝑘 (x) and 𝑈 𝑘 (x) are exact and approximate solutions of Eq. (1), respectively.

1.1. The main aim of this paper

𝑖

(10)

where 𝜆 = 2∕𝛿𝑡, 𝑓 (𝑢) = |𝑢|2 𝑢 and |𝑅𝑘∗+1 | < 𝐶∗ 𝛿𝑡2 , for a positive constant C∗ . By omitting the truncation local error 𝑅𝑘∗+1 we obtain the following form

and three-dimensional Schrödinger equations 𝑖𝑢𝑡 + 𝑢𝑥𝑥 + 𝑢𝑦𝑦 + 𝑢𝑧𝑧 + 𝛾|𝑢|2 𝑢 = 0.

(9)

where 𝑒𝑘+1 (x) = 𝑈 𝑘+1 (x) − 𝑈̂ 𝑘+1 (x) whereas 𝑈 𝑘+1 (x) and 𝑈̂ 𝑘+1 (x) are exact and approximate solutions of Eq. (11), respectively. Multiplying Eq. (12) by 𝑒𝑘+1 (x) and integrating on Ω, we obtain ⟨𝑒𝑘+1 (x), 𝑒𝑘+1 (x)⟩ − 𝑖𝜆−1 𝛼⟨Δ𝑒𝑘+1 (x), 𝑒𝑘+1 (x)⟩ − 𝑖𝜆−1 ⟨𝑤(x)𝑒𝑘+1 (x), 𝑒𝑘+1 (x)⟩ = ⟨𝑒𝑘 (x), 𝑒𝑘+1 (x)⟩ + 𝑖𝜆−1 𝛼⟨Δ𝑒𝑘 (x), 𝑒𝑘+1 (x)⟩ + 2𝑖𝜆−1 𝛽⟨(𝑓 (𝑈 𝑘 ) − 𝑓 (𝑈̂ 𝑘 )), 𝑒𝑘+1 (x)⟩ + 𝑖𝜆−1 ⟨𝑤(x)𝑒𝑘 (x), 𝑒𝑘+1 (x)⟩.

(8)

in which 𝑢𝑘 (x) = 𝑢(x, 𝑘𝛿𝑡), 𝛿𝑡 = 𝑡𝑘+1 − 𝑡𝑘 and 𝑅𝑘+1 is truncation error of time discretization such that |𝑅𝑘+1 | < 𝐶𝛿𝑡, for a positive constant C. By simplification we have

(13)

From 𝑒𝑘+1 ∈ 𝐿20 (Ω) = {𝑢 ∈ 𝐿2 (Ω) ∶ 𝑢|𝜕Ω = 0} and using Green’s formula, we have 75

E. Shivanian, A. Jafarabadi

‖𝑒𝑘+1 (x)‖2 2

𝐿 (Ω)

Engineering Analysis with Boundary Elements 83 (2017) 74–86

+ 𝑖𝜆−1 𝛼‖∇𝑒𝑘+1 (x)‖2 2

𝐿 (Ω)

− 𝑖𝜆−1 ⟨𝑤(x)𝑒𝑘+1 (x), 𝑒𝑘+1 (x)⟩

= ⟨𝑒𝑘 (x), 𝑒𝑘+1 (x)⟩ + 𝑖𝜆−1 𝛼⟨Δ𝑒𝑘 (x), 𝑒𝑘+1 (x)⟩ + 2𝑖𝜆−1 𝛽⟨(𝑓 (𝑈 𝑘 ) − 𝑓 (𝑈̂ 𝑘 )), 𝑒𝑘+1 (x)⟩ + 𝑖𝜆−1 ⟨𝑤(x)𝑒𝑘 (x), 𝑒𝑘+1 (x)⟩.

Similar to the proof of Theorem 1, since 𝑤(x) is a real function and by definition of inner product, we conclude that the term ⟨𝑤(x)𝜉 𝑘+1 (x), 𝜉 𝑘+1 (x)⟩ is a real value. Also, taking the absolute value of both sides using the definition of absolute value in complex numbers and triangle inequality, implies

(14)

Since 𝑤(x) is a real function, we conclude that term ⟨𝑤(x)𝑒𝑘+1 (x), 𝑒𝑘+1 (x)⟩ is a real value. By taking the absolute value of both sides, the definition of absolute value in complex numbers and triangle inequality, we have ‖𝑒𝑘+1 (x)‖2 2

𝐿 (Ω)

‖𝜉 𝑘+1 (x)‖2 2

𝐿 (Ω)

+ 2𝜆−1 |𝛽||⟨(𝑓 (𝑈 𝑘 ) − 𝑓 (𝑈̂ 𝑘 )), 𝜉 𝑘+1 (x)⟩|

≤ |⟨𝑒𝑘 (x), 𝑒𝑘+1 (x)⟩| + 𝜆−1 |𝛼||⟨Δ𝑒𝑘 (x), 𝑒𝑘+1 (x)⟩|

+ 𝜆−1 |⟨𝑤(x)𝜉 𝑘 (x), 𝜉 𝑘+1 (x)⟩| + |⟨𝑅𝑘∗+1 , 𝜉 𝑘+1 (x)⟩|.

+ 2𝜆−1 |𝛽||⟨(𝑓 (𝑈 𝑘 ) − 𝑓 (𝑈̂ 𝑘 )), 𝑒𝑘+1 (x)⟩| + 𝜆−1 |⟨𝑤(x)𝑒𝑘 (x), 𝑒𝑘+1 (x)⟩|.

(15)

‖𝑒

(x)‖2 2 𝐿 (Ω)

𝑘

𝐿 (Ω)

≤ {‖𝑒 (x)‖𝐿2 (Ω) + 𝜆 |𝛼|‖Δ𝑒 (x)‖𝐿2 (Ω) 𝑘

𝑘

+ 2𝜆 |𝛽|‖𝑒 (x)‖𝐿2 (Ω) + 𝜆 1 ‖𝑒 (x)‖𝐿2 (Ω) }‖𝑒 −1

‖𝜉 𝑘+1 (x)‖2 2

𝑘

−1

−1

𝑘+1

≤ {‖𝜉 𝑘 (x)‖𝐿2 (Ω) + 𝜆−1 |𝛼|‖Δ𝜉 𝑘 (x)‖𝐿2 (Ω) + 2𝜆−1 |𝛽|‖𝜉 𝑘 (x)‖𝐿2 (Ω) + 𝜆−1 1 ‖𝜉 𝑘 (x)‖𝐿2 (Ω)

(x)‖𝐿2 (Ω) ,

+ ‖𝑅𝑘∗+1 ‖𝐿2 (Ω) }‖𝜉 𝑘+1 (x)‖𝐿2 (Ω) ,

(16)

where |𝑤(x)| ≤ 1 . From maximum modulus principle for analytic functions [52], we can set ‖Δ𝑒𝑘 (x)‖𝐿2 (Ω) ≤ 2 ‖𝑒𝑘 ‖𝐿2 (Ω) for a positive number such as 2 . Therefore, we have ( ) ‖𝑒𝑘+1 (x)‖𝐿2 (Ω) ≤ 1 + 2𝜆−1 |𝛽| + 𝜆−1 1 + 𝜆−1 |𝛼|2 ‖𝑒𝑘 (x)‖𝐿2 (Ω) ( )2 ≤ 1 + 2𝜆−1 |𝛽| + 𝜆−1 1 + 𝜆−1 |𝛼|2 ‖𝑒𝑘−1 (x)‖𝐿2 (Ω)

(23)

where |𝑤(x)| ≤ 1 . By definition of 2 in Theorem 1, we have ‖𝜉 𝑘+1 (x)‖𝐿2 (Ω) ≤ (1 + 2𝜆−1 |𝛽| + 𝜆−1 1 + 𝜆−1 |𝛼|2 )‖𝜉 𝑘 (x)‖𝐿2 (Ω) + 𝐶∗ 𝛿𝑡2 ) 2( ≤ (1 + 2𝜆−1 |𝛽| + 𝜆−1 1 + 𝜆−1 |𝛼|2 ) ‖𝜉 𝑘−1 (x)‖𝐿2 (Ω) + 2𝐶∗ 𝛿𝑡2 ≤⋯

≤⋯ ( )𝑘+1 0 ≤ 1 + 2𝜆−1 |𝛽| + 𝜆−1 1 + 𝜆−1 |𝛼|2 ‖𝑒 (x)‖𝐿2 (Ω) . (17)

𝑘+1 (

≤ (1 + 2𝜆−1 |𝛽| + 𝜆−1 1 + 𝜆−1 |𝛼|2 )

) ‖𝜉 0 (x)‖𝐿2 (Ω) + (𝑘 + 1)𝐶∗ 𝛿𝑡2 . (24)

𝑘+1

The sequence 𝑣𝑘 = (1 + 2𝜆−1 |𝛽| + 𝜆−1 1 + 𝜆−1 |𝛼|2 ) converges to exp( 𝑇2 (2|𝛽| + 1 + |𝛼|2 )) because (1 + 2𝜆−1 |𝛽| + 𝜆−1 1 + 𝜆−1 |𝛼|2 )𝑘+1 is equivalent to exp((𝑘 + 1)𝜆−1 (2|𝛽| + 1 + |𝛼|2 )) and lim𝑘→∞ (𝑘 + 1)𝛿𝑡 = 𝑇 . Then (vk )k is bounded and so ∃C∗ ∗ > 0, ∀k, vk ≤ C∗ ∗ (vk is positive for all k) and then we get ‖𝑒𝑘+1 (x)‖𝐿2 (Ω) ≤ 𝐶∗∗ ‖𝑒0 (x)‖𝐿2 (Ω) ,

Clearly, taking into account 𝜉 0 (x) = 0 and definition of 𝛿t, similar to Theorem 1, leads to ‖𝜉 𝑘+1 (x)‖𝐿2 (Ω) ≤ 𝐶∗∗ 𝐶∗ 𝑇 𝛿𝑡, which completes the proof.

(18)

which means the stability of the scheme.

(22)

Using Schwarz inequality, Lipschitz condition and boundedness of 𝑤(x), we can write

Now, using Schwarz inequality, Lipschitz condition and boundedness of 𝑤(x), we can write 𝑘+1

≤ |⟨𝜉 𝑘 (x), 𝜉 𝑘+1 (x)⟩| + 𝜆−1 |𝛼||⟨Δ𝜉 𝑘 (x), 𝜉 𝑘+1 (x)⟩|

(25) □

3. Spectral meshless radial point interpolation scheme (SMRPI)

□

Let 𝑢(x, 𝑡𝑘+1 ) = 𝑢𝑘+1 ∈ 𝐿2 (Ω) be the exact solution of Eq. (1) with initial and boundary conditions (2) and (3). Also, 𝑈 (x, 𝑡𝑘+1 ) = 𝑈 𝑘+1 ∈ 𝐿2 (Ω) be approximate solution with same initial and boundary condī × [0, 𝑇 ], 𝑢(x, 𝑡) is an analytical functions. Also, we assume that over Ω tion. Now, the following theorem holds for the convergence of time discretization:

Consider a continuous function 𝑢(x) defined in a domain Ω, which is represented by a set of field nodes. The 𝑢(x) at a point of interest x is approximated in the form of

Theorem 2. The semi time-discrete solution (11) is L2 -convergent with the convergence order (𝛿𝑡).

where 𝑅𝑖 (x) is a radial basis function (RBF), n is the number of RBFs, 𝑃𝑗 (x) is monomial in the space coordinate x, and m is the number of polynomial basis functions. Coefficients ai and bj are unknown which should be determined. In order to determine ai and bj in Eq. (26), 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 ). Coefficients ai and bj can be determined by enforcing Eq. (26) to be satisfied at these n nodes surrounding the point of interest x. Therefore, by the idea of interpolation Eq. (26) is converted to the following form:

𝑢(x) =

𝑖=1

Proof. Subtracting Eqs. (11) from (10), yields 𝜉 𝑘+1 (x) − 𝑖𝜆−1 𝛼Δ𝜉 𝑘+1 (x) − 𝑖𝜆−1 𝑤(x)𝜉 𝑘+1 (x) = 𝜉 𝑘 (x) + 𝑖𝜆−1 𝛼Δ𝜉 𝑘 (x) + 2𝑖𝜆−1 𝛽(𝑓 (𝑈 𝑘 ) − 𝑓 (𝑈̂ 𝑘 )) + 𝑖𝜆−1 𝑤(x)𝜉 𝑘 (x) + 𝑅𝑘∗+1 ,

(19)

where 𝜉 𝑘+1 (x) = 𝑢𝑘+1 (x) − 𝑈 𝑘+1 (x). Multiplying both sides of the above equation by 𝜉 𝑘+1 (x) and integrating on Ω, we obtain ⟨𝜉

𝑘+1

( x) , 𝜉

𝑘+1

𝑘

= ⟨𝜉 (x), 𝜉

(x)⟩ − 𝑖𝜆 𝛼⟨Δ𝜉

𝑘+1

−1

𝑘+1 𝑘

( x) , 𝜉

𝑘+1

(x)⟩ − 𝑖𝜆 ⟨𝑤(x)𝜉 −1

𝑘+1

(x), 𝜉

𝑘+1

(x)⟩

𝑘+1

(20)

Since 𝜉 𝑘+1 (x) ∈ 𝐿20 (Ω), using Greens’ formula implies

𝑗=1

𝑛 ∑ 𝑖=1

𝑃𝑗 (x)𝑏𝑗 = R𝑡𝑟 (x)a + P𝑡𝑟 (x)b,

𝜙𝑖 (x)𝑢𝑖 .

(26)

(27)

This is because the RPIM shape functions are created to pass through nodal values. Moreover, the shape functions are the partitions of unity, i.e.

‖𝜉 𝑘+1 (x)‖𝐿2 (Ω) + 𝑖𝜆−1 𝛼‖∇𝜉 𝑘+1 (x)‖𝐿2 (Ω) − 𝑖𝜆−1 ⟨𝑤(x)𝜉 𝑘+1 (x), 𝜉 𝑘+1 (x)⟩ = ⟨𝜉 𝑘 (x), 𝜉 𝑘+1 (x)⟩ + 𝑖𝜆−1 𝛼⟨Δ𝜉 𝑘 (x), 𝜉 𝑘+1 (x)⟩ + 2𝑖𝜆−1 𝛽⟨(𝑓 (𝑈 𝑘 ) − 𝑓 (𝑈̂ 𝑘 )), 𝜉 𝑘+1 (x)⟩ + 𝑖𝜆−1 ⟨𝑤(x)𝜉 𝑘 (x), 𝜉 𝑘+1 (x)⟩ + ⟨𝑅𝑘∗+1 , 𝜉 𝑘+1 (x)⟩.

𝑚 ∑

where 𝜙𝑖 (x)’s are called the RPIM shape functions which have the Kronecker delta function property, that is { 1 for 𝑖 = 𝑗, 𝑗 = 1, 2, … , 𝑛 , 𝜙𝑖 (x𝑗 ) = (28) 0 for 𝑖 ≠ 𝑗, 𝑖, 𝑗 = 1, 2, … , 𝑛.

−1

+ 𝑖𝜆−1 ⟨𝑤(x)𝜉 𝑘 (x), 𝜉 𝑘+1 (x)⟩ + ⟨𝑅𝑘∗+1 , 𝜉 𝑘+1 (x)⟩.

𝑅 𝑖 ( x) 𝑎 𝑖 +

𝑢(x) = 𝚽𝑡𝑟 (x)U𝑠 =

(x)⟩ + 𝑖𝜆 𝛼⟨Δ𝜉 (x), 𝜉 (x)⟩ + 2𝑖𝜆 𝛽⟨(𝑓 (𝑈 𝑘 ) − 𝑓 (𝑈̂ 𝑘 )), 𝜉 𝑘+1 (x)⟩ −1

𝑛 ∑

𝑛 ∑

(21)

𝑖=1

76

𝜙𝑖 (x) = 1.

(29)

E. Shivanian, A. Jafarabadi

Engineering Analysis with Boundary Elements 83 (2017) 74–86

4. Discretization and numerical implementation for SMRPI method

For more details about radial point interpolation (RPI) shape functions and the way they are constructed, one can see [49]. Now, we construct operational matrices that are essential tools of present approach. Operational matrices make the technique more appropriate to handle partial differential equations with high derivatives. Suppose that the number ̄ = Ω ⋃ 𝜕Ω is of total nodes covering the domain of the problem i.e. Ω N. On the other hand, we know that n is depend on point of interest x (so, after that we call it 𝑛x ) in Eq. (27) which is the 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 ). Therefore, we have 𝑛x ≤ 𝑁 and Eq. (27) can be modified as 𝑢(x) = 𝚽𝑡𝑟 (x)U𝑠 =

𝑁 ∑ 𝑗=1

𝜙𝑗 (x)𝑢𝑗 .

In this section, we consider Eq. (9) and explain how to implement SMRPI method to obtain discrete equations. Eq. (9) can be rewritten as ( 2 𝑘+1 ) 𝜕 𝑢 (x) 𝜕 2 𝑢𝑘+1 (x) 𝑖𝜆𝑢𝑘+1 (x) + 𝛼 + + 𝑤(x)𝑢𝑘+1 (x) 𝜕𝑥2 𝜕𝑦2 ) ( 2 𝑘 𝜕 𝑢 ( x) 𝜕 2 𝑢 𝑘 ( x) = 𝑖𝜆𝑢𝑘 (x) − 𝛼 + − 2𝛽𝑓 (𝑢̃ 𝑘 ) − 𝑤(x)𝑢𝑘 (x), (41) 𝜕𝑥2 𝜕𝑦2 where where 𝑢̃ is the latest available approximation of u. Consider N nodes with arbitrary distribution on the boundary and domain of the problem. Assuming that 𝑢(x𝑖 , 𝑘𝛿𝑡), 𝑖 = 1, 2, … , 𝑁 are known, our aim is to compute 𝑢(x𝑖 , (𝑘 + 1)𝛿𝑡), 𝑖 = 1, 2, … , 𝑁. So, we have N unknown values and to compute these unknown values 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 x𝑖 ∈ Ω, to obtain the discrete equations from Eq. (41), substituting approximation formulas (30) and (33) into Eq. (41) yields (𝑁 ) 𝑁 𝑁 ∑ ∑ 𝜕 2 𝜙𝑗 (x) ∑ 𝜕 2 𝜙𝑗 (x) 𝑘+1 𝑖𝜆 𝜙𝑗 (x)𝑢𝑘𝑗 +1 + 𝛼 𝑢𝑘𝑗 +1 + 𝑢𝑗 𝜕𝑥2 𝜕𝑦2 𝑗=1 𝑗=1 𝑗=1

(30)

Since corresponding to node x𝑗 there is a shape function 𝜙𝑗 (x), 𝑗 = 1, 2, … , 𝑁, if we define Ω𝑐x = {x𝑗 ∶ x𝑗 ∉ Ωx }, then obviously from Eq. (28) we obtain ∀x𝑗 ∈ Ω𝑐x ∶ 𝜙𝑗 (x) = 0.

(31)

The derivatives of 𝑢(x) with respect to x and y are 𝑁 𝜕𝑢(x) ∑ 𝜕𝜙𝑗 (x) = 𝑢 , 𝜕𝑥 𝜕𝑥 𝑗 𝑗=1

𝑁 𝜕𝑢(x) ∑ 𝜕𝜙𝑗 (x) = 𝑢𝑗 , 𝜕𝑦 𝜕𝑦 𝑗=1

(32)

𝑗=1

and for high derivatives of 𝑢(x) are 𝜕 𝑠 𝑢(x) = 𝜕𝑥𝑠 where

𝑁 ∑

𝜕 𝑠 𝜙𝑗 (x)

𝑗=1

𝜕𝑥𝑠

𝜕 𝑠 (.) 𝜕𝑥𝑠

and

𝜕 𝑠 𝑢(x) = 𝜕𝑦𝑠

𝑢𝑗 ,

𝜕 𝑠 (.)

𝑁 ∑

𝜕 𝑠 𝜙𝑗 (x)

𝑗=1

𝜕𝑦𝑠

= 𝑖𝜆 𝑢𝑗 ,

𝜙𝑗 (x)𝑢𝑘𝑗 +1

𝜙𝑗 (x)𝑢𝑘𝑗 − 𝛼

− 2𝛽𝑓 (𝑢̃ 𝑘 ) − 𝑤(x)

⎛ 𝜕𝜙1 (x1 ) ⎛ 𝑢′𝑥1 ⎞ ⎜ 𝜕𝑥 ⎜ 𝑢′ ⎟ ⎜⎜ 𝜕𝜙1 (x2 ) ⎜ 𝑥2 ⎟ = ⎜ ⎜ ⋮ ⎟ ⎜ 𝜕𝑥 ⋮ ⎜ ′ ⎟ ⎜ ⎝𝑢 𝑥 𝑁 ⎠ 𝜕𝜙 (x ) ⎜ 1 𝑁 ⎝ 𝜕𝑥

𝜕𝜙2 (x1 ) 𝜕𝑥 𝜕𝜙2 (x2 ) 𝜕𝑥 ⋮ 𝜕𝜙2 (x𝑁 ) 𝜕𝑥

⎛ 𝜕𝜙1 (x1 ) ⎞ ⎜ 𝜕𝑦 ⎟ ⎜⎜ 𝜕𝜙1 (x2 ) ⎟= ⎜ ⋮ ⎟ ⎜⎜ 𝜕𝑦 ⋮ ⎜ ′ ⎟ ⎜ ⎝𝑢 𝑦 𝑁 ⎠ 𝜕𝜙 (x ) ⎜ 1 𝑁 ⎝ 𝜕𝑦

𝜕𝜙2 (x1 ) 𝜕𝑦 𝜕𝜙2 (x2 ) 𝜕𝑦 ⋮ 𝜕𝜙2 (x𝑁 ) 𝜕𝑦

⋯ ⋯ ⋱ ⋯ ⋯ ⋯ ⋱ ⋯

𝜕𝜙𝑁 (x1 ) ⎞ ⎟⎛ ⎞ 𝜕𝑥 ⎟ 𝑢1 𝜕𝜙𝑁 (x2 ) ⎟⎜ 𝑢 ⎟ 2 ⎟⎜⎜ ⋮ ⎟⎟, 𝜕𝑥 ⎟⎜ ⎟ ⋮ 𝜕𝜙𝑁 (x𝑁 ) ⎟⎝𝑢𝑁 ⎠ ⎟ ⎠ 𝜕𝑥 𝜕𝜙𝑁 (x1 ) ⎞ ⎟⎛ ⎞ 𝜕𝑦 𝑢 𝜕𝜙𝑁 (x2 ) ⎟⎜ 1 ⎟ ⎟⎜ 𝑢 2 ⎟ 𝜕𝑦 ⎟⎜ ⋮ ⎟. ⎟⎜ ⎟ ⋮ 𝜕𝜙𝑁 (x𝑁 ) ⎟⎝𝑢𝑁 ⎠ ⎟ ⎠ 𝜕𝑦

𝜕𝑥2

𝑗=1

𝑁 ∑

𝑢𝑘𝑗 +

𝑁 ∑ 𝜕 2 𝜙𝑗 (x) 𝑗=1

𝜕𝑦2

) 𝑢𝑘𝑗

𝜙𝑗 (x)𝑢𝑘𝑗 .

(42)

Now, setting x = x𝑖 , 𝑖 = 1, 2, … , 𝑁Ω (NΩ is the number of nodes in Ω) in the above equation and applying notations (36)–(40) implies 𝑖𝜆𝑢𝑘𝑖 +1 + 𝛼

𝑁 ( ) ∑ 𝐷𝑥(2) + 𝐷𝑦(2) 𝑢𝑘𝑗 +1 + 𝑤(x𝑖 )𝑢𝑘𝑖 +1 𝑗=1

= 𝑖𝜆𝑢𝑘𝑖 − 𝛼

(34)

𝑖𝑗

𝑖𝑗

𝑁 ( ) ∑ 𝐷𝑥(2) + 𝐷𝑦(2) 𝑢𝑘𝑗 − 2𝛽𝑓 (𝑢̃ 𝑘 ) − 𝑤(x𝑖 )𝑢𝑘𝑖 , 𝑗=1

𝑖𝑗

𝑖𝑗

𝑖 = 1, 2, … , 𝑁Ω . (43)

For nodes which are located on the boundary 𝜕Ω using Eq. (3) we have the following relation 𝑢𝑘𝑖 +1 = ℎ(x𝑖 , (𝑘 + 1)𝛿𝑡),

𝑖 = 1, 2, … , 𝑁𝜕Ω . (44) ̄ i.e. 𝑁 = 𝑁Ω + Suppose that N is the total number of nodes covering Ω, 𝑁𝜕Ω . Therefore, we have 𝑁Ω + 𝑁𝜕Ω , i.e. N, equations and N unknowns. If we show the sparse system as the following matrix form

(35)

̃ 𝑘 ) + b, AU𝑘+1 = BU𝑘 + F(U

(45)

then matrixes A, B, U and b are defined as follows: ( ) ⎧𝑖𝜆𝛿𝑖𝑗 + 𝐷𝑥(2)𝑖𝑗 + 𝐷𝑦(2) 𝑖𝑗 ⎪ A𝑖𝑗 = ⎨ + 𝑑𝑖𝑎𝑔𝑖 (𝑤(x𝑖 )), 𝑖 = 1, 2, … , 𝑁Ω , for ⎪ 𝑖 = 𝑁Ω + 1, … , 𝑁 ⎩𝛿𝑖𝑗 ,

We denote the above coefficients matrices as Dx, Dy and also for high order derivatives we have the following matrix-vector multiplications: 𝑈𝑦(𝑠) = 𝐷(𝑠) 𝑦𝑈 ,

(𝑁 ∑ 𝜕 2 𝜙𝑗 (x)

𝑗=1

are s’th derivatives with respect to x and y. Denot-

𝜕𝑦𝑠 𝜕 𝑠 (.) 𝜕 𝑠 (.) ing = , 𝑢(𝑦𝑠) (.) = and setting x = x𝑖 in Eq. (32) then the 𝜕𝑥𝑠 𝜕𝑦𝑠 following matrix form is given

𝑈𝑥(𝑠) = 𝐷(𝑠) 𝑥𝑈 ,

𝑁 ∑ 𝑗=1

(33)

𝑢(𝑥𝑠) (.)

⎛ 𝑢′𝑦1 ⎜ 𝑢′ ⎜ 𝑦2

𝑁 ∑

+ 𝑤 ( x)

(36)

𝑗 = 1, 2, … , 𝑁, (46)

where (

𝑈𝑥(𝑠) = 𝑢(𝑥𝑠) , 𝑢(𝑥𝑠) , … , 𝑢(𝑥𝑠) 1

𝐷𝑥(𝑠) = 𝑖𝑗

𝐷𝑦(𝑠) = 𝑖𝑗

𝜕 𝑠 𝜙𝑗 (x𝑖 ) 𝜕𝑥𝑠

𝑁

2

,

)𝑡𝑟

,

𝑈𝑦(𝑠)

( )𝑡𝑟 = 𝑢(𝑦𝑠) , 𝑢(𝑦𝑠) , … , 𝑢(𝑦𝑠) , 1

2

𝑁

( ) ⎧𝑖𝜆𝛿𝑖𝑗 − 𝐷𝑥(2) + 𝐷𝑦(2) 𝑖𝑗 𝑖𝑗 ⎪ B𝑖𝑗 = ⎨ − 𝑑𝑖𝑎𝑔𝑖 (𝑤(x𝑖 )), 𝑖 = 1, 2, … , 𝑁Ω , ⎪0 , 𝑖 = 𝑁Ω + 1, … , 𝑁 ⎩

(37)

(38)

for 𝑗 = 1, 2, … , 𝑁, (47)

𝜕 𝑠 𝜙𝑗 (x𝑖 ) , 𝜕𝑦𝑠

b𝑖 = 0, for 𝑖 = 1, 2, … , 𝑁Ω ,

(39)

𝑡𝑟

𝑈 = (𝑢1 , 𝑢2 , … , 𝑢𝑁 ) .

b𝑖 = ℎ

𝑘+1

(x𝑖 ), for 𝑖 = 𝑁Ω + 1, … , 𝑁, (48)

( )𝑡𝑟 U𝑘 = 𝑢𝑘1 , 𝑢𝑘2 , … , 𝑢𝑘𝑁 . In Eqs. (46) and (47), 𝛿 ij is the Kronecker delta function, i.e.

(40) 77

(49)

E. Shivanian, A. Jafarabadi

Engineering Analysis with Boundary Elements 83 (2017) 74–86

Fig. 1. Considered domains in the current paper.

Table 1 Numerical results of absolute errors and relative errors for real and imaginary parts with different h, 𝛿t and 𝑇 = 1 on [0, 1]2 for Example 1. Real part

ℎ = 𝛿𝑡 = 1∕10 ℎ = 1∕20, 𝛿𝑡 = 1∕50 ℎ = 1∕24, 𝛿𝑡 = 1∕100 ℎ = 1∕32, 𝛿𝑡 = 1∕1000

(where l denotes the number of iterations in each time level). We put 𝑘

̃ )= F(U

Imaginary part

𝜀∞ (u)

𝜀R (u)

𝜀∞ (u)

𝜀R (u)

1.0000𝑒−004 8.2238𝑒−007 2.3506𝑒−007 4.1033𝑒−008

8.0579𝑒−005 4.0418𝑒−007 1.0445𝑒−007 1.8468𝑒−008

4.6741𝑒−005 9.7089𝑒−007 1.9839𝑒−007 6.1393𝑒−008

5.9660𝑒−005 1.1873𝑒−006 2.4121𝑒−007 5.6992𝑒−008

Real part 𝜀R (u)

Imaginary part 𝜀R (u)

CPU time (s)

1 2 3 4

5.0727𝑒−008 5.2717𝑒−008 7.1667𝑒−008 4.7897𝑒−008

4.8855𝑒−008 1.7078𝑒−007 1.0034𝑒−007 1.2685𝑒−007

64.5981 68.1999 72.2288 72.5101

𝛿𝑖𝑗 =

{ 1, 𝑖 = 𝑗 , 0, 𝑖 ≠ 𝑗

+1 ‖U𝑘𝑙 +1 − U𝑘𝑙−1 ‖ < 𝜖,

(53)

where 𝜖 is a fixed number, in this case is equal to 10−15 . Finally we put U𝑘+1 = U𝑘𝑙 +1 . After this, we can go to the next time level. Therefore the nonlinearity of equation is eliminated by the described iterative (P-C) method. This process is iterated until reaching to the desirable time T. 5. Numerical results In this section, we present the numerical results of the proposed method on three test problems. We test the accuracy and stability of the method described in this paper by performing the mentioned scheme for different values of final time T, h (the distance between the nodes in x or y direction) and 𝛿t. To show the accuracy and convergence of the method, two kinds of error measures, maximum absolute error 𝜀∞ and relative error 𝜀R defined by:

(50)

𝜀∞ (𝑢) = ‖𝑢𝑒𝑥𝑎𝑐𝑡 − 𝑢𝑎𝑝𝑝𝑟𝑜𝑥 ‖∞

and also diagi (·) denotes the i’th diagonal element of diagonal matrix diag(·).

= max{|𝑢𝑒𝑥𝑎𝑐𝑡 (x𝑖 , 𝑡) − 𝑢𝑎𝑝𝑝𝑟𝑜𝑥 (x𝑖 , 𝑡)| ∶ 𝑖 = 1, 2, … , 𝑁},

Remark 1. For dealing with the nonlinearity we use the following prõ 0 ) = −2𝛽|U0 |2 ∗ U0 (𝑘 = 0) where “∗ ” decedure [47]. First, we put F(U notes the Hadamard product. Now, Eq. (45) can be solved as a system of linear algebraic equations for the unknown U1 , then we use the following Crank–Nicolson scheme: ̃ 0 ) = 1 (−2𝛽|U0 |2 ∗ U0 − 2𝛽|U1 |2 ∗ U1 ). F(U 2

(52)

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

Table 2 Numerical results of relative errors and CPU times used for real and imaginary parts with ℎ = 1∕32, 𝛿𝑡 = 0.01 and different time level on domain Ω1 for Example 1. Time T

1 (−2𝛽|U𝑘 |2 ∗ U𝑘 − 2𝛽|U𝑘𝑙 +1 |2 ∗ U𝑘𝑙 +1 ). 2

√∑ √ 𝑁 √ 𝑖=1 (𝑢𝑒𝑥𝑎𝑐𝑡 (x𝑖 , 𝑡) − 𝑢𝑎𝑝𝑝𝑟𝑜𝑥 (x𝑖 , 𝑡))2 𝜀𝑅 (𝑢) = √ , ∑𝑁 2 𝑖=1 (𝑢𝑒𝑥𝑎𝑐𝑡 (x𝑖 , 𝑡))

(54)

(55)

are used, where uexact and uapprox are exact and approximate solutions, respectively. There are a number of types of RBFs, and the characteristics of RBFs have been widely investigated. In the current work, we have chosen the thin plate spline (TPS) as radial basis functions in Eq. (26). This RBF is defined as follows:

(51)

̃ ), for unknown U1 in iteration Again Eq. (45) is solved using the new F(U (𝑙 = 2). Generally, at the time level (𝑘 + 1), we iterate between calculat̃ 𝑘 ) and computing the approximation value of the unknown U𝑘+1 ing F(U 𝑙 0

𝑅(x) = 𝑟2𝛽 ln(𝑟), 78

𝛽 = 1, 2, …

(56)

E. Shivanian, A. Jafarabadi

Engineering Analysis with Boundary Elements 83 (2017) 74–86

Fig. 2. Surface plots of the numerical and analytical solutions u(x, y, T) with 𝑁 = 441, 𝛿𝑡 = 0.01 and 𝑇 = 1 on [0, 1]2 for Example 1.

Fig. 3. Graphs of approximate and exact solutions u(x, 0, t) at various time levels with 𝑁 = 1089, 𝛿𝑡 = 0.01 and 𝑇 = 1 on −5 ≤ 𝑥 ≤ 5 for Example 1.

For the second-order partial differential equation (1), 𝛽 = 2 is used for thin plate splines and also in test problems, we let 𝑟𝑠 = 4.2ℎ. In Eq. (26) set 𝑚 = 28 that results polynomial basis functions as follows:

𝑢(𝑥, 𝑦, 𝑡) =

𝑖 exp(𝑖𝑡) . cosh 𝑥 cosh 𝑦

(57)

We solve this problem by the present method with several values of h, 𝛿t and final time T on domains [0, 1]2 and Ω1 . Numerical results of 𝜀∞ and 𝜀R on [0, 1]2 for real and imaginary parts with different h, 𝛿t and 𝑇 = 1 have been listed in Table 1. Also in Table 2, we have given the numerical results of 𝜀R for real and imaginary parts of u(x, y, t) with ℎ = 1∕32, 𝛿𝑡 = 0.01 and different time level on domain Ω1 . In Fig. 2, it is shown the surface plots of the numerical and analytical solutions for real and imaginary parts with 𝑁 = 441, 𝛿𝑡 = 0.01 and 𝑇 = 1 on [0, 1]2 . Fig. 3 shows the graphs of approximate and exact solutions u(x, 0, t) at various time levels with 𝑁 = 1089, 𝛿𝑡 = 0.01 and 𝑇 = 1 on 𝑥 ∈ [−5, 5]. The absolute errors of real and imaginary parts with 𝑁 = 1681, 𝛿𝑡 = 0.01 at several time levels T on [0, 1]2 , have been plotted in Fig. 4. Graphs of approximate solutions and absolute errors for real and imaginary parts of u(x, y, T) with ℎ = 1∕40, 𝛿𝑡 = 0.01 and 𝑇 = 1 on Ω1 , have been shown in Fig. 5. From figures and tables, one can see that the SMRPI provides an accurate and stable approximate solution to the cubic nonlinear Schrödinger equation.

𝑃 𝑇 (x) = {1, 𝑥, 𝑦, 𝑥2 , 𝑥𝑦, 𝑦2 , 𝑥3 , 𝑥2 𝑦, 𝑥𝑦2 , 𝑦3 , 𝑥4 , 𝑥3 𝑦, 𝑥2 𝑦2 , 𝑥𝑦3 , 𝑦4 , 𝑥5 , 𝑥4 𝑦, 𝑥3 𝑦2 , 𝑥2 𝑦3 , 𝑥𝑦4 , 𝑦5 , 𝑥6 , 𝑥5 𝑦, 𝑥4 𝑦2 , 𝑥3 𝑦3 , 𝑥2 𝑦4 , 𝑥𝑦5 , 𝑦6 }. Fig. 1 presents the considered irregular domains in test problems which are defined as follows: Ω1 = {(𝑥, 𝑦) ∈ [0, 1]2 } ⧵ {(𝑥, 𝑦) ∈ (1∕4, 3∕4)2 }, Ω2 = {(𝑥, 𝑦) ∈ [−1, 1]2 } ⧵ {(𝑥, 𝑦) ∶ 𝑥2 + 𝑦2 < 9∕16}, and also, the circumference of Ω3 is 𝑟 = 1 − 1∕4 cos 4𝜃, where 0 ≤ 𝜃 ≤ 2𝜋. Finally, Ω4 presents the irregular distribution of collocation nodes covering the domain [0, 1]2 (using the MATLAB routine ‘haltonset’). Example 1. Consider Eqs. (1)–(3) with 𝛼 = 1, 𝛽 = 0 and potential function 𝑤(𝑥, 𝑦) = 3 − 2 tanh2 𝑥 − 2 tanh2 𝑦. The exact solution of this example is [14,15,18,19] 79

E. Shivanian, A. Jafarabadi

Engineering Analysis with Boundary Elements 83 (2017) 74–86

Fig. 4. Absolute error of real and imaginary parts u(x, y, T) with 𝑁 = 1681, 𝛿𝑡 = 0.01 at several time levels T, on [0, 1]2 for Example 1. Table 3 Numerical results of absolute errors and relative errors for real and imaginary parts with different h, 𝛿t and 𝑇 = 1 on [0, 1]2 for Example 2. Real part

ℎ = 𝛿𝑡 = 1∕10 ℎ = 1∕20, 𝛿𝑡 = 1∕100 ℎ = 1∕26, 𝛿𝑡 = 1∕500 ℎ = 1∕34, 𝛿𝑡 = 1∕1000

Table 4 Numerical results of relative errors and CPU times used for real and imaginary parts with ℎ = 1∕16, 𝛿𝑡 = 0.1 and different time level on domain Ω2 for Example 2.

Imaginary part

𝜀∞ (u)

𝜀R (u)

𝜀∞ (u)

𝜀R (u)

1.9015𝑒−004 3.7862𝑒−007 2.2082𝑒−008 5.7132𝑒−009

2.3292𝑒−004 3.4666𝑒−007 2.2237𝑒−008 5.5729𝑒−009

8.3516𝑒−005 2.7738𝑒−006 9.8194𝑒−008 2.1764𝑒−008

9.5293𝑒−005 3.5986𝑒−006 1.2969𝑒−007 3.1479𝑒−008

Example 2. Consider Eqs. (1)–(3) with 𝛼 = 1 and potential function w(x, y) ≡ 0. The exact solution of this example is [22,23] (( )) 𝑢(𝑥, 𝑦, 𝑡) = 𝐴 exp 𝑖 𝑘1 𝑥 + 𝑘2 𝑦 − 𝐵𝑡 , 𝐵 = 𝑘21 + 𝑘22 − 𝛽|𝐴|2 . (58)

Time T

Real part 𝜀R (u)

Imaginary part 𝜀R (u)

CPU time (s)

1 2 3 4

1.0888𝑒−004 5.6762𝑒−005 9.2350𝑒−005 8.7799𝑒−005

1.1407𝑒−004 5.3798𝑒−005 9.2308𝑒−005 8.9011𝑒−005

45.5005 45.6087 46.8991 47.4177

[18]. We solve this problem by the present method with several values of h, 𝛿t and final time T on domains [0, 1]2 and Ω2 . Table 3 presents the numerical results of 𝜀∞ and 𝜀R on [0, 1]2 for real and imaginary parts with different h, 𝛿t and 𝑇 = 1. In Table 4, the numerical results of 𝜀R for real and imaginary parts u(x, y, t) with ℎ = 1∕16, 𝛿𝑡 = 0.1 and different time level on domain Ω2 have been listed. The main aim of

The equation describes a progressive plane wave [23]. For a case study, the current model is investigated by letting 𝐴 = 1∕2 and 𝑘1 = 𝑘2 = 𝛽 = 1

Table 5 Numerical results of absolute errors and relative errors for real and imaginary parts with ℎ = 1∕10, 𝛿𝑡 = 0.01 and 𝑇 = 1 on [0, 1]2 for Example 2. Real part

𝑚=6 𝑚 = 10 𝑚 = 15 𝑚 = 28

Imaginary part

CPU time(s)

𝜀∞ (u)

𝜀R (u)

𝜀∞ (u)

𝜀R (u)

5.3618𝑒−005 1.6161𝑒−005 1.4156𝑒−006 4.4859𝑒−007

5.3298𝑒−005 1.8871𝑒−005 1.0276𝑒−006 3.4649𝑒−007

7.7048𝑒−005 6.5307𝑒−005 3.4319𝑒−006 2.5875𝑒−006

9.5919𝑒−005 8.6061𝑒−005 4.4537𝑒−006 3.2484𝑒−006

80

4.237939 4.267952 4.286870 4.420315

E. Shivanian, A. Jafarabadi

Engineering Analysis with Boundary Elements 83 (2017) 74–86

Fig. 5. Graphs of approximate solution and absolute error for real and imaginary parts u(x, y, T) with ℎ = 1∕40, 𝛿𝑡 = 0.01 and 𝑇 = 1, on domain Ω1 for Example 1.

Fig. 6. Surface plots of the numerical and analytical solutions u(x, y, T) with 𝑁 = 121, 𝛿𝑡 = 0.01 and 𝑇 = 1 on [0, 1]2 for Example 2.

Example 3. Consider Eqs. (1)–(3) with 𝛼 = 1, 𝛽 = 2𝜋 2 − 1 and poten( ) tial function 𝑤(𝑥, 𝑦) = (2𝜋 2 − 1) 1 − cos2 𝜋𝑥 cos2 𝜋𝑦 . The exact solution of this example is [15,18]

Table 5 is to explain the influence of the number of polynomial basis functions on the errors. It is seen that by increasing the value of “m” the errors are improved. The surface plots of the numerical and analytical solutions for real and imaginary parts with 𝑁 = 121, 𝛿𝑡 = 0.01 and 𝑇 = 1 on [0, 1]2 , have been shown in Fig. 6. In Fig. 7, the absolute error of real and imaginary parts with 𝑁 = 1681, 𝛿𝑡 = 0.01 at several time levels T on [0, 1]2 is plotted. Fig. 8 shows the graphs of approximate solutions and absolute errors for real and imaginary parts of u(x, y, T) with ℎ = 1∕17, 𝛿𝑡 = 0.1 and 𝑇 = 1 on Ω2 .

𝑢(𝑥, 𝑦, 𝑡) = exp(−𝑖𝑡) cos 𝜋𝑥 cos 𝜋𝑦.

(59)

We solve this problem by the current method with several values of h, 𝛿t and final time T on domains [0, 1]2 , Ω3 and Ω4 . Table 6 lists the numerical results of 𝜀∞ and 𝜀R on [0, 1]2 for real and imaginary parts with different h, 𝛿t and 𝑇 = 1. In Table 7, the numerical results of 𝜀R for 81

E. Shivanian, A. Jafarabadi

Engineering Analysis with Boundary Elements 83 (2017) 74–86

Fig. 7. Absolute error of real and imaginary parts u(x, y, T) with 𝑁 = 1681, 𝛿𝑡 = 0.01 at several time levels T, on [0, 1]2 for Example 2.

Fig. 8. Graphs of approximate solution and absolute error for real and imaginary parts u(x, y, T) with ℎ = 1∕17, 𝛿𝑡 = 0.1 and 𝑇 = 1, on domain Ω2 for Example 2. 82

E. Shivanian, A. Jafarabadi

Engineering Analysis with Boundary Elements 83 (2017) 74–86

Fig. 9. Surface plots of the numerical and analytical solutions u(x, y, T) with 𝑁 = 625, 𝛿𝑡 = 0.01 and 𝑇 = 1 on [0, 1]2 for Example 3.

Fig. 10. Absolute error of real and imaginary parts u(x, y, T) with 𝑁 = 1089, 𝛿𝑡 = 0.01 at several time levels T, on [0, 1]2 for Example 3.

83

E. Shivanian, A. Jafarabadi

Engineering Analysis with Boundary Elements 83 (2017) 74–86

Fig. 11. Graphs of approximate solution and absolute error for real and imaginary parts u(x, y, T) with ℎ = 1∕17, 𝛿𝑡 = 0.01 and 𝑇 = 1, on domain Ω3 for Example 3.

real and imaginary parts u(x, y, t) with ℎ = 1∕17, 𝛿𝑡 = 0.01 and different time level on domain Ω3 are listed. Fig. 9 shows the surface plots of the numerical and analytical solutions for real and imaginary parts with 𝑁 = 625, 𝛿𝑡 = 0.01 and 𝑇 = 1 on [0, 1]2 . In Fig. 10, the absolute error of real and imaginary parts with 𝑁 = 1089, 𝛿𝑡 = 0.01 at several time levels T on [0, 1]2 is shown. Fig. 11 presents the graphs of approximate solutions and absolute errors for real and imaginary parts of u(x, y, T) with ℎ = 1∕17, 𝛿𝑡 = 0.01 and 𝑇 = 1 on Ω3 . Graphs of approximate solution and absolute error for real and imaginary parts u(x, y, T) with 𝑁 = 440, 𝛿𝑡 = 0.01 and 𝑇 = 1, on domain Ω4 have been demonstrated in Fig. 12. It is clear that the accuracy of the present method do not depend on the distribution of collocation nodes.

6. Conclusion In this article, the spectral meshless radial point interpolation (SMRPI) has been applied to the numerical solutions of the two-dimensional cubic nonlinear Schrödinger equations on regular and irregular domains. For discretization, firstly, we discretized the time derivative using a finite difference formula and obtained a time discrete scheme. We proved that the time semi-discrete scheme is unconditionally stable and convergent using the energy method. 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

Fig. 12. Graphs of approximate solution and absolute error for real and imaginary parts u(x, y, T) with 𝑁 = 440, 𝛿𝑡 = 0.01 and 𝑇 = 1, on domain Ω4 for Example 3. 84

E. Shivanian, A. Jafarabadi

Engineering Analysis with Boundary Elements 83 (2017) 74–86

Table 6 Numerical results of absolute errors and relative errors for real and imaginary parts with different h, 𝛿t and 𝑇 = 1 on [0, 1]2 for Example 3. Real part

ℎ = 𝛿𝑡 = 1∕14 ℎ = 1∕20, 𝛿𝑡 = 1∕100 ℎ = 1∕26, 𝛿𝑡 = 1∕1000 ℎ = 1∕34, 𝛿𝑡 = 1∕2000

[15] Dehghan M, Mirzaei D. The meshless local Petrov–Galerkin (MLPG) method for the generalized two-dimensional non-linear Schrödinger equation. Eng Anal Bound Elem 2008;32(9):747–56. [16] Dehghan M, Taleei A. A Chebyshev pseudospectral multidomain method for the soliton solution of coupled nonlinear Schrödinger equations. Comput Phys Commun 2011;182(12):2519–29. [17] Mohebbi A, Dehghan M. The use of compact boundary value method for the solution of two-dimensional Schrödinger equation. J Comput Appl Math 2009;225(1):124–34. [18] Abbasbandy S, Ghehsareh HR, Hashim I. A meshfree method for the solution of two-dimensional cubic nonlinear Schrödinger equation. Eng Anal Bound Elem 2013;37(6):885–98. [19] Abdur R, Ismail AIBM. Numerical studies on two-dimensional schrodinger equation by Chebyshev spectral collocation method. UPB Sci Bull Series A 2011;73(1):101–10. [20] Dereli Y. The meshless kernel-based method of lines for the numerical solution of the nonlinear Schrödinger equation. Eng Anal Bound Elem 2012;36(9):1416–23. [21] Tian ZF, Yu P. High-order compact adi (hoc-adi) method for solving unsteady 2d Schrödinger equation. Comput Phys Commun 2010;181(5):861–8. [22] Xu Y, Zhang L. Alternating direction implicit method for solving two-dimensional cubic nonlinear Schrödinger equation. Comput Phys Commun 2012;183(5):1082–93. [23] Gao Z, Xie S. Fourth-order alternating direction implicit compact finite difference schemes for two-dimensional Schrödinger equations. Appl Numer Math 2011;61(4):593–614. [24] Manafian J. Optical soliton solutions for Schrödinger type nonlinear evolution equations by the tan (𝜙 (𝜉)/2)-expansion method. Optik-Int J Light Electron Optics 2016;127(10):4222–45. [25] Arbabi S, Najafi M. Exact solitary wave solutions of the complex nonlinear Schrödinger equations. Optik-Int J Light Electron Optics 2016;127(11):4682–8. [26] Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int J Numer Methods Eng 1994;37(2):229–56. [27] 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. [28] Duan Y, Tan Y-J. A meshless Galerkin method for Dirichlet problems using radial basis functions. J Comput Appl Math 2006;196(2):394–401. [29] 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. [30] 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. [31] Peng M, Li D, Cheng Y. The complex variable element-free Galerkin (CVEFG) method for elasto-plasticity problems. Eng Struct 2011;33(1):127–35. [32] Dai B, Cheng Y. An improved local boundary integral equation method for two-dimensional potential problems. Int J Appl Mech 2010;2(02):421–36. [33] Fu-Nong B, Dong-Ming L, Jian-Fei W, Yu-Min C. An improved complex variable element-free Galerkin method for two-dimensional elasticity problems. Chin Phys B 2012;21(2):020204. [34] Shen SNAS. The meshless local Petrov-Galerkin (MLPG) method: a simple & less– costly alternative to the finite element and boundary element methods. Comput Model Eng Sci 2002;3:11–51. [35] Nayroles B, Touzot G, Villon P. Generalizing the finite element method: diffuse approximation and diffuse elements. Comput Mech 1992;10(5):307–18. [36] Liu G, Yan L, Wang J, Gu Y. Point interpolation method based on local residual formulation using radial basis functions. Struct Eng Mech 2002;14(6):713–32. [37] 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. [38] 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(12):1693–702. [39] Atluri SN, Zhu T. A new meshless local Petrov-Galerkin (MLPG) approach in computational mechanics. Comput Mech 1998;22(2):117–27. [40] Dehghan M, Mirzaei D. Meshless local Petrov–Galerkin (MLPG) method for the unsteady magnetohydrodynamic (MHD) flow through pipe with arbitrary wall conductivity. Appl Numer Math 2009;59(5):1043–58. [41] Gu Y, Liu G-R. A meshless local Petrov-Galerkin (MLPG) method for free and forced vibration analyses for solids. Comput Mech 2001;27(3):188–98. [42] Shivanian E. Meshless local Petrov–Galerkin (MLPG) method for three-dimensional nonlinear wave equations via moving least squares approximation. Eng Anal Bound Elem 2015a;50:249–57. [43] 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:98–105. [44] Hosseini VR, Shivanian E, Chen W. Local integration of 2-d fractional telegraph equation via local radial point interpolant approximation. Eur Phys J Plus 2015;130(2):1–21. [45] Aslefallah M, Shivanian E. Nonlinear fractional integro-differential reaction-diffusion equation via radial basis functions. Eur Phys J Plus 2015;130(3):1–9. [46] Shivanian E. On the convergence analysis, stability, and implementation of meshless local radial point interpolation on a class of three-dimensional wave equations. Int J. Numer Methods Eng 2016;105(2):83–110. [47] 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.

Imaginary part

𝜀∞ (u)

𝜀R (u)

𝜀∞ (u)

𝜀R (u)

9.1304𝑒−006 6.3766𝑒−006 1.0660𝑒−006 4.1518𝑒−007

1.6899𝑒−005 1.1802𝑒−005 1.9730𝑒−006 7.6843𝑒−007

2.4102𝑒−005 5.5303𝑒−006 3.1469𝑒−006 1.1395𝑒−006

2.8643𝑒−005 6.5722𝑒−006 3.7397𝑒−006 1.3542𝑒−006

Table 7 Numerical results of relative errors and CPU times used for real and imaginary parts with ℎ = 1∕17, 𝛿𝑡 = 0.01 and different time level on domain Ω3 for Example 3. Time T

Real part 𝜀R (u)

Imaginarypart 𝜀R (u)

CPU time (s)

1 2 3 4

5.7460𝑒−005 8.2930𝑒−005 7.3746𝑒−006 5.8837𝑒−005

2.5700𝑒−005 2.8893𝑒−005 5.8006𝑒−005 3.5377𝑒−005

111.8982 131.1241 139.2076 152.3902

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 is considerable that the simple implementation of SMRPI is the main advantage of the present method. The presented results through the figures and tables show the accuracy and efficiency of the method. The numerical examples of nonlinear Schrödinger equations have been worked out which show that the SMRPI method is simple and useful to solve the the two-dimensional cubic nonlinear Schrödinger equations. Acknowledgments The authors are grateful to anonymous reviewers for carefully reading this paper and for their comments and suggestions that have improved the paper. References [1] Arnold A. Numerically absorbing boundary conditions for quantum evolution equations. VLSI Des 1998;6(1–4):313–19. [2] Agrawal GP. Nonlinear fiber optics. Amsterdam: Academic Press; 2007. [3] Saleh BEA, Teich MC. Fundamentals of Photonics. New York, USA: John Wiley & Sons, Inc.; 1991. [4] Levy M. Parabolic equation methods for electromagnetic wave propagation. London: Inst. of Electr. Eng.; 2000. [5] Wazwaz A-M. A study on linear and nonlinear schrodinger equations by the variational iteration method. Chaos Solitons Fractals 2008;37(4):1136–42. [6] Biswas A, Milovic D. Bright and dark solitons of the generalized nonlinear Schrödingers equation. Commun Nonlinear Sci Numer Simul 2010;15(6):1473–84. [7] Zhang Z, Liu Z, Miao X, Chen Y. Qualitative analysis and traveling wave solutions for the perturbed nonlinear Schrödinger’s equation with kerr law nonlinearity. Phys Lett A 2011;375(10):1275–80. [8] Kivshar YS, Agrawal G. Optical solitons: from fibers to photonic crystals. Amsterdam: Academic press; 2003. [9] Remoissenet M. Waves called solitons: concepts and experiments. Berlin: Springer Science & Business Media; 2013. [10] Subaşi M. On the finite-differences schemes for the numerical solution of two dimensional Schrödinger equation. Numer Methods Partial Differ. Equ. 2002;18(6):752–8. [11] Antoine X, Besse C, Mouysset V. Numerical schemes for the simulation of the two-dimensional Schrödinger equation using non-reflecting boundary conditions. Math Comput 2004;73(248):1779–99. [12] Delfour M, Fortin M, Payr G. Finite-difference solutions of a non-linear Schrödinger equation. J. Comput Phys 1981;44(2):277–88. [13] Xie S-S, Li G-X, Yi S. Compact finite difference schemes with high accuracy for one-dimensional nonlinear Schrödinger equation. Comput Methods Appl Mech Eng 2009;198(9):1052–60. [14] Dehghan M, Shokri A. A numerical method for two-dimensional Schrödinger equation using collocation and radial basis functions. Comput Math Appl 2007;54(1):136–46. 85

E. Shivanian, A. Jafarabadi

Engineering Analysis with Boundary Elements 83 (2017) 74–86

[48] Tadeu A, Chen C, António J, Simoes N. A boundary meshless method for solving heat transfer problems using the fourier transform. Adv Appl Math Mech 2011;3(05):572–85. [49] Shivanian E. A new spectral meshless radial point interpolation (SMRPI) method: a well-behaved alternative to the meshless weak forms. Eng Anal Bound Elem 2015b;54:1–12.

[50] Shivanian E, Jafarabadi A. Numerical solution of two-dimensional inverse force function in the wave equation with nonlocal boundary conditions. Invers Probl Sci Eng 2017:1–25. [51] Folland GB. Real analysis: modern techniques and their applications. New York: John Wiley & Sons; 2013. [52] Brown JW, Churchill RV, Lapidus M. Complex variables and applications, vol. 7. New York: McGraw-Hill; 1996.

86