Meshless spectral method for solution of time-fractional coupled KdV equations

Applied Mathematics and Computation 341 (2019) 321–334

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Meshless spectral method for solution of time-fractional coupled KdV equations Manzoor Hussain a, Sirajul Haq a,∗, Abdul Ghafoor a Faculty of Engineering Sciences, GIK Institute, Topi 23640, KPK, Pakistan

a r t i c l e

i n f o

Keywords: Shallow water waves Coupled KdV equations Meshless spectral method Radial basis functions Caputo fractional derivative

a b s t r a c t In this article, an efficient and accurate meshless spectral interpolation method is formulated for the numerical solution of time-fractional coupled KdV equations that govern shallow water waves. Meshless shape functions constructed via radial basis functions (RBFs) and point interpolation are used for discretization of the spatial operator. Approximation   of fractional temporal derivative is obtained via finite differences of order O τ 2−α and a quadrature formula. The formulated method is applied to various test problems available in the literature for its validation. Approximation quality and efficiency of the method is measured via discrete error norms E2 , E∞ and Erms . Convergence analysis of the proposed method in space and time is numerically determined by varying nodal points M and time step-size τ respectively. Stability of the proposed method is discussed and affirmed computationally, which is an important ingredient of the current study. © 2018 Elsevier Inc. All rights reserved.

1. Introduction Shallow water waves describe the flow below a pressure surface in a fluid usually observed in lake shores and beaches. Mathematical models describing dynamics of this phenomena offers a variety of applications in engineering sciences and specially in ocean engineering. Nonlinear partial differential equations (PDEs) are used for the description of these model phenomenons and many others [9,32,33]. Time-dependent coupled KdV equations got much researchers attention in this respect [20]. For instance, they are used to model interaction of two long waves having different dispersion relation [13]. Other applications include, the Fermi–Pasta–Ulam problem in continuum limit, the evolution of long, one-dimensional waves in shallow water having weakly non-linear restoring forces, long internal waves in a density stratified ocean, ion acoustic waves in a plasma and acoustic waves on a crystal lattice. With the advent of fractional calculus of arbitrary order derivatives, fractional PDEs (FPDEs) are believed to have more accurate description of real world problems [10,16,21,35]. They have applications in many fields, e.g. engineering, physics, acoustics, viscoelasticity, fluid dynamics and even in finance. FPDEs containing time-fractional derivatives received importance as of their infinitesimal generating nature in evolution processes, when considering long time scale limit, stability of states and equilibrium [6,11,12]. The Caputo sense of fractional derivative provides modeling of real-world phenomenons more realistically and it simplifies initial conditions in formulation of the problem like the integer case [3]. FPDEs always appear as a challenging task when it comes to obtaining their analytical solutions. Therefore, numerical techniques are always of interest and preferable for approximate solutions of such complex equations. ∗

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

https://doi.org/10.1016/j.amc.2018.09.001 0 096-30 03/© 2018 Elsevier Inc. All rights reserved.

322

M. Hussain et al. / Applied Mathematics and Computation 341 (2019) 321–334

In this work, a robust and efficient meshless spectral radial point interpolation method (MSRPIM) is formulated to solve numerically a general class of shallow water waves equations by letting time derivative of arbitrary order. In this regard, we have the time-fractional coupled KdV (TFCKdV) equations of fractional order 0 < α , β ≤ 1, read as:

∂tα U − p∂xxxU − p1U ∂xU − p2V ∂xV = Y (x, t ),

(1.1)

∂tβ V − q∂xxxV + q1U ∂xV = Z (x, t ),

(1.2)

subject to initial conditions

U (x, 0 ) = U0 (x ),

V (x, 0 ) = V0 (x ),

x ∈ [a, b],

(1.3)

and boundary conditions

U (l, t ) = χl (t )

and V (l, t ) = γl (t ),

l = a, b, t > 0,

(1.4)

where p, q < 0, p1 , p2 , q1 are constants and Y, Z are known source functions. It is observed that, time-fractional derivative in coupled KdV equations has a sound physical impact. As integer order KdV equations consider only the instant of time (local property), thus the obtained soliton solutions may not show the solitary waves very well. Therefore, TFCKdV equations may estimate the effect of higher order dispersion of the regular KdV equation to increase the amplitude of the soliton [4,5]. TFCKdV equations have been solved by various methods, e.g. Homotopy decomposition method (HDM) [2], generalized Differential transform method (gDTM) [18], and very recently by spectral collocation method using shifted Chebyshev polynomials [1]. Amongst a variety of existing numerical methods in literature like finite differences and finite element based method, meshfree (meshless) methods using radial basis functions (RBFs) established their own importance over the usual meshdependent methods (e.g. FDM and FEM). The beauty of meshfree methods is their pair-wise distance dependence of nodal points, ease implementation and high accuracy. There are a number of meshfree methods developed by researchers both in global and local form, e.g. global RBFs approximation method [8,14,34] as well meshless local Petrov–Galerkin (MLPG) method [7,23] for solution of potential engineering problems. However, a serious drawback of these methods is either the severe ill-conditioning of the system matrix in global form or the costly background integration schemes required in local form. To resolve these issues, a new technique called meshfree spectral radial point interpolation method (MSRPIM) is developed that based on spectral collocation approach and meshfree shape functions. RBFs together with point interpolation are utilized for the generation of meshfree shape functions which are then used for building explicit differentiation matrices for approximation of spatial part. An important feature of these shape functions is their Kronecker delta functions property, which ease the implementation of boundary conditions. Other features are their partition of unity and patch test property, (see e.g. [17]). It is remarked that the current methodology is quite different from the usual RBFs approximation method [22]. MSRPIM has been applied recently to solve various important time-fractional PDEs, e.g. 2D convection– diffusion–reaction equations [24], fractional integral equations [25], Shrödinger equations [26], cable equations [27], inverse source and diffusion problems [28,29] and reaction-subdiffusion equations [30]. Our objective in this paper is to showcase application of MSRPIM for obtaining accurate solutions of TFCKdV equations. Rest of the paper is organized in the following manner: Section 2 gives a brief description of solution methodology for the considered problem whereas in Section 3 stability analysis of the proposed method is carried out. Computational results are presented in Section 4 and finally we end the paper with a concise conclusion given in Section 5. 2. Proposed solution methodology This part briefly outlines the formulation of proposed method for the TFCKdV Eqs. (1.1)–(1.4). Let the time interval [0, T] is discretized using grid points {tn = nτ |n = 0, . . . , N } where τ = T /N is time-step length and T is final time. Adopting the notation for solutions at nth time level as U n (x ) = U (x, tn ), V n (x ) = V (x, tn ) and same for the others. Then approximation to the temporal fractional derivative at the (n + 1 )th time level is obtained as [8]:

∂tα U n+1 (x ) =

⎧ n n+1−m

  α (x ) − U n−m (x ) + O τ 2−α , ⎪ ⎨α m=0 m U

0 < α < 1,

n ≥ 0,

n+1 n ⎪ ⎩ U (x ) − U (x ) + O (τ ), τ

α = 1,

n ≥ 0,

where

α =

τ −α , αm = (m + 1 )1−α − (m )1−α .

(2 − α )

The same procedure is applied for approximation of temporal derivative of order β . Using the above relation while waiving the error term, Eqs. (1.1) and (1.2) in semi-discrete implicit form becomes

α

n  m=0



αm U n+1−m (x ) − U n−m (x ) − p∂xxxU n+1 (x ) − p1U n+1 (x )∂xU n+1 (x ) = p2V n (x )∂xV n (x ) + Y n+1 (x ),

(2.1)

M. Hussain et al. / Applied Mathematics and Computation 341 (2019) 321–334

β

n 



βm V n+1−m (x ) − V n−m (x ) − q∂xxxV n+1 (x ) = −q1U n (x )∂xV n (x ) + Z n+1 (x ).

323

(2.2)

m=0

We linearized the nonlinear involved term as:

U n+1 (x )∂xU n+1 (x ) = U n (x )∂xU n+1 (x ) + U n+1 ∂xU n (x ) − U n (x )∂xU n (x ). Exploiting above relation into Eqs. (2.1) and (2.2) and then simplification lead us to



α U n+1 (x ) − p∂xxxU n+1 (x ) − p1 cn (x )∂xU n+1 (x ) + dn (x )U n+1 (x ) = [α − p1 dn (x )]U n (x ) + Y n+1 (x ),

(2.3)

β V n+1 (x ) − q∂xxxV n+1 (x ) = β V n (x ) − q1 cn (x )∂xV n (x ) + Z n+1 (x ),

(2.4)

where cn (x ) = U n (x ), dn (x ) = ∂xU n (x ) and

Y n+1 (x ) = Y n+1 (x ) + p2V n (x )∂xV n (x ) − Z1n (x ), Z1n (x ) = α

n 

Z n+1 (x ) = Z n+1 (x ) − Z2n (x ),



αm U n+1−m (x ) − U n−m (x ) ,

Z2n (x ) = β

m=1

n 



βm V n+1−m (x ) − V n−m (x ) .

m=1

Zrn (x )

Note that, = 0, (r = 1, 2 ) whenever n = 0. For approximation of spatial part, we now employ meshfree shape functions. RBFs and point interpolation method (PIM) are used to obtain these shape functions. Now according to PIM, assume U(x) and V(x) are given smooth functions which are approximated at the point of interest x as: M 

U (x ) =

M    μ j g x − x j + μM+1 = μ j g j (x ) + μM+1 ,

j=1

(2.5)

j=1

and

V (x ) =

M 

M    ν j h x − x j + νM+1 = ν j h j (x ) + νM+1 ,

j=1

(2.6)

j=1

where fj (x), gj (x) are given RBFs with  ·  as the Euclidean norm and μj , ν j are expansion coefficients. In the matrix-vector form above equations gives:

U = gϒ

⇒

ϒ = g−1 U and V = h

= h−1 V,

⇒

where

U = [U1 , . . . , UM , 0]T , V = [V1 , . . . , VM , 0]T , with U j = U (x j ), V j = V (x j ), and ϒ = [μ1 , . . . , μM+1 ]T , = [ν1 , . . . , νM+1 ]T are the vectors of expansion coefficients. Distinct collocation points always results in non-singularity of g and h [19]. The matrices g, h entries are listed as:

g=



{gi j }1≤i, j≤M

0 , 1

0T

h=

{hi j }1≤i, j≤M 0T



0 , 1

where 0T = [0, . . . , 0]M×1 . Now one can rewrite Eqs. (2.5) and (2.6) in the form:

U ( x ) = U =

M+1 

ω j (x )U j , V (x ) = V =

j=1

M+1 

χ j (x )V j ,

(2.7)

j=1

where UM+1 = VM+1 = 0. These shape functions possess Kronecker delta function property, i.e.



ω j (xi ) = δiωj =

1,



i = j, and

0,

i = j,

χ j (xi ) = δiχj =

1,

i = j,

0,

i = j.

(2.8)

The Kronceker delta function property makes the implementation of boundary conditions easy. Also, these satisfy partition of unity property, viz. M+1  j=1

ω j (x ) = 1 and

M+1 

χ j ( x ) = 1,

j=1

where it is not necessary 0 ≤ ωj (x), χ j (x) ≤ 1 (see [17]).

324

M. Hussain et al. / Applied Mathematics and Computation 341 (2019) 321–334

Moreover, approximation to the space derivatives is obtained as:

∂xl U (x ) =

M+1 

∂xl ω j (x )U j = Dl,u ( x )U, j

(2.9)

∂xl χ j (x )V j = Dl,j v (x )V,

(2.10)

j=1

and

∂xl V (x ) =

M+1  j=1

l where ∂xl = ∂ l and Dl,u , Dl,j v are the lth order explicit differentiation matrices corresponds to U, V components respectively. j ∂x Utilizing Eqs. (2.5)–(2.10) into Eqs. (2.3) and (2.4), there comes

M+1 





α ω j (x ) − p1 dn (x )ω j (x ) − pD3j ,u (x ) − p1 cn (x )D1j ,u (x ) U jn+1

j=1

=

M+1 



α ω j (x ) − p1 dn (x )ω j (x ) U jn + Y n+1 (x ),

(2.11)

j=1 M+1 



M+1 



β χ j (x ) − qD3j ,v (x ) V jn+1 = β χ j (x ) − q1 cn (x )D1j ,v (x ) V jn + Z n+1 (x ).

j=1

(2.12)

j=1

For collocation, let {xi = a + (i − 1 ) × h}M be the given collocation points such that x1 = a, xM = b and h is the separation i=1 distance between two points. Then collocating the above equations at xi ’s and using Eq. (2.8) together with boundary conditions, we have

(α − p1 dn (xi ) )Uin+1 −

M+1 









pD3i j,u + p1 cn (xi ) ∗ D1i j,u U jn+1 = α δiωj − p1 dn (x ) Uin + Rn+1 (xi ),

(2.13)

j=1

β Vin+1 − q

M+1 

D3i j,vV jn+1 = β Vin − q1

j=1

where

∗

R

M+1 

cn (xi ) ∗ D1i j,vV jn + S n+1 (xi ),

(2.14)

j=1

denotes element-wise multiplication of a vector and matrix, and n+1

(xi ) = [χan+1 , Y n+1 (x2 ), . . . , Y n+1 (xM ), χbn+1 , 0]T , S n+1 (xi ) = [γan+1 , Z n+1 (x2 ), . . . , Z n+1 (xM ), γbn+1 , 0]T .

By rewriting Eqs. (2.13) and (2.14) in compact matrix-vector form

Au Un+1 = Bu Un + Rn+1 ,

Av Vn+1 = Bv Vn + Sn+1 ,

n = 0, 1, . . . , N,

(2.15)

where Ad , Bd , (d = u, v ) are the matrices corresponds to components U, V. From Eqs. (2.15) the solution vectors are advanced from time t = tn to t = tn+1 in iterative manner. The iteration is started by providing the initial solution vector via initial conditions as:

U0 = {U0 (xi )}M i=1 ,

V0 = {V0 (xi )}M i=1

3. Stability Analysis −1 In this part, we analyze stability of the proposed computational method. Assume A−1 u , Av ’s exist then Eqs. (2.13) gives n −1 n+1 Un+1 = A−1 u Bu U + Au R

and

n −1 n+1 Vn+1 = A−1 , v Bv V + Av S

∀n ≥ 0.

(3.1)

Assume Uexact and Vexact be the exact solutions, then error vectors are defined as

Eu = |Uexact − U| and

Ev = |Vexact − V|

(3.2)

From Eqs. (3.1) and (3.2), one can write

Enu+1 = MEnu

and

Env +1 = NEnv ,

∀n ≥ 0.

(3.3)

−1 Now from Eqs. (3.3) one can easily infer that the method is stable if the amplification matrices M = A−1 u Bu , N = Av Bv satisfy the Lax–Richtmyer criterion of stability, i.e.

M. Hussain et al. / Applied Mathematics and Computation 341 (2019) 321–334

325

M ≤ 1 and N ≤ 1. For the amplification matrices to be normal, M = ρ (M ), N = ρ (N ) holds, where ρ denotes spectral radius of a matrix. However, the inequalities

ρ (M ) ≤ M and ρ (N ) ≤ N always remain true. From the above two inequalities it is obvious that

ρ (M ) ≤ 1 and ρ (N ) ≤ 1. Thus for the numerical scheme (3.1) to be stable spectral radii of matrices M, N should be less than or equal to unity [31]. The elements of these matrices mainly depend on the ratio [34]

hk

τ

,

where k being the order of highest spatial derivative in a given PDE. Thus to keep this ratio constant, for h to sufficiently small one should have τ → 0. In the next section, we shall computationally explore the stability criterion of the proposed method for shape parameter c dependent RBFs by keeping h and τ fixed. 4. Computational results This part shows the method application for the problem considered in Section 1. Approximation quality of the obtained solutions are measured via discrete error norms, read as:



E∞ (L ) = max |L (xi ) − L(xi )|,

E2 ( L ) =

∗

i

M 

(L∗ (xi ) − L(xi ) )2 ,

i=1

 Erms (L ) =

h

(1/M )

M 

(L∗ (xi ) − L(xi ) )2 (L = U or V ),

i=1

where L∗ , L denotes the exact and computed solutions respectively. The RBFs used in this work, are defined as follow: Multiquadric (MQ):

gi j =



ri2j + c12 ,

hi j =



ri2j + c22 ,

Shifted Thin Plate Splines (STPS):



gi j = ri2j + c12 Soliton: 2

gi j = sech



2







log ri2j + c12 ,



c1 ri j ,

hi j = ri2j + c22 2

hi j = sech



2





log ri2j + c22 ,



c2 ri j ,

where ri j = {xi − x j }M , and c1,2 ∈ R+ are the shape parameters. Accuracy and numerical stability of the method mainly i, j=1 depends on this parameter. To explore stability of the proposed scheme, we let

ci = εi × h,

i = 1, 2

where ε i ’s are given positive constants. For verification of stability of the proposed method, spectral radii ρ ’s are computed for given h and τ . The intervals of stability and corresponding error norms have been reported in Figs. 1 and 2 for MQ and STPS RBFs, while taking h = 0.1, τ = 0.05/4 and α = β = 0.3, 0.5, 0.95. From these figures one can see that MSRPIM remains stable i.e. the condition of stability ρ ≤ 1 is satisfied. Hence produce convergent solution whenever ε i ’s lie within the stable interval. Also, from these figure one can easily observe that matrix N corresponds to component V has relatively large stability interval. Anyhow, for simplicity, we have taken ε i ’s the same for both U, V such that the matrices M, N have stable spectrum. For all the test problems, we computed the solutions in the interval [0,1] for both spatial and temporal part, unless stated otherwise. Problem 1. Consider the TFCKdV Eqs. (1.1)–(1.4) with p = q = −1, p1 = −6, p2 = q1 = 3, and

Y (x, t ) = 3xt 4 +

2xt 2−α ,

(3 − α )

Z (x, t ) = 3xt 4 +

2xt 2−β .

(3 − β )

Initial and boundary conditions are adjusted according to the exact solution U (x, t ) = V (x, t ) = xt 2 . This problem is solved

326

M. Hussain et al. / Applied Mathematics and Computation 341 (2019) 321–334

Fig. 1. Graphs of spectral radii’s and corresponding error norms for problem 1 using MQ RBF.

for h = 0.1, τ = 0.05/4 using MQ (ε1 = ε2 = 27.4; ε1 = ε2 = 24.8; ε1 = ε2 = 20.5 ) and STPS (ε1 = ε2 = 24.6; ε1 = ε2 = 28.8; ε1 = ε2 = 15 ) for α = β = 0.3, 0.5, 0.95 respectively. Computer simulation is run up to time T = 1 and error norms are recorded in Table 1 at various time-levels. Moreover, in Table 2 maximum errors are reported using MQ RBF for time-step lengths τ = 0.05/2, 0.05/4, 0.05/6. From the given tables one can observed that MQ yields comparatively better accuracy in contrast to STPS RBF for the same set of parameters. Also, in Figs. 1 and 2 stability plots are displayed for various values of α , β . Fig. 3 shows solution profiles and absolute errors in space-time, which shows that the proposed method gives

M. Hussain et al. / Applied Mathematics and Computation 341 (2019) 321–334

327

Fig. 2. Graphs of spectral radii’s and corresponding error norms for problem 1 using STPS RBF.

satisfactory accuracy and also produce almost zero error at boundaries (see Fig. 3: (c) and (d)). Finally for α = 0.3, β = 0.5 and t = 1, Tables 3 and 4 shows time and space convergence of MSRPIM respectively. By varying τ ∈ {0.05, 0.05/2, 0.05/3, 0.05/4, 0.05/5, 0.05/6}, we have documented error norms in Table 3 showing decreasing time-step size results in solutions accuracy improvement. Similarly on varying M ∈ {2, 4, 6, 8, 10} error norms are shown in Table 4, clearly indicates that increasing number of nodal points yield better accuracy. It is clear from these tables that the proposed method gives convergence towards exact solution and yields acceptable accuracy for all values of τ and M.

328

M. Hussain et al. / Applied Mathematics and Computation 341 (2019) 321–334 Table 1 Error norms at various time-levels corresponds to problem 1. RBFs

Time t

MQ U

0.1 0.4 0.7 1.0 0.1 0.4 0.7 1.0 0.1 0.4 0.7 1.0 0.1 0.4 0.7 1.0

V

STPS U

V

α = β = 0.3

α = β = 0.5

α = β = 0.95

E∞

Erms

E∞

Erms

E∞

Erms

3.311E − 05 1.191E − 03 1.153E − 03 3.591E − 03 5.112E − 05 5.281E − 04 6.839E − 04 3.138E − 03 1.342E − 05 2.196E − 03 2.251E − 02 7.280E − 02 8.055E − 06 2.449E − 04 8.988E − 03 1.113E − 01

2.190E − 05 8.297E − 04 7.898E − 04 2.402E − 03 3.464E − 05 3.584E − 04 4.052E − 04 2.102E − 03 7.724E − 06 1.528E − 03 1.563E − 02 5.061E − 02 5.132E − 06 1.499E − 04 6.228E − 03 7.714E − 02

1.420E − 04 2.987E − 03 6.043E − 03 4.050E − 03 1.440E − 04 3.224E − 03 5.438E − 03 3.332E − 03 2.270E − 05 2.173E − 03 1.193E − 02 5.177E − 03 2.099E − 05 1.450E − 03 7.797E − 04 4.556E − 03

9.706E − 05 2.085E − 03 4.198E − 03 2.785E − 03 9.826E − 05 2.257E − 03 3.790E − 03 2.141E − 03 1.281E − 05 1.519E − 03 8.279E − 03 3.478E − 03 1.286E − 05 1.003E − 03 4.465E − 04 2.907E − 03

2.931E − 04 4.499E − 03 1.232E − 02 2.008E − 02 2.951E − 04 4.906E − 03 1.304E − 02 1.607E − 02 1.400E − 04 6.569E − 03 2.867E − 02 7.539E − 02 1.425E − 04 6.282E − 03 2.493E − 02 3.892E − 02

1.919E − 04 3.140E − 03 8.561E − 03 1.388E − 02 1.929E − 04 3.432E − 03 9.101E − 03 1.123E − 02 8.369E − 05 4.585E − 03 1.993E − 02 5.226E − 02 8.540E − 05 4.394E − 03 1.737E − 02 2.711E − 02

Table 2 Maximum error norms using MQ RBF for different values of τ corresponds to problem 1. Time

τ

α = β = 0.3

α = β = 0.5

α = β = 0.95

E∞ (U)

E ∞ (V)

E∞ (U)

E ∞ (V)

E∞ (U)

E ∞ (V)

t = 0.5

0.05/2 0.05/4 0.05/6

8.339E − 04 1.318E − 03 1.579E − 03

5.452E − 04 8.219E − 04 1.257E − 03

4.446E − 03 4.200E − 03 3.738E − 03

4.920E − 03 4.368E − 03 4.323E − 03

7.606E − 03 6.973E − 03 6.768E − 03

9.167E − 03 7.710E − 03 7.251E − 03

Table 3 Time-convergence for problem 1 using different τ ’s when t = 1, h = 0.2, α = 0.3, β = 0.5.

τ

U

0.05 0.05/2 0.05/3 0.05/4 0.05/5 0.05/6

V

E∞

E2

Erms

E∞

E2

Erms

3.592E − 02 2.567E − 02 1.739E − 02 1.267E − 02 9.717E − 03 7.736E − 03

2.599E − 02 1.894E − 02 1.275E − 02 9.246E − 03 7.060E − 03 5.601E − 03

2.372E − 02 1.729E − 02 1.164E − 02 8.440E − 03 6.445E − 03 5.113E − 03

4.553E − 01 4.767E − 02 4.005E − 02 3.428E − 02 2.989E − 02 2.650E − 02

3.377E − 01 3.552E − 02 2.993E − 02 2.565E − 02 2.239E − 02 1.986E − 02

3.083E − 01 3.242E − 02 2.733E − 02 2.342E − 02 2.044E − 02 1.813E − 02

Table 4 Space-convergence for problem 1 using different M’s when t = 1, τ = 0.05/6, α = 0.3, β = 0.5. M

U

V

E∞

E2

Erms

E∞

E2

Erms

2.0 4.0 6.0 8.0 10.0

3.774E − 02 2.915E − 02 1.614E − 02 1.338E − 02 7.469E − 04

2.668E − 02 2.115E − 02 1.170E − 02 9.758E − 03 4.712E − 04

2.179E − 02 1.892E − 02 1.083E − 02 9.199E − 03 4.493E − 04

6.585E − 02 3.384E − 02 1.406E − 02 1.074E − 03 8.493E − 04

4.656E − 02 2.463E − 02 1.029E − 02 6.548E − 04 5.432E − 04

3.802E − 02 2.203E − 02 9.528E − 03 6.173E − 04 5.180E − 04

Problem 2. As a second example, we reconsider the TFCKdV Eqs. (1.1)–(1.4) with p = q = −1, p1 = −6, p2 = q1 = 3, and

Y (x, t ) = 3xt 5 + xt ( 2 −α ) 5

( 72 ) ,

( 72 − α )

Z (x, t ) = 3xt 5 + xt ( 2 −β ) 5

( 72 ) .

( 72 − β )

√ Initial and boundary conditions are obtained from exact solutions U (x, t ) = V (x, t ) = x t 5 . On setting h = 0.1, τ = 0.01 and α = β = 0.3, 0.5, 0.9, the solutions are computed via MQ (ε1 = ε2 = 27.4, ε1 = ε2 = 24.8, ε1 = ε2 = 20.5 ) and STPS (ε1 = ε2 = 24.4, ε1 = ε2 = 23.5, ε1 = ε2 = 15 ) for the mentioned α , β ’s respectively. The computed error norms are listed in Table 5 at times-levels t = 0.1, 0.5, 1, showing that the current method produce reasonable accuracy whilst using limited number of nodal points M = 10. Also, in Table 6 E∞ , E rms are recorded using MQ RBF for α = β = 0.7 along with error rates

M. Hussain et al. / Applied Mathematics and Computation 341 (2019) 321–334

329

Fig. 3. Solution profiles and absolute errors for problem 1 using MQ RBF. Table 5 Error norms at various time-levels corresponds to problem 2. RBFs

MQ U V

STPS U V

Time t

0.1 0.5 1.0 0.1 0.5 1.0 0.1 0.4 1.0 0.1 0.4 1.0

α = β = 0.3

α = β = 0.5

α = β = 0.95

E∞

Erms

E∞

Erms

E∞

Erms

1.697E − 05 1.266E − 03 5.279E − 03 1.906E − 05 5.523E − 04 3.087E − 03 1.011E − 04 4.223E − 03 7.430E − 03 1.016E − 04 4.579E − 03 9.040E − 03

1.154E − 05 8.815E − 04 3.559E − 03 1.303E − 05 3.733E − 04 2.015E − 03 7.058E − 05 2.939E − 03 5.119E − 03 7.099E − 05 3.194E − 03 6.219E − 03

4.130E − 05 3.242E − 03 3.171E − 03 4.145E − 05 3.483E − 03 2.810E − 03 4.929E − 05 3.152E − 03 3.868E − 03 4.948E − 05 2.862E − 03 3.613E − 03

2.796E − 05 2.262E − 03 2.145E − 03 2.805E − 05 2.439E − 03 1.671E − 03 3.364E − 05 2.200E − 03 2.655E − 03 3.376E − 05 2.008E − 03 2.247E − 03

9.879E − 05 5.041E − 03 2.074E − 02 9.901E − 05 5.513E − 03 1.618E − 02 8.872E − 05 7.650E − 03 7.326E − 02 8.900E − 05 7.297E − 03 3.492E − 02

5.991E − 05 3.518E − 03 1.433E − 02 6.004E − 05 3.857E − 03 1.126E − 02 5.656E − 05 5.338E − 03 5.079E − 02 5.674E − 05 5.102E − 03 2.436E − 02

Table 6 Error norms and Error rates for problem 2 when α = β = 0.7. Time

t = 0.5

τ 0.05 0.05/2 0.05/4

U

V

E∞

Erms

E.R.

E∞

Erms

E.R.

4.835E − 03 4.710E − 03 4.604E − 03

3.375E − 03 3.290E − 03 3.217E − 03

– 1.026 1.023

7.040E − 03 5.778E − 03 5.118E − 03

4.939E − 03 4.050E − 03 3.585E − 03

– 1.218 1.128

330

M. Hussain et al. / Applied Mathematics and Computation 341 (2019) 321–334 Table 7 Comparison of error norms for problem 3.

α=β

t

U E2

m = 4 [1]

E2

m = 4 [1]

1.0

0.1 0.5 1.0

2.46235373E − 03 8.96035397E − 03 5.07910497E − 03

8.44695952E − 03 – –

2.16027970E − 03 5.76274881E − 03 1.23437273E − 03

8.44695952E − 03 – –

V

Fig. 4. Graphs of spectral radii’s for problem 3 using MQ RBF.

(E.R.) in time computed using the formula:

E.R. =

E∞ (τi ) E∞ (τi+1 )

where E∞ (τ i ) is the maximum error observed at time-step τ i . From the given tables one can observed that MQ yields comparatively better accuracy in contrast to STPS RBF for the same set of parameters. Problem 3. [1] Consider again the TFCKdV Eqs. (1.1)–(1.4) with p1 = −6, p2 = q1 = 3, and Y (x, t ) = Z (x, t ) = 0. The initial conditions are set as:

U (x, 0 ) = V (x, 0 ) =

4c2 ecx

(1 + ecx )2

,

(c = p, q ),

whereas boundary conditions are taken as:

4c2 e−c t , (1 + e−c3 t )2 3

U ( 0, t ) = V ( 0, t ) =

3

U ( 1, t ) = V ( 1, t ) =

4c2 ec−c t , (1 + ec−c3 t )2

(c = p, q ).

Exact solution is known only for the case α = β = 1 and given as:

4 p2 e px−p t , (1 + e px−p3 t )2 3

U (x, t ) =

4q2 eqx−q t . (1 + eqx−q3 t )2 3

V (x, t ) =

We solved this problem for p = q = 1 and τ = 0.01, h = 0.1. In Table 7, E2 error norms for U, V are computed and compared with earlier work at times t = 0.1, 0.5, 1.0 when α = β = 1 using MQ RBF ε1 = ε2 = 20.5. The table shows the current method produce relatively better accuracy in comparison to that of [1]. Also, as the exact solution is not known for various values of α and β . Thus to compute the solutions for α = β = 0.3, 0.5, we go for the stability criterion of the proposed method. As shown in Fig. 4, the method remains stable whenever 18.4 < ε 1 , ε 2 < 19.5 for MQ RBF. Hence the obtained solutions are convergent and trustful for any value of ε 1 , ε 2 chosen from the stability interval. Measure of goodness of the computed solutions are considered as absolute difference of the consecutive time-levels defined as:

|U n+1 (x ) − U n (x )| and |V n+1 (x ) − V n (x )|. The solutions are computed for α = β = 0.3, 0.5(ε1 = ε2 = 19 ) using h = 0.1, τ = 0.1 and T = 0.5. The 3D graphs of solutions components U, V as well their absolute differences at consecutive time-levels are pictured in Fig. 5.

M. Hussain et al. / Applied Mathematics and Computation 341 (2019) 321–334

331

Fig. 5. 3D solutions graphs and corresponding absolute errors for problem 3. Table 8 Comparison of maximum errors at various time-levels corresponds to problem 4 when α = β = 1. MQ-RBF

Time t

M = 5, h = 20 Present

[15]

Present

[15]

Present

[15]

U

0.1 0.5 1.0 1.5 2.0 0.1 1.0 2.0

8.968E − 07 4.487E − 06 8.981E − 06 1.348E − 05 1.799E − 05 1.118E − 05 1.121E − 04 2.246E − 04

– 6.22E − 05 1.25E − 04 1.87E − 04 2.49E − 04 – – –

2.191E − 06 1.099E − 05 2.208E − 05 3.327E − 05 4.454E − 05 1.124E − 05 1.127E − 04 2.260E − 04

– 6.23E − 05 1.25E − 04 1.87E − 04 2.49E − 04 – – –

2.205E − 06 1.107E − 05 2.226E − 05 3.355E − 05 4.496E − 05 1.205E − 05 1.205E − 04 2.409E − 04

– 6.01E − 05 1.20E − 04 1.81E − 04 2.41E − 04 – – –

V

M = 10, h = 10

M = 20, h = 5

Problem 4. [15] Now consider the homogenous TFCKdV Eqs. (1.1)–(1.4) together with p1 = 6 p, p2 = 2ν and q1 = 3. The initial conditions are set as:

U0 (x ) =

−c2 (1 + p) 4c2 ecx + , 3 + 6p (1 + ecx )2

V0 (x ) =

decx

(1 + ecx )2

,

whereas boundary conditions are taken as:

χl (t ) =

−c2 (1 + p) 4c2 ec(l+d1 t ) + , 3 + 6p (1 + ec(l+d1 t ) )2

γl (t ) =

dec(l+d1 t ) , (1 + ec(l+d1 t ) )2

(l = a, b).

Exact solution is known only for the case α = β = 1 and given as:

U (x, t ) =

−c2 (1 + p) 4c 2 ec ( x+d1 t ) + , 3 + 6p (x + ec(l+d1 t ) )2

V (x, t ) =

dec(x+d1 t )

( 1 + ec ( x+d1 t ) )2

,

where c is a constant, and



d = −c

2

−24 p

ν

,

d1 =

−pc . 1 + 2p

For comparison purpose, the constants are taken as ν = c = 0.1, p = −1.5, q = −1, q1 = 3, τ = 0.1 and α = β = 1. To validate effectiveness of the proposed method, this problem is solved over a large space interval x ∈ [−50, 50] for various values of spatial step-size h. The obtained maximum error is presented in Table 8 in comparison to the work of Khater et al. [15]. In Table 9 error norms are computed and compared to [15] for M = 15 while keeping other parameters fixed. From these tables,

332

M. Hussain et al. / Applied Mathematics and Computation 341 (2019) 321–334 Table 9 Comparison of error norms at various time-levels corresponds to problem 4 when α = β = 1, M = 15. U/V

Time t

STPS E∞

Erms

E∞

Erms

E∞

Erms

U

0.5 1.0 1.5 2.0 0.5 1.0 1.5 2.0

1.079E − 05 2.169E − 05 3.269E − 05 4.379E − 05 5.852E − 05 1.171E − 04 1.756E − 04 2.342E − 04

5.218E − 06 1.043E − 05 1.565E − 05 2.087E − 05 3.577E − 05 7.154E − 05 1.073E − 04 1.431E − 04

1.110E − 05 2.232E − 05 3.365E − 05 4.509E − 05 5.824E − 05 1.165E − 04 1.748E − 04 2.331E − 04

5.340E − 06 1.068E − 05 1.602E − 05 2.136E − 05 3.543E − 05 7.087E − 05 1.063E − 04 1.417E − 04

6.10E − 05 1.22E − 04 1.83E − 04 2.42E − 04 – – – –

– – – – – – – –

V

MQ

[15]

Fig. 6. Graphs of spectral radii’s and corresponding solutions profiles for problem 4 using STPS RBF.

we see that the current MSRPIM produce very good accuracy in contrast to the cited work. Also, solutions are computed using STPS RBF for h = 5, τ = 0.1, T = 2 when α = β = 0.5, 0.9. These solutions are obtained once again on the basis of stability criterion of the proposed numerical scheme which remains stable for 17 < ε 1, 2 < 18.4 (see Fig. 6: (a),(d)). In the same figure solutions profiles are pictured for α = β = 0.5, 0.9 using STPS (ε1 = ε2 = 17.2), showing standing solitary wave prorogation as time marches towards final time. Problem 5. [2] As a final test problem, consider the homogenous TFCKdV Eqs. (1.1)–(1.4) together with p = −6μ, p1 = 2ν, p2 = −μ and q = −3ν, q1 = ν . The initial and boundary conditions are set as:

U0 (x ) =

d

μ

 

2

sech

1 2



d

μ

x ,

V0 (x ) =

d



2μ

 

2

sech

1 2



d

μ

x ,

and

χl (t ) = 0, γl (t ) = 0, (l = a, b). where d is a constant. Exact solution is unknown for this problem. The initial conditions represent solitary waves that travels with fixed amplitudes Ampu = dp and Ampv = √d . For computation purpose, the parameters are taken as ν = 1, μ = 2p

1, d = 1. The solution is obtained for α = 0.75, β = 0.45 via Soliton RBF using ε1 = ε2 = 0.0 0 0 04 and τ = 0.1, h = 0.5 over the space interval x ∈ [−10, 10]. Computer simulations are run up to t = 2 and obtained solution profiles are pictured in Fig. 7. This√figure shows prorogation of solitary wave as time marches. The true values of the amplitudes are Ampu = 1 and Ampv = 1/ 2. From the given figure we see that they travel with fixed amplitude. For verification, amplitudes of the two

M. Hussain et al. / Applied Mathematics and Computation 341 (2019) 321–334

333

Fig. 7. Graphs of solution profiles for problem 5 using Soliton RBF. Table 10 Amplitudes of solitons at various time levels using α = 0.75, β = 0.45 for problem 5. Time t

0.0

0.1

0.3

0.5

0.7

0.9

1.0

1.1

1.3

1.5

1.7

1.9

2.0

Ampu Ampv

1.0 0 0 0.707

0.999 0.707

0.999 0.706

1.0 0 0 0.706

1.0 0 0 0.706

1.001 0.707

1.001 0.707

1.002 0.707

1.003 0.707

1.004 0.707

1.006 0.707

1.008 0.708

1.009 0.708

solitons are computed at various time levels and documented in Table 10. The tabulated values of the amplitudes remain constant as the time advances. This confirm capability of the proposed method to produce solution of a problem having no exact solution. 5. Conclusion In this paper, MSRPIM was applied successfully for the numerical solutions of a class of TFCKdV equations. The numerically computed solutions were matched with exact solutions as well to some earlier works showing very good agreement for various α , β ’s values. Numerical convergence of the current scheme in space and time was also studied through variation of time step-size τ and number of nodal points M. Also, stability of the proposed method is established theocratically and justified computationally. The outcomes of the present methodology shows potential applicability in obtaining accurate solutions of various other important time-dependent PDEs and FPDEs. Acknowledgments The authors are thankful for the comments of the worthy reviewers which improved quality of the manuscript. The first author acknowledges the financial support of GIK Institute for his Ph.D. studies. References [1] B. Albuohimad, H. Adibi, S. Kazem, A numerical solution of time-fractional coupled Korteweg-de vries equation by using spectral collocation method, Ain Shams Eng. J. (2017), doi:10.1016/j.asej.2016.10.010 (in press). [2] A. Atangana, A. Secer, in: The time-fractional coupled-Korteweg-de-vries equations, 2013, 2013, p. 947986. [3] M. Caputo, Linear models of dissipation whose q is almost frequency independent, part II, J. R. Astral. Soc. 13 (1967) 529–539. [4] S.A. El-Wakil, E.M. Abulwafa, E.K. El-Shewy, A.A. Mahmoud, Time-fractional KDV equation for plasma of two different temperature electrons and stationary ion, Phys. Plasmas. 18 (2011a) 092116. [5] S.A. El-Wakil, E.M. Abulwafa, M.A. Zahran, A.A. Mahmoud, Time-fractional KDV equation: formulation and solution using variational methods, Nonlinear Dyn. 65 (2011b) 55–63. [6] Y. Fujita, Cauchy problems of fractional order and stable processes, Jpn. J. Appl. Math. 7 (3) (1990) 459–476. [7] Y. Gu, G. Liu, A meshless local Petrov-Galerkin (MLPG) method for free and forced vibration analyses for solids, Comput. Mech. 27 (2011) 188–198. [8] S. Haq, M. Hussain, Selection of shape parameter in radial basis functions for solution of time-fractional black-scholes models, Appl. Math. Comput. 335 (2018) 248–263. [9] S. Haq, M. Ishaq, Solution of coupled Whitham-Broer-Kaup equations using optimal homotopy asymptotic method, Ocean Eng. 84 (2014) 81–88. [10] R. Hilfer, Applications of Fractional Calculus in Physics, World Scientific Publishing, River Edge, USA., 20 0 0a. [11] R. Hilfer, Fractional diffusion based on Riemann–Liouville fractional derivative, J. Phys. Chem 104 (20 0 0b) 3914–3917. [12] R. Hilfer, Foundations of fractional dynamics, Fractals 3 (3) (1995) 549–556. [13] R. Hirota, J. Satsuma, Soliton solutions of a coupled Kortewegde vries equation, Phys. Lett. A 85 (1981) 407–408. [14] E.J. Kansa, Multiquadrics–a scattered data approximation scheme with application to computation fluid dynamics, II: Solutions to hyperbolic, parabolic, and elliptic partial differential equations, Comput. Math. Appl. 19 (1990) 149–161.

334

M. Hussain et al. / Applied Mathematics and Computation 341 (2019) 321–334

[15] A.H. Khatera, R.S. Temsah, D.K. Callebaut, Numerical solutions for some coupled nonlinear evolution equations by using spectral collocation method, Math. Comput. Model. 48 (2008) 1237–1253. [16] A. Kilbas, H.M. Srivastava, J.J. Trujillo, Theory and applications of fractional differential equations, in: North-Holland Mathematics Studies, Elsevier Science, Amsterdam, the Netherlands, 2006. [17] G.R. Liu, T.Y. Gu, An Introduction to Meshfree Methods and Their Programming, Springer Press, Berlin, 2005. [18] J.C. Liu, G.L. Hou, New approximate solution for time-fractional coupled KDV equations by generalised differential transform method, Chin. Phys. B 19 (11) (2010) 110203. [19] C.A. Micchelli, Interpolation of scattered data: distance matrix and conditionally positive definite functions, Constr. Approx. 2 (1986) 11–22. [20] L. Ostrovsky, A. Stepanyants-Yu, Do internal solutions exist in the ocean? Rev. Geophys. 27 (1989) 293–310. [21] I. Podlubny, U. San Diego, Fractional Differential Equations, Academic Press, 1999. [22] E. Shivanian, A new spectral meshless radial point interpolation (SMRPI) method: a well-behaved alternative to the meshless weak forms, Eng. Anal. Bound. Elem. 54 (2015) 1–12. [23] E. Shivanian, Analysis of the time factional 2-d diffusion-wave equation via moving least square (MLS) approximation, Int. J. Appl. Comput. Math. 3 (2017a) 2447–2466. [24] E. Shivanian, Local radial basis function interpolation method to simulate 2d fractional-time convection-diffusion-reaction equations with error analysis, Numer. Methods Partial Differ. Equ. 33 (2017b) 974–994. [25] E. Shivanian, A. Jafarabadi, An improved spectral meshless radial point interpolation for a class of time-dependent fractional integral equations: 2d fractional evolution equation, J. Comput. Appl. Math. 325 (2017a) 18–33. [26] E. Shivanian, A. Jafarabadi, Error and stability analysis of numerical solution for the time fractional nonlinear schrȵdinger equation on scattered data of general-shaped domains, Numer. Methods Partial Differ. Equ. 33 (2017b) 1043–1069. [27] E. Shivanian, Jafarabadi, An improved meshless algorithm for a kind of fractional cable problem with error estimate, Chaos Solitons Fractals 110 (2018a) 138–151. [28] E. Shivanian, A. Jafarabadi, The spectral meshless radial point interpolation method for solving an inverse source problem of the time-fractional diffusion equation, Appl. Numer. Math. 129 (2018b) 1–25. [29] E. Shivanian, A. Jafarabadi, The numerical solution for the time-fractional inverse problem of diffusion equation, Eng. Anal. Bound. Elem. 91 (2018c) 50–59. [30] E. Shivanian, A. Jafarabadi, Analysis of the spectral meshless radial point interpolation for solving fractional reaction-subdiffusion equation, J. Comput. Appl. Math. 336 (2018d) 98–113. [31] E. Shivanian, A. Jafarabadi, Capillary formation in tumor angiogenesis through meshless weak and strong local radial point interpolation, Eng. Comput. 34 (2018e) 603–619. [32] H. Triki, T. Ak, S. Moshokoa, A. Biswas, Soliton solutions to KDV equation with spatio-temporal dispersion, Ocean Eng. 114 (2016) 192–203. [33] H. Triki, D. Milovic, A. Biswas, Solitary waves and shock waves of the KDV6 equation, Ocean Eng. 73 (2013) 119–125. [34] M. Uddin, S. Haq, RBF approximation method for time fractional partial differential equations, Commun. Nonlinear Sci. Numer. Simul. 16 (11) (2011) 4208–4214. [35] B.J. West, M. Bologna, P. Grigolini, Physics of Fractal Operators, Springer, New York, USA., 2003.