An implicit RBF meshless approach for solving the time fractional nonlinear sine-Gordon and Klein–Gordon equations

Engineering Analysis with Boundary Elements 50 (2015) 412–434

Contents lists available at ScienceDirect

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

An implicit RBF meshless approach for solving the time fractional nonlinear sine-Gordon and Klein–Gordon equations Mehdi Dehghan a,n, Mostafa Abbaszadeh a, Akbar Mohebbi b a Department of Applied Mathematics, Faculty of Mathematics and Computer Science, Amirkabir University of Technology, No. 424, Hafez Ave., 15914, Tehran, Iran b Department of Applied Mathematics, Faculty of Mathematical Science, University of Kashan, Kashan, Iran

art ic l e i nf o

a b s t r a c t

Article history: Received 15 February 2014 Received in revised form 19 May 2014 Accepted 10 September 2014

In this paper, we propose a numerical method for the solution of time fractional nonlinear sine-Gordon equation that appears extensively in classical lattice dynamics in the continuum media limit and Klein–Gordon equation which arises in physics. In this method we first approximate the time fractional derivative of the mentioned equations by a scheme of order Oðτ3  α Þ; 1 o α o 2 then we will use the Kansa approach to approximate the spatial derivatives. We solve the two-dimensional version of these equations using the method presented in this paper on different domains such as rectangular and nonrectangular domains. Also, we prove the unconditional stability and convergence of the time discrete scheme. We show that convergence order of the time discrete scheme is OðτÞ. We solve these fractional PDEs on different non-rectangular domains. The aim of this paper is to show that the meshless method based on the radial basis functions and collocation approach is also suitable for the treatment of the nonlinear time fractional PDEs. The results of numerical experiments are compared with analytical solutions to confirm the accuracy and efficiency of the presented scheme. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Time fractional sine-Gordon equation Time fractional Klein–Gordon equation Kansa method Radial basis functions (RBFs) Unconditional stability Convergence Energy method Non-rectangular domains

1. Introduction In recent years there has been a growing interest in the field of fractional calculus [2,3,21,24,31,53,54,58,66,69,70]. Fractional differential equations have attracted increasing attention because they have applications in various fields of science and engineering. Many phenomena in fluid mechanics, viscoelasticity, chemistry, physics, finance and other sciences can be described very successfully by models using mathematical tools from fractional calculus, i.e., the theory of derivatives and integrals of fractional order. Some of the most applications are given in the book of Oldham and Spanier [68], the book of Podlubny [70] and the papers of Metzler and Klafter [56], Bagley and Trovik [4]. Many considerable works on the theoretical analysis [20,89] have been carried on, but analytic solutions of most fractional differential equations cannot be obtained explicitly. So many authors have resorted to numerical solution strategies based on convergence and stability analysis [7,8,12,23,44,52,57,64,65,62,73,76,80,93,94]. There are several definitions of a fractional derivative of order α 4 0 [68,69]. The two most commonly used are the Riemann–Liouville and Caputo. The difference between the two definitions is in the order of evaluation [67]. We start with recalling the essentials of the fractional calculus. The fractional calculus is a name for the theory of integrals and derivatives of arbitrary order, which unifies and generalizes the notions of integer-order differentiation and n-fold integration. Now, we give some basic definitions and properties of the fractional calculus theory. Definition 1. For μ A R and x 4 0, a real function f(x), is said to be in the space C μ if there exists a real number p 4 μ such that f ðxÞ ¼ xp f 1 ðxÞ, m where f 1 ðxÞ A Cð0; 1Þ, and for m A N it is said to be in the space C m μ if f A C μ . Definition 2. The Riemann–Liouville fractional integral operator of order α 4 0 for a function f ðxÞ A C μ ; μ Z  1, is defined as J α f ðxÞ ¼

n

1 Γ ðα Þ

Z

x 0

ðx  tÞα  1 f ðtÞ dt;

α 4 0;

x 4 0;

J 0 f ðxÞ ¼ f ðxÞ:

Corresponding author. E-mail addresses: [email protected], [email protected] (M. Dehghan), [email protected] (M. Abbaszadeh), [email protected] (A. Mohebbi).

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

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

413

Also we have the following properties:  J α J β f ðxÞ ¼ J α þ β f ðxÞ;  J α J β f ðxÞ ¼ J β J α f ðxÞ; Γ ðγ þ1Þ α þ γ x :  J α xγ ¼ Γ ðα þ γ þ 1Þ Definition 3. If m be the smallest integer that exceeds α, the Caputo time fractional derivative operator of order α 4 0 is defined as 8 Z t m > 1 ∂ uðx; sÞ > > ðt  sÞm  α  1 ds; m 1 o α o m; m A N; < α ∂ uðx; tÞ Γ ðm  αÞ 0 ∂sm C α ð1:1Þ D uðx; tÞ ¼ ¼ 0 t > ∂t α ∂m uðx; tÞ > > ; m ¼ α: : m ∂t

In this paper, we consider the time fractional nonlinear equation ∂α uðx; tÞ ¼ ∇2 uðx; tÞ  Fðuðx; tÞÞ þ f ðx; tÞ; ∂t α

x A Ω;

t A ½0; T;

ð1:2Þ

with the initial conditions uðx; 0Þ ¼ ωðxÞ; ∂uðx; tÞ  ¼ ψ ðxÞ;  ∂t t¼0

x A Ω;

ð1:3Þ

and boundary condition uðx; tÞ ¼ hðx; tÞ;

x A ∂Ω;

t 40;

ð1:4Þ

where uðx; tÞ represents the solute concentration, f ðx; tÞ is source term, hðx; tÞ is boundary solute concentration, ωðxÞ, ψ ðxÞ are initial solute concentration, α is the order of time derivative, ∂α uðx; tÞ=∂t α is the Caputo fractional derivative and is defined as follows: Z t 2 ∂α uðx; tÞ 1 ∂ uðx; sÞ ds ¼ ; α A ð1; 2Þ; ð1:5Þ α ∂t Γ ð2  αÞ 0 ∂s2 ðt  sÞα  1 and assume Fðuðx; tÞÞ satisfies the Lipschitz condition with respect to u   FðuÞ  FðuÞ ~  r Lju  uj; ~ ~ 8 u; u:

ð1:6Þ

Now we consider two important cases of Eq. (1.2):

 Firstly, in Eq. (1.2), if we put FðuÞ ¼ sin ðuÞ then we obtain time fractional nonlinear sine-Gordon equation. The main aim of [45] is to



present new developments on a 1D lattice model with fractional power-law interatomic interaction defined by fractional values of the parameter s and nonlinear self-interaction potential for classical and quantum mechanical consideration. Also the authors of [45] showed that the dynamics on the 1D lattice can be equivalently presented by the corresponding fractional nonlinear equation in the long-wave limit and they concentrated on the conditions of such equivalence, type of the equations with fractional derivatives and some related properties. As examples of their developments, fractional sine-Gordon and wave-Hilbert nonlinear equations have been found for classical lattice dynamics. Also fractional nonlinear Schrödinger and nonlinear Hilbert–Schrodinger equations have been obtained for quantum lattice dynamics in the long-wave limit. The study of [83] shows that the list of possible applications of fractional equations can be naturally expanded to include non-chaotic and non-random dynamics as well. Authors of [32] applied the homotopy perturbation method for finding approximate analytical solutions of the fractional nonlinear Klein–Gordon equation as a numerical algorithm. The main aim of [39] is to use a fractional subequation method to construct the exact analytical solutions of the space-time fractional Cahn–Hilliard and the fractional nonlinear Klein–Gordon equations. Secondly, in Eq. (1.2), if we set FðuÞ ¼  r 1 u  r 2 u2  r 3 u3 then we get time fractional nonlinear klein–Gordon equation. As mentioned in [28] similar results about fractional Klein–Gordon-type equations have been recently discussed in [29], where an application to the fractional telegraph-type processes has been investigated. The main aim of [28] is to present some exact results related to the fractional Klein–Gordon equation involving fractional powers of the D'Alembert operator also, the authors of the mentioned paper found an exact analytic solution by using the McBride theory of fractional powers of hyper-Bessel operators. A wavelet method for a class of spacetime fractional Klein–Gordon equations with constant coefficients is proposed in [36], by combining the Haar wavelet and operational matrix together and efficaciously dispersing the coefficients. The homotopy analysis method is used to construct an approximate solution for the nonlinear space-time fractional derivatives Klein–Gordon equation in [30]. Kurulay [43] proposed the homotopy analysis method to obtain the solution of nonlinear fractional Klein–Gordon equation [61]. Authors of [85] presented a high order finite difference scheme for solving the two-dimensional nonlinear fractional Klein–Gordon equation subject to Neumann boundary conditions. Authors of [12] applied the homotopy analysis method to solve linear fractional problems such as fractional wave, Burgers, Korteweg-de Vries (KdV), KdV-Burgers, and Klein-Gordon equations with initial conditions.

1.1. A brief review of the meshless method In recent years radial basis functions (RBFs) have been extensively used in different context and emerged as a potential alternative in the field of numerical solution of PDEs. The use of RBFs in the numerical solution of partial differential equations (PDEs) has gained

414

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

popularity [11] in engineering and science community as it is meshless and can readily be extended to multi-dimensional problems. The key idea of the meshless methods is that they can obtain accurate and stable solution of integral equations or partial differential equations with various boundary conditions with a set of particles without using any mesh [59]. The most important advantages of meshless methods compared to finite element methods are the following: their high-order continuous shape functions, simpler incorporation of h- and p-adaptivity and certain advantages in crack problems. Wendland [87] combined the theory of radial basis functions with the field of Galerkin methods to solve partial differential equations. Also, he showed convergence and derived error estimates for smooth problems in arbitrary dimensions. The main aim of [27] is to study meshless collocation methods using radial basis functions to approximate regular solutions of systems of equations with linear differential or integral operators. Author of [88] considered error estimates for interpolation by a special class of compactly supported radial basis functions. Wu and Shaback [90] introduced a suitable variational formulation for the local error of scattered data interpolation by radial basis functions ΦðrÞ, the error can be bounded by a term depending on the Fourier transform of the interpolated function f and a certain ‘Kriging function’, which allows a formulation as an integral involving the Fourier transform of Φ. Authors of [22] combined the domain decomposition method with the so-called meshless Galerkin method for solving PDE using RBF with convergence analysis. Cai [6] used the theory of radial basis functions together with Galerkin methods to deal with the partial differential equations with Dirichlet boundary conditions. Also we refer the interested reader to [14,17,77,78,82] for applications of the meshless method of radial basis functions for the numerical solution of various partial differential equations. Recently, many fractional partial differential equations are solved using meshless approach based on the radial basis functions and moving least squares (MLS) technique. In [33] authors presented an implicit meshless collocation technique for solving time-fractional diffusion equation. Also, the stability and convergence of this meshless method are investigated theoretically and numerically. Authors of [34] presented an implicit meshless technique based on the radial basis functions for the numerical simulation of the anomalous subdiffusion equation. Also, they have discussed the stability and convergence of their method. Authors of [50] presented an implicit meshless method based on the radial basis functions for the numerical simulation of time-fractional diffusion equation. Authors of [95] presented an implicit meshless approach based on the moving least squares (MLS) approximation for the numerical simulation of fractional advectiondiffusion equation. Authors of [60] proposed a numerical method for the solution of the time-fractional nonlinear Schrödinger equation in one and two dimensions which appears in quantum mechanics. Liu et al. [48] for solving the time fractional mobile/immobile transport model proposed a finite difference method to discretize the time variable and obtained a semi-discrete scheme. Then they discussed its stability and convergence and finally considered a meshless method based on radial basis functions (RBFs) to discretize the space variable. The main aim of [46] is to study the high order finite difference method for the reaction and anomalous-diffusion equations. Some new numerical approaches based on piecewise interpolation for fractional calculus, and some new improved approaches based on the Simpson's rule for numerical integration for solving the fractional differential equations are proposed in [47]. The meshless method has already proved successful in standard quantum mechanics as well as for several other engineering and physical problems [1,14–17,51,82]. The aim of the current paper is to show that the meshless method based on the radial basis functions and collocation approach is also suitable for the treatment of the nonlinear partial differential equations. We refer the reader to [15,16,18,25,74,75,81] for other type of meshless techniques and their applications in solving partial differential equations or to [59] for integral equations. 1.2. The main aim of this paper In this paper we apply the meshless collocation method based on radial basis functions for solving Eq. (1.2). Firstly, we obtain a time semi-discretization scheme then full discretization is introduced using radial basis functions (RBFs) collocation method. The outline of this paper is as follows. In Section 2 we express some basis concepts for radial basis functions. In Section 3, we obtain a time discrete scheme to the mentioned equation. In this section, we also prove the unconditional stability and convergence of the time discrete scheme using the energy method and prove that computational order of time discrete scheme is OðτÞ. The proposed numerical scheme for the one and two dimensional time fractional sine-Gordon equations based on meshless method is given in Section 4. In Section 5 we report the numerical experiments of solving the mentioned equation with the proposed method for some test problems. Finally a conclusion is given in Section 6.

2. Basic concepts for RBFs approximation method As mentioned in [49] the definition of a meshfree method is given as follows: A meshfree method is a method used to establish system of algebraic equations for the whole problem domain without the use of a predefined mesh for the domain discretization. Also, as said in [49] meshfree methods use a set of nodes scattered within the problem domain as well as sets of nodes scattered on the boundaries of the domain to represent (not discretize) the problem domain and its boundaries. These sets of scattered nodes are called field nodes, and they do not form a mesh, meaning it does not require any a priori information on the relationship between the nodes for the interpolation or information of the unknown functions of field variables. In this paper, we use the meshfree method based on RBFs collocation approach. The reason we use the RBFs collocation method is that it works for arbitrary geometry with high dimensions and it does not require a mesh at all. The meshfree method using RBFs is the so-called Kansa's method [40–42], where the RBFs are directly implemented for the approximation of the solution of PDEs. Kansa's method was developed in 1990, in which the concept of solving PDEs by using RBFs, especially MQ, was initiated. As mentioned in [84], the MQ approximation scheme was first introduced by Hardy [35] who successfully applied this method for approximating surface and bodies from field data. In this section we introduce the basis definitions of radial basis functions in general case and we express some basic theorems for the interpolation error using radial basis functions. Definition 4 (Wendland [86]). A symmetric function ϕ A Rd  Rd ⟶R is strictly conditionally positive definite of order m, if for all sets X ¼ fx1 ; …; xN g  Rd of distinct points, and all vectors λ A Rd satisfying ∑N i ¼ 1 λi pðxi Þ ¼ 0 for any polynomial p of degree at most m 1, the

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

415

quadratic form N

N

λT Aλ ¼ ∑ ∑ λi λj ϕðxi  xj Þ; i¼1j¼1

is positive, whenever λ a 0. We interpolate a continuous function f : Rd ⟶R on a set X ¼ fx1 ; …; xN g with choosing the radial basis function ϕ : Rd ⟶R that is radial in the sense that ϕðxÞ ¼ Ψ ð J x J Þ, where J : J is the usual Euclidean norm on Rd as we will explain in the next section. Now, we assume ϕ to be strictly conditionally definite of order m, then the interpolation function has the following form [55]: N

l

i¼1

j¼1

I f ðxÞ ¼ ∑ λi ϕðx  xi Þ þ ∑ γ j pj ðxÞ; d þ m  1

where l ¼ coefficients

and fp1 ; p2 ; …; pl g is a basis for polynomial space of degree at most m i.e. ⊓dm . The basic problem is to find N þl unknown and γj in which N interpolation conditions are given in the following form:

m1

λi

I f ðxi Þ ¼ f i ;

i ¼ 1; …; N;

and for l remaining conditions we use the following equations: N

∑ λi pj ðxi Þ ¼ 0;

i¼1

1 r j rl:

Also, as mentioned in [63], the exponential convergence proof in applying RBFs in Sobolov space was given by Yoon [91,92], spectral convergence of the method in the limit of flat RBFs was given by Fornberg et al. [26]. Some popular choices of RBFs [78] are listed in the following 1. Linear radial basis function:

φðrÞ ¼ r; 2. Cubic radial basis function:

φðrÞ ¼ r 3 ; 3. Multiquadratics (MQ) radial basis function: pffiffiffiffiffiffiffiffiffiffiffiffiffiffi φðrÞ ¼ r 2 þ c2 ; 4. Gaussian (GS) radial basis function:

φðrÞ ¼ expð  cr2 Þ; 5. Polyharmonic splines radial basis function:

φðrÞ ¼ r 2n lnðrÞ;

where the free parameter c is called the shape parameter [77] of the RBFs. As mentioned in [72] the accuracy of many schemes for interpolating scattered data with radial basis functions depends on a shape parameter, c, of the radial basis function. Author of [72] showed, numerically, that the optimal value of c depends on the number and distribution of data points, on the data vector, and on the precision of the computation . Also this author presented an algorithm for selecting a good value for c that implicitly takes all the above considerations into account. Authors of [37] showed, numerically, that RBF in fact performs better than polynomials, as the optimal shape parameter exists not in the limit, but at a finite value. Also we refer the interested reader to [9,10,19,38] for various research works on RBFs.

3. Time discrete scheme Firstly, we define the functional spaces endowed with standard norms and inner products. 3.1. Preliminary of functional analysis Let Ω denote a bounded and open domain in Rd, for d ¼2 or 3 and let dx be the Lebesgue measure on Rd. For p o þ 1, we denote by R Lp ðΩÞ the space of the measurable functions u : Ω-R such that Ω juðxÞjp dx r þ 1. It is a Banach space for the norm [5] Z 1=p J u J Lp ðΩÞ ¼ juðxÞjp dx : Ω

416

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

The space L2 ðΩÞ is a Hilbert space with the inner product Z ðu; vÞ ¼ uðxÞvðxÞ dx; Ω

which induces the norm Z 1=2 juðxÞj2 dx : J u J L2 ðΩÞ ¼ Ω

Ω is an open domain in Rd, α ¼ ðα1 ; …; αd Þ is a d-tuple of nonnegative integers and jαj ¼ ∑pi¼ 1 αi , and also we put

Also we assume that Dα v ¼

∂ jα j v α : ∂x1 ∂xα2 2 …∂xd p α1

We define

n o H m ðΩÞ ¼ v A L2 ðΩÞ; Dα v A L2 ðΩÞ for all jαjr m :

This is a Hilbert space for the inner product Z ðu; vÞm ¼ ∑ Dα uðxÞDα vðxÞ dx; jα j r m Ω

which induces the norm J u J Hm ðΩÞ ¼

!

∑ ‖Dα u‖2L2 ðΩÞ :

jα j r m

The Sobolev space W 1;p ðIÞ is defined to be  Z Z W 1;p ðIÞ ¼ u A Lp ðIÞ; ( g A Lp ðIÞ : uφ0 ¼ g φ0 ; 8 φ A C 1 ðIÞ : I

I

In the following corollary we express the Poincare's inequality. Corollary 1 (Brezis [5]). Suppose that 1 r p o 1 and that J u J Lp ðΩÞ r C Ω J ∇u J Lp ðΩÞ ;

Ω is a bounded open set. Then there exists a constant C Ω (depending on Ω and p) such

8 u A W 1;p 0 ðΩÞ:

3.2. Analysis of the time discrete scheme For discretization of time variable, we need some preliminary. We define t k ¼ kτ ;

k ¼ 0; 1; …; N;

where τ ¼T/N is the step size of time variable. We introduce the following notations: un  ð1=2Þ ¼

 1 n u þun  1 ; 2

δt un  ð1=2Þ ¼

1

τ

 un  un  1 ;

where un ¼ uðx; t n Þ. In this section, we discrete the time variable using the following lemmas. Lemma 1 (Sun and Wu [79]). Suppose gðtÞ A C 2 ½0; t n  and 1 o α o 2 then  " # " # Z tn 0 n1   g ðtÞ 1 1 2  α 23  α 1α   þ dt  a gðt Þ  ∑ ða  a Þgðt Þ  a gðt Þ r  ð1 þ 2 Þ maxjg″ðtÞj τ3  α ; n1 0   α1 τ 0 n k ¼ 1 nk1 nk k 3α 0 r t r tn  0 ðt n  tÞ  2  α 12 where ak ¼

τ2  α h

2α

ðk þ 1Þ2  α  k

2α

i

:

Lemma 2 (Sun and Wu [79]). Let 0 o β o 1 and am ¼ ðm þ1Þ1  β  m1  β ; m ¼ 0; 1; …, then 1 ¼ a0 4 a1 4 ⋯ 4 am -0

as m-1:

Let vðx; tÞ ¼

∂uðx; tÞ ; ∂t

ð3:1Þ

and wðx; tÞ ¼

1 Γ ð2  αÞ

Z 0

t

∂vðx; sÞ ds : ∂s ðt sÞα  1

ð3:2Þ

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

417

Using Taylor expansion, it follows from (3.1) and (1.2) that vn  ð1=2Þ ¼ δt un  ð1=2Þ þ ðr 1 Þn  ð1=2Þ ; wn  ð1=2Þ ¼ ∇2 un  ð1=2Þ þ f

n  ð1=2Þ

 Fðun  1 Þ þ ðr 2 Þn  ð1=2Þ ;

ð3:3Þ x A Ω;

n Z 1;

ð3:4Þ

and there exists a constant c1 such that jðr 1 Þn  ð1=2Þ jr c1 τ2 ;

jðr 2 Þn  ð1=2Þ j r c1 ðτÞ:

ð3:5Þ

We define the operator

 n1 P un  ð1=2Þ ; q ¼ a0 un  ð1=2Þ  ∑ ðan  k  1  an  k Þuk  ð1=2Þ  an  1 q; k¼1

so from Lemma 1 we have Z tn 1 ∂vðx; tÞ dt 1 1  n  wn ¼ ¼ P v ; ψ þ Oðτ3  α Þ: ∂t ðt n  tÞα  1 Γ ð2  αÞ τ Γ ð2  αÞ 0 From Lemma 1 and v0 ¼ vðx; 0Þ ¼ ψ ðxÞ ¼ ψ , the above relation can be written as

  1 1 P vn  ð1=2Þ ; ψ þ Oðτ3  α Þ: wn  ð1=2Þ ¼ wn þ wn  1 ¼ 2 Γ ð2  αÞτ Consequently we can write

 1 P vn  ð1=2Þ ; ψ þ ðr 3 Þn  ð1=2Þ ; wn  ð1=2Þ ¼ Γ ð2  αÞτ

ð3:6Þ

and there exists a constant c2 such that  jðr 3 Þn  ð1=2Þ  rc2 τ3  α :

ð3:7Þ

If we substitute (3.3) into (3.6), we obtain



 1 1 P δt un  ð1=2Þ ; ψ þ P ðr 1 Þn  ð1=2Þ ; 0 þ ðr 3 Þn  ð1=2Þ : wn  ð1=2Þ ¼ Γ ð2  αÞτ Γ ð2  αÞτ Substituting the above result into (3.4) yields

 1 n  ð1=2Þ P δt un  ð1=2Þ ; ψ ¼ ∇2 un  ð1=2Þ þ f  Fðun  1 Þ þ Rn  ð1=2Þ ; Γ ð2  αÞτ where Rn  ð1=2Þ ¼



1 r ir M  1;

n Z 1;

ð3:8Þ



 1 P ðr 1 Þn  ð1=2Þ ; 0 þ ðr 3 Þn  ð1=2Þ þ ðr 2 Þn  ð1=2Þ : Γ ð2  αÞτ

According to Lemma 1, 2 and inequalities (3.5) we can write ( " # ) n1 1 a0 c1 τ2 þ ∑ ðan  k  1 an  k Þc1 τ2 þ c2 τ3  α þ c1 τ jRn  ð1=2Þ jr Γ ð2  αÞτ k¼1 

1 r a0 c1 τ2 þ ða0  an  1 Þc1 τ2 þ c2 τ3  α þc1 τ Γ ð2  αÞτ    2α 

1 1 2τ c1 τ 2 þ c2 τ 3  α þ c1 τ r a0 c1 τ2 þa0 c1 τ2 þc2 τ3  α þ c1 τ r Γ ð2  αÞτ Γ ð2  αÞτ ð2  αÞ 

r

  2c1 τ 3  α 2c1 2c1 þ c 2 τ 3  α þ c1 τ r þ c2 τ 3  α þ c1 τ r þ c2 þ c1 τ ; ð2  αÞΓ ð2  αÞ ð2  αÞΓ ð2  αÞ ð2  αÞΓ ð2  αÞ

then we have jRn  ð1=2Þ j r C τ; where C¼



ð3:9Þ

2c1 þ c2 þ c1 : ð2  αÞΓ ð2  αÞ

Now by omitting the small term Rn  ð1=2Þ , we can construct the following difference scheme for the equation (3.8): ( ) 8 n1 1 > 1 n  ð1=2Þ n  ð1=2Þ k  > 2 >  ∑ ðan  k  1  an  k Þδt U  an  1 ψ ¼ ∇2 U n  ð1=2Þ þ f  FðU n  1 Þ; > < Γ ð2  αÞτ a0 δt U k¼1

> x A Ω; n Z 1; > > > : n U ¼ hðx; t n Þ; x A ∂Ω:

ð3:10Þ

418

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

Lemma 3 (Sun and Wu [79]). For any G ¼ fG1 ; G2 ; …g and q, we obtain " # N n1 t1  α N t2  α 2 q ; ∑ a0 Gn  ∑ ðan  k  1  an  k Þ Gk an  1 q Gn Z N τ ∑ G2n  N 2 n¼1 2ð2  αÞ n¼1 k¼1 N ¼ 1; 2; …; in which al for all l ¼ 0; 1; 2; … are defined in Lemma 2. Proposition 1 (Quarteroni and Valli [71]). Assume that an is a non-negative sequence, and that the sequence bn satisfies 8 > < b0 r c0 ; n1

n1

k¼0

k¼0

> : bn r c0 þ ∑ dk þ ∑ ak bk ; then bn satisfies 8 b r c0 ð1 þ a0 Þ þd0 ; > < 1

n Z 1;

n1

n2

n1

k¼0

k¼0

r ¼ kþ1

> : bn r c0 ∏ ð1 þak Þ þ ∑ dk ∏ ð1 þ ar Þ þ dn  1 ;

n Z 2:

Moreover, if c0 Z 0 and d0 Z 0, it follows ! ! n1

bn r c0 þ ∑ dk

n1

exp

k¼0

∑ ak ;

k¼0

n Z1:

n n Theorem 5. Suppose Un and U~ be the exact and approximate solutions of Eq. (3.10), respectively, where as U n ; U~ A H 10 . Then the time discrete scheme (3.10) is unconditionally stable and we have the following inequality:

J en J L2 ðΩÞ r C J ∇e0 J L2 ðΩÞ ; in which en ¼ U n  U~ . n

Proof. Firstly, we obtain the following roundoff error equation: ( ) h i n1 1 n1 a0 δt en  ð1=2Þ  ∑ ðan  k  1  an  k Þ δt ek  ð1=2Þ ¼ ∇2 en  ð1=2Þ  FðU n  1 Þ  FðU~ Þ : τΓ ð2  αÞ k¼1 Multiplying the above equation by δt en  ð1=2Þ and integrating on Ω, result ( )

 n1

 1 a0 δt en  ð1=2Þ ; δt en  ð1=2Þ  ∑ ðan  k  1  an  k Þ δt ek  ð1=2Þ ; δt en  ð1=2Þ τΓ ð2  αÞ k¼1

 h  i n1 2 n  ð1=2Þ n  ð1=2Þ  FðU n  1 Þ  FðU~ ¼ ∇ e ; δt e Þ ; δt en  ð1=2Þ : Then we can write ( ) n1 1 a0 ‖δt en  ð1=2Þ ‖2L2 ðΩÞ  ∑ ðan  k  1  an  k Þ‖δt ek  ð1=2Þ ‖L2 ðΩÞ ‖δt en  ð1=2Þ ‖L2 ðΩÞ τΓ ð2  αÞ k¼1

 h  i n1 n  ð1=2Þ n  ð1=2Þ  FðU n  1 Þ  FðU~ ; ∇δ t e Þ ; δt en  ð1=2Þ : r  ∇e

ð3:11Þ

Also, we have

 Z



 Z ∇en þ ∇en  1 ∇en  ∇en  1  ∇en  ð1=2Þ ; ∇δt en  ð1=2Þ ¼ ∇en  ð1=2Þ ∇δt en  ð1=2Þ ¼ 2 τ Ω Ω Z h i  n 2  n  1 2 i 1 h 1 n 2 n1 2 ¼ ∇e  ∇e ‖∇e ‖L2 ðΩÞ  ‖∇e ‖L2 ðΩÞ : ¼ 2τ Ω 2τ Summing (3.11) from n¼ 1 to m and using the above equality, we obtain the following relation: ( ) m n1 1 n  ð1=2Þ k  12 ∑ a ‖δ e ‖L2 ðΩÞ  ∑ ðan  k  1  an  k Þ‖δt e ‖L2 ðΩÞ ‖δt en  ð1=2Þ ‖L2 ðΩÞ τΓ ð2  αÞ n ¼ 1 0 t k¼1 i i  m h 1h n1 þ ‖∇em ‖2L2 ðΩÞ  ‖∇e0 ‖2L2 ðΩÞ r  ∑ FðU n  1 Þ  FðU~ Þ ; δt en  ð1=2Þ : 2τ n¼1 We have m

 ∑

n¼1





n1   F U n  1 ; x; t n  1 F U~ ; x; t n  1 ; δt en  ð1=2Þ m

rL ∑

n¼1

     

 m   n  1 ~ n  1   n  ð1=2Þ      U U ; δt e  ¼ L ∑ en  1 ; δt en  ð1=2Þ  n¼1

m L2 Γ ð2  αÞ m t 1m α ∑ ‖δt en  ð1=2Þ ‖2L2 ðΩÞ : ∑ ‖en  1 ‖2L2 ðΩÞ þ r 2Γ ð2  α Þ n ¼ 1 2t 1m α n ¼ 1

ð3:12Þ

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

419

Using Lemma 3 and the above relations we can obtain i m t 1m α 1h ∑ ‖δt en  ð1=2Þ ‖2L2 ðΩÞ þ ‖∇em ‖2L2 ðΩÞ  ‖∇e0 ‖2L2 ðΩÞ 2τ 2Γ ð2  αÞ n ¼ 1

r

m L2 Γ ð2  αÞ m t 1m α ∑ ‖δt en  ð1=2Þ ‖2L2 ðΩÞ : ∑ ‖en  1 ‖2L2 ðΩÞ þ 1α 2Γ ð2  α Þ n ¼ 1 2t m n¼1

Simplifying, changing index from m to n and multiplying 2τ on both sides of the above relation, give n1

‖∇en ‖2L2 ðΩÞ r ‖∇e0 ‖2L2 ðΩÞ þ τL2 t αn  1 Γ ð2  αÞ ∑ ‖ej ‖2L2 ðΩÞ :

ð3:13Þ

j¼0

Using the Poincare inequality [71] ‖en ‖L2 ðΩÞ r C Ω ‖∇en ‖L2 ðΩÞ ; we obtain n1

‖en ‖2L2 ðΩÞ r C 2Ω ‖∇e0 ‖2L2 ðΩÞ þ τC 2Ω L2 t αn  1 Γ ð2  αÞ ∑ ‖ej ‖2L2 ðΩÞ :

ð3:14Þ

j¼0

Considering the above relation and Proposition 1 with dk ¼ 0;

c0 ¼ C 2Ω ‖∇e0 ‖2L2 ðΩÞ ;

ak ¼ τC 2Ω L2 t αn  1 Γ ð2  αÞ;

bk ¼ ‖ek ‖2L2 ðΩÞ ;

we can write 8 ‖e0 ‖22 rC 2Ω ‖∇e0 ‖2L2 ðΩÞ ; > > < n1

2 2 2 α1 k 2 n 2 0 2 > > : ‖e ‖L2 ðΩÞ r C Ω ‖∇e ‖L2 ðΩÞ þ τC Ω L t n Γ ð2  αÞ ∑ ‖e ‖L2 ðΩÞ ; k¼0

so from Proposition 1 we can conclude n1

‖en ‖22 rC 2Ω ‖∇e0 ‖2L2 ðΩÞ ∏

r C Ω exp nτ 2

k¼0



 n1 1 þ τC 2Ω L2 t αn  1 Γ ð2  αÞ r C 2Ω ‖∇e0 ‖2L2 ðΩÞ ∏ exp τC 2Ω L2 t αn  1 Γ ð2  αÞ

C Ω L2 t αn  1 2

Γ ð2  α Þ



k¼0

‖∇e0 ‖2L2 ðΩÞ r C 2Ω



 exp TC 2Ω L2 t αn  1 Γ ð2  αÞ ‖∇e0 ‖2L2 ðΩÞ rC 2Ω exp T 2 C 2Ω L2 Γ ð2  αÞ ‖∇e0 ‖2L2 ðΩÞ :

Finally, we have ‖en ‖L2 ðΩÞ r C‖∇e0 ‖L2 ðΩÞ ;

 where C ¼ C Ω exp ðT 2 C 2Ω L2 Γ ð2  αÞÞ=2 ; which completes the proof.

□

Theorem 6. Let un ¼ uðx; t n Þ A H 10 be the exact solution of Eq. (3.8) and Un be the approximate solution obtained by scheme (3.10), then the time discrete scheme (3.10) is convergent with the convergence order OðτÞ. Proof. We set

ρn ¼ un U n ;

n Z 1;

where ρ ¼ 0. Subtracting (3.10) from (3.8), we obtain the following roundoff error equation: ( ) h i n1 1 n  ð1=2Þ k  1=2 a δρ  ∑ ðan  k  1  an  k Þ δt ρ ¼ ∇2 ρn  ð1=2Þ  Fðun  1 Þ  FðU n  1 Þ þ Rn  ð1=2Þ : τΓ ð2  αÞ 0 t k¼1 0

Multiplying the above equation by δt ρn  ð1=2Þ and integrating on Ω, result ( )

 n1

 1 a0 δt ρn  ð1=2Þ ; δt ρn  ð1=2Þ  ∑ ðan  k  1  an  k Þ δt ρk  ð1=2Þ ; δt ρn  ð1=2Þ τΓ ð2  αÞ k¼1





 ¼ ∇2 ρn  ð1=2Þ ; δt ρn  ð1=2Þ  Fðun  1 Þ  FðU n  1 Þ; δt ρn  ð1=2Þ þ Rn  ð1=2Þ ; δt ρn  ð1=2Þ ; i.e.

( ) n1 1 n  ð1=2Þ k  ð1=2Þ a ‖δ ρ ‖L2 ðΩÞ  ∑ ðan  k  1 an  k Þ‖δt ρ ‖L2 ðΩÞ ‖δt ρn  ð1=2Þ ‖L2 ðΩÞ τΓ ð2  αÞ 0 t k¼1



 þ ∇ρn  ð1=2Þ ; ∇δt ρn  ð1=2Þ r  Fðun  1 Þ  FðU n  1 Þ; δt ρn  ð1=2Þ þ ‖Rn  ð1=2Þ ‖L2 ðΩÞ ‖δt ρn  ð1=2Þ ‖L2 ðΩÞ :

Summing the above relation from n ¼1 to m we obtain ( ) m n1 1 ∑ a0 ‖δt ρn  ð1=2Þ ‖L2 ðΩÞ  ∑ ðan  k  1  an  k Þ‖δt ρk  ð1=2Þ ‖L2 ðΩÞ ‖δt ρn  ð1=2Þ ‖L2 ðΩÞ τΓ ð2  αÞ n ¼ 1 k¼1  m

m 1 þ ‖∇ρm ‖2L2 ðΩÞ r  ∑ Fðun  1 Þ  FðU n  1 Þ; δt ρn  ð1=2Þ þ ∑ ‖Rn  ð1=2Þ ‖L2 ðΩÞ ‖δt ρn  ð1=2Þ ‖L2 ðΩÞ : 2τ n¼1 n¼1

420

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

Using Lemma 3 and a relation similar to inequality (3.12) yield m t 1m α 1 L2 Γ ð2  αÞ m ∑ ‖δt ρn  ð1=2Þ ‖2L2 ðΩÞ þ ‖∇ρm ‖2L2 ðΩÞ r ∑ ‖ρn  1 ‖2L2 ðΩÞ 2τ 2Γ ð2  αÞ n ¼ 1 t 1m α n¼1

þ

m m t 1m α t 1m α Γ ð2  αÞ m ∑ ‖δt ρn  ð1=2Þ ‖2L2 ðΩÞ þ ∑ ‖δt ρn  ð1=2Þ ‖2L2 ðΩÞ þ 1  α ∑ ‖Rn  ð1=2Þ ‖2L2 ðΩÞ : 4Γ ð2  α Þ n ¼ 1 4Γ ð2  αÞ n ¼ 1 tm n¼1

So simplifying, changing index from m to n, multiplying both sides of the above relation by 2τ and using inequality (3.9) we can write n1

n

j¼0

j¼1

‖∇ρn ‖2L2 ðΩÞ r 2τL2 t αn  1 Γ ð2  αÞ ∑ ‖ρ j ‖2L2 ðΩÞ þ 2τt αn  1 Γ ð2  αÞ ∑ ‖Rj  ð1=2Þ ‖2L2 ðΩÞ :

ð3:15Þ

Using the Poincare inequality [71] ‖ρn ‖L2 ðΩÞ r C Ω ‖∇ρn ‖L2 ðΩÞ ; then inequality (3.15) yields n1

‖ρn ‖2L2 ðΩÞ r τC 2Ω L2 t αn  1 Γ ð2  αÞ ∑ ‖ρ j ‖2L2 ðΩÞ þ TC 2 C 2Ω t αn  1 Γ ð2  αÞτ2 :

ð3:16Þ

j¼0

Now, from Proposition 1 we have ‖ρ

n 2 ‖L2 ðΩÞ r TC 2 C 2Ω t αn  1

n1

Γ ð2  αÞτ exp ∑ τ 2

j¼0

! C Ω L2 t αn  1 2

Γ ð2  αÞ

 r TC 2 C 2Ω t αn  1 Γ ð2  αÞτ2 exp TC 2Ω L2 t αn  1 Γ ð2  αÞ ¼ CðT; α; C Ω Þτ2 ;

as a result of the above inequality we have ‖ρn ‖L2 ðΩÞ r CðT; α; C Ω Þτ; which completes the proof.

□

4. Numerical scheme Now we present our numerical technique for both one and two-dimensional versions of Eqs. (1.2)–(1.4). 4.1. The one-dimensional case In this case, we consider the approximate expansion of uðxi ; t n Þ which is as follows: M

uðxi ; t n Þ ¼ ∑ cnj φðr ij Þ þcnM þ 1 xi þ cnM þ 2 ;

ð4:1Þ

j¼1

where

φðr ij Þ ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxi  xj Þ þ c2 ¼ r 2 þc2 :

Besides M equations resulting from collocating (4.1) at M points, additional two equations are required by the following regularization conditions: M

M

j¼1

j¼1

∑ cnj ¼ ∑ cnj xj ¼ 0:

For positive integer numbers M and N, let τ ¼ T=N denote the step size of time variable, t, and h ¼ ðb aÞ=M denote the step size of space variable. So we define xj ¼ a þ jh; t k ¼ kτ ;

j ¼ 0; 1; 2; …; M;

k ¼ 0; 1; 2; …; N:

The exact and approximate solutions at the point ðxj ; t k Þ are denoted by ukj and U kj ; respectively. In this case, the Ω domain is ½a; b and Eqs. (1.2)–(1.4) change to the following equations: ∂α uðx; tÞ ∂2 uðx; tÞ ¼  Fðuðx; tÞÞ þ f ðx; tÞ; ∂t α ∂x2

x A ½a; b;

t A ½0; T;

ð4:2Þ

with the initial conditions uðx; 0Þ ¼ ωðxÞ;  ∂uðx; tÞ  ¼ ψ ðxÞ; ∂t t ¼ 0

x A ½a; b;

ð4:3Þ

and boundary conditions uða; tÞ ¼ h1 ða; tÞ;

uðb; tÞ ¼ h2 ðb; tÞ;

t 4 0:

ð4:4Þ

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

Considering the initial and boundary conditions, (4.3) and (4.4), respectively, we have 8 2 1 2 > > > μ U 1  1 ∂ U ¼ μ ωðxÞ þ 1 ∂ ωðxÞ þ μb0 ψ ðxÞ FðωðxÞÞ þ f 1=2 ; n ¼ 1; > > 2 > 2 ∂x 2 ∂x2 τ τ > > > 2 n > μ 1 ∂ U μ 1 ∂2 U n  1 μ n  1 > > < Un  ¼ Un  1 þ þ ∑ ðbn  k  1  bn  k ÞU k 2 2 ∂x 2 ∂x2 τ τ τ k¼1 > > n1 > μ μ k1 > > > >  τ ðbn  2  bn  1 ÞωðxÞ  τ ∑ ðbn  k  1  bn  k ÞU > k¼2 > > > > n  ð1=2Þ : þ μbn  1 ψ ðxÞ  FðU n  1 Þ þf ; n ¼ 2; 3; …; N;

421

ð4:5Þ

where

μ¼

a0

Γ ð2  αÞτ

;

α

bk ¼ ðk þ 1Þ2  α k :

Also substituting Eq. (4.1) in (4.5) we obtain the matrix form Acn ¼ Bn ; in which 2

φðr1;1 Þ

6 6 Lðφðr 2;1 ÞÞ 6 6 6 ⋮ 6 6 A ¼ 6 Lðφðr M  1;1 ÞÞ 6 6 6 φðr M;1 Þ 6 6 x1 4 1

φðr 1;2 Þ

⋯

φðr 1;M  1 Þ

φðr 1;M Þ

Lðφðr 2;2 ÞÞ

⋯

Lðφðr 2;M  1 ÞÞ

Lðφðr 2;M ÞÞ

⋮

⋱

Lðφðr M  1;2 ÞÞ

⋯

Lðφðr M  1;M  1 ÞÞ

Lðφðr M  1;M ÞÞ

φðrM;2 Þ

⋯

φðrM;M  1 Þ

φðrM;M Þ

x2 1

⋯ 1

xM  1 1

xM 1

x1

μ x τ 2

1

μ τ

μ μ x τ M1 τ xM

1

0 0

0 0

3 7 7 7 7 7 7 7 7; 7 7 7 7 7 5

where the operator L is

μ τ

Lðφðr ij ÞÞ ¼ φðr ij Þ 

1 d φðr ij Þ ; 2 dx2 2

and

T cn ¼ cn1 ; cn2 ; …; cnM ; cnM þ 1 ; cnM þ 2 ;

n n

n Bn ¼ b1 ; b2 ; …; bM ; 0; 0 ;

in which 8 μ 1 ∂2 ωðxi Þ > ð1=2Þ > þ μb0 ψ ðxi Þ  Fðωðxi ÞÞ þ f i ; n ¼ 1; 1 oi o M; > ωðxi Þ þ > > 2 ∂x2 τ > > > M 2 > M M n  1 ∂ φj ðxi Þ μ > μ 1 > n1 n1 > þ ∑ ∑ ðbn  k  1  bn  k Þckj φj ðxi Þ > ∑ cj φj ðxi Þ þ ∑ cj > 2j¼1 τ j¼1 τ j ¼ 1k ¼ 1 ∂x2 > > > > > > < μ μ M n1 n ∑ ðbn  k  1 bn  k Þckj  1 φj ðxi Þ bi ¼  τ ðbn  2  bn  1 Þωðxi Þ  τ j ∑ ¼ 1k ¼ 2 > > ! > > > M > > > þ μbn  1 ψ ðxi Þ  F ∑ cnj  1 φj ðxi Þ þ f n  ð1=2Þ ; n Z2; 1 o i o M; > i > > j¼1 > > > > > > h1 ðt n  ð1=2Þ Þ; n Z 1; and i ¼ 1; > > > : h ðt n Z 1; and i ¼ M: 2 n  ð1=2Þ Þ;

4.2. The two-dimensional case In this case, we assume the approximate expansion of uðxi ; t n Þ is as follows: M

uðxi ; t n Þ ¼ ∑ cnj φðr ij Þ;

ð4:6Þ

j¼1

in which

φðr ij Þ ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxi  xj Þ2 þ c2 ¼ r 2 þc2 :

In the case, we consider the following equations: ∂α uðx; tÞ ¼ Δuðx; tÞ  Fðuðx; tÞÞ þ f ðx; tÞ; ∂t α with the initial conditions  ∂uðx; tÞ uðx; 0Þ ¼ ωðxÞ; ¼ ψ ðxÞ; ∂t t ¼ 0

x A Ω;

x A Ω;

0 r t r T;

ð4:7Þ

ð4:8Þ

422

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

and boundary conditions uðx; tÞ ¼ hðx; tÞ;

x A ∂Ω:

ð4:9Þ

Now, in the two-dimensional case we assume x ¼ ðx; yÞ then for the use of the Kansa method, we let

fðxi ; yi ÞgM i¼1

be M collocation points in

M I Ω in which fðxi ; yi ÞgM i ¼ 1 are boundary points and fðxi ; yi Þgi ¼ MI þ 1 are interior points. We consider only Hardy's MQ. For each point fðxi ; yi Þg,

let us denote qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi φj ðx; yÞ ¼ ðx  xj Þ2 þ ðy  yj Þ þ c2 : By considering initial conditions (4.8) we have 8 μ 1 1 1 μ ð1=2Þ > > U  ΔU 1 ¼ ΔωðxÞ þ ωðxÞ þ μb0 ψ ðxÞ  FðωðxÞÞ þ f ; > > 2 2 τ τ > > > > n  1 > μ n 1 1 μ μ > n n1 > þ U n  1 þ ∑ ðbn  k  1  bn  k ÞU k < U  ΔU ¼ Δ U 2 2 τ τ τ k¼1 > n1 > μ μ > >  ðb  bn  1 Þ ωðxÞ  ∑ ðbn  k  1 bn  k ÞU k  1 > > > τ n2 τ k¼2 > > > > : þ μb ψ ðxÞ  FðU n  1 Þ þf n  ð1=2Þ ; n ¼ 2; 3; …; N:

n ¼ 1;

ð4:10Þ

n1

Also by substituting Eq. (4.6) in (4.10) we obtain the matrix form Acn ¼ Bn ; in which 2

φðr1;1 Þ φðr2;1 Þ

6 6 6 6 ⋮ 6 6 φðr 6 M I ;1 Þ A¼6 6 Lðφðr M þ 1;1 ÞÞ I 6 6 6 Lðφðr MI þ 2;1 ÞÞ 6 6 ⋮ 4 Lðφðr M;1 ÞÞ

φðr 1;2 Þ φðr 2;2 Þ ⋮

… ⋱

φðr1;MI þ 1 Þ φðr2;MI þ 1 Þ

⋯ ⋯ ⋯

φðr MI ;2 Þ ⋯ ⋱ φðrMI ;MI þ 1 Þ ⋯ Lðφðr MI þ 1;2 ÞÞ ⋯ Lðφðr MI þ 1;MI þ 1 ÞÞ ⋯ Lðφðr MI þ 2;2 ÞÞ ⋯ Lðφðr MI þ 2;MI þ 1 ÞÞ ⋱ ⋯ ⋮

⋯

Lðφðr M;2 ÞÞ

⋯

 1 Lðφðr ij ÞÞ ¼ φj ðxÞ  Δφj ðxÞx ¼ xi ; 2

T cn ¼ cn1 ; cn2 ; …; cnM ;

⋯

⋱ Lðφðr M;MI þ 1 ÞÞ

M I þ 1 ri r M;

⋯

φðr 1;M Þ φðr 2;M Þ

3

7 7 7 7 ⋮ 7 φðr MI ;M Þ 7 7 7; Lðφðr MI þ 1;M ÞÞ 7 7 7 Lðφðr MI þ 2;M ÞÞ 7 7 7 ⋮ 5 Lðφðr M;M ÞÞ

1 rj r M;

n n n Bn ¼ b1 ; b2 ; …; bM ;

where

8 μ 1 > > ωðxi Þ þ Δωðxi Þ þ μb0 ψ ðxi Þ  Fðωðxi ÞÞ þ f 1=2 ; n ¼ 1; M I þ 1 r ir M; > > i > 2 τ > > > M M M n  1 >μ 1 μ > > > ∑ cnj  1 φj ðxi Þ þ ∑ cnj  1 Δφj ðxi Þ þ ∑ ∑ ðbn  k  1  bn  k Þckj φj ðxi Þ > >τ j¼1 2j¼1 τ j ¼ 1k ¼ 1 > > > < M n1 n μ μ bi ¼  ðbn  2  bn  1 Þωðxi Þ  ∑ ∑ ðbn  k  1  bn  k Þck  1 φ ðxi Þ; j j > > τ τ j ¼ 1k ¼ 2 > > >   > > M > > > þ μbn  1 ψ ðxi Þ  F ∑ cnj  1 φj ðxi Þ þ f n  1=2 ; n Z2; M I þ 1 r i r M; > i > > j¼1 > > > > : hðxi ; t n  1=2 Þ; n Z 1 and 1 r i rM I :

5. Numerical results In this section we present the numerical results of the new method on several test problems. We tested the accuracy and stability of the new method described in this paper by performing the mentioned method for different values of h and τ. We performed our computations using Matlab 7 software on a Pentium IV, 2800 MHz CPU machine with 2 GB of memory. To illustrate the accuracy of the method, we compute the error norms L1 and RMS. vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u M  2 u u ∑ uðxj ; TÞ Uðxj ; TÞ tj ¼ 1 : L1 ¼ max jUðxj ; TÞ uðxj ; TÞj; RMS ¼ M 1rjrM1 Also we calculated the computational orders of the method presented in this paper (denoted by C-order) with the following formula [13]:

 log EE12

; log hh12 in which E1 and E2 are errors correspond to grids with mesh size h1 and h2, respectively.

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

423

The criterion of irregular domains O1 and O2 is r ¼ n12 ½1 þ 2n þ n2  ðn þ 1Þ cos ðnyÞ, as n ¼4 for O1 and n ¼8 for O2. Also, the criterion of rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi irregular domains O3 and O4 is r ¼0.3

cos ð2yÞ þ

2

1:1  sin ð2yÞ and r ¼ 1+cos(y)sin2(y), respectively.

5.1. Test problem 1 We consider the time fractional sine-Gordon with the following form: 8 α ∂ uðx; tÞ ∂2 uðx; tÞ > > > ¼  sin ðuðx; tÞÞ þ f ðx; tÞ; > > ∂t α ∂x2 > > > < φðxÞ ¼ 0; ψ ðxÞ ¼ 0; 0 o x o 1;

ϕ1 ðtÞ ¼ 0; ϕ2 ðtÞ ¼ t 2 sin ð1Þ; 0 r t r1; > > >   > > 2t 2  α > 2 > > sin ðxÞ þ sin ðt 2 sin ðxÞÞ; : f ðx; tÞ ¼ Γ ð2  αÞ þ t

The exact solution of this test problem is uðx; tÞ ¼ t 2 sin ðxÞ: We solve this problem with the method presented in this paper with several values of h, τ, α, c for a ¼0 and b¼1 at final time T ¼1. The L1 and RMS errors, C-order are shown in Tables 1 and 2. Table 1 shows that the proposed method produces results with good accuracy. Computational orders in Table 1 confirm the first order of accuracy in the time component. Fig. 1 presents the irregular domains considered. Fig. 2 shows the graphs of exact and approximate solutions (right panel) and absolute error (left panel) obtained with α ¼ 1.75, M ¼ 33, c ¼ 0.01, and N ¼ 100 in final time T ¼ 0.25 for Test problem 1. 5.2. Test problem 2 In this test, we consider two-dimensional time fractional sine-Gordon equation with the following form ∂α uðx; y; tÞ ∂2 uðx; y; tÞ ∂2 uðx; y; tÞ 2t 2  α sin ðx þ yÞ þ 2t 2 sin ðx þ yÞ þ sin ðt 2 sin ðx þ yÞÞ; ¼ þ  sin ðuðx; y; tÞ þ α 2 2 ∂t Γ ð3  αÞ ∂x ∂y with the initial conditions uðx; y; 0Þ ¼ 0;

∂uðx; y; 0Þ ¼ 0; ∂t

x; y A Ω;

and boundary condition uðx; y; tÞ ¼ t 2 sin ðx þ yÞ;

x; y A ∂Ω:

The exact solution of the above problem is uðx; y; tÞ ¼ t 2 sin ðx þ yÞ;

; x; y A Ω:

Now, we suppose that Ω ¼ ½0; 1  ½0; 1. We solve this problem with the method presented in this paper with several values of h, c, τ, α for L ¼1 at final time T ¼1. Table 3 shows errors and computational orders obtained for Test problem 2 with h ¼ 1=5 and c¼ 0.5 on interval Ω ¼ ½  1; 1  ½  1; 1. Table 1 Errors and computational orders obtained for Test problem 1 with h ¼ 1/50 and c¼ 0.1. τ

α ¼1.15

α¼ 1.85

L1

C-order

L1

C-order

1=10

1:3757  10  2

–

9:8260  10  3

–

1=20

6:9463  10  3

0.9858

4:5919  10  3

1.0975

1=40

3:4855  10  3

0.9949

2:2413  10  3

1.0347

1=80

1:7414  10  3

1.0011

1:1127  10  3

1.0103

1=160

8:6602  10  4

1.0077

5:5212  10  4

1.0110

Table 2 Errors obtained for Test problem 1 with τ ¼ 1/512 and c ¼ 0.1. h

α ¼ 1:5

α ¼ 1:75

RMS

L1

RMS

L1

1=4

1:3102  10  2

1:8898  10  2

1:3059  10  2

1:8849  10  2

1=8

3:1704  10  3

4:7116  10  3

3:0629  10  3

4:6463  10  3

4

4

4

9:3782  10  4

7:3824  10  5

1:1580  10  4

1=16

4:0916  10

1=32

1:0785  10  4

9:3517  10

1:6654  10  4

4:1832  10

424

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

Ω

1

1.5

1

1

0.5

0.5

0

0

y

y

Ω 1.5

−0.5

−0.5

−1

−1

−1.5 −1.5

−1

−0.5

0

0.5

1

−1.5 −1.5

1.5

−1

−0.5

0

x

Ω 0.2

2

0.5

1

1.5

x

Ω

3

1.5

0.15

4

1

0.1

0.5

0

y

y

0.05

−0.05

0 −0.5

−0.1

−1

−0.15 −0.2 −0.5

0

0.5

−1.5 −1

−0.5

x

0

0.5

1

1.5

x Fig. 1. The irregular domains considered.

Fig. 2. Graphs of exact and approximate solutions(right panel) and absolute error (left panel) obtained with α ¼ 1:75, M¼ 33, c¼ 0.01, and N ¼ 100 in final time T ¼ 0.25 for Test problem 1.

Fig. 3 demonstrates the graph of approximate solution and contour of resulting error with α ¼ 1:75, τ ¼ 1=320, c¼ 0.6 and h ¼ 1=10 for Test problem 2. Graph of approximate solution and contour of resulting error with α ¼ 1:50, τ ¼ 1=640, c ¼0.8 and h ¼ 1=5 are shown in Fig. 4 for Test problem 2. Fig. 5 presents graph of approximate solution and contour of resulting error on Ω1 domain with α ¼ 1:25,

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

425

Table 3 Errors and computational orders obtained for Test problem 2 with h ¼ 1/5 and c¼ 0.5 on interval Ω ¼ ½  1; 1  ½  1; 1. τ

α ¼ 1:25

α ¼ 1:75 L2

L1

L1

L2

1=10

1:6803  10  2

8:6063  10  2

1:1295  10  2

5:9453  10  2

1=20

8:4912  10

3

2

3

2:7347  10  2

1=40

4:2901  10  3

2:3370  10  2

2:2178  10  3

1:4099  10  2

3

3

3

8:7794  10  3

4:4257  10

4:8061  10

1=160

2:1780  10

1=160

1:6190  10  3

8:6396  10  3

1:4745  10  3

6:7635  10  3

1=320

3

3

3

6:0290  10  3

1:4288  10

6:7554  10

1:7179  10 1:3831  10

2e−005

e−0 05 1e−0 05

0.06

5e−006

−00

−

6

00

6

e−

−1e

−1

−0.04

− 5e

5

0 −5e

0

0

−00

1e−005 0 5e−006 0 −0.02

5e

−0

0.02

1e−005

06

1.5

0.04

y

1:3239  10

−005 −1.5e 5 −00 −2e

00 5

−0.06 −0.2

−0.1

0 x

0.1

0.2

Fig. 3. Graph of approximate solution and contour of resulting error with α¼ 1.75, τ ¼ 1/320, c ¼0.6 and h ¼1/10 for Test problem 2.

Fig. 4. Graph of approximate solution and contour of resulting error with α ¼1.50, τ¼ 1/640, c¼ 0.8 and h ¼ 1/5 for Test problem 2.

τ ¼ 1=320, c¼ 0.8 and h ¼ 1=5 for Test problem 2. Also, graph of approximate solution and contour of resulting error on Ω2 domain with α ¼ 1:80, τ ¼ 1=160, c¼ 0.8 and h ¼ 1=5 for Test problem 2 are graphed in Fig. 6. In Tables 4–7 errors and computational orders are reported for Test problem 2 on different domains. In these tables we conclude that the computational order of the new method is approximately OðτÞ in time variable. 5.3. Test problem 3 In this test, we consider two-dimensional nonlinear time fractional Klein–Gordon equation ∂α uðx; y; tÞ ∂2 uðx; y; tÞ ∂2 uðx; y; tÞ ¼ þ þ u2 ðx; y; tÞ þ u3 ðx; y; tÞ þf ðx; y; tÞ; ∂t α ∂x2 ∂y2

426

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

Fig. 5. Graphs of approximate solution and resulting error with α ¼1.25, τ ¼ 1/320, c ¼0.8 and h ¼1/5 for Test problem 2.

Fig. 6. Graphs of approximate solution and resulting error with α¼ 1.80, τ ¼1/160, c¼ 0.8 and h ¼ 1/5 for Test problem 2.

Table 4 Errors and computational orders obtained for Test problem 2 on Ω1 domain with h ¼ 1/5. τ

α ¼ 1:15

α ¼ 1:85

c

L1

C-order

L1

C-order

1=10

1:3694  10  2

–

1:3163  10  2

–

0.1

1=20

6:9271  10  3

0.9832

6:9308  10  3

0.9253

0.2

1=40

3:5123  10  3

0.9798

3:5133  10  3

0.9801

0.3

1=160

1:7612  10  3

0.9958

1:7721  10  3

0.9873

0.4

1=160

8:7695  10  4

1.0060

8:6004  10  4

1.0429

0.6

1=320

4:2924  10  4

1.0370

4:1727  10  4

1.0434

0.8

with the initial condition uðx; y; 0Þ ¼ 0; ∂uðx; y; 0Þ ¼ 0; ∂t

x; y A Ω;

and boundary conditions

 2 2 uðx; y; tÞ ¼ t 2 sec h ðx þ y  rÞ þ sec h ðx þ y þ rÞ ;

x; y A ∂Ω;

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

427

Table 5 Errors and computational orders obtained for Test problem 2 on Ω2 domain with h ¼ 1=10. τ

α ¼ 1:75

α ¼ 1:25

L1

c

C-order

L1

C-order

1=10

1:2595  10  2

–

1:2252  10  2

–

0.1

1=20

4:8950  10

3

1.4035

4:8959  10  3

1.3234

0.2

1=40

1:9909  10  3

1.2979

1:9873  10  3

1.3008

0.4

1=160

6:3005  10

4

1.6599

8:4773  10  4

1.2291

0.8

1=160

3:3702  10  4

0.9026

4:0400  10  4

1.0692

0.8

Table 6 Errors and computational orders obtained for Test problem 2 on Ω3 domain with h ¼ 1=5 and α ¼ 1:15. τ

α ¼ 1:25

α ¼ 1:75

c

L1

C-order

L1

C-order

1=10

1:4732  10  3

–

1:4671  10  3

–

0.1

1=20

8:9139  10  4

0.7248

8:9001  10  4

0.7211

0.2

1=40

3:2327  10  4

1.4633

3:2323  10  4

1.4613

0.3

1=160

1:2067  10  4

1.4217

1:2106  10  4

1.4168

0.4

1=160

5:1664  10  5

1.2238

5:1964  10  5

1.2201

0.5

1=320

2:4756  10  5

1.0614

2:4909  10  5

1.0608

0.6

Table 7 Errors and computational orders obtained for Test problem 2 on Ω4 domain with h ¼ 1=5. τ

α ¼ 1:25

α ¼ 1:50

L1

c

C-order

L1

C-order

1=10

6:1422  10  3

–

6:0519  10  3

–

0.5

1=20

3:1708  10

3

0.9539

3:2425  10  3

0.8963

0.5

1=40

1:6140  10  3

0.9742

1:6826  10  3

0.9464

0.5

1=160

8:1485  10

4

0.9860

8:5999  10  4

0.9683

0.5

1=160

4:0244  10  4

1.0177

4:2977  10  4

1.0007

0.8

1=320

1:9940  10  4

1.0131

2:1517  10  4

0.9981

0.8

Table 8 Errors and computational orders obtained for Test problem 3 with h ¼ 1=10, r ¼0 and c¼ 0.5 on interval Ω ¼ ½0; 1  ½0; 1. τ

α ¼ 1:25

α ¼ 1:75

CPU time (s)

L1

C-order

L1

C-order

1=10

6:6533  10  2

–

6:1693  10  2

–

0.1256

1=20

3:8488  10  2

0.7897

3:4433  10  2

0.8413

2104

1=40

2:0823  10  2

0.8862

1:8224  10  2

0.9179

0.3981

1=160

1:0813  10  2

0.9454

9:3480  10  3

0.9631

0.8167

1=320

5:4659  10  3

0.9842

4:6927  10  3

0.9942

2.0600

1=640

2:6996  10  3

1.0177

2:3071  10  3

1.0243

5.5539

1=1280

1:2923  10  3

1.0628

1:1205  10  3

1.0419

17.5537

1=2560

6:1201  10  4

1.0783

5:4035  10  4

1.0522

62.6263

where f¼

2t 2  α

Γ ð3  αÞ



2 2 2 2 sec h ðx þ y rÞ þ sec h ðx þ y þ rÞ 4t 2 3 sec h ðx þ y rÞ tan h ðx þ y rÞ

 2 2 2 2  sec h ðx þ y rÞ þ 3 sec h ðx þ y þ rÞ tan h ðx þy þ rÞ  sec h ðx þy þ rÞ

2

3 2 2 2 2  t 4 sec h ðx þy  rÞ þ sec h ðx þ y þ rÞ  t 6 sec h ðx þ y  rÞ þ sec h ðx þ y þ rÞ :

428

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

Then the exact solution of this test problem is

 2 2 uðx; y; tÞ ¼ t 2 sec h ðx þ y  rÞ þ sec h ðx þ y þ rÞ ;

x; y A Ω:

where r is a real constant. We solve this problem with the method presented in this paper with several values of h, τ, α, c on different intervals ½a; b  ½a; b at final time T ¼1. Table 8 shows errors and computational orders obtained for Test problem 3 with h ¼ 1=10, r ¼0 and c ¼0.5 on interval Ω ¼ ½0; 1  ½0; 1. Also, Table 9 presents errors and computational orders obtained for Test problem 3 with h ¼ 1=5, r ¼0, T ¼2 and c ¼0.5 on interval Ω ¼ ½0; 1  ½0; 1. The graphs of approximation solution using the present method with h ¼ 1=5, τ ¼ 1=1000, α ¼ 1:45, T¼ 1 and c¼ 0.5 on rectangular domain Ω ¼ ½0; 1  ½0; 1 with different vales of r for Test problem 3 are shown in Fig. 7. Table 9 Errors and computational orders obtained for Test problem 3 with h ¼ 1=5, r¼ 0, T ¼ 2 and c ¼ 0.5 on interval Ω ¼ ½0; 1  ½0; 1. τ

α ¼ 1:65

α ¼ 1:95

L1

CPU time (s)

C-order

L1

C-order

3:4785  100

–

2:0794  100

–

1=200

2:6588  10

0

0.3877

1:2494  100

0.7349

0.3811

1=400

1:7207  100

0.6278

6:5373  10  1

0.9344

0.7838

1=800

8:7591  10  1

0.9741

2:8782  10  1

1.1835

2.1341

1=1600

2:5909  10  1

1.7573

8:3381  10  2

1.7919

6.6096

1=3200

1:2688  10  1

1.0300

4:2251  10  2

0.9807

22.8121

2.5

1.4

2

1.2

1.5

U(x,y,1)

U(x,y,1)

1=100

1 0.5 0 1

0.8 0.6

0.5

0

−0.5

−1 1

0.5

0

−1

0.4 1

−1

x

−0.5

0.5

y

0 x

0 −0.5

0.5

−0.5

y

−1 1

0.08

U(x,y,1)

1.5

U(x,y,1)

1

1 0.5

0.06 0.04 0.02

0 1

−1 −0.5

0.5

0

0

x

−0.5

0.5 −1 1

y

0 1

−1 −0.5

0.5 0

0 x

−0.5

0.5

y

−1 1

Fig. 7. Graphs of approximation solution using the present method with h ¼1/5, τ ¼ 1/1000, α¼ 1.45, T ¼1 and c ¼0.5 on rectangular domain Ω ¼ ½0; 1  ½0; 1 with different vales of r for Test problem 3.

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

429

5.4. Test problem 4 In this test, we consider two-dimensional nonlinear time fractional Klein–Gordon equation: ∂α uðx; y; tÞ ∂2 uðx; y; tÞ ∂2 uðx; y; tÞ ¼ þ  u3 ðx; y; tÞ þ f ðx; y; tÞ; ∂t α ∂x2 ∂y2 with the initial condition uðx; y; 0Þ ¼ 0; ∂uðx; y; 0Þ ¼ 0; ∂t

x; y A Ω;

and boundary conditions uðx; y; tÞ ¼ t 2 cos ðπ xÞ cos ðπ yÞ; where



f ðx; y; tÞ ¼

2t 2  α þ2π 2 Γ ð2  αÞ



x; y A ∂Ω;

cos ðπ xÞ cos ðπ yÞ þt 6 cos 3 ðπ xÞ cos 3 ðπ yÞ:

Then the exact solution of this test problem is uðx; y; tÞ ¼ t 2 cos ðπ xÞ cos ðπ yÞ;

x; y A Ω:

We solve this problem with the method presented in this paper with several values of h, τ, α, c on different intervals ½a; b  ½a; b at final time T¼ 1. Table 10 shows errors and computational orders obtained for Test problem 4 with h ¼ 1=10 and T ¼1 on interval Ω ¼ ½0; 1  ½0; 1. Fig. 8 shows the graphs of approximation solution using the present method with h ¼ 1=10, τ ¼ 1=160, α ¼ 1:15, T ¼1 and c¼0.05 on rectangular domain Ω ¼ ½  1; 1  ½ 1; 1 for Test problem 4. The graphs of approximation solution using the present method with h ¼ 1=10, τ ¼ 1=160, α ¼ 1:15, T ¼1 and c ¼0.05 on rectangular domain Ω ¼ ½0; 1  ½0; 1 using irregular points for Test problem 4 are presented in Fig. 9. Fig. 10 demonstrates the graphs of approximation solution using the present method with h ¼ 1=5, τ ¼ 1=160, α ¼ 1:15, T ¼1 and c¼ 0.5 on rectangular domain Ω ¼ ½  2; 2  ½  2; 2 using irregular points for Test problem 4. Table 10 Errors and computational orders obtained for Test problem 4 with h ¼ 1=10 and T ¼ 1 on interval Ω ¼ ½0; 1  ½0; 1. τ

α ¼ 1:45

α ¼ 1:85

c

L1

C-order

L1

C-order

1=10

3:2987  10  3

–

6:1693  10  3

–

0.1

1=20

8:9043  10  4

1.8893

3:4433  10  4

2.2025

0.2

1=40

3:8768  10  4

1.1996

1:8224  10  4

0.4465

0.4

1=40

1:7125  10  4

1.1787

1:8224  10  4

0.9776

0.6

1=160

6:0272  10  5

1.5065

9:3480  10  5

1.5004

0.8

−3

x 10 1.5

2 1

1

0

0.5 −1

0 1

−2 1 0.5 0

−0.5 −1 −1

−0.5

0

0.5

1

0.5 0 −0.5 −1 −1

−0.5

0

0.5

1

Fig. 8. Graphs of approximation solution using the present method with h ¼1/10, τ¼ 1/160, α¼ 1.15, T ¼ 1 and c ¼0.05 on rectangular domain Ω ¼ ½  1; 1  ½  1; 1 for Test problem 4.

430

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

1

0.015

0.01

0

Error

U(x,y,1)

0.5

−0.5

0.005

−1 1

y

0 1 0.5 0

0.4

0.2

0

0.8

0.6

0.5

1

y

x

0

0

0.2

0.8

0.6

0.4

1

x

Fig. 9. Graphs of approximation solution using the present method with h ¼ 1/10, τ ¼1/160, α ¼1.15, T ¼1 and c ¼0.05 on rectangular domain Ω ¼ ½0; 1  ½0; 1 using irregular points for Test problem 4.

is Numerical solution is Exact solution

−3

x 10 6

2

4

Error

U(x,y,1)

1 0

2

−1

0 2

−2 2 1 0

y

−1 −2

0

−1

−2

1

2

x

1

y

0 −1 −2

−2

−1

0

1

2

x

Fig. 10. Graphs of approximation solution using the present method with h ¼1/5, τ ¼ 1/160, α¼ 1.15, T ¼ 1 and c¼ 0.5 on rectangular domain Ω ¼ ½  2; 2  ½  2; 2 using irregular points for Test problem 4.

5.5. Test problem 5 In this test, we consider two-dimensional nonlinear time fractional sine-Gordon equation: ∂α uðx; y; z; tÞ ∂2 uðx; y; z; tÞ ∂2 uðx; y; z; tÞ ∂2 uðx; y; z; tÞ ¼ þ þ  sin ðuðx; y; z; tÞÞ þ f ðx; y; z; tÞ; ∂t α ∂x2 ∂y2 ∂z2 with the initial condition uðx; y; 0Þ ¼ 0; ∂uðx; y; 0Þ ¼ 0; ∂t

x; y A Ω;

and boundary condition uðx; y; z; tÞ ¼ t 2 exp  and the source term is

ðx  r Þ2

β



ðy  r Þ2

β

(



ðz  r Þ2

β

! ;

ðx; y; zÞ A ∂Ω;

!! ðx  r Þ2 ðy r Þ2 ðz  r Þ2   Γ ð3  αÞβ2 r r r Γ ð3  αÞβ2 !    ) ðx  r Þ2 ðy  r Þ2 ðz  r Þ2 1 1 2 2 2 1 1 1 2   :  12 exp   t 2  α β þ t 2 Γ ð3  αÞ  β þ r 2 þ  z  x  y r þ x2 þ y2 þ z2 6 2 3 3 3 3 3 3 r r r

f ðx; y; z; tÞ ¼

1

sin exp 

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

431

Then the exact solution of this test problem is 2

uðx; y; z; tÞ ¼ t exp 

ðx r Þ2

β



ðy  r Þ2

β



ðz  r Þ2

!

β

ðx; y; zÞ A Ω  R3 :

;

We solve this problem with the method presented in this paper with several values of h, τ, spherical domains at final time T ¼1.

α, c on different domains such as cubical and

Table 11 Errors and computational orders obtained for Test problem 5 on cubical domain with h ¼ 1=5, c¼ 0.9, β ¼ 1=10000, r¼ 0.5 and T ¼ 1. τ

α ¼ 1:50

α ¼ 1:25 C-order

L1

CPU time (s)

L1

C-order

1=10

8:7799  10

–

9:8251  10

–

0.1437

1=20

4:4262  10  3

0.9881

4:8564  10  3

1.0166

0.2151

1=40

2:2048  10  3

1.0054

2:3943  10  3

1.0203

0.3587

1=80

1:0821  10  3

1.0268

1:1711  10  3

1.0317

0.8219

1=160

5:1875  10  4

1.0607

5:6187  10  4

1.0596

2.3152

1=320

2:3679  10  4

1.1314

2:5787  10  4

1.1236

7.5786

1=640

9:5777  10  5

1.3058

1:0603  10  4

1.2822

27.1448

1=1280

3:4025  10  5

1.4931

3:1724  10  5

1.7408

102.9832

3

3

−4

x 10

1.001

6

U(x,y,z=r,T)

Absolute Error

8

4 2

1.0005 1 0.9995 0.999 1

0 1 0.5 y 0

0.2

0

0.4

0.6

1

0.8

y

x

0.5 0

0

0.4

0.2

0.6 x

0.8

1

Fig. 11. Graphs of approximation solution and absolute error using the present method with h ¼1/5, τ ¼ 1/1280, α¼ 1.45, r ¼ 0.5 and β¼ 1/1000, T ¼ 1 and c¼ 0.9 on cubical domain in the plate z ¼r for Test problem 5.

−4

x 10 16

0.9998 1

12

0.8

0.6

10

0.6

0.4

8

0.4

0.2

6

0.2

0 1

4

0 1

1

0.9996 0.9994 0.9992

z

0.8

z

14

0.5

y

0 0

0.2

0.4

0.6

x

0.8

1

0.999 0.9988

2 0

y

0.5 0

0

0.2

0.4

0.6

x

0.8

1

0.9986 0.9984

Fig. 12. Approximation solution and absolute error using the present method with h ¼ 1/15, τ¼ 1/100, α¼ 1.85, T ¼ 1 and c ¼0.001, r¼ 0.5 and β ¼1/1000 on cubical domain in the plate z¼ r for Test problem 5.

432

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434 −4

−4

x 10 0

x 10 0

1

−2

0.5

−4

−2

0.6 0.4

−4

0.2

0

−6

−0.5

y

z

0.8 1

−6

0 −0.2

−8

−8

−0.4

−1 1 0.5 0 −0.5

y

−1

−1

−0.5

0

0.5

1

−10

−0.8

−12

x

−10

−0.6

−1 −1

−12 −0.5

0

0.5

1

x 1

1

0.8

1.0005

1

0.6

1

0.4

0.9995

0.9995 0.999

0

0.9985

−0.5 −1 1 0.5

y

0 −0.5 −1 −1

−0.5

0.5

0

x

1

0.2

y

z

0.5

1.0005

0.999

0

0.9985

−0.2

0.998

0.998

−0.4

0.9975

−0.6

0.9975

0.997

−0.8

0.997

−1 −1

−0.5

0

0.5

1

x

Fig. 13. Graphs of approximation solutions and absolute error using the present method with τ ¼1/20, α ¼ 1.25, r¼ 0.5, β¼ 1/1000 and T ¼ 1 on spherical domain for Test problem 5.

Table 11 shows errors and computational orders obtained for Test problem 5 on cubical domain with h ¼ 1=5, c¼ 0.9, β ¼ 1=10; 000, r ¼0.5 and T¼ 1. Table 5 presents that the computational orders of time discrete are close to the theoretical orders, i.e the order of method is OðτÞ in the time variable. Fig. 11 shows graphs of approximation solution and absolute error using the present method with h ¼ 1=5, τ ¼ 1=1280, α ¼ 1:45, T¼ 1 and c¼ 0.9 on cubical domain in the plate z ¼r for Test problem 5. Approximation solution and absolute error using the present method with h ¼ 1=15, τ ¼ 1=100, α ¼ 1:85, T ¼1 and c ¼0.001, r ¼0.5 and β ¼ 1=1000 on cubical domain in the plate z ¼r for Test problem 5 are presented in Fig. 12. Also, Fig. 13 presents graphs of approximation solutions and absolute error using the present method with τ ¼ 1=20, α ¼ 1:25, r ¼ 0.5, β ¼ 1=1000 and T¼ 1 on spherical domain for Test problem 5.

6. Conclusion In this paper, we solved the nonlinear time fractional sine-Gordon and Klein–Gordon equations using the meshless method of radial basis functions. Firstly, we discretized the time derivative using a scheme with order Oðτ3  α Þ and obtained a time semi-discrete scheme. Also, we obtained a full discrete scheme using meshless method based on radial basis functions and Kansa's approach. Also we showed that the time discrete scheme is unconditionally stable and its convergence order is OðτÞ. We used the energy method for proving the stability and convergence of time discrete scheme. Since the coefficient matrix is ill-conditioned, therefore, we used the LU decomposition method for solving the linear system of algebraic equations which arises from the process of collocating points. We solved the mentioned equation on different irregular domains. Numerical results showed that the computational orders of time discrete are close to the theoretical convergence orders and confirm the efficiency of the new method developed in the current paper.

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

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

433

References [1] Abbasbandy S, Roohani Ghehsareh H, Hashim I. A meshfree method for the solution of two-dimensional cubic nonlinear Schrödinger equation. Eng Anal Bound Elem 2013;37:885–98. [2] Agrawal OP. A general solution for the fourth-order fractional diffusion-wave equation. Fract Calculus Appl Anal 2000;3(1):1–12. [3] Agrawal OP. Ageneral solution for the fourth-order fractional diffusion-wave equation defined in bounded domain. Comput Struct 2001;79:1497–501. [4] Bagley R, Torvik P. A theoretical basis for the application of fractional calculus to viscoelasticity. J Rheol 1983:201–10. [5] Brezis H. Functional analysis, Sobolev spaces and partial differential equations. Dordrecht, Heidelberg, London: Springer, New York; 2011. [6] Cai Z. Convergence and error estimates for meshless Galerkin methods. Appl Math Comput 2007;184:908–16. [7] Chen S, Liu F, Zhuang P, Anh V. Finite difference approximations for the fractional Fokker–Planck equation. Appl Math Model 2009;33:256–73. [8] Chen Y, Wu Y, Cui Y, Wang Z, Jin D. Wavelet method for a class of fractional convection-diffusion equation with variable coefficients. J Comput Sci 2010;1:146–9. [9] Cheng AHD, Golberg MA, Kansa EJ, Zammito G. Exponential convergence and h–c multiquadric collocation method for partial differential equations. Numer Methods Partial Diff Eqs 2003;19:571–94. [10] Cheng AHD. Multiquadric and its shape parameter–a numerical investigation of error estimate, condition number, and round-off error by arbitrary precision computation. Eng Anal Bound Elem 2012;36:220–39. [11] Dehghan M, Mohammadi V. The numerical solution of Fokker–Planck equation with radial basis functions (RBFs) based on the meshless technique of Kansa's approach and Galerkin method. Eng Anal Bound Elem 2014;47:38–63. [12] Dehghan M, Manafian J, Saadatmandi A. The solution of the linear fractional partial differential equations using the homotopy analysis method. Z Naturfor—Sect A 2010;65:935–49. [13] Dehghan M, Mohebbi A. High-order compact boundary value method for the solution of unsteady convection-diffusion problems. Math Comput Simul 2008;79: 683–99. [14] Dehghan M, Shokri A. A numerical method for two-dimensional Schrodinger equation using collocation and radial basis functions. Comput Math Appl 2007;54:136–46. [15] Dehghan M, Mirzaei D. The meshless local Petrov–Galerkin (MLPG) method for the generalized two-dimensional non-linear Schrodinger equation. Eng Anal Bound Elem 2008;32:747–56. [16] Dehghan M, Salehi R. The solitary wave solution of the two-dimensional regularized long-wave equation in fluids and plasmas. Comput Phys Commun 2011;182:2540–9. [17] Dehghan M, Tatari M. Use of radial basis functions for solving the second-order equation with nonlocal boundary conditions. Numer Methods Partial Diff Eqs 2008;24:924–38. [18] Dehghan M, Ghesmati A. Combination of meshless local weak and strong (MLWS) forms to solve the two dimensional hyperbolic telegraph equation. Eng Anal Bound Elem 2010;34:324–36. [19] Dehghan M, Nikpour A. Numerical solution of the system of second-order boundary value problems using the local radial basis functions based differential quadrature collocation method. Appl Math Model 2013;37:8578–99. [20] Diethelm K, Ford NJ. Analysis of fractional differential equations. J Math Anal Appl 2002;265:229–48. [21] Du R, Cao WR, Sun ZZ. A compact difference scheme for the fractional diffusion-wave equation. Appl Math Model 2010;34:2998–3007. [22] Duan Y, Tan YJ. Meshless Galerkin method based on regions partitioned into subdomains. Appl Math Comput 2005;162:317–27. [23] Esmaeili S, Shamsi M. A pseudo-spectral scheme for the approximate solution of a family of fractional differential equations. Commun Nonlinear Sci Numer Simul 2011;16:3646–54. [24] Esmaeili S, Shamsi M, Luchko Y. Numerical solution of fractional differential equations with a collocation method based on Müntz polynomials. Comput Math Appl 2011;62:918–29. [25] Fasshauer GE. Meshfree approximation methods with MATLAB. USA: World scientific; 2007. [26] Fornberg B, Wright G, Larsson E. Some observations regarding interpolants in the limit of flat radial basis functions. Comput Math Appl 2004;47:37–55. [27] Franke C, Schaback R. Convergence order estimates of meshless collocation methods using radial basis functions. Adv Comput Math 1998;8:381–99. [28] Garra R, Polito F. Fractional Klein–Gordon equation for linear dispersive phenomena: analytical methods and applications; 2014, submitted for publication. [29] Garra R, Orsingher E, Polito F. Fractional Klein–Gordon equations and related stochastic processes; 2014, submitted for publication. [30] Gepreel KA, Mohamed MS. Analytical approximate solution for nonlinear space-time fractional Klein–Gordon equation. Chin Phys B 2013;22:010201-1. [31] Ginoa M, Cerbelli S, Roman HE. Fractional diffusion equation and relaxation in complex viscoelastic materials. Physica A 1992;191:449–53. [32] Golmankhaneh Alireza K, Golmankhaneh Ali K, Baleanu D. On nonlinear fractional Klein–Gordon equation. Signal Process 2011;91:446–51. [33] Gu YT, Zhuang P, Liu Q. An advanced meshless method for time fractional diffusion equation. Int J Comput Methods 2011;8:653–65. [34] Gu YT, Zhuang P, Liu F. An advanced implicit meshless approach for the non-linear anomalous subdiffusion equation. Comput Model Eng Sci 2010;56:303–34. [35] Hardy RL. Multiquadric equations of topography and other irregular surfaces. J Geophys Res 1971;76:1705–915. [36] Hariharan G. Wavelet method for a class of fractional Klein–Gordon equations. J Comput Nonlinear Dyn 2013;8:021008-1. [37] Huang CS, Yen HD, Cheng AHD. On the increasingly flat radial basis function and optimal shape parameter for the solution of elliptic PDEs. Eng Anal Bound Elem 2010;34:802–9. [38] Huang CS, Lee CF, Cheng AHD. Error estimate, optimal shape factor, and high precision computation of multiquadric collocation method. Eng Anal Bound Elem 2007;31:614–23. [39] Jafari H, Tajadodi H, Kadkhoda N, Baleanu D, Fractional subequation method for Cahn–Hilliard and Klein–Gordon equations. Abstract and applied analysis, vol. 2013. Article ID 587179. 5p. [40] Kansa EJ. Multiquadrics A scattered data approximation scheme with applications to computational fluid-dynamics. I. Comput Math Appl 1990;19:127–45. [41] Kansa EJ. Multiquadrics A scattered data approximation scheme with applications to computational fluid dynamics. II. Comput Math Appl 1990;19:147–61. [42] Kansa EJ, Aldredge RC, Ling L. Numerical simulation of two-dimensional combustion using mesh-free methods. Eng Anal Bound Elem 2009;33:940–50. [43] Kurulay M. Solving the fractional nonlinear Klein–Gordon equation by means of the homotopy analysis method. Adv Diff Eqs 2012;187:1–8. [44] Langlands TAM, Henry BI. The accuracy and stability of an implicit solution method for the fractional diffusion equation. J Comput Phys 2005;205:719–36. [45] Laskin N, Zaslavsky G. Nonlinear fractional dynamics on a lattice with long range interactions. Physica A 2006;368:38–54. [46] Li CP, Ding H. Higher order finite difference method for the reaction and anomalous-diffusion equation. Appl Math Model; in Press, http://www.dx.doi.org.10.1016/j. apm.2013.12.002. [47] Li CP, Chen A, Ye J. Numerical approaches to fractional calculus and fractional ordinary differential equation. J Comput Phys 2011;230:3352–68. [48] Liu Q, Liu F, Turner I, Anh V, Gu YT. A RBF meshless approach for modeling a fractal mobile/immobile transport model. Appl Math Comput 2014;226:336–47. [49] Liu GR, Gu YT. An introduction to meshfree methods and their programming. Berlin, Heidelberg, New York: Springer, Dordrecht; 2005. [50] Liu Q, Gu Y, Zhuang P, Liu F, Nie Y. An implicit RBF meshless approach for time fractional diffusion equations. Comput Mech 2011;48:1–12. [51] Liu GR, Gu YT. Meshfree methods: moving beyond the finite element method. Second ed. Boca Raton: CRC Press; 2002. [52] Liu F, Zhuang P, Burrage K. Numerical methods and analysis for a class of fractional advection-dispersion models. Comput Math Appl 2012;64:2990–3007. [53] Mainardi F. Fractional relaxation-oscillation and fractional diffusion-wave phenomena. Chaos Solitons Fractals 1996;7(9):1461–77. [54] Mainardi F. The fundamental solutions for the fractional diffusion-wave equation. Appl Math Lett 1996;9(6):23–8. [55] Madych WR, Nelson SA. Multivariate interpolation and conditionally positive definite functions. Approx Theory Appl 1988;4:77–89. [56] Metzler R, Klafter J. The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics. J Phys A 2004;37:R161–208. [57] Metzler R, Klafter J. Boundary value problems for fractional diffusion equations. Physica A 2000;278:107–25. [58] Miller KS, Ross B. An introductional the fractional calculus and fractional differential equations. New York, London: Academic Press; 1974. [59] Mirzaei D, Dehghan M. A meshless based method for solution of integral equations. Appl Numer Math 2010;60:245–62. [60] Mohebbi A, Abbaszadeh M, Dehghan M. The use of a meshless technique based on collocation and radial basis functions for solving the time fractional nonlinear Schrödinger equation arising in quantum mechanics. Eng Anal Bound Elem 2013;37:475–85. [61] Mohebbi A, Abbaszadeh M, Dehghan M. High-order difference scheme for the solution of linear time fractional Klein–Gordon equations. Numer Methods Partial Diff Eqs 2014;30:1234–53. [62] Murillo JQ, Yuste SB. An explicit difference method for solving fractional diffusion and diffusion-wave equations in the Caputo form. J Comput Nonlinear Dyn 2011;6 [article no. 021014].

434

M. Dehghan et al. / Engineering Analysis with Boundary Elements 50 (2015) 412–434

[63] Mohyud-Din ST, Negahdary E, Usman M. A meshless numerical solution of the family of generalized fifth-order Korteweg-de Vries equations. Int J Numer Methods Heat Fluid Flow 2012;22:641–58. [64] Momani S. An algorithm for solving the fractional convection-diffusion equation with nonlinear source term. Commun Non Linear Sci Numer Simulat 2007;12:1283–90. [65] Momani S, Yildirim A. Analytical approximate solutions of the fractional convection-diffusion equation with nonlinear source term by He's homotopy perturbation method. Int J Comput Math 2010;87:1057–65. [66] Momani S, Odibat ZM. Fractional Green function for linear time-fractional inhomogeneous partial differential equations in fluid mechanics. J Appl Math Comput 2007;24:167–78. [67] Odibat ZM. Computational algorithms for computing the fractional derivatives of functions. Math Comput Simul 2009;79:2013–20. [68] Oldham KB, Spanier J. The fractional calculus: theory and application of differentiation and integration to arbitrary order. Academic Press; 1974. [69] Oldham KB, Spanier J. The fractional calculus. New York, London: Academic Press; 1974. [70] Podulbny I. Fractional differential equations. New York: Academic Press; 1999. [71] Quarteroni A, Valli A. Numerical approximation of partial differential equations. New York: Springer-Verlag; 1997. [72] Rippa S. An algorithm for selecting a good value for the parameter c in radial basis function interpolation. Adv Comput Math 1999;11:193–210. [73] Saadatmandi A, Dehghan M, Azizi MR. The Sinc-Legendre collocation method for a class of fractional convection-diffusion equations with variable coefficients. Commun Non Linear Sci Numer Simulat 2012;17:4125–36. [74] Salehi R, Dehghan M. A meshless local Petrov–Galerkin method for the time-dependent Maxwell equations. J Comput Appl Math 2014;268:93–110. [75] Salehi R, Dehghan M. A moving least square reproducing polynomial meshless method. Appl Numer Math 2013;69:34–58. [76] Schneider WR, Wess W. Fractional diffusion and wave equations. J Math Phys 1989;30(1):134–44. [77] Shokri A, Dehghan M. Meshless method using radial basis functions for the numerical solution of two-dimensional complex Ginzburg–Landau equation. Comput Model Eng Sci 2012;34:333–58. [78] Shokri A, Dehghan M. A Not-a-Knot meshless method using radial basis functions and predictor-corrector scheme to the numerical solution of improved Boussinesq equation. Comput Phys Commun 2010;181:1990–2000. [79] Sun ZZ, Wu XN. A fully discrete difference scheme for a diffusion-wave system. Appl Numer Math 2006;56:193–209. [80] Tadjeran C, Meerschaert MM, Scheffler HP. A second-order accurate numerical approximation for the fractional diffusion equation. J Comput Phys 2006;213:205–13. [81] Taleei A, Dehghan M. Direct meshless local Petrov–Galerkin method for elliptic interface problems with applications in electrostatic and elastostatic. Comput Methods Appl Mech Eng 2014;278:479–98. [82] Tatari M, Dehghan M. A method for solving partial differential equations via radial basis functions: application to the heat equation. Eng Anal Bound Elem 2010;34:206–12. [83] Tarasova VE, Zaslavsky GM. Fractional dynamics of systems with long-range space interaction and temporal memory. Physica A: Stat Mech Appl 2007;383:291–308. [84] Vanani SK, Aminataei A. On the numerical solution of neural delay differential equations using multiquadric approximation scheme. Bull Korean Math Soc 2008;45:663–70. [85] Vong S, Wang Z. A compact difference scheme for a two dimensional fractional Klein–Gordon equation with Neumann boundary conditions. J Comput Phys 2014;274:268–82. [86] Wendland H. Scattered data approximation. In: Cambridge monograph on applied and computational mathematics. Cambridge University Press; 2005. [87] Wendland H. Meshless Galerkin methods using radial basis functions. Math Comput 1999;68:1521–31. [88] Wendland H. Error estimates for interpolation by compactly supported radial basis functions of minimal degree. J Approx Theory 1998;93:258–72. [89] Wess W. The fractional diffusion equation. J Math Phys 1996;27:2782–5. [90] Wu Z, Shaback R. Local error estimates for radial basis function interpolation of scattered data. IMA J Numer Anal 1993;13:13–27. [91] Yoon J. Spectral approximation orders of radial basis function interpolation on the Sobolov space. SIAM J Math Anal 1999;33:946–58. [92] Yoon J. Lp-error estimates for “shifted” surface spline interpolation on Sobolev space. Math Comp 2003;72(243):1349–67. [93] Yuste SB. Weighted average finite difference methods for fractional diffusion equations. J Comput Phys 2006;216:264–74. [94] Yuste SB, Acedo L. An explicit finite difference method and a new Von Neumann-type stability analysis for fractional diffusion equations. SIAM J Numer Anal 2005;42 (5):1862–74. [95] Zhuang P, Gu YT, Liu F, Turner I, Yarlagadda PKDV. Time-dependent fractional advection-diffusion equations by an implicit MLS meshless method. Int J Numer Methods Eng 2011;88:1346–62.