Numerical solution of time-fractional coupled Korteweg–de Vries and Klein–Gordon equations by local meshless method

(1)

Numerical solution of time-fractional coupled Korteweg–de Vries and Klein–Gordon equations by local meshless method

MUHAMMAD NAWAZ KHAN¹, IMTIAZ AHMAD², ALI AKGÜL³ ^,∗, HIJAZ AHMAD¹and PHATIPHAT THOUNTHONG⁴

1Department of Basic Sciences, University of Engineering and Technology Peshawar, Peshawar, Pakistan

2Department of Mathematics, University of Swabi, Swabi, Pakistan

3Department of Mathematics, Art and Science Faculty, Siirt University, Siirt, Turkey

4Renewable Energy Research Centre, Department of Teacher Training in Electrical Engineering, Faculty of Technical Education, King Mongkut’s University of Technology North Bangkok, 1518 Pracharat 1 Road, Bangsue, Bangkok10800, Thailand

∗Corresponding author. E-mail: aliakgul00727@gmail.com

MS received 16 June 2020; revised 25 July 2020; accepted 18 August 2020

Abstract. This article provides numerical simulations of the time-fractional coupled Korteweg–de Vries and Klein–Gordon equations via the local meshless collocation method (LMCM) utilising the radial basis functions. The recommended local meshless technique is utilised for the space derivatives of the models whereas Caputo fractional definition is used for time-fractional derivative. Numerical experiments are performed for one-dimensional coupled Korteweg–de Vries and two-dimensional Klein–Gordon equations. In order to verify the efficiency and accuracy of the proposed meshless method, numerical results are compared with exact and numerical techniques reported in recent literature which reveals that the method is computationally attractive and produces better results.

Keywords. Meshless collocation method; radial basis functions; coupled Korteweg–de Vries equations; Klein–

Gordon equation; irregular domain.

PACS Nos 02.60.−x; 02.30.Jr; 02.90.+p

1. Introduction

The theory of nonlinear evolution equations (NLEEs) has been applied in numerous areas of theoretical physics, applied mathematics and engineering [1].

In this way, these equations need to be studied in more depth. So far, these NLEEs have been examined for integer-order evolution or integral-order temporal derivative. Now is the time to summarise these evo- lutionary terms in a fractional order and to derive the evolution equations of a fractional order. Fractional cal- culus [2] is an extension of derivatives and integrals of fractional order. Fractional partial differential equations (FPDEs) have recently been applied in various areas of materials science, engineering, chemistry, dynam- ics systems and statistical mechanics [2–10]. Taking time limits into account, the FPDEs containing time- fractional derivatives are very important because they

have infinitesimal generating nature in evolution pro- cesses, state stability and equilibrium [11]. Caputo’s sense of partial derivatives provides a more realistic modelling of phenomena in the real world and simpli- fies the initial conditions in the formulation of problems like the integer case [12]. The FPDEs consistently show up as a difficult assignment with regard to getting their analytical solutions. Therefore, numerical methods are always interesting and desirable tools that can be utilised to approximate solution to these complex equations. This article presents some of the NLEEs, par- ticularly coupled Korteweg–de Vries (KdV) equation and Klein–Gordon equation, that occur in engineering and theoretical physics in fractional evolution, which are considered for numerical solution by the proposed mechanism.

The one-dimensional time-fractional coupled KdV equation [13]

0123456789().: V,-vol

(2)

∂^αw1(x,t)

∂t^α =γ∂³w1

∂x³ +μw1∂w1

∂x +νw2∂w2

∂x + f₁(x,t),

∂^βw2(x,t)

∂t^β =γ∂³w2

∂x³ −μw1∂w2

∂x + f2(x,t), x ∈, t ≥0, αandβ ∈(0,1),

(1)

with initial and boundary conditions w1(x,0)=φ1(x), w2(x,0)=φ2(x), w1(x,t)=ψ1(x), w2(x,t)=ψ2(x),

x ∈∂, t ≥0.

The two-dimensional time-fractional Klein–Gordon equation [14]

∂^αw(x,y,t)

∂t^α

=w−w³+ f, (x,y)∈, t ≥0, α∈(1,2), (2) with initial and boundary conditions

w(x,y,0)=φ(x,y), ∂w(x,y,0)

∂t =ψ(x,y),

∂w

∂n(x,y,t)=0, (x,y)∈∂, t ≥0,

whereis the Laplacian differential operator,is the domain and∂is the boundary.

Much effort has recently been put into developing algorithm for the exact and numerical solutions of the FPDEs. In this work, several methods have been applied to the FPDE solution, including variational iteration method [15–18], homotopy analysis method [19,20], spectral method [21], Haar wavelet hybrid method [22], finite difference method [23], expansion method [24,25], method of the Laplace transform [26], Riccati transformation method [27] and meshless methods [28–

30].

Meshless techniques are often considered for the efficient and accurate numerical solution of almost all types of partial differential equation models. Partic- ularly, the meshless techniques utilising radial basis functions (RBFs) are the most attractive methods among these techniques. These meshless techniques utilise only the Euclidian distance between two points and do not need mesh in the computational domain contrast to the mesh-based method, for example, finite-element and finite-difference methods. Based on this fact, the meshless technique is a useful very flexible mathematical tool and can be applied to high-dimensional models with irregular and complex domains in various types of real- world problems [31–33].

Meshless method based on RBFs can be found in both local and global versions. The standard meshless

method based on globally radial functions is an effective technique but it leads to ill-conditioning dense system matrix and also the execution time increase with the increase in data points. To overcome these lim- itations, local version of the meshless technique has been recommended by the researchers [34,35]. The local technique uses only neighbouring data points and produces sparse matrix system which can be solved effi- ciently. Furthermore, the execution time of the local version is considerably less than the execution time of the global meshless technique. Recently, researchers utilised local meshless techniques for the efficient and accurate numerical simulations of complex PDE models [36–39].

Current research is devoted to the use of explicit time discretisation scheme combined with local meshless method (LMM) for the numerical investigation of the model equations (1) and (2). Both rectangular and non-rectangular computational domains with uniform, random, Chebyshev and Halton nodes are considered in numerical experiments.

2. Discretisation of space derivatives using the local meshless collocation method

The local meshless collocation method (LMCM) [40]

is utilised for time-fractional models described in §1.

The derivatives ofw1(¯z,t), w2(¯z,t)andw(¯z,t)at the centresz¯h are approximated by the neighbourhood of

¯

zh,{¯zh1,z¯h2,z¯h3, . . . ,¯zhn_h} ⊂ {¯z1,z¯2, . . . ,z¯Nⁿ},nh Nⁿ, whereh =1,2, . . . ,Nⁿ. In the case of one dimen- sion (1D), ¯z = x and for two dimensions (2D),z¯ = (x,y).

Now for the 1D case w^(m)₁ (xh)≈

n_h

k=1

λ^(m)_k w1(xhk),

w⁽₂^m⁾(xh)≈

nh

k=1

η⁽_k^m⁾w2(xhk), h =1,2, . . . ,N. (3) Substituting RBFψ(x −xp)in (3),

ψ⁽^m⁾(xh −xp)=

n_h

k=1

λ⁽_hk^m⁾ψ(xhk−xp),

ψ⁽^m⁾(xh −xp)=

n_h

k=1

η⁽_hk^m⁾ψ(xhk−xp),

p=h1,h2, . . . ,hn_h, (4) where

ψ(xhk−xp)=

1+(cxhk −xp)²,

(3)

ψ(xhk−xp)=exp(−xhk−xp²/c²) in the case of multiquadric (MQ) RBF.

Matrix form of (4) is

⎡

⎢⎢

⎣

ψ_h1⁽^m⁾(xh) ψ_h2⁽^m⁾(x_h)

...

ψ_hn⁽^m_h⁾(xh)

⎤

⎥⎥

⎦

ψ⁽_nh^m⁾

=

⎡

⎢⎢

⎣

ψh1(xh1) ψh2(xh1) · · · ψhn_h(xh1) ψh1(xh2) ψh2(xh2) · · · ψhn_h(xh2)

... ... ... ...

ψh1(x hn_h) ψh2(xhn_h)· · ·ψhn_h(x hn_h)

⎤

⎥⎥

⎦

A_nh

⎡

⎢⎢

⎣ λ⁽_h1^m⁾ λ⁽_h2^m⁾ ...

λ⁽_hn^m_h⁾

⎤

⎥⎥

⎦

λ⁽_nh^m⁾

,

⎡

⎢⎢

⎣

ψ_h1⁽^m⁾(xh) ψ_h2⁽^m⁾(xh)

...

ψ_hn⁽^m_h⁾(xh)

⎤

⎥⎥

⎦

ψ^(m)_nh

=

⎡

⎢⎢

⎣

ψh1(xh1) ψh2(xh1) · · · ψhn_h(xh1) ψh1(xh2) ψh2(xh2) · · · ψhn_h(xh2)

... ... ... ...

ψh1(x h_n_h) ψh2(x_hn_h)· · ·ψhn_h(x h_n_h)

⎤

⎥⎥

⎦

A_nh

⎡

⎢⎢

⎣ η⁽_h1^m⁾ η⁽_h2^m⁾ ...

η_hn⁽^m_h⁾

⎤

⎥⎥

⎦

η⁽_nh^m⁾

,

(5) where

ψp(xk)=ψ(xk−xp), p =h1,h2, . . . ,hnh, for eachk=i1,h2, . . . ,hnh. Equation (5) can be writ- ten as

ψ⁽n^m_h⁾=An_hλ⁽n^m_h⁾, ψ⁽n^m_h⁾=An_hη⁽n^m_h⁾, (6) the invertibility of the matrixAn_his guaranteed by [41].

From (6), we obtain

λ⁽_n^m_h⁾=A⁻_n_h¹ψ⁽_n^m_h⁾, η⁽_n^m_h⁾=A⁻_n_h¹ψ⁽_n^m_h⁾. (7) Equations (3) and (7) imply

w₁⁽^m⁾(x_h)=(λ⁽_n^m_h⁾)^T(w₁)nh, w₂⁽^m⁾(xh)=(η⁽n^m_h⁾)^T(w2)n_h,

where (w1)n_h =

w1(xh1), w1(xh2), . . . , w1(xhn_h)T

, (w2)n_h =

w2(xh1), w2(xh2), . . . , w2(xhn_h)T

. The above technique can be rehashed for two- and three- dimensional case.

2.1 Time-fractional discretisation

Several definitions of fractional derivatives are available in the literature including Riemann–Liouville derivative [42], Atangana–Baleanu–Caputo–Fabrizio [43], He and Li [44] and Caputo [12]. But in this paper we utilised Caputo definition [12] forα ∈(0,1)for time-fractional derivative∂^αw(¯z,t)/∂t^α is as follows:

∂^αw(¯z,t)

∂t^α

=

⎧⎪

⎨

⎪⎩ 1 (1−α)

_t

0

∂w(¯z, ϑ)

∂ϑ (t−ϑ)^−αdϑ, 0< α <1

∂w(¯z,t)

∂t , α=1.

Considert0,t1, . . . ,tQ as equally spaced time intervals for[0,T]wheret_q =qτ,q = 0,1,2, . . . ,Q,τ is the time-step size and the time-fractional derivative (q + 1)th time level is approximated as

∂^αw(¯z,t_q+1)

∂t^α

= 1 (1−α)

t_q₊₁ 0

∂w(¯z, ϑ)

∂ϑ (tq+1−ϑ)^−αdϑ

= 1 (1−α)

q s=0

₍s+1)τ sτ

∂w(¯z, ϑ)

∂ϑ (ts+1−ϑ)^−αdϑ

≈ 1 (1−α)

q s=0

₍_s₊₁_)τ

sτ

∂w(¯z, ϑs)

∂ϑ (ts+1−ϑ)^−αdϑ.

The term∂w(¯z, ϑs)/∂ϑ is approximated as

∂w(¯z, ϑs)

∂ϑ = w(¯z, ϑs+1)−w(¯z, ϑs)

ϑ +O(τ).

Then

(4)

∂^αw(¯z,tq+1)

∂t^α ≈ 1

(1−α) q s=0

w(¯z,ts+1)−w(¯z,ts) τ

×

₍s+1)τ

sτ (ts+1−ϑ)^−αdϑ

= 1 (1−α)

q s=0

w(¯z,tq+1−s)−w(¯z,tq−s) τ

×

₍_s₊₁_)τ

sτ (ts+1−ϑ)^−αdϑ

=

⎧⎪

⎪⎪

⎨

⎪⎪

⎪⎩ τ^−α

(2−α)(w^q⁺¹−w^q)+ τ^−α (2−α)

×^q

s=1

(w^q⁺¹⁻^s−w^q⁻^s)[(s+1)¹^−α−s¹^−α],q ≥1, τ^−α

(2−α)(w¹−w⁰). q =0. Leta0=[τ^−α/(2−α)]andbs=(s+1)^1−α−s^1−α,s= 0,1, . . . ,q, where in simplified form, we have

∂^αw(¯z,t_q₊₁)

∂t^α ≈

⎧⎨

⎩

a0(w^q⁺¹−w^q)+a0

q s=1

bs(w^q⁺¹⁻^s−w^q⁻^s), q ≥1,

a0(w¹−w⁰), q =0,

(8)

wherew=w1orw2in the case of (1).

The Caputo fractional derivative forα ∈(1,2)is

∂^αw(¯z,t)

∂t^α =

⎧⎪

⎪⎨

⎪⎪

⎩ 1 (2−α)

_t

0

∂²w(¯z, ϑ)

∂ϑ² (t−ϑ)¹^−αdϑ,1< α <2,

∂w²(¯z,t)

∂t² , α=2.

Like the previous case

∂^αw(¯z,tq+1)

∂t^α

= 1 (2−α)

t_q₊₁ 0

∂²w(¯z, ϑ)

∂ϑ² (tq+1−ϑ)¹^−αdϑ, which leads to

∂^αw(¯z,tq+1)

∂t^α ≈ 1

(2−α) q s=0

_(s+1)τ

sτ

∂²w(¯z, ϑs)

∂ϑ²

×(ts+1−ϑ)^1−αdϑ.

The term∂²w(¯z, ϑs)/∂ϑ²can be approximated as

∂²w(¯z, ϑs)

∂ϑ²

= w(¯z, ϑs+1)−2w(¯z, ϑs)+w(¯z, ϑs−1)

ϑ² +O(τ²).

The fractional derivative in simplified form is

∂^αw(¯z,tq+1)

∂t^α ≈a0(w^q+1−2w^q +w^q−1) +a0

q s=1

bs(w^q⁺¹⁻^s−2w^q⁻^s+w^q⁻^s⁻¹), q ≥0. Leta₀= [τ^−α/(3−α)]andb_s =(s+1)²^−α−s²^−α, s =0,1, . . . ,q, in more compact form, we have

∂^αw(¯z,t_q₊₁)

∂t^α

≈

⎧⎪

⎪⎪

⎪⎨

⎪⎪

⎩

a0(w^q⁺¹−2w^q+w^q⁻¹)+a0 q−1 s=1

bs(w^q⁺¹⁺^s

−2w^q⁻^s+w^q⁺¹⁻^s) +2a0bq

w¹−w⁰−τ∂w(¯z,0)

∂t

, q ≥1

2a0

w¹−w⁰−τ∂w(¯z,0)

∂t

, q =0.

3. Results and discussions

The LMCM is validated for finding numerical results of 1D time-fractional coupled KdV and 2D time-fractional Klein–Gordon model equations. The numerical and graphical results of the suggested LMCM are accurate and significant in both regular and irregular domains.

In all numerical simulations, multiquadric radial basis function with shape parameter c = 0.001 and c = 1 is used in coupled KdV and Klein–Gordon model equations respectively with time step sizeτ =0.0001 unless mentioned explicitly. Furthermore, MATLAB R2015a (8.5.0.197613) software is utilised for all computations.

To test the accuracy, we have utilised maximum error norm (L_∞) and root mean square error (RMS) norm which are defined as

L_∞=max|wex−w|,

(5)

Table 1. Comparison of RMS for Test Problem1.

α=β=0.3 α=β =0.5 α=β =0.95

t [13] LMCM [13] LMCM [13] LMCM

w1 0.1 2.190e−05 1.2763e−07 9.706e−05 1.0434e−07 1.919e−04 1.9188e−06

0.4 8.297e−04 1.1428e−05 2.085e−03 7.8759e−06 3.140e−03 1.0545e−05

0.7 7.898e−04 4.0790e−05 4.198e−03 3.7380e−05 8.561e−03 3.8685e−05

1.0 2.402e−03 7.0326e−05 2.785e−03 6.9251e−05 1.388e−02 7.4748e−05

w2 0.1 3.464e−05 6.9304e−08 9.826e−05 7.6654e−08 1.929e−04 1.9137e−06

0.4 3.584e−04 5.1676e−06 2.257e−03 3.7734e−06 3.432e−03 9.2408e−06

0.7 4.052e−04 1.5837e−05 3.790e−03 1.5247e−05 9.101e−03 2.5418e−05

1.0 2.102e−03 2.4720e−05 2.141e−03 2.4929e−05 1.123e−02 3.4857e−05

RMS= 1

N N i=1

(wex−w)²,

wherewexandwdenote exact and numerical solutions respectively. In the case of (1),w=w1orw2.

Test Problem 1. First, consider the time-fractional coupled KdV equation

∂^αw1(x,t)

∂t^α = −∂³w1

∂x³ −6w1∂w1

∂x +3w2∂w2

∂x + f1,

∂^βw2(x,t)

∂t^β = −∂³w2

∂x³ −3w1∂w2

∂x + f2,

x ∈(0,1), t ≥0, αandβ ∈(0,1), (9) where

f1(x,t)=3xt⁴+ 2xt²^−α (3−α), f₂(x,t)=3xt⁴+ 2xt²^−β

(3−β),

and the exact solution isw1(x,t)=w2(x,t)=xt². The numerical results of the suggested LMCM are reported in table 1 corresponding to Test Problem 1.

These results are calculated using nodal points N = 10. The tabulated results revealed the better accuracy of the LMCM compared to the method given in [13].

Numerical solutions and absolute errors ofw1(x,t)and w2(x,t)forα = β = 0.5 andt = 1 are visualised in figures1and2respectively, which show better accuracy of the LMCM.

Test Problem 2. Consider the following time-fractional coupled KdV equation:

∂^αw1(x,t)

∂t^α = −∂³w1

∂x³ −6w1∂w1

∂x +3w2∂w2

∂x + f1,

∂^βw2(x,t)

∂t^β = −∂³w2

∂x³ −3w1∂w2

∂x + f₂,

x ∈(0,1), t≥0, αandβ ∈(0,1), (10) where

f1(x,t)=3xt⁵+ x(⁷₂) (⁷₂ −α)t⁵²^−α, f₂(x,t)=3xt⁵+ x(⁷₂)

(⁷₂ −β)t⁵²^−β,

and the exact solution isw1(x,t)=w2(x,t)=x√ t⁵. Table2lists the numerical results for Test Problem2 for various values of α, β and t when N = 10. It is observed that the results of LMCM are better than the results from the method given in [13]. Absolute errors ofw1(x,t)andw2(x,t)are appeared in figures3and4 respectively.

Test Problem 3. Consider the homogeneous time-fracti onal coupled KdV equation

∂^αw1(x,t)

∂t^α = p

∂³w1

∂x³ +6w1∂w1

∂x

+2qw2∂w2

∂x ,

∂^βw2(x,t)

∂t^β = −∂³w2

∂x³ −3w1∂w2

∂x ,

x ∈(−50,50), t ≥0, αandβ ∈(0,1), (11) and the exact solution is given as

w1(x,t)= −c₁²(1+p)

3+6p + 4c²₁e^c¹⁽^x⁺^d¹^t⁾ (1+e^c¹⁽^x⁺^d¹^t⁾)², w2(x,t)= de^c¹⁽^x⁺^d¹^t⁾

(1+e^c¹⁽^x⁺^d¹^t⁾)², where

d = −c₁²

−24p

q , d1 = −pc1

1+2p.

For the sake of comparison, the parameters are considered as p = −1.5,q = 0.1, c1 = 0.1 and α = β =1.

(6)

0 0.2 0.4 0.6 0.8 1 0

0.5 10

0.2 0.4 0.6 0.8 1

t x w1(x,t)

0 0.2 0.4 0.6 0.8 1

0 0.5

10 1 2 3

x 10⁻⁴

t x Absolute error w1

Figure 1. Numerical solution (left) and the corresponding absolute error (right) ofw1(x,t)whenα=β =0.5 andN=10 for Test Problem1.

0 0.2 0.4 0.6 0.8 1

0 0.5

10 0.2 0.4 0.6 0.8 1

t x w2(x,t)

0 0.2 0.4 0.6 0.8 1

0 0.5

10 2 4 6 8

x 10⁻⁵

Figure 2. Numerical solution (left) and the corresponding absolute error (right) ofw2(x,t)whenα=β =0.5 andN=10 for Test Problem1.

Table 2. Comparison of RMS for Test Problem2.

α=β=0.3 α=β =0.5 α=β =0.95

t [13] LMCM [13] LMCM [13] LMCM

w1 0.1 1.154e−05 1.8235e−08 2.796e−05 2.8350e−08 5.991e−05 7.6503e−07

0.5 8.815e−04 1.3565e−05 2.262e−03 9.4946e−06 3.518e−03 1.2065e−05

1.0 3.559e−03 8.6262e−05 2.145e−03 8.4211e−05 1.433e−02 8.9128e−05

w2 0.1 1.303e−05 1.1556e−08 2.805e−05 2.5296e−08 6.004e−05 7.6452e−07

0.5 3.733e−04 6.0734e−06 2.439e−03 4.5114e−06 3.857e−03 1.0461e−05

1.0 2.015e−03 3.0225e−05 1.671e−03 3.0729e−05 1.126e−02 4.4670e−05

0 0.2 0.4 0.6 0.8 1

0 0.5

10 1 2 3

x 10⁻⁴

0 0.2 0.4 0.6 0.8 1

0 0.5

10 0.2 0.4 0.6 0.8 1

x 10⁻⁴

Figure 3. Absolute error whenα=β =0.3 andN =10 for Test Problem2.

(7)

0 0.2 0.4 0.6 0.8 1 0

0.5 10

1 2 3

x 10⁻⁴

0 0.2 0.4 0.6 0.8 1

0 0.5

10 0.5 1 1.5

x 10⁻⁴

Figure 4. Absolute error whenα=β=0.95 andN =10 for Test Problem2.

Table 3. Comparison of RMS whenα=β =1 for Test Problem3.

N =5 N =10 N =20

t LMCM [13] [45] LMCM [13] [45] LMCM [13] [45]

w1 0.1 6.20e−07 8.96e−07 – 4.05e−07 2.19e−06 – 2.97e−07 2.24e−06 –

0.5 3.10e−06 4.48e−06 6.22e−05 2.02e−06 1.09e−05 6.23e−05 1.49e−06 1.10e−05 6.01e−05 1.0 6.22e−06 8.98e−06 1.25e−04 4.06e−06 2.20e−05 1.25e−04 2.98e−06 2.22e−05 1.20e−04 1.5 9.34e−06 1.34e−05 1.87e−04 6.09e−06 3.32e−05 1.87e−04 4.47e−06 3.35e−05 1.81e−04 2.0 1.24e−05 1.79e−05 2.49e−04 8.12e−06 4.45e−05 2.49e−04 5.97e−06 4.49e−05 2.41e−04

w2 0.1 3.54e−06 1.11e−05 – 3.45e−06 1.12e−05 – 2.37e−06 1.20e−05 –

1.0 3.54e−05 1.12e−04 – 3.46e−05 1.12e−04 – 2.38e−05 1.20e−04 –

2.0 7.11e−05 2.24e−04 – 6.92e−05 2.26e−04 – 4.76e−05 2.40e−04 –

−50

0

50 0

0.5 1.5 1

2 0 2 4 6 8 10

x 10⁻³

t x w1(x,t)

−50

0

50 0

0.5 1.5 1

20 1 2 3 4

x 10⁻⁵

Figure 5. Approximate solution (left) and the corresponding absolute error (right) ofw1(x,t)whenα=β =0.5 andN=20 for Test Problem3.

To show the applicability of the LMCM in comparison with the methods given in [13,45], RMS are reported in table3for different values oft andN. These results ensure that the LMCM is more accurate and reliable.

Furthermore, the accuracy increases with increasing nodal points N. Figures5and6show numerical solutions and absolute errors respectively, while absolute error when α =β = 0.3, N =10 andα =β =0.7, N =20 are shown in figures7and8.

Test Problem 4. The two-dimensional time-fractional Klein–Gordon equation (2) have the exact solution

w(x,y,t)

=4000e^x⁺^yx⁶(1−x)⁶y⁶(1−y)⁶t^α+³,0≤t ≤3, where

f =

2000(4+α) 3

e^x+yx⁶(1−x)⁶y⁶(1−y)⁶t³

(8)

−50

0

50 0

0.5 1 1.5

−0.052

−0.04

−0.03

−0.02

−0.01 0

t x w2(x,t)

−50

0

50 0

0.5 1 1.5 20 1 2 3

x 10⁻⁴

Figure 6. Approximate solution (left) and the corresponding absolute error (right) ofw2(x,t)whenα=β =0.5 andN=20 for Test Problem3.

−50

0

50 0

0.5 1 1.5 20 1 2 3 4

x 10⁻⁵

−50

0

50 0

0.5 1 1.5 20 1 2 3

x 10⁻⁴

Figure 7. Absolute error ofw1(x,t)(left) andw2(x,t)(right) whenα=0.3, β=0.5 andN =10 for Test Problem3.

−50

0

50 0

0.5 1 1.5 20 1 2 3 4

x 10⁻⁵

−50

0

50 0

0.5 1 1.5 20 1 2 3

x 10⁻⁴

Figure 8. Absolute error ofw1(x,t)(left) andw2(x,t)(right) whenα=0.5, β=0.7 andN =20 for Test Problem3.

−4000e^x⁺^y(30x⁴−240x⁵+757x⁶−1206x⁷ +1005x⁸−380x⁹+15x¹⁰+18x¹¹+x¹²)

×y⁶(1−y)⁶t^3+α

−4000e^x⁺^yx⁶(1−x)⁶(30y⁴−240y⁵ +757y⁶−1206y⁷+1005y⁸−380y⁹ +15y¹⁰+18y¹¹+y¹²)t^3+α

(9)

Table 4. Comparison ofL_∞and RMS for Test Problem4.

α=1.3 α=1.5 α=1.7

t L_∞ RMS L_∞ RMS L_∞ RMS

0.5 3.8520e−07 1.8343e−08 3.1614e−07 1.5054e−08 2.5005e−07 1.1907e−08

1 7.8163e−06 3.7220e−07 7.6986e−06 3.6660e−07 7.5282e−06 3.5848e−07

2 1.5535e−04 7.3979e−06 2.0402e−04 9.7154e−06 2.3507e−04 1.1194e−05

3 8.8948e−04 4.2356e−05 1.1073e−03 5.2732e−05 1.3782e−03 6.5632e−05

0 0.2 0.4 0.6 0.8 1

0 0.5

10 1 2 3 4 5 6

x 10⁻⁴

y x

w(x,y,t)

0 0.2 0.4 0.6 0.8 1

0 0.5

10 1 2 3 4 5 6

x 10⁻⁴

y x

w(x,y,t)

Figure 9. Exact solution (left) and approximate solution (right) whenα=1.7 andt=1 for Test Problem4.

0 0.2 0.4 0.6 0.8 1

0 0.5

10 0.01 0.02 0.03 0.04 0.05 0.06

y x

w(x,y,t)

0 0.2 0.4 0.6 0.8 1

0 0.5

10 0.01 0.02 0.03 0.04 0.05 0.06

y x

w(x,y,t)

Figure 10. Exact solution (left) and approximate solution (right) whenα=1.3 andt =3 for Test Problem4.

+4000³e^3x⁺^3yx¹⁸(1−x)¹⁸y¹⁸(1−y)¹⁸t³^α+⁹. Table4 shows the computed results of the LMCM for various values ofαat different time levels up to final time t =3 using spatial domain[0,1]andN =20 for Test Problem4, which are in excellent agreement with the exact solutions. The simulation results are also shown graphically in figures9and10which show the behaviour of approximate and exact solutions whereas figure 11 shows absolute error at α = 1.3 (left) and α = 1.7 (right). The suggested LMCM is also tested on non- uniform nodes including Halton nodes, random nodes

and Chebyshev nodes, and the results are visualised in figure12in the form of absolute error. It is revealed also that in these cases the recommended algorithm produced better accuracy.

Test Problem 5. Finally, the two-dimensional time- fractional Klein–Gordon equation (2) has exact solution w(x,y,t)=cos(πx)cos(πy)t^α+³, 0≤t ≤1,