Novel operational matrices-based method for solving fractional-order delay differential equations via shifted Gegenbauer polynomials

Applied Mathematics and Computation 372 (2020) 124985

Contents lists available at ScienceDirect

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

Novel operational matrices-based method for solving fractional-order delay differential equations via shifted Gegenbauer polynomials M. Usman a,b,c,∗, M. Hamid d, T. Zubair d, R.U. Haq e, W. Wang d, M.B. Liu a,b,c,∗ a

BIC-ESAT, College of Engineering, Peking University, Beijing 100871, PR China State Key Laboratory for Turbulence and Complex Systems, Department of Mechanics and Engineering Science, Peking University, Beijing 100871, PR China c Institute of Ocean Research, Peking University, Beijing 100871, PR China d School of Mathematical Sciences, Peking University, Beijing 100871, PR China e Department of Electrical Engineering, Bahria University Islamabad, 44000, Pakistan b

a r t i c l e

i n f o

Article history: Received 31 January 2019 Revised 4 October 2019 Accepted 15 December 2019

Keywords: Shifted Gegenbauer polynomials Delay differential equation Operational matrix Fractional calculus

a b s t r a c t Accurate solutions of nonlinear multi-dimensional delay problems of fractional-order arising in mathematical physics and engineering recently have been found to be a challenging task for the research community. This paper witnesses that an efficient fully spectral operational matrices-based scheme is developed and successfully applied for stable solutions of time-fractional delay differential equations (DDEs). Monomials are introduced in order to proposed the novel operational matrices for fractional-order integration Iν and derivative Dν by means of shifted Gegenbauer polynomials. Some ordinary and partial delay differential equations of fractional-order are considered to show reliability, efficiency and appropriateness of the proposed method. In order to approximate the delay term in DDEs a novel delay operational matrix ab is introduced with the help of shifted Gegenbauer polynomials. The proposed algorithm transform the problem understudy into a system of algebraic equations which are easier to tackle. Analytical solutions of the mentioned problem are effectively obtained, and an inclusive comparative study is reported which reveals that the proposed computational scheme is effective, accurate and well-matched to investigate the solutions of aforementioned problems. Error bound analysis is enclosed in our investigation to reveal the consistency and support the mathematical formulation of the algorithm. This proposed scheme can be extended to explore the solution of more dervisfy problem of physical nature in complex geometry. © 2019 Elsevier Inc. All rights reserved.

1. Introduction Fractional calculus (FC) and differential equations with fractional-order derivative have now established as a renowned and another revolutionary division of fractional calculus because it has extensive range of applications in engineering, stochastic, biological and industrial mechanisms [1,2]. This antique mathematical theme is approved from a conversant dialogue among two recognized mathematicians (Hospital and Leibniz). Recently, it is also utilized to simulate the complex ∗

Corresponding author at: BIC-ESAT, College of Engineering, Peking University, Beijing 100871, PR China. E-mail addresses: [email protected] (M. Usman), [email protected] (M.B. Liu).

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

2

M. Usman, M. Hamid and T. Zubair et al. / Applied Mathematics and Computation 372 (2020) 124985

and multifaceted systems with the properties of viscoelastic materials. Several experts have offered comprehensive definition to calculate the fractional-order derivative [3,4]. From the comparison between integer and fractional-order differential equations (FDEs), it can be established that FDEs have huge range of benefits in other regions of disciplines like neurology, economics, signal dispensation, filtering, electrical circuits, electro-analytical chemistry, solid, and statistical mechanics, which can be seen in refs. [5–11]. Fractional differential equations can also be observed in physical-chemistry because of local chemical symmetry by mass transmission which is appropriate for solubility concentration from contributor stage into acceptor stage. Above-mentioned prospective applications in the field of science and engineering motivated the community of scientist to utilize this concept in the modeling of complex marvels of nature. In the last ten years, number of problems have been simulated for different sorts of complex physical marvels. After development of these models, next and most significant task was to evaluate the appropriate solution of formulated problems. Although it was most difficult task for FDEs but numerous analytical and numerical methods have widely been implemented to appraise the solution of existing problem [12–16]. Ul Hassan [17] established a procedure grounded on soliton concept to examine the result of equation (biological population) of fractional-order. Shah et al. [18] inspected the heat transmission for fluid of second grade over a wavering perpendicular sheet. Recently, a detail study on fractional and variable order DDEs is clarified by Aguilar et al. [19]. In this study, they appraised the solution of the mentioned problem with the help of artificial neural networks. The generalized form of delay differential [20,21] equation with fractional-order (FDDE) derivative is given as, C ν 0 Dt y

(t ) = F (y(t )) + β1 y(at − b) + β2 y(t ) + G(t ),

(1)

In Eq. (1), y(t) is dependent variable, F is the nonlinear operator and ν is fractional-order derivative with 0 < ν ≤ 1. Furthermore, β i , i = 1, 2, a = 0, b are arbitrary constants which is different for different kinds of delay differential problems. y(at − b) and G(t) are delay function and nonhomogeneous function of the expression (1) respectively. The suitable initial condition is given as,

y(t ) = φ (t ), t ≥ 0.

(2)

Solutions of delay differential equations are challenging task among the research community due to its huge involvement in all field of science and technology including chemistry, finance and physics [22]. Daftardar-Gejji and Jafari developed [23] proposed an efficient new iterative scheme to search numerical solution of the problem presented in (1)-(2). A new method based on predictor-corrector is reported by Daftardar-Gejji et al. [24] to explore the numerical solution of fractionalorder delay differential equation. In this paper, they shows that the proposed method is more accurate and time efficient as compare to other methods. Recently, it is observed that the time-fractional delay partial differential equations arising in the modeling of dynamics of populations structured w.r.t. size of cell, specimen age, maturation level etc. [25,26]. The general form of time-fractional delay partial differential equations defined in 0 ≤ x ≤ L, 0 ≤ t ≤ T and 0 < ν ≤ 1 is given as, C ν 0 Dt u

(x, t ) + μ1C0 Dtν u(x, a1t + b1 ) ∂ 2u ∂2 = μ2 2 + μ3 2 u(x, a2 t + b2 ) + μ4 u(x, a3 t + b3 ) + μ5 u(a4 x + b4 , t ) + μ6 u(a5 x + b5 , a6 t + b6 ) ∂x ∂x + μ7 F (u(x, t ) ) + f (x, t ).

(3)

Conditions subject to above problems are given as,

u(x, 0 ) = u ( 0, t ) =

α (x ), β1 (t ), u(L, t ) = β2 (t ).

(4)

where μi , aj , bj for i = 1, 2, …, 7, j = 1, 2, …, 6 are constants, F is nonlinear operator and f(x, t) is source term. Solodushkina et al. [27] are suggested an efficient scheme based on finite difference method to investigate the stable solution of first-order delay partial differential equation which used frequently in modeling the structured cell populations dynamics when maturity and age level are considered. Neutral delay parabolic differential equations is discussed by Zhang and Zhang [28] via a novel fourth-order stable implicit compact difference algorithm. In order to improve the accuracy of temporal variable different type of Richardson extrapolation techniques are utilized in this paper. Niu et al. [29] proposed a novel collocation reproducing kernel method to investigate the solution of delay heat conduction model. They shows the efficiency and reliability by considering two test problems. The procedure of orthogonal basis polynomials such as Haar, Chebyshev, Laguerre, Gegenbauer, Legendre, Laurent, and few others have opened new inquires in numerous areas of science [30–36]. Heydari et al. [37] presented the numerical results of fractional-order partial differential equations. They implemented the polynomial based scheme with Dirichlet boundary conditions. A modified version of wavelet process grounded with Chebyshev polynomials for fractional-order delay equation which can be perceived in refs. [38]. Recently, Gegenbauer polynomial based techniques offered by the Usman et al. [39,40] have provided better solutions because it improved the accuracy and reduced the computational work at a concrete level. With help of these techniques, they clarified the solution of two difficult problems i.e. one is fractional differential equations with variable order [39] and other is fluid problem [40] which is highly nonlinear coupled differential equations. Some other methods can be found in [41–44]. The concepts of origin shifted are further utilized for Gegenbauer polynomials and formulated shifted Gegenbauer polynomials. By means of shifted Gegenbauer polynomials novel operational matrices for fractional order integration Iν and

M. Usman, M. Hamid and T. Zubair et al. / Applied Mathematics and Computation 372 (2020) 124985

3

derivative Dν developed with the help of monomials. In order to make the method more accurate and efficient operational matrix is proposed ab to approximate the delay terms present in the DDEs. With the help of these operational matrices, a scheme is designed for the fractional-order delay differential equation which is actually main finding of this article. Development and implementations of this scheme in this context cannot be seen in literature. With the help designed scheme, delay differential equation is further convert into system of algebraic equation by means of proposed methodology. To proof the accuracy of this method comparison among evaluated solution obtained from the suggested method, exact and published results of the problem are also offered. It is detected that it is an accurate, effective and compatible methodology to treat the linear and nonlinear delay differential equation of fractional-order and it can be further be extended for the major finding of nonlinear problems of fractional order. This paper is divided into different sections as follows: Some basic definitions, important properties of fraction calculus and shifted Gegenbauer polynomials can be seen in the beginning part of the article (Section 2). In the next section, operational matrices for derivative, integration and delay terms are deliberated with the support of some important theorems and lemmas (Section 3). In section 4, error bounds is discussed. Five steps of proposed methodology can be seen Section 5. In order to proof the compatible and accuracy of the results of the suggested scheme, test problems along with error table and graphs are explained in Section 6. Conclusion and references are given in last section of the article. 2. Preliminaries and basic definitions 2.1. Fractional calculus The basic definitions regarding Riemann-Liouville and the Caputo’s fractional derivative are explained in this section. Definition 1. Fractional derivative of order δ in Riemann-Liouville sense [45] is given as follow: RL δ 0 Dt

1 dγ (γ − δ ) dt γ

f (t ) =



t

1

(t − s )δ−γ +1

a

f (s )ds, for γ − 1 ≤ δ < γ ∈ Z + .

(5)

Definition 2. The fractional derivative of n − 1 < δ < n in Caputo’s sense [45] is given as bellow: C δ 0 Dt

f (t ) =

1 (n − δ )



t

1

0+

(t − s )δ

f (n ) (s )ds, n ∈ N.

(6)

The operator C0 Dtδ satisfy the following properties:

1. C0 Dtδ (λ f (t ) + γ g(t )) = λC0 Dtδ f (t ) + γ C0 Dtδ g(t ),

⎧ ⎨ (β + 1 ) t β −δ , C δ β 2. 0 Dt t = (β − δ + 1 ) ⎩ 0,

otherwise,

(7)

δ ∈ N0 , β < δ.

3. C0 Dtδ λ = 0, where γ and λ are the constants. Definition 3. The fractional integral of order δ > 0 in Riemann-Liouville sense [45] is given as bellow: RL δ 0 Jt

f (t ) =

1

(δ )



t

1

0+

(t − s )1−δ

The above integral operator

RL J δ 0 t

f (s )ds.

satisfy the following properties.

α +β f (t ), (t )]=RL 0 Jt ( β + 1 ) β +α α β 2. RL t . 0 Jt t = (β + α + 1 )

1.

(8)

RL α RL β 0 Jt 0 Jt [ f

(9)

2.2. Shifted Gegenbauer polynomials and their properties The well-known Gegenbauer polynomials Cqν (t ) which are defined on the interval [ − 1, 1] can be computed by means of the following explicit formula:

Cqν (t ) =

q  2

k=0

(−1 )k (q − k + ν ) (2t )q−2k . k!(ν )(q − 2k )!

In order to use these polynomials on the interval [0, 1], we need to defined the shifted Gegenbauer polynomials [46] by introducing the new variable t = 2t − 1. The qth-order shifted Gegenbauer polynomials can be found with the help of an

4

M. Usman, M. Hamid and T. Zubair et al. / Applied Mathematics and Computation 372 (2020) 124985

explicit formula given as:

(λ + 1/2)  (−1 )q−k (q + k + 2λ) k t . (2 λ ) (k + λ + 1/2)(q − k )!k! q

Gλq (t ) =

(10)

k=0

The shifted Gegenbauer polynomials Gλq (t ) given in Eq. (10) are orthogonal [46] with respect to L2 -space on the interval [0, 1].



1

0



Gλp (t )Gλq (t )ϑ λ (t )dt =

0, for p = q, hλq , for p = q,

where ϑλ (t) and hλq are denotes the weight function and normalizing factor respectively and given as.

 λ− 1 ϑ λ (t ) = t − t 2 2 , hλq =

  2  λ + 12 (q + 2 λ) ((2 λ ) )2 (2q + 2 λ)q!

.

Connection of the shifted Gegenbauer polynomials (8) [46] with other polynomials are given as:

Lq (t ) = G1q /2 (t ), Tq (t ) = G0q (t ), Uq∗ (t ) = G1q (t ). In above expression Tq (t ), Uq∗ (t ), and Lq (t) are represents the first kind shifted Chebyshev, second kind shifted Chebyshev and shifted Legendre polynomials respectively. 2.3. Function approximation Let f(t) be a function from L2 (R )-space can be expanded with the help of truncated shifted Gegenbauer polynomials

Gλq (t ) as.

f (t ) =

M−1 

aq+1 Gλq (t ) = C T (t ),

(11)

q=0

where aq+1 =

1 0

f (t )Gλq (t )ϑ λ (t )dt. Moreover, C and (t) are the matrices of order M × 1 given as.



C = [a1 , a2 , a3 , . . . , aM ]T , (t ) = Gλ0 , Gλ1 , Gλ2 , . . . , GλM−1

T

.

(12)

Similarly, a function of two variables f(x, t) can be expanded with the help of truncated shifted Gegenbauer polynomials Gλq (t ) as:

f (x, t ) =

M−1  M−1 

a p+1,q+1 Gλp (x )Gλq (t ) = T (x )C (t ).

p=0 q=0

where (x) and (t) are the vectors of order M × 1 whereas C is a square matrix given below of order M × M as,

⎡

a1,1 ⎢ a2,1 C=⎢ . ⎣ .. aM,1

a1,2 a2,2 .. . aM,2

··· ··· .. . ···

⎤

a1,M  a2,M ⎥ (x ) = Gλ (x ), Gλ (x ), Gλ (x ), . . . , Gλ (x ) T , 0 1 2 M−1 ⎥ ,  .. ⎦ T (t ) = Gλ0 (t ), Gλ1 (t ), Gλ2 (t ), . . . , GλM−1 (t ) . . aM,M

(13)

3. Operational matrices of derivatives In the section, the operational matrices of integration and derivative of fractional-order is developed. Also operational matrix of delay term is constructed in this section. In order to construct these matrices first we defined the following family of monomials defined on [0, 1] as.

ωq = t q ,

(14)

here q = 0, 1, 2, …, M − 1. The qth-monomial functions can be stated as a vector form:

 = [ω0 , ω1 , ω2 , . . . , ωM−1 ].

(15)

Theorem 1. Consider (t) be a shifted Gegenbauer polynomials vector present in Eq. (12) then the delay expression (at + b) must satisfied the following relation.

(at + b) = ab (t ),

(16)

M. Usman, M. Hamid and T. Zubair et al. / Applied Mathematics and Computation 372 (2020) 124985

5

where ab is an operational matrix of order M × M which is stated as.

⎡

ζ00 ⎢ ζ01 ⎢ ζ2 ab = ⎢ ⎢ .0 ⎣ .. ζ0M−1

ζ10 ζ11 ζ12

ζ20 ζ21 ζ22

ζ1M−1

ζ2M−1

.. .

··· ··· ··· .. . ···

.. .

⎤ 0 ζM−1 1 ζM−1 ⎥ ⎥ 2 ζM−1 ⎥. .. ⎥ . ⎦ M−1 ζM−1

Then the following relation must satisfied:

⎛ ⎜

ζ ji =

 

⎛

(λ + 1/2 )(l + s + λ + 1/2 )  k−l ⎟⎟ l ⎟⎟ b ⎟⎟. ⎟⎟ (l + 2λ + s + 1 ) a ⎠⎠

⎜

j ⎜ i ⎜ k i−k j−s k   ⎜ ⎜ (2 j + 2λ ) j!(−1 ) (i + k + 2λ )(−1 ) ( j + s + 2λ )a ⎜ ⎜ ( j + 2λ )(k + λ + 1/2 )(i − k )!k!(s + λ + 1/2 )( j − s )!s! k=0

⎞⎞

k

⎝ s=0 ⎝

l=0

Proof. By means of function approximation, we obtain the subsequent relation as:

ζ ji = =



1

ω (t )Gλi (t ) Gλj (t )dt,

q 0 j i  k=0 s=0

2( j + λ ) j!(−1 )i−k (i + k + 2λ ) (−1 ) j−s (s + k + 2λ ) ( j + 2λ )(k + λ + 1/2 )(i − k )!k! (s + λ + 1/2 )( j − s )!s!



1 0

k



t s (at + b) t − t 2

λ−1/2

dt.

Integrating the integral term by parts and after a lot of simplification, we acquire the following required expression:

ζ ji =

  k−l

j i  

k 2( j + λ ) j!(−1 )i−k (i + k + 2λ ) (−1 ) j−s (s + j + 2λ )ak  k ( j + 2λ )(k + λ + 1/2 )(i − k )!k! (s + λ + 1/2 )( j − s )!s! l l=0

k=0 s=0

b a

(λ + 1/2 )(l + s + λ + 1/2 ) . (l + 2λ + s + 1 )

Remark 1. When b = 0 then ζ ji must satisfied the following relation:

ζ ji =

j i   k=0 s=0

2( j + λ ) j!(−1 )i−k (i + k + 2λ ) (−1 ) j−s (s + j + 2λ )ak ak (ν + 1/2 )(s + k + ν + 1/2 ) . ( j + 2λ )(k + λ + 1/2 )(i − k )!k! (s + λ + 1/2 )( j − s )!s! (k + 2ν + s + 1 )

Theorem 2. Consider (t) be a vector of monomials presented in Eq. (15) and

(t ) =  (t ),

(17)

where  is a matrix having order M × M which is stated as:

⎡

ξ00 ⎢ ξ01 ⎢ ξ2 =⎢ ⎢ .0 ⎣ .. ξ0M−1

ξ10 ξ11 ξ12

ξ20 ξ21 ξ22

ξ1M−1

ξ2M−1

.. .

.. .

··· ··· ··· .. . ···

⎤ 0 ξM−1 1 ξM−1 ⎥ ⎥ 2 ξM−1 ⎥. .. ⎥ . ⎦ M−1 ξM−1

Then the following relation must satisfied.

ξ = q j

(18)



j  2(−1 ) j−i (2ν )( j + ν ) j!(i + j + 2ν ) i + q + ν +

 ( j − i )!i! i + ν + 12 ( j + 2ν )(i + q + 2ν + 1)

i=0

1 2

.

Proof. Since that

(t ) =  (t ). By means of function approximation, we obtain

ξ = q j



1

0

ωq (t )Gλj (t ) ϑ λ (t )dt =

 j  λ− 12  λ + 12  (−1 ) j−i (i + j + 2λ) 1 i+q   t t − t2 dt. 1 h j (2λ )  i + λ + 2 ( j − i )! i ! 0 i=0

Integrating the integral term by parts and after simplification, we attained the following required expression:



ξ jq =

j  2(−1 ) j−i (2ν )( j + ν ) j!(i + j + 2ν ) i + q + ν +



i=0

( j − i )!i! i + ν +

1 2

( j + 2ν )(i + q + 2ν + 1)

1 2



.

(19)

6

M. Usman, M. Hamid and T. Zubair et al. / Applied Mathematics and Computation 372 (2020) 124985

For M = 4 the matrix  takes the following form.

⎡

⎢ =⎢ ⎣

1

0

1 2 3+2λ 8(1+λ ) 5+2λ 16(1+λ )

1 4λ 1 4λ 3(5+2λ ) 32λ(2+λ )

0 0

⎤

0 0 0

1 8λ(1+λ ) 3 16λ(1+λ )

⎥ ⎥. ⎦

3 32λ(1+λ )(2+λ )

Lemma 1. The differentiation of fractional order ν > 0; γ − 1 < ν < γ of Eq. (14) is defined as follows: C ν 0 Dt

ωq (t ) =

q!

(q − ν + 1 )

t q −ν , q = γ , γ + 1 , . . . , M − 1 .

Proof. The proof is very simple by using the Eq. (7). Lemma 2. The fractional differentiation of order ν of Eq. (15) in Caputo sense is defined as follow, where γ − 1 < ν < γ be a positive function: C ν 0 Dt [

(t )] = Pν (t ).

Here Pν is square having order M × M defined as:

⎡

0 ⎢ .. ⎢.

⎢ ⎢ ⎢0 ⎢ ⎢ ν −ν ⎢0 P =t ⎢ ⎢ .. ⎢. ⎢ ⎢ ⎢0 ⎣

··· .. .

0 .. .

···

γ! (γ −ν +1 )

···

0 .. . .. . .. .

··· ··· ···

0

···

0 .. .

···

0 .. .

0

···

0

(γ +1 )! (γ −ν +2 )

··· .. .

0 .. .

···

(M−2 )! (M−ν −1 )

0

···

0

(M−1 )! (M−ν )

.. . .. . .. .

0 .. . .. . .. . .. .

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

Proof. The proof is very simple by using the Lemma 1. Lemma 3. The integration of fractional order ν > 0; γ − 1 < ν < γ of expression (15) is defined as follow: RL ν 0 Jt

ωq (t ) =

q!

(q + ν + 1)

t q+ν , q = γ , γ + 1, . . . , M + 1.

Proof. The proof is very simple by using the Eq. (9). Lemma 4. The integration of fractional order ν > 0; γ − 1 < ν < γ of Eq. (15) in Riemann–Liouville sense is defined as follow: RL ν 0 Jt [

(t )] = Qν (t ).

Here Qν is square having order M × M defined as:

⎡

⎢ ⎢ ⎢ Qν = t ν ⎢ ⎢ ⎣

γ! (γ +ν +1)

0 0 .. . 0

0 (γ +1 )! (γ +ν +2)

0 .. . 0

0

0

0

0

0

0

0 .. .

0 .. .

(γ +2 )! (γ +ν +3)

.. . 0

0

⎤ ⎥ ⎥ ⎥ ⎥. ⎥ ⎦

(M−1 )! (M+ν )

Proof. The proof is very simple by means of Lemma 3. Theorem 3. The differentiation of fractional order ν > 0; γ − 1 < ν < γ of Eq. (12) in Caputo sense is defined as follow: C ν 0 Dt [

 (t )] = Dν (t ) = −1 Pν  (t ).

In above the matrices  and Pν are presented in above and Dν is known as the operational matrix for differentiation of fractional order ν for the Gegenbauer polynomials. Proof. By means of the relation (17), we obtain the subsequent form:

(t ) = −1 (t ). Differentiating both side of above expression w.r.t. t of order δ and using the Lemma 2, we achieved the following form: C ν 0 Dt [

 (t )] = −1C0 Dtν (t ) = −1 Pν  (t ),

M. Usman, M. Hamid and T. Zubair et al. / Applied Mathematics and Computation 372 (2020) 124985

7

which completes the proof. For λ = 1/2, M = 3 and 0 < ν < 1 the matrix Dν is given as below:

⎡

0

Dν = ⎣ (2−ν ) − (23−ν ) + (34−ν )

0

⎤

1

0 0

− (23−ν ) + (36−ν )

2 (3−ν )

1

(2−ν )

⎦.

Theorem 4. The integration of fractional order ν > 0; γ − 1 < ν < γ of Eq. (12) in Riemann-Liouville sense is defined as follow. RL ν 0 Jt [

 (t )] = Iν (t ) = −1 Qν  (t ),

here the matrices  and Qν are presented in above and Iν is known as the operational matrix for integration of fractional order ν for the Gegenbauer polynomials. Proof. Proof is similar as Theorem 3. For λ = 1, M = 3 and 0 < ν the matrix Iν is given as below:

⎡

1

0

(1+ν )

Iν = ⎣

2 2 − (1+ ν ) + (2+ν ) 8 10 − (2+ν ) + (3+ ν)

3 (1+ν )

⎤

0 0

1 (2+ν ) 4 8 − (2+ν ) + (3+ ν)

⎦.

2 (3+ν )

4. Error bound analysis In this section, error bound analysis of the shifted Gegenbauer polynomials method are discussed comprehensively. Theorem 5. Consider that u(t) ∈ Cm + 1 [0, 1] and Yq = span{Gλ0 (t ), Gλ1 (t ), Gλ2 (t ), . . . , GλN (t )}. Let Hm (),  = [0, 1] be Sobolve space that contain all functions u(t) on  = [0, 1] such that u(t) and its all weak derivative up to of order m are in L2 [0,1]. If the best approximation from ϒq of u(t) is uN (t) = CT (t) then the error bound is given as follows.

u(t ) − uN (t ) 2 ≤

β

(m + 1 )!



2α 2m+3 , 2m + 3

where β = max |u(m+1) (t )|, α = max(1 − t0 , t0 ). t∈[0,1 )

Proof. First, we consider the Taylor’s formula

u¯ (t ) = u(t0 ) + u p (t0 )(t − t0 ) + · · · + um p (t0 )

(t − t0 )m m!

.

Also, we know that

  (t − t0 )m+1 , ∃ξ ∈ (0, 1 ). |u(t ) − u¯ (t )| ≤ u(m+1) (ξ ) (m + 1 )!

(20)

Since, the best approximation of u(t) ∈ ϒ is CT  (t). Therefore, by means of the definition of best approximation, we have

   (m+1)  (t − t0 )m+1 2 u dt , (ξ ) (m + 1 )! 0 0  2  1  2 2m+3 2α β β ≤ . (t − t0 )2m+2 dt ≤ m + 1 ! m + 1 ! 2m + 3 ( ) ( ) 0

u(t ) − uN (t ) 22 ≤ u(t ) − u¯ (t ) 22 =



1

(u(t ) − u¯ (t ))2 dt =



1

After take the square root of the above expressions we achieved the required result.

5. Proposed methodology This section dedicated to the development of novel algorithms by utilizing the operational matrices for derivative Dν , integration Iν and delay term ab to investigate the solutions of delay differential equations of fractional-order. The proposed method have the following steps:

8

M. Usman, M. Hamid and T. Zubair et al. / Applied Mathematics and Computation 372 (2020) 124985

5.1. Working steps for ODEs Step 1. Let us consider the delay differential equation [21] given in (1)-(2) as, C ν 0 Dt u

(t ) = μ1 u(t ) + μ2 u(at + b) + μ3 F (u(t ) ) + f (t ), ν ≤ 1, 0 < t ≤ 1.

(21)

Associated with the following Dirichlet boundary conditions

u (0 ) = α .

(22)

Step 2. For the solution of problem (21), (22) via proposed method, we consider the following trial solution as:

u˜ (t ) = C T (t ).

(23)

Here, C = [a1 , a2 , a3 , …, aM ], is a matrix having order M × 1, which needs to be computed and is a vector present in Eq. (12). Approximate the each term in Eq. (21) with the help of operational matrices as, C ν 0 Dt u

(t ) = C T Dν (t ), u(at + b) = C T ab (t ).

(24)

Step 3. Problem presented in (21) takes the following matrix form after incorporating the above expressions:





C T Dν (t ) − μ1C T (t ) − μ2C T ab (t ) − μ3 F C T (t ) − f (t ) ∼ = (t ).

(25)

and the initial condition takes the following form:

 = C T (0 ) α .

(26)

Step 4. In order to find the matrix C, we must need M algebraic equations. Therefore, we collocate matrix Eq. (25) as:

(ti ) = 0, for i = 1, 2, 3, . . . , M − 1, i where ti = M−1 , From the above expression we achieve M − 1 equations, to calculate the values of M unknowns we need one more equation which can be obtain by using the initial condition as:

(0 ) = α . Step 5. Finally, the required unknown matrix C is achieved after solving the M system of algebraic equations with the aid of Maple 2015. Approximate solution of problem (21), (22) is obtained after inserting the matrix C into trial solution.

5.2. Working steps for PDEs Step 1. Consider the following delay partial differential equation given in (3), (4) as, C ν 0 Dt u

(x, t ) + μ1C0 Dtν u(x, a1t + b1 ) ∂ 2u ∂2 = μ2 2 + μ3 2 u(x, a2 t + b2 ) + μ4 u(x, a3 t + b3 ) + μ5 u(a4 x + b4 , t ) + μ6 u(a5 x + b5 , a6 t + b6 ) ∂x ∂x + μ7 F (u(x, t ) ) + f (x, t ), 0 < ν ≤ 1.

(27)

Conditions subject to above problems are given as,

u(x, 0 ) = u ( 0, t ) =

α (x ), β1 (t ), u(L, t ) = β2 (t ).

(28)

where 0 ≤ x ≤ L, 0 ≤ t ≤ T, μi , aj , bj for i = 1, 2, …, 7, j = 1, 2, …, 6 are constant, F is nonlinear operator and f(x, t) is source term. Step 2. In order to investigate the approximate solution via suggested method of above problem, we consider the following trial solution:

u˜ (x, t ) = T (x )C (t ).

(29)

Here, C = [ap,q ]M × M , is a square matrix of order M × M, which needs to be calculated and are the vectors present in above (13). Let us approximate each term in above delay partial differential equations by means of derivative, integration and delay operational matrices as,

M. Usman, M. Hamid and T. Zubair et al. / Applied Mathematics and Computation 372 (2020) 124985

9

Fig. 1. Flow chart of proposed method.

C ν 0 Dt u

(x, t ) = T (x )CDν (t ), C0 Dtν u(x, a1t + b1 ) = T (x )C ab11 Dν (t ),  T  T ∂ 2u ∂2 = T (x ) D2 C (t ), u(x, a2 t + b2 ) = T (x ) D2 C ab2 (t ) 2 2 2 ∂x ∂x  a 4 T a3 T T u(x, a3 t + b3 ) = (x )C b (t ), u(a4 x + b4 , t ) = (x ) b C (t ), 3 4  a 5 T a6 T u(a5 x + b5 , a6 t + b6 ) = (x ) b C b (t ) 5 6

(30)

Step 3. The following matrix equation attained after incorporating the above relations (30) and trial solution (29) into discussed problem (27):



(x ) T

 T

CDν + μ1 C ab1 Dν − μ2 D2 1

−μ4 C 

a3 b3



a4 T

− μ 5 b

4

 T

C − μ3 D 2



a5 T

C − μ 6 b

5

C

C ab2

2

a6 b6

 (t ) − μ7 F T (x )C (t ) − f (x, t ) ∼ = R(x, t ).

(31)

And associated conditions (28) takes the following form:

1 (x ) = T (x )C (0 ) − α (x ) 0. 2 (t ) = T (0 )C (t ) − β1 (t ) 0, 3 (t ) = T (L )C (t ) − β2 (t ) 0. Step 4. In order to find the square matrix C, we need



R xi =

i j L, t j = T M M



M2

(32)

algebraic equations. For this, we collocate matrix Eq. (31) as:

= 0, for i = 2, 3, 4, . . . M − 1, j = 2, 3, 4, . . . , M.

(33)

Above residual vector gives M − 3M + 2 algebraic equations. In order to construct rest 3M − 2 equations we must collocate initial and boundary conditions as,

  i 1 xi = L = 0, for i = 1, 2, 3, . . . , M,

(34)

    j j 2 t j = T = 3 t j = T = 0, for j = 2, 3, 4, . . . , M.

(35)

M

M

M

Step 5. System of algebraic equations attained from (33)–35) can be solved by suitable approach. After computing the matrix C, approximate solution of ((27), (28) obtained by using C in trial solution (29). 6. Test problems In this section some numerical problems are taken to validate the efficiency of the shifted Gegenbauer polynomials based method. For numeric computation we used Maple 2015. Absolute-error, L2 , L∞ , and RMS norms can be computed by means of the following relations.

|E(xi , t = 0.5)| = |u(xi , t = 0.5) − u˜ (xi , t = 0.5)|, L∞ = max (|E(xi , t = 0.5)| ), i   M M ! ! 2 1 L2 = |u(xi , t = 0.5) − u˜ (xi , t = 0.5)| , RMS = M+1 |u(xi , t = 0.5) − u˜ (xi , t = 0.5)|2 . i=1

i=1

10

M. Usman, M. Hamid and T. Zubair et al. / Applied Mathematics and Computation 372 (2020) 124985

Problem 1. Consider the linear delay differential equation of fractional order is given as [21]: c 0.9 0 Dη u

(η ) =

2 η 1. 1 η 0.1 − + u(η − 0.1 ) − u(η ) + 0.2η − 0.11, (2.1) (1.1)

(36)

along with the initial condition is given as follow [21]:

u ( η ) = 0,

η ≤ 0.

(37)

η2

Exact solution of above problem is u(η) = − η. Step 2. Consider the following approximate solution M = 3 to find the solution of the above delay differential equation:

u˜ (η ) =

2 

ak+1 Gλk (η ) = C T (η ),

(38)

k=0

where C = [a1 ,a2 ,a3 ]T , is a matrix need to be determine. Step 3. Inserting the trial solution (38), and delay relations (16) for a = 1 and b = −0.1 into the given delay differential equation, we achieve the following matrix equations as:

C T K0.9 (η ) =

2 η 1. 1 η 0.1 − + C T 1−0.1 (η ) − C T (η ) + 0.2η − 0.11, (2.1) (1.1)

(39)

and initial condition takes the following form:

C T (0 ) = 0. where

K0.9 ,C,

(40)

(η) and

1−0.1

the matrices are given as:

 9 0  10 9 0.9 − 10  π η × 20 K = 120 π csc 10 11

0 10

0 0

360 11

200 11

,

   1 τ1 25 1 1 2 ( 2η − 1 ) −10 C = τ2 , (η ) = , −0.1 = 25 τ3 4 16η2 − 16η + 3

0 25 −20

0 0 25

.

Step 4. In order to find C collocate the matrices Eqs. (39), (40) as follow to construct system of algebraic equations.

a1 − 2 a2 + 3 a3 = 0, 1 1100π 1 1100π









 

9 π 9 π 10 −80 0 02 10  10  9 a3 sin π10 + 50 02  10 sin 10 9 +220 0 02 10  10 a2 sin 10 − 176a3 π + 11π + 440a2 π 9







9

 

9 π π 9 1440 0 0 10  π a3 9sin 10 − 90 0 0 sin 10  10 +440 0 0 sin 10  10 a2 + 1584a3 π − 99π + 440a2 π

= 0,

 = 0.

Step 5. The matrix C is obtained as following after solving the above system of equation with the help of Maple 2015.

"

C= −

3 16

1 16

0

#T

.

Exact solution is achieved here as bellow after incorporating the value of C into the trial solution (38).

"

3 u˜ (η ) = − 16

0

1 16

#



1 2 ( 2η − 1 ) 16η2 − 16η + 3

= η2 − η.

The analytical solution is obtained within 6 seconds for M = 3 which is the beauty of this method whereas an approximate solution is attained in the ref. [21] via two approaches which took 88.296875 seconds to compuate the results. It is clearly shown through the problem understudy that the proposed method is much effective, less time consuming and accurate. Problem 2. Consider the following nonlinear delay differential equation [21]. c ν 0 Dη u

(η ) =

$ 2 1− ν [u(η )] 2 + u(η − τ ) − u(η ) + 2τ u(η ) − τ 2 , (3 − ν )

(41)

Initial condition associated with problem (41) is given as.

u ( η ) = 0,

η ≤ 0.

(42)

M. Usman, M. Hamid and T. Zubair et al. / Applied Mathematics and Computation 372 (2020) 124985

11

1

Exact solution Approximate solution

0.8

u(η)

0.6

0.4

0.2

0

0

0.2

0.4

η

0.6

0.8

1

Fig. 2. Comparison of exact and approximate solutions for M = 3. Table 1 Comparison of absolute error of the proposed method with [21] of Problem 2.

η

[21]

[21]

Proposed method

0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

0.078120 0.129928 0.190687 0.248601 0.307649 0.366427 0.425322 0.484208 0.543110 0.602019

0.078155 0.129978 0.190760 0.248694 0.307763 0.366563 0.425479 0.484387 0.543312 0.602243

3.00 1.00 1.00 0.00 1.00 0.00 1.00 1.00 2.00 1.00

× × × × × × × × × ×

10−11 10−10 10−10 10−00 10−10 10−00 10−09 10−09 10−09 10−09

Given delay problem (41)–(42) have an exact solution u(η) = η2 . Now, following the above methodology by considering M = 3. Here, we got three nonlinear system of equations and proceeding as before Maple 2015 is used for solving these equations for unknown vector C which is given bellow.



C = 0.31250 0 0 0 02

0.25

T

0.06249999993 .

Consequently, we got the following approximate solution as:

u(η ) = 0.9999999989η2 + 1.1 × 10−9 η. Fig. 2 demonstrate comparison of approximate solution of (41), (42) with the exact solution u(η) = η2 . It is observed that the proposed method very accurate and efficient to investigate the solution of the given nonlinear delay differential equation. Table 1 shows the comparison of the absolute error with the method given in ref. [21]. Again this table is witnesses that the proposed method is well-match to solve such type of physical problem. Moreover, exact solution of this problem can be attained for M > 3. Problem 3. Consider the system of linear delay differential equations as [20]:

%η& η + u(η ) − v(η ), ( η ) = e −η − e 2 + u 2 % & η η C ν η 2 − u (η ) − v (η ) 0 Dη v ( η ) = e + e − u C ν 0 Dη u

2

with the following initial conditions:

u ( 0 ) = v ( 0 ) = 1. The exact solutions of the above system when ν = 1 is,

u ( η ) = e η , v ( η ) = e −η . Proceeding as before for M. Table 2 shows the error analysis for different values of ν which shows that as ν → 1 our solution converges to the exact solution for ν = 1. Table 3 construct to show that as the order of approximation M increase our approximate solution converges uniformly to the analytical solution. Fig. 3 is plotted to exhibit the behavior of L2 , L∞ ,

12

M. Usman, M. Hamid and T. Zubair et al. / Applied Mathematics and Computation 372 (2020) 124985 Table 2 Comparison between the exact and approximate solution for different values of ν when M = 5.

ηi

u(η )

0.2 0.4 0.6 0.8 1.0

v(η )

ν = 0.95

ν = 0.99

ν = 1.0

ν = 0.95

ν = 0.99

ν = 1.0

2.34E−02 5.19E−02 8.83E−02 1.38E−01 2.05E−01

3.67E−03 8.59E−03 1.47E−02 2.28E−02 3.40E−02

8.67E−04 1.30E−03 1.98E−03 2.97E−03 4.28E−03

1.36E−02 2.38E−02 3.58E−02 5.37E−02 8.11E−02

2.26E−03 3.91E−03 5.72E−03 8.47E−03 1.28E−02

3.80E−04 6.33E−04 9.89E−04 1.48E−03 2.12E−03

Table 3 Comparison between the exact and approximate solution for different values of M when ν = 1.

ηi

u(η )

0.2 0.4 0.6 0.8 1.0

M = 10

M = 20

M = 30

M = 10

M = 20

M = 30

1.86E−10 2.76E−10 4.09E−10 6.04E−10 8.83E−10

1.50E−26 2.22E−26 3.30E−26 4.87E−26 7.11E−26

8.58E−43 1.36E−42 2.07E−42 3.10E−42 4.54E−42

1.38E−11 8.40E−11 1.68E−10 2.74E−10 4.13E−10

1.19E−27 6.82E−27 1.36E−26 2.22E−26 3.34E−26

3.39E−43 1.08E−43 5.95E−43 1.17E−42 1.90E−42

6

1

10

-4

10

-9

5

Times (s)

Nroms Error

10

v(η )

-14

10

-19

10

u: L∞ u: L2 u: RMS v: L∞ v: L2 v: RMS

-24

10

10-29 3

6

9

4 3 2 1

12

(a)

M

15

18

21

24

3

6

9

12

15

18

M (b)

21

24

27

30

Fig. 3. 2D behavior of (a) Norms against M and (b) computational time (s) against M.

RMS norms and computational time (s) against the convergence control parameter M. It is important to highlight that all error norms are stable and converges to 0 as convergence control parameter enhancing gradually. Also from this figure we can observe that highly accurate results can be attained in very short time like 10−30 accuracy can be achieved via proposed method no more than 5 seconds. Therefore, the proposed method very efficient tool for such problems. Problem 4. Now consider the heat conduction partial differential equation with delay against μi = 0, for i = 1, 3, 5, 6 and aj = bj = 0, for j = 1, 2, 4, 5, 6 as given bellow [29]: C ν 0 Dt u

(x, t ) = (sin (t ) + 2 )2

% & 1 ∂ 2u + cos (t )u x, t − + f (x, t ), 2 2 ∂x

(43)

initial and boundary conditions are given as,

u(x, 0 ) = x(x − 1 ), u(0, t ) = u(1, t ) = 0.

(44)

Analytic solution of above problem is u(x, t) = x(x − 1)exp (t). Source term f(x, t) can be computed correspondingly. In order to show the detail of solution procedure via proposed method, we select M = 3 and ν = 1 then problem (43) reduces

M. Usman, M. Hamid and T. Zubair et al. / Applied Mathematics and Computation 372 (2020) 124985

13

to the following matrix equation,

" #  T T (x ) CD1 − (sin t + 2 )2 D2 C − cos t C 1−0.5 (t ) − f (x, t ) ∼ = R(x, t ),

where the matrices are given as,



(x ) =

1 4x − 2 , (t ) = 16x2 − 16x + 3



0 2 D2 = 2 0 x 5

0 0 4

0 0 , 1−0.5 = 1





1 4t − 2 ,C= 16t 2 − 16t + 3

1 −2 4

0 1 −4



a1,1 a2,1 a3,1

a1,2 a2,2 a3,2

(45)



a1,3 0 1 a2,3 , D1 = 2 t a3,3 2

0 1 4

0 0 , 2

0 0 . 1

and associated conditions takes the following form:

1 (x ) = T (x )C (0 ) − x(x − 1 ) 0. 2 (t ) = T (0 )C (t ) 0, 3 (t ) = T (1 )C (t ) 0.

(46)

In order to compute unknown matrix C, we must have 9 algebraic equations. Therefore, collocate (45) i.e., R(xi = i/3, tj = j/3) = 0, for i = 2, j = 2, 3, we get

% &

980 16 8 32 2320 2 a1,3 + a2,2 + a2,3 − 160a3,1 − a3,2 + a3,3 − 128 sin a3,1 3 3 9 9 27 3 % & % & % & % & % & 256 7 2 640 2 2 4 2 2 − sin a3,2 + sin a3,3 − cos a1,1 + cos a1,2 − cos a1,3 3 3 9 3 3 3 3 9 3 %2& % & % & % & % & 14 20 8 2 2 5 2 2 − 2/3 cos a2,1 + cos a2,2 − cos a2,3 + cos a3,1 − cos a3,2 3 9 3 27 3 9 3 27 3 4 a1,2 +

+

% &

%

% &&2

35 2 2 2 cos a3,3 − 2e 3 cos 81 3 3

%

% &&2

%

2

+ 8e 3 sin

% &&

%2& 3

%

+ 32 cos

% &

% 2 &&2 3

a3,1

160 2 2 2 2 1 92 2 cos a3,3 − cos e6 + e 3 = 0, 9 3 9 3 9 2900 4400 8 32 92 4a1,2 + 16a1,3 + a2,2 + a2,3 − 160a3,1 − a3,2 − a3,3 + e − 128 sin (1 )a3,1 3 3 9 9 9 2 2 − 256 sin (1 )a3,2 − 384 sin (1 )a3,3 − cos (1 )a1,1 + cos (1 )a1,3 − cos (1 )a2,1 + cos (1 )a2,3 3 3 5 5 2 2 + cos (1 )a3,1 − cos (1 )a3,3 − 2e(cos (1 ) ) + 8e sin (1 ) + 32(cos (1 ) ) a3,1 9 9 2 2 + 64(cos (1 ) ) a3,2 + 96(cos (1 ) )2 a3,3 − cos (1 )e1/2 = 0, 9 +

64 2 cos 3 3

a3,2 −

From above matrix Eq. (45) we attained only 2 algebraic equations. Rest 7 algebraic equations can be obtained by collocating initial 1 (xi = i/3) = xi (xi − 1), for i = 1, 2, 3 as,

2 5 5 4 10 2 a2,1 − a3,1 − 2 a1,2 + a2,2 + a3,2 + 3 a1,3 − 2 a2,3 − a3,3 + = 0 , 3 9 3 9 3 9 5 4 5 2 10 2 a1,1 + a2,1 − a3,1 − 2 a1,2 − a2,2 + a3,2 + 3 a1,3 + 2 a2,3 − a3,3 + = 0 , 3 9 3 9 3 9 a1,1 + 2 a2,1 + 3 a3,1 − 2 a1,2 − 4 a2,2 − 6 a3,2 + 3 a1,3 + 6 a2,3 + 9 a3,3 = 0 , a1,1 −

and boundary conditions 2 (tj = j/3) = 3 (tj = j/3) = 0, for j = 2, 3 as,

4 5 5 2 10 a1,2 − a2,2 + 2 a3,2 − a1,3 + a2,3 − a3,3 = 0 , 3 3 9 9 3 4 5 10 5 2 a1,1 + 2 a2,1 + 3 a3,1 + a1,2 − a2,2 + 2 a3,2 − a1,3 − a2,3 − a3,3 = 0 , 3 3 9 9 3 a1,1 − 2 a2,1 + 3 a3,1 + 2 a1,2 − 4 a2,2 + 6 a3,2 + 3 a1,3 − 6 a2,3 + 9 a3,3 = 0 , a1,1 − 2 a2,1 + 3 a3,1 +

a1,1 + 2 a2,1 + 3 a3,1 + 2 a1,2 + 4 a2,2 + 6 a3,2 + 3 a1,3 + 6 a2,3 + 9 a3,3 = 0 . Step 5. After solving above system of 9 algebraic equations we get the following unknown matrix C as,



C=

−0.3170595647 −1.889934961 × 10−13 0.1056865216

−0.08060316892 0 0.02686772297

−0.01054892447 −3.40188293 × 10−13 0.003516308146

14

M. Usman, M. Hamid and T. Zubair et al. / Applied Mathematics and Computation 372 (2020) 124985 Table 4 Absolute error comparison of attained results with published results [29] against numerous values of M, x and t. N=M

x = 0.25

x = 0.5

[29] t = 0.25 4 8 16 32 t = 0.5 4 8 16 32

Present

x = 0.75

[29]

Present

[29]

Present

5.54 2.43 6.27 1.56

× × × ×

10−5 10−5 10−6 10−6

3.52 2.78 6.65 4.24

× × × ×

10−4 10−9 10−15 10−19

8.97 3.49 8.83 2.20

× × × ×

10−5 10−5 10−6 10−6

4.64 3.84 9.47 5.33

× × × ×

10−4 10−9 10−15 10−19

8.29 2.62 6.41 1.57

× × × ×

10−5 10−5 10−6 10-6

3.45 2.78 6.65 3.15

× × × ×

10−4 10−9 10−15 10−19

5.41 5.88 1.49 3.75

× × × ×

10−4 10−5 10−5 10−6

2.20 4.85 1.58 6.97

× × × ×

10−5 10−11 10−16 10−21

7.77 8.40 2.10 5.28

× × × ×

10−4 10−5 10−5 10−6

2.70 6.84 2.67 9.91

× × × ×

10−5 10−11 10−16 10−21

6.24 6.27 1.52 3.78

× × × ×

10−4 10−5 10−5 10−6

1.85 4.84 1.58 7.05

× × × ×

10−5 10−11 10−16 10−21

Table 5 Analysis of L2 , L∞ and RMS norms against ν when M = 10.

ν

L∞

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

2.95 3.04 3.18 3.38 3.64 3.95 4.11 3.41 4.26 1.47

L2 × × × × × × × × × ×

10−4 10−4 10−4 10−4 10−4 10−4 10−4 10−4 10−4 10−4

RMS

2.46 2.61 2.79 2.99 3.22 3.37 3.05 1.75 2.61 7.73

× × × × × × × × × ×

10−3 10−3 10−3 10−3 10−3 10−3 10−3 10−3 10−3 10−3

2.44 2.58 2.76 2.96 3.19 3.34 3.02 1.73 2.59 7.66

× × × × × × × × × ×

10−5 10−5 10−5 10−5 10−5 10−5 10−5 10−5 10−5 10−5

10-12

Nroms Error

Absolute Error

1.5E-13 10-13

1E-13

-14

10

5E-14

0

ν = 0.8 ν = 0.9 ν = 1.0 0

0.1

0.2

0.3

0.4

0.5

x

0.6

L∞ L2 RMS

10-15

0.7

0.8

0.9

1

0.1

0.2

(a)

0.3

0.4

0.5

ν

0.6

0.7

0.8

0.9

1

(b)

Fig. 4. Impact of fractional parameter ν on (a) absolute error and (b) L2 , L∞ and RMS norms.

Approximate solution of (43)-(44) obtained by using C in trial solution as,

u˜ (x, t ) = −x + x2 + 0.9001749x2t 2 − 0.9001749xt 2 + 0.81935939x2t − 0.81935939xt. In the similar way simulations are performed for various values of convergence control parameter M and fractional parameter ν . Table 4 is construct to show the accuracy of the proposed method. This table is witnesses that our attained solution is more accurate than the published results and converges uniformly to the exact solution as M enhancing. Fig. 4 are plotted to see the influence of fractional parameter ν on absolute error and L2 , L∞ , RMS norms. Approximate solutions converges to exact solution and absolute error decreases gradually as fractional parameter ν → 1. Table 4 and Fig. 4 evident that proposed method very accurate and stable against M and ν .

M. Usman, M. Hamid and T. Zubair et al. / Applied Mathematics and Computation 372 (2020) 124985

100

10

Absolute Error

Nroms Error

10

10

-1

10-3 -4

10

-5

10

-6

-4

10-5

-2

10

15

L∞ L2 RMS 6

8

10

12

M

10

-6

10

-7

10

-8

14

M = 10 M = 12 M = 14 0

0.5

1

1.5

x

(a)

2

(b)

Fig. 5. Effect of M when ν = 0.5 on (a) L2 , L∞ and RMS norms and (b) absolute error.

10-2 10-4

Nroms Error

Nroms Error

10

10

-7

-6

10-8

10

10-10 10-12 10-14

L∞ L2 RMS

-16

10

4

6

8

-8

10-9

10

M

12

14

16

18

0.5

L∞ L2 RMS 0.6

0.7

(a)

ν

0.8

0.9

1

(b)

Fig. 6. Behavior of L2 , L∞ and RMS norms against (a) M and (b) ν .

Problem 5. Consider the following first order partial differential equation with delay defined in  = [0, 2] × [0, 2π ] as [27]:

%x ∂u π& x + 2C0 Dtν u(x, t ) = u(x, t ) − e 2 u ,t − + f (x, t ), ∂x 2 2 initial and boundary conditions are given as,

u(x, 0 ) = sin t, t ∈ [0, 2π ], u(x, t ) = exp (x ) sin t, x ∈ [0, 2], t ∈ [−π /2, 0]. Source term f(x, t) can be computed correspondingly. Proceeding as before Table 5 is created to see the behavior of L2 , L∞ and RMS norms as varying ν when M = 10. It is observed that all norms decreases gradually as fractional parameter approaches to unity. Also 10−4 accuracy of L∞ norm is achieved by solving 10 × 10 algebraic equations whereas the method in the ref. [27] attain 10−2 accuracy of L∞ norm against h = 1/40 and ࢞ = π /200, which clearly shows that proposed schemes gives better results than found in the method given in ref. [27] against less computational cost. Fig. 5 illustrate that approximate solution and L2 , L∞ and RMS norms is significantly affected by M. In other words, accuracy enhance as enhancing the convergence control parameter M.

16

M. Usman, M. Hamid and T. Zubair et al. / Applied Mathematics and Computation 372 (2020) 124985

10

Absolute Error

10

0

-5

10-10 -15

10

M=5 M = 10 M = 15 M = 20

-20

10

0

0.2

0.4

x

0.6

0.8

1

(a)

(b)

Fig. 7. Comparison of absolute error attained via (a) proposed method and (b) published results [28].

Problem 6. Consider the following neutral delay parabolic partial differential equation defined in  = [0, 1] × [0, 2] as [28]: C ν 0 Dt u

(x, t ) +

& 1 ∂ 2u 1 ∂ 2 % & 1 1 ∂ν % u x, t − = + u x, t − + f (x, t ), ∂tν 10 2 ∂ x2 2 ∂ x2 10

Following are the initial and boundary conditions,

u(x, 0 ) = 0, u(0, t ) = u(1, t ) = t 2 . Select f(x, t) so that given problems satisfied the exact solutions u(x, t) = t2 cos π x. By means of proposed scheme Fig. 6 is plotted to analyze the behavior of all norms for numerous values of fractional parameter ν and convergence control parameter M. Again same behavior of ν and M is observed. All discussed norms decreases gradually and varying value of M. Proposed algorithm to investigate the solution of discussed problem is stable w.r.t. M and ν . Fig. 7 contain the comparison study of attained and published results for absolute error. Very interesting and significant results attained here, highly accurate results (10−15 accuracy) obtained against M = 20 whereas the method give in the ref. [28] cannot attain this level of accuracy after a huge computational work h = 1/80 and τ = 1/6400, which clearly demonstrate the effectiveness and appropriateness of solve such kind of problems. 7. Conclusion In the current work, we developed an algorithm based on operational matrices of derivatives and integrations. The matrices has been genertared via a new approach with the help of shifted Gegenbauer polynomials. The validatity of the method is checked through different fractional-order differential equations mainly delay equations. The proposed method successfully converts the delay problems into system of algebraic equations assisted by developed operational matrices. Obtained solutions of delay problems by means of suggested method are compared with existing results. The outcomes found via suggested method are endorsing the accuracy and effectiveness of the proposed technique. The error-bound is analyzed in our investigation to show the consistency and support the mathematical formulation of the algorithm. The norms L2 and L∞ as well as RMS and absolute error has been presented in tabular form while some graphical plots are being made to show the credibility of the scheme. It is noted that it is an accurate and efficient mathematical tool to tackle the nonlinear fractional order problems of complex nature and it can be further extended for the nonlinear dynamical problems arisning in applied sciences. Acknowledgment Dr. M. Usman is grateful for the support of Peking University through the Boya Post-doctoral Fellowship. References [1] F. Mainardi, Fractional Calculus and Waves in Linear Visco-Elasticity: An Introduction to Mathematical Model, Imperial college press, London, 2010. [2] M. Morgado, N. Frod, P Lima, Analysis and numerical methods for fractional differential equations with delay, J. Comput. Appl. Math. (2013) 159–168. [3] E. Oliveria, J. Machado, A review of definitions for fractional derivatives and integral, Math. Probl. Eng. 2014 (2014) 238459.

M. Usman, M. Hamid and T. Zubair et al. / Applied Mathematics and Computation 372 (2020) 124985

17

[4] K. Miller, B. Ross, An Introduction to the Fractional Calculus and Fractional Differential Equations, John Wiley and Sons, 1993. [5] D. Tien, Fractional stochastic differential equations wfith applications to finance, J. Math. Anal. Appl. 397 (1) (2013) 334–348. [6] R. Matusu, Fractional order calculus in control theory, in: Proceedings of the 13th WSEAS International Conference on Automatic Control, Modelling and Simulation, 2011. [7] E. Barkai, R. Metzler, J. Klafter, From continuous time random walks to the fractional Fokker-Planck equation, Phys. Rev. E 61 (1) (20 0 0) 132. [8] S. Yuste, L. Acedo, K. Lindenberg, Reaction front in an A+ B→ C reaction-subdiffusion process, Phys. Rev. E 69 (3) (2004) 036126. [9] D. Benson, S. Wheatcraft, M.M. Meerschaert, Application of a fractional advection-dispersion equation, Water Resour. Res. 36 (6) (20 0 0) 12–1403. [10] R. Gorenflo, F. Mainardi, E. Scalas, M. Raberto, Fractional calculus and continuous-time finance III: the diffusion limit, Math. Finance 20 01 (20 01) 171–180. [11] Y. Chen, Y. Wei, D. Liu, H. Yu, Numerical solution for a class of nonlinear variable order fractional differential equations with Legendre wavelets, Appl. Math. Lett. 46 (2015) 83–88. [12] O. Mohamad, A. Khlaif, Homotopy analysis method for solving delay differential equations of fractional order, Math. Theory Model. 4 (14) (2014) 48–56. [13] U. Saeed, Hermite wavelet method for fractional delay differential equations, J. Differ. Equ. (2014). [14] S. Davaeifar, J. Rashidinia, Solution of a system of delay differential equations of multi pantograph type, J. Taibah Univ. Sci. 11 (2017) 1141–1157. [15] A. Ardjouni, A. Djoudi, Existence of positive periodic solutions for a nonlinear neutral differential equation with variable delay, Appl. Math. E-Notes 12 (2012) 94–101. [16] K. Moaddy, A. Freihat, M. Smadi, E. Abuteen, I. Hashim, Numerical investigation for handling fractional-order Rabinovich–Fabrikant model using the multistep approach, Soft Comput. 22 (3) (2018) 773–782. [17] Q. Ul Hassan, S. Mohyud-Din, Investigating biological population model using exp-function method, Int. J. Biomath. 9 (02) (2016) 1650026. [18] N. Shah, I. Khan, Heat transfer analysis in a second grade fluid over and oscillating vertical plate using fractional Caputo–Fabrizio derivatives, Eur. Phys. J. C 76 (7) (2016) 362. [19] C. Aguilar, A. Escamilla, J. Aguilar, V. Martinez, H. Ugalde, New numerical approximation for solving fraction delay differential equation of veriable order using artificial neural networks, Eur. Phys. J. Plus 133 (2) (2018) 75. [20] S. Edeki, G. Akinlabi, Analytical solution of a time fractional system of proportional delay differential equations, in: Proceedings of the Second International Conference on Knowledge Engineering and Application, 2017. [21] V. Gejji, Y. Sukale, S. Bhalekar, Solving fractional delay differential equations: a novel approach, Int. J. Theory Appl. 18 (2015) 400–418. [22] W. Zhen, H. Xia, S. Guodong, Analysis of nonlinear dynamics and chaos in a fractional order financial system with time delay, Comput. Math. Appl. 62 (2011) 1531–1539. [23] V. Daftardar-Gejji, H. Jafari, An iterative method for solving non linear functional equations, J. Math. Anal. Appl. 316 (2006) 753–763. [24] V. Daftardar-Gejji, Y. Sukale, S. Bhalekar, A new predictor-corrector method for fractional differential equations, Appl. Math. Computations 244 (2014) 158–182. [25] M. Gyllenberg, J. Heijmans, An abstract delay differential equation modeling size dependent cell growth and division, SIAM J. Math. Anal. 18 (1987) 74–88. [26] M.C. Mackey, R. Rudnicki, A new criterion for the global stability of simultaneous cell replication and maturation processes, J. Math. Biol. 38 (1999) 195–219. [27] S.I. Solodushkina, I.F. Yumanovaa, R.H. De Staelen, First order partial differential equations with time delay and retardation of a state variable, J. Comput. Appl. Math. 289 (2015) 322–330. [28] Q. Zhang, C. Zhang, A compact difference scheme combined with extrapolation techniques for solving a class of neutral delay parabolic differential equations, Appl. Math. Lett. 26 (2013) 306–312. [29] J. Niu, L. Sun, M. Xu, J. Hou, A reproducing kernel method for solving heat conduction equations with delay, Appl. Math. Lett. 100 (2020) 106036. [30] B. Esmail, F. Fattahzadeh, Numerical solution of differential equations by using Chebyshev wavelet operational matrix of integration, Appl. Math.Comput. 188 (1) (2007) 417–426. [31] F. Mohammadi, M. Hosseini, S. Mohyud-Din, Legendre wavelet galerkin method for solving ordinary differential equations with non-analytic solution, Int. J. Syst. Sci. 42 (4) (2011) 579–585. [32] M. Khader, A. Hendy, The approximate and exact solution of the fractional order delay differential equations using Legendre seudo-spectral method, Int. J. Pure Appl. Math. 3 (74) (2012) 287–297. [33] P. Rahimkhani, Y. Ordokhani, A numerical scheme based on Bernoulli wavelets and collocation method for solving fractional partial differential equations with Dirichlet boundary conditions, Numer. Methods Part. Differ. Equ. 35 (1) (2019) 34–59. [34] B. Jin, B. Li, Z. Zhou, Subdiffusion with a time-dependent coefficient: Analysis and numerical solution, Math. Comput. 88 (319) (2019) 2157–2186. [35] D. Li, C. Wu, Z. Zhang, Linearized Galerkin FEMs for nonlinear time fractional parabolic problems with non-smooth solutions in time direction, J. Sci. Comput. 2019 (2019) 1–17. [36] L. Li, D. Li, Exact solutions and numerical study of time fractional Burgers’ equations, Appl. Math. Lett. 100 (2020) 106011. [37] M. Heydari, M. Hooshmandasl, F. Mohammadi, Legendre wavelets method for solving fractional partial differential equations with Dirichlet boundary conditions, Appl. Math. Comput. 234 (2014) 267–276. [38] U. Saeed, M. Rehman, M. Iqbal, Modified Chebyshev wavelet methods for fractional delay-type equations, Appl. Math. Comput. 264 (2015) 431–442. [39] M. Usman, M. Hamid, R.U. Haq, W. Wang, An efficient algorithm based on Gegenbauer wavelets for numerical solutions of fractional differential equations of variable-order, Eur. Phys. J. Plus 133 (8) (2018) 327. [40] M. Usman, M. Hamid, M.M. Rashidi, Gegenbauer wavelets collocation based scheme to explore the solution of free bio-convection of nanofluid in 3D nearby stagnation point, Neural Comput. Appl. 2018 (2018) 1–17. [41] A. Wakif, Z. Boulahia, R. Sehaqui, Numerical analysis of the onset of longitudinal convective rolls in a porous medium saturated by an electrically conducting nanofluid in the presence of an external magnetic field, Results Phys. 7 (2017) 2134–2152. [42] A. Wakif, Z. Boulahia, R. Sehaqui, A semi-analytical analysis of electro-thermo-hydrodynamic stability in dielectric nanofluids using Buongiorno’s mathematical model together with more realistic boundary conditions, Results Phys. 9 (2018) 1438–1454. [43] A. Wakif, Z. Boulahia, R. Sehaqui, Numerical study of the onset of convection in a Newtonian nanofluid layer with spatially uniform and non uniform internal heating, J. Nanofluids 6 (1) (2017) 136–148. [44] A. Wakif, Z. Boulahia, S.R. Mishra, M.M. Rashidi, Influence of a uniform transverse magnetic field on the thermo-hydrodynamic stability in water-based nanofluids with metallic nanoparticles using the generalized Buongiorno’s mathematical model, Eur. Phys. J. Plus 133 (5) (2018) 181. [45] Z. Odibat, Approximations of fractional integrals and Caputo fractional derivatives, Appl. Math. Comput. 178 (2) (2006) 527–533. [46] E.H. Doha, A.H. Bhrawy, S.S. Ezz-Eldien, A new Jacobi operational matrix: an application for solving fractional differential equations, Appl. Math. Modell. 36 (10) (2012) 4931–4943.