BDF-type shifted Chebyshev approximation scheme for fractional functional differential equations with delay and its error analysis

Applied Numerical Mathematics 118 (2017) 266–276

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

BDF-type shifted Chebyshev approximation scheme for fractional functional differential equations with delay and its error analysis V.G. Pimenov a,c , A.S. Hendy b,c,∗ a

Institute of Mathematics and Mechanics, Ural Branch of the RAS, 620000, Yekaterinburg, 16 st. S. Kovalevskoy, Russia Department of Mathematics, Faculty of Science, Benha University, Benha, 13511, Egypt Department of Computational Mathematics and Computer Science, Institute of Natural sciences and Mathematics, Ural Federal University, Yekaterinburg 620002, Russia b c

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 13 November 2014 Received in revised form 22 February 2017 Accepted 5 March 2017 Available online 28 March 2017

A numerical method for fractional order differential equations (FDEs) and constant or timevarying delayed fractional differential equations (FDDEs) is constructed. This method is of BDF-type which is based on the interval approximation of the true solution by truncated shifted Chebyshev series. This approach can be reformulated in an equivalent way as a Runge–Kutta method and its Butcher tableau is given. A detailed local and global truncating errors analysis is deduced for the numerical solutions of FDEs and FDDEs. Illustrative examples are included to demonstrate the validity and applicability of the proposed approach. © 2017 IMACS. Published by Elsevier B.V. All rights reserved.

Keywords: Fractional differential equations with delay BDF-type methods Shifted Chebyshev polynomials Runge–Kutta methods Local and global truncating errors

1. Introduction This paper is devoted to present a new BDF-type method based on approximations by truncated shifted Chebyshev series, appropriate for Initial value problems of FDEs and FDDEs.

D (β) y (t ) = f (t , y (t ), y (t − τ )),

t ∈ [0, L ],

0 < β ≤ 1,

(1)

with initial condition

y (t ) = φ(t ),

t ∈ [−τ , 0],

(2)

such that τ > 0 is the delay term which may be constant or variable with respect to t. The Caputo fractional derivative operator D (β) of order β is defined in the following form [12]:

D

(β)

f (x) =

1

(m − β)

x 0

f (m) (ξ )

(x − ξ )ν −m+1

dξ,

ν > 0,

where m − 1 < β ≤ m, m ∈ N, x > 0.

*

Corresponding author. E-mail addresses: [email protected] (V.G. Pimenov), [email protected] (A.S. Hendy).

http://dx.doi.org/10.1016/j.apnum.2017.03.013 0168-9274/© 2017 IMACS. Published by Elsevier B.V. All rights reserved.

V.G. Pimenov, A.S. Hendy / Applied Numerical Mathematics 118 (2017) 266–276

267

It is worth mentioning that FDDEs (1)–(2) will be reduced to FDEs when τ = 0. Fractional differential equations have been an active area of research owing to the intensive development of the theory of fractional calculus itself and their applications in various fields of science, such as physics [23], chemistry [13], engineering [21,26], automatic control [27,7,18], etc. The insertion of a time delay on the right side of the fractional differential equation balances the equation by also adding a degree of system memory on the right side of the equation. By including both fractional derivatives and finite time delays in some models (with memory and heritage properties) [4], more complete and realistic models can be established. Machado [28] proposed the calculation of fractional algorithms based on time-delay systems. The study analyzed the memory properties of fractional operators and their relation with time delay. The effect of delay on the chaotic behavior of fractional order Liu system has been investigated in [3]. A delayed fractional order financial system was proposed and the complex dynamical behaviors of such a system are discussed by numerical simulations [29]. In the past few years, many scholars studied the existence problems of solutions for fractional differential equations. Existence, uniqueness, and structural stability of solutions of nonlinear differential equations of fractional order were studied by Diethelm and Ford [12]. Existence and uniqueness theorems for the initial value problem for the system of fractional differential equations were proposed by Gejji and Babakhani in [9]. Analysis of solutions for FDDEs has been carried in some approaches since time delay is unavoidable in the real systems. By the Banach fixed point theorem and the nonlinear alternative of Leray–Schauder type, Benchohra et al. [2] first discussed the existence of solutions for initial value problems of fractional order 0 < α < 1 with infinite delay and Riemann–Liouville fractional derivative. They also discussed operator type neutral equations in the same way. Deng and Qu [11] discussed the same problem in [2] and obtained some new uniqueness results. By fixed point theorem, Zhou et al. [32] considered the initial value problem of fractional neutral differential equations with infinite delay and Caputo fractional derivative. Agarwal et al. [1] considered the similar problems of Caputo fractional neutral equations with bounded delay. More recently, Yang et al. [31] gave a global existence and uniqueness of initial value problems for nonlinear higher fractional equations with delay by fixed point theory. By fixed point theory the nonlinear alternative of Leray–Schauder type, and the properties of absolutely continuous functions space, the authors obtained some new results involving local and global solutions. Numerical methods for solving such problems which contain both fractional derivatives and time delays are considerably less developed. In this regard, Wang [30] lately approximated the delayed fractional order differential equation by combining the general Adams–Bashforth–Moulton method with the linear interpolation method. Morgado et al. [22] studied numerically a special form of initial value problem for a linear fractional differential equation with finite delay by adaptation of a fractional backward difference method. A new predictor–corrector method had been developed to solve fractional delay differential equations [10]. Related error analysis of the method had been presented and the error bounds also obtained. Bahrawy proposed an accurate and robust approach to approximate the solution of functional Dirichlet boundary value problem with a type of variable order Caputo fractional derivative [5]. A method which is a generalization from finite difference method, has been provided to solve a class of fractional delay differential equations [20]. In this approach, we present a BDF difference scheme based on approximations of Clenshaw and Curtis type [8]. We do not use an approximation of the function f in the right hand side of the fractional differential equation but we use a truncated Chebyshev series instead. This kind of series results to be the interpolating polynomial that passes through certain intermediate points known as interpolation points and impose as in the BDF methods that the derivatives at these points coincide with the values of the function f . The structure of this paper is arranged as follows: we present a definition of shifted Chebyshev polynomials and their analytic formula in the next section. In the third section, a derivation of the method is obtained. The proposed method is represented as one step recurrence formula, local truncating error and global truncating error are clarified in the fourth section. The fifth section is devoted to introduce some numerical examples to ensure the accuracy of the presented method. Finally, paper ends with a brief conclusion. 2. Shifted Chebyshev polynomials The shifted Chebyshev polynomials T n∗ (x) of degree n which are defined on the interval [0, L ] such that T 0∗ (x) = 1,

T ∗ (x) = 2x 1

L

− 1, have the following analytic form [19]:

T n∗ (x) = n

n  (n + k − 1)!22k k (−1)n−k x , (n − k)!(2k)! Lk

(3)

k =0

where, T n∗ (0) = (−1)n , T n∗ ( L ) = 1. The orthogonality condition of these polynomials is

L

T ∗j (x) T k∗ (x) w (x)dx = δ jk hk ,

(4)

0

where, the weight function w (x) = 

1 Lx−x2

, hk =

bk 2

π , with b0 = 2, bk = 1, k ≥ 1.

The function y (x) which belongs to the space of square integrable in [0, L ], may be expressed in terms of shifted Chebyshev polynomials as

268

V.G. Pimenov, A.S. Hendy / Applied Numerical Mathematics 118 (2017) 266–276

y (x) =

∞ 

cn T n∗ (x),

n =0

where the coefficients cn are given by:

L

1

cn =

hn

y (x) T n∗ (x) w (x)dx,

n = 0, 1, 2, ... .

(5)

0

Clenshaw and Curtis [8] introduced an approximation of the function y (x). It is used in terms of shifted Chebyshev polynomials as follows, N 

y N (x) =



an T n∗ (x),

N 2   y (xr ) T n∗ (xr ). N

an =

n =0

(6)

r =0

The summation symbol with double primes denotes a sum with both first and last terms halved. The grid (interpolation) points are chosen to be the Chebyshev–Gauss–Lobatto points associated with the interval [0, L ], xr = 2L − 2L cos( πNr ), r = 0, 1, ..., N. 3. Derivation of the method Let us consider the initial value problem for arbitrary order fractional differential equation with delay (1)–(2). The numerical solution is taken at the point t s+1 = t s + h such that h is the step size. Rewrite the solution y (t ) in the interval [t s , t s+1 ] in terms of a new variable α defined by

t = ts +

hα L

α ∈ [0, L ].

,

Then, the solution of FDDEs will be expressed in the form



y (t ) = y t s +

hα  L

= y¯ (α ),

(7)

by a finite sum approximation of Clenshaw and Curtis type (6)

y¯ (α ) =

N 



an T n∗ (α ),

an =

n =0

N 2   y (t s + ζr h) T n∗ (αr ), N

(8)

r =0

such that

αr

ζr = (β)

Dt

L

αr =

,

y¯ (α ) =

L 2

−

L 2

cos(

πr N

).

N N 2     (β) y (t s + ζr h) T n∗ (αr ) D t T n∗ (α (t )). N n =0

(9)

r =0

Using the analytic form of shifted Chebyshev polynomials (3) and the properties of Caputo fractional derivatives, we get (β) ∗

Dt

T n (α (t )) =

n  n(−1)(n−k) 2(2k) (n + k − 1)! k =0

=

(n − k)!(2k)! Lk

(β)

Dt

 L h

t−

L k  , ts h

n  n(−1)(n−k) 2(2k) (n + k − 1)!(k + 1) L β k−β α . (n − k)!(2k)! Lk (k + 1 − β)hβ

(10)

(11)

k=β

Now,

αk−β can be approximated in terms of shifted Chebyshev series

αk−β ∼ =

N 

ckj T ∗j (α ),

(12)

j =0

where ckj is obtained from (5) as clarified in [14],

ckj =

⎧ ⎨ ⎩

L k−β (k−β+ 12 ) (k−β+1) , k−β j (−1) j−r ( j +r −1)!22r +1 (k+r −β+ 12 ) jL √ , r =0 ( j −r )!(2r )!(k+r −β+1) π

√1

π

for

j = 0;

for

j = 1, 2, ... .

(13)

V.G. Pimenov, A.S. Hendy / Applied Numerical Mathematics 118 (2017) 266–276

269

From Eqs. (11)–(13), we get ∞ 1  T n (α (t )) = β ωβ (n, j ) T ∗j (α ), h

(β) ∗

Dt where

n = β , β + 1, ... ,

(14)

j =0



ωβ (n, j ) = nk=β ψn, j,k , ⎧ 1 n−k 2k ⎨ n(−1) (n+k−1√)!2 k!(k−β+ 2 ) , β 2 L (n−k)!(2k)! π (k−β+1) ψn, j ,k = 1 (−1) j−r ( j +k−1)!22r (k+r −β+ 12 ) k! j ⎩ (−1)n−k nj (n+k−1)!22k+√ , L β (n−k)!(2k)!(k−β+1)

π

r =0

( j −r )!(2r )!(k+r −β+1)

for

j = 0;

for

j = 1, 2, ... .

(15)

Then, ψn, j ,k can be put in the following form

ψn, j ,k =

(−1)n−k 2n(n + k − 1)!(k − β + 12 ) b j L β (k + 12 )(n − k)!(k − β − j + 1)(k + j − β + 1)

,

j = 0, 1, ..., .

(16)

Then, n N (−1)(n−k) 2n(n + k − 1)!(k − β + 12 ) 1   T n (α (t )) = β T ∗j (α ) 1 h b ( k + )( n − k )!( k − β + j + 1 )( k − β − j + 1 ) j 2 j =0 k=β

(β) ∗

Dt

(17)

From Eqs. (8) and (17), we obtain the approximate equation as follows

f (t , y (t ), y (t − τ ))

=

4 Nhβ

N 



n=β

N 

(18) 

N 

n 

(−1)(n−k)n(n + k − 1)!(k − β

b j (k + j =0 k=β

r =0

Evaluating the formula (18) at

1 )(n 2

+

1 ) y (t s 2

+ ζr

h) T ∗ (α n

r)

− k)!(k − β + j + 1)(k − β − j + 1)

T ∗j (α ).

α = 2L − 2L cos( πN ),  = 1, 2, ..., N gives an implicit system of algebraic equations

f (t s + ζ h, y (t s + ζ h), y (t s + ζ h − τ ))

=

4 Nhβ

N  n=β



N 



N 

r =0

n 

(19)

(−1)(n−k)n(n + k − 1)!(k − β

b j (k + j =0 k=β

1 )(n 2

+

1 ) y (t s 2

+ ζr

h) T ∗ (α n

r)

− k)!(k − β + j + 1)(k − β − j + 1)

T ∗j (α ).

Note that for any delayed term τ , (t s + ζ h − τ ) may not be a grid point t  for any . So, we establish an approximation for the delayed function y (t − τ ) as follows

d = y (t  − τ (t  )) = y (t s + ζ h − τ (t s + ζ h)). We define the delayed term as τ = (ms + δ )h such that 0 ≤ δ < 1 and ms is a positive integer, then d = y (t 0 + (s − ms )h + (ζ − δ )h). Two different cases will be considered: The 1st case (δ = 0):



y (t s−ms + ζ h),  ,

for for

s > ms ; s ≤ ms .

(20)

⎧ ⎨ y (t 0 + ζ h),  = y (t 0 ), ⎩ φ(t  − τ (t  )),

for for for

s = ms ; s + ζ = m s ; s + ζ < m s .

(21)

d = Such that

The 2nd case (0 < δ < 1): In this case, y (t − τ ) is not a grid point. So, we interpolate it by the two nearest points,

d = y (t 0 + (s − ms − 1)h + ζ h + (1 − δ )h),

  = y (t s−ms −1 + ζ h) + (1 − δ )h ,

= y (t s−ms −1 + ζ h) + (1 − δ )hy  (t s−ms −1 + ζ h) + O (h2 ), y (t

s−ms −1 + ζ h + h ) − y (t s−ms −1 + ζ h ) y (t s−ms −1 + ζ h) + (1 − δ )h , (1 − δ ) y (t s−ms + ζ h) + δ y (t s−ms −1 + ζ h).

h

270

V.G. Pimenov, A.S. Hendy / Applied Numerical Mathematics 118 (2017) 266–276

Then, we have

⎧ ⎨ (1 − δ ) y (t s−ms + ζ h) + δ y (t s−ms −1 + ζ h), d = (1 − δ ) y (t 0 + ζ h) + δ φ(t  − τ (t  )), ⎩ h¯  ,

Such that



h¯  =

(1 − δ ) y (t 0 ) + δ φ(t  − τ (t  )), φ(t  − τ (t  )),

for for

for for for

s > ms ; s = ms ; s < ms .

(22)

s + ζ = m s ; s + ζ < m s .

(23)

With the aid of Eqs. (19)–(23) and the abbreviation f := f (t s + ζ h, y (t s + ζ h), d ), we get the final form of our method

hβ f =

n N N N 1 4       (−1)(n−k) n(n + k − 1)!(k − β + 2 ) y (t s + ζr h) T n∗ (αr ) ∗ T j (α ), N b j (k + 12 )(n − k)!(k − β + j + 1)(k − β − j + 1) n=β r =0 j =0 k=β

(24)

where the unknowns are the values of the solution at the intermediate points y (t s + ζ h). Explicitly, the algebraic system (24) at (β = 1) and N = 4 is given by the following four equations:

hf1 = −(2 +

√

2) y (t s ) +

√

√

√

2 y (t s + ζ1 h) + 2 2 y (t s + ζ2 h) −

√

2 y (t s + ζ3 h) + (2 −

√

hf2 = 1 y (t s ) − 2 2 y (t s + ζ1 h) + 0 y (t s + ζ2 h) + 2 2 y (t s + ζ3 h) − 1 y (t s + ζ4 h), hf3 = −(2 −

√

2) y (t s ) +

√

√

√

2 y (t s + ζ1 h) − 2 2 y (t s + ζ2 h) −

√

2 y (t s + ζ3 h) + (2 +

√

√ √

2) y (t s + ζ4 h), 2) y (t s + ζ4 h),

hf4 = 1 y (t s ) − (8 − 4 2) y (t s + ζ1 h) + 4 y (t s + ζ2 h) − (8 + 4 2) y (t s + ζ3 h) + 11 y (t s + ζ4 h). Also, the algebraic system (24) at (β = 12 ) and N = 4 is given by the following four equations:

√

√

√

√

√

h1/2 f1 = −( 2(335 + 182 2)a) y (t s ) + ((280 + 320 2)a) y (t s + ζ1 h) − (12 2(1 − 7 2)a) y (t s + ζ2 h)

√ √ − ((56 + 16 2)a) y (t s + ζ3 h) + ((−28 + 43 2)a) y (t s + ζ4 h), √ √ √ h1/2 f2 = −(71 2/d) y (t s ) − ((168 + 40 2)/d) y (t s + ζ1 h) + (180 2/d) y (t s + ζ2 h) √ √ + ((168 − 40 2)/d) y (t s + ζ3 h) − ((29 2)/d)8 y (t s + ζ4 h), √ √ √ √ √ h1/2 f3 = −( 2(335 − 182 2)b) y (t s ) + ((56 − 16 2)b) y (t s + ζ1 h) − (12 2(1 + 7 2)b) y (t s + ζ2 h) √ √ + (−280 + 320 2)by (t s + ζ3 h) + ((28 + 43 2)b) y (t s + ζ4 h), √ √ √ h1/2 f4 = −(43 2c ) y (t s ) + ((168 − 152 2)c ) y (t s + ζ1 h) + (12 2c ) y (t s + ζ2 h) √ √ − ((168 + 152 2)c ) y (t s + ζ3 h) + 335 2c y (t s + ζ4 h), such that

a=



1 105

2

π

sin

π 

,

8

b=

1 105



2

π

cos

π  8

,

c=

1 105



2

π

d = 105

,

√

π.

Solving this system, we obtain in particular the required value at the final point on the interval y (t s +ζ4 h) = y (t s + h) = y s+1 . If we repeat this procedure along the integration interval [0, L ], a discrete solution for the problem in (1)–(2) will be obtained. 4. Detailed error analysis of the proposed method In this section, a local truncating error and global truncating error are analyzed. Let us rewrite the method in (24) at N = 4 by one step recurrence formula

C Y s = B Y s −1 + h β F s ,

(25)

where



F s = f1 , f2 , f3 , f4





T ,

T

Y s = y (t s + ζ1 h), y (t s + ζ2 h), y (t s + ζ3 h), y (t s + ζ4 h)

T

Y s−1 = y (t s−1 + ζ1 h), y (t s−1 + ζ2 h), y (t s−1 + ζ3 h), y (t s−1 + ζ4 h)

.

,

V.G. Pimenov, A.S. Hendy / Applied Numerical Mathematics 118 (2017) 266–276

271

Remark 1. It is noted that

y (t s − ζ4− h) = y (t s−1 + ζ h),

 = 1, 2, 3, 4 ,



(26)

T Y s−1 = y (t s + ζ3 h), y (t s + ζ2 h), y (t s + ζ1 h), y (t s + ζ0 h) . At (β = 1), the constant matrices B and C will be given by

⎛

0 ⎜0 ⎜ B =⎝ 0 0

√ ⎞

0 0 2+ 2 0 0 −1√ ⎟ ⎟, 0 0 2− 2⎠ 0 0 −1

√ √ ⎞ √ √ 2 2 2 −√ 2 2− 2 √ ⎜ −2 2 0√ 2 √2 −1√ ⎟ ⎟. √ C =⎜ ⎝ 2√ −2 2 − 2√ 2 + 2 ⎠ −8 + 4 2 4 −8 − 4 2 11 ⎛

Which have a complete agreement with the results obtained by Ramos and Vigo-Aguiar [25]. At (β = 1/2), the constant matrices B and C can be written as

⎛

0 ⎜0 B =⎜ ⎝0 0

0 0 0 0

√

√

⎞

0 2(335 √ + 182 2)a ⎟ 0 √ 71 2/d √ ⎟, ⎠ 2(335 − 0 √182 2)b 0 43 2c

√ √ √ √ √ ⎞ 43 2)a (280 + 320√ 2)a −12 2(1√− 7 2)a −(56 + 16√ 2)a (−28 +√ ⎜ (−168 − 40 2)/d 180 2/d√ (168 − 40 2√)/d (−29 2√)/d ⎟ ⎟. √ √ C =⎜ ⎠ ⎝ (56 − 16 2)b −12 2(1√+ 7 2)b (−280 + 320 √2)b (28 + 43 √ √ 2)b (168 − 152 2)c 12 2c (−168 − 152 2)c 335 2c ⎛

As the matrix C is non-singular, multiplying the formula (25) by the inverse matrix C −1 yields

Y s = C −1 B Y s −1 + h β C −1 F s .

(27)

⎛

0 ⎜0 − 1 ⎜ Remark 2. The product C B is a square matrix of dimension 4 on the following form ⎝ 0 0

0 0 0 0

0 0 0 0

⎞

1 1⎟ ⎟. 1⎠ 1

Remark 3. The method (27) at (β = 1) can be considered as a Runge–Kutta method where its Butcher tableau [6] is given by α1 α2 α3 α4

√

√

√

1 (88 − 4 2) 384 √ 1 (64 + 48 2) 384 √ 1 (88 + 28 2) 384 1 3

1 (40 − 64 2) 384 1 6 √ 1 (40 + 64 2) 384 1 3

1 (88 − 28 2) 384 √ 1 (64 − 48 2) 384 √ 1 (88 + 4 2) 384 1 3

1 3

1 3

1 3

1 − 16

0 1 − 16

0 0

which ensures a complete agreement with the tableau obtained by Ramos and Aguiar [25]. Consider the difference operators associated with each stage in our method as follows

L i ( y (t ), h) = y (t s + ζi h) − y (t s ) − hβ

4 

c i j D (β) y (t s + ζi h),

i = 1, 2, 3, 4,

(28)

j =1

L i ( y (t ), h) = y (t s + ζi h) − y (t s ) − hβ

4 



c i j f (t s + ζi h, y (t s + ζi h), di ) + O (h2 ) .

(29)

j =1

After expanding y (t s + ζi h) and D (β) y (t s + ζi h) in fractional Taylor series about t s

y (t s + ζi h) = y (t s ) + D (β) y (t s )

(ζi h)β (ζi h)2β + D (2β) y (t s ) + ... (β + 1) (2β + 1)

+ D (5β) y (t s + θζi h)

(ζi h)5β , (5β + 1)

and due to the two different cases of the approximation of the delay term, we obtain:

(30)

272

V.G. Pimenov, A.S. Hendy / Applied Numerical Mathematics 118 (2017) 266–276

The 1st case (δ = 0):

L i ( y (t ), h) = O (h5β ).

(31)

The 2nd case (0 < δ < 1):

L i ( y (t ), h) = O (hq ),

q = min(5β, β + 2).

(32)

Now, we are going to compute the global truncating error. Let us start with the following theorem which will be used later in the proof of global error theorem. Theorem 1. Suppose that y (t ) is the exact solution of the system of FDDEs (1)–(2) and ν (t ) is a polygon formed from the output β discrete points resulted by the formula (27) so that D t ν (t ) = f (t s + ζ h, y (t s + ζ h), d ), t s < t < t s+1 . If the following inequalities are held

ν (t 0 ) − y (t 0 ) ≤ , (β)

Dt

t 0 ∈ [−τ , 0],

(33)

ν (t ) − f (t , ν (t ), ν (t − τ )) ≤ ε,

(34)

f (t , ν (t ), ν (t − τ )) − f (t , y (t ), ν (t − τ )) ≤ l1 ν (t ) − y (t ) ,

(35)

f (t , y (t ), ν (t − τ )) − f (t , y (t ), y (t − τ )) ≤ l2 ν (t − τ ) − y (t − τ ) ,

(36)

such that l1 , l2 are two positive constants, then for t ≥ t 0 , we have the following error estimate

t (t − s)β−1 E β,β (l(t − s)β )ds.

ν (t ) − y (t ) ≤  E β,1 (l(t − t 0 ) ) + ε β

(37)

t0

Proof. For any chosen norm, we investigate the error estimate by

m(t ) = ν (t ) − y (t ) ,

mτ (t ) = ν (t − τ ) − y (t − τ ) ,

(38)

and now, we estimate its natural growth. (β) and using the Operate on both sides of the 1st equation in (38) by time Caputo fractional derivative operator D t triangle inequality yields (β)

Dt

(β)

ν (t ) − f (t , ν (t ), ν (t − τ )) + f (t , ν (t ), ν (t − τ )) − f (t , y (t ), y (t − τ ))

(β)

ν (t ) − f (t , ν (t ), ν (t − τ )) + f (t , ν (t ), ν (t − τ )) − f (t , y (t ), ν (t − τ ))

m(t ) ≤ D t

≤ Dt

(39)

+ f (t , y (t ), ν (t − τ )) − f (t , y (t ), y (t − τ )) . (β)

Define δ(t ) = D t ν (t ) − f (t , ν (t ), ν (t − τ )) which is called the defect of the approximate solution and use Lipchitz conditions (35)–(36), we get (β)

Dt

m(t ) ≤ δ(t ) + l1 m(t ) + l2 mτ (t ),

m(t ) = 

δ(t ) ≤ ε ,

ν (t ) such that δ(t ) ≤ ε (40)

∀t ∈ [−τ , 0].

(41)

Assume that m(t ) ≤ u (t ), mτ (t ) ≤ u (t ) and define l := l1 + l2 . This leads to solve the following fractional differential equation instead of (41) (β)

Dt

u (t ) = l u (t ) + ε ,

u (t 0 ) = ,

(42)

the solution of eq. (42) exists and is unique [24]. The explicit solution of (42) is given by:

t (t − s)β−1 E β,β (l(t − s)β )ds,

u (t ) =  E β,1 (l(t − t 0 )β ) + ε

(43)

t0

such that E β1 ,β2 (l(t − t 0 )β ) is the two parameter Mittag-Leffler function which defined by:

E β1 ,β2 (l(t − t 0 )β ) =

∞  (l(t − t 0 )β )k , (β1k + β2 )

β1 , β2 > 0,

E 1,1 (l(t − t 0 )β ) = exp(l(t − t 0 )β ).

k =0

A combination of m(t ) ≤ u (t ) and (43) yields the desired result (37).

2

V.G. Pimenov, A.S. Hendy / Applied Numerical Mathematics 118 (2017) 266–276

273

Fig. 1. Lady Windermere fan [17].

Theorem 2. Let U be a neighborhood of (t , y (t )|t 0 ≤ t ≤ t N ) where y (t ) is the exact solution of (1)–(2). Suppose that achieved in U and also local error estimates (31)–(32) are valid in U . Then the global error E can be bounded by

E ≤ C  h4β , 

E ≤ C (h

4β

(δ = 0), 2

+ h ),

h = max h s−1 , 1< s ≤ N





C =

≤ l is

(44)

(0 < δ < 1),

where

df dy

(45)

for l ≥ 0; for l < 0.

c 1 , c 2 ,

(46)

Where c 1 , c 2 > 0 and h is small enough for the numerical solution to remain in U . Proof. The global error E = y (t N ) − y N is found along N − s steps of the numerical method (see Fig. 1). Estimate the local 5β 5β β+2 errors e s with the help of (31)–(32) to obtain e s ≤ C h s−1 , (δ = 0) and e s ≤ C (h s−1 + h s−1 ), (0 < δ < 1). From Theorem 1 with ε = 0, we obtain

E s ≤  E β,1 (l(t − t 0 )β ).

(47)

By use of Lady Windermere fan (Fig. 1), it is obvious that

E ≤

N  s =1

With the help of

Es .

(48)



5β

β

β+2

β



h = max h s−1 , ⇒ h s−1 ≤ h4β h s−1 , h s−1 ≤ h2 h s−1 , Eqs. (47)–(48) and due to the two different cases of 1< s ≤ N

the approximation of the delay term, we obtain: The 1st case (δ = 0):

E ≤ C  h4β

N  s =1

β

h s−1 E β,1 (l(t N − t s )β ) = C

 h 4β I .

(49)

The 2nd case (0 < δ < 1):

E ≤ C  (h4β + h2 ) t I≤

N  s =1

N

tt0N t0

β

h s−1 E β,1 (l(t N − t s )β ) = C

 (h4β + h2 ) I ,

E β,1 (l(t N − t )β )dt ≤ c 1 ,

for l ≥ 0;

E β,1 (l(t N − (t + h))β )dt ≤ c 2 ,

for l < 0.

A combination of Eqs. (49)–(51) yields the desired result (44).

(50)

(51)

2

Remark 4. Since the global truncating error E (44) tends to zero as h tends to zero then the proposed method (24) is convergent, i.e. | E | → 0 as h → 0.

274

V.G. Pimenov, A.S. Hendy / Applied Numerical Mathematics 118 (2017) 266–276

Table 1 Maximum absolute error and numerical convergence order are given respectively for the 1st model (52)–(53) for τ = 1,  = 0.1. n=8

n = 16

n = 32

n = 64

2.2 × 10−3

1.1 × 10−3

5.6 × 10−4

2.8 × 10−4

0.98

0.988

0.997

β = 0.50

5.3 × 10−3 1.978

1.35 × 10−3 1.988

3.4 × 10−4 1.999

8.5 × 10−5

β = 0.75

1.6 × 10−4 2.995

2.01 × 10−5 2.998

2.5 × 10−6 2.999

3.1 × 10−7

β = 0.25

Table 2 Maximum absolute error and numerical convergence order are given respectively for the 1st model (52)–(53) for τ = 0.3,  = 0.1. n=8

n = 16

n = 32

n = 64

β = 0.25

8.5 × 10−2 0.97

4.3393 × 10−2 0.985

2.19233 × 10−2 0.993

1.10149 × 10−2

β = 0.50

2.5 × 10−3 1.96

6.4 × 10−4 1.977

1.6 × 10−4 1.986

4.1 × 10−5

β = 0.75

7.8 × 10−3 1.988

1.96 × 10−3 1.993

4.9 × 10−4 1.998

1.2 × 10−4

5. Numerical verification In this section, we illustrate the effectiveness of theoretical results which are obtained in the previous sections. We will consider some real life models in Caputo sense. The first model is based on the effect of noise on light which is reflected from laser to mirror. It has been introduced by Pieroux [20].

D β y (t ) =

−1

x(t ) = 0.9,



y (t ) +

1



y (t ) y (t − τ ),

t ∈ [0, 1],

(52)

t ∈ [−τ , 0].

(53)

The second model is Hutchinson model which is related to the rate of population growth. The rate of population growth in each specific time is dependent on some unique relationships on that time, and its equation is given as [16]

D β y (t ) = r y (t ) − y (t ) = 0.5,

r k

y (t ) y (t − τ ),

t ∈ [0, 1],

(54)

t ∈ [−τ , 0].

(55)

The main aim of this section is to ensure the order of convergence for our method. But the exact solutions for the considered examples are not available, then the maximum absolute error E (h) and convergence order R (h) for the numerical solution of each of the considered examples are calculated using the double mesh principle [15]:

E (h) = max | y s (h) − y s (2h) |, s

R (h) = log2

 E (2h)  E (h)

,

such that the time step is taken as h = 21n . The calculated maximum absolute error and numerical convergence order for the numerical solution of the considered models are arranged in the form of Tables 1–4. 6. Conclusion and remarks In this paper, the construction of an implicit method appropriate for solving fractional differential equations with nonlinear delay is done. This method is of BDF-type which is based on the interval approximation of the true solution by truncated shifted Chebyshev polynomials. Local and global truncating errors are deduced for the numerical solutions of FDDEs. It is proved that the convergence order of the proposed method for FDDEs depending on the two introduced cases of delay term approximation. It is clear that the convergence order of FDEs is equal to 4β . The convergence order of FDDEs is equal to 4β if the delay term is approximated according to the 1st case (δ = 0) and equal to min(4β, 2) if the delay term is approximated according to the 2nd case (0 < δ < 1). Theoretical results coincide with the results obtained from numerical examples in the previous section.

V.G. Pimenov, A.S. Hendy / Applied Numerical Mathematics 118 (2017) 266–276

275

Table 3 Maximum absolute error and numerical convergence order are given respectively for the 2nd model (54)–(55) for τ = 2, r = 0.15, k = 1. n=8

n = 16

n = 32

n = 64

5.3 × 10−4

2.7 × 10−4

1.4 × 10−4

6.97 × 10−5

0.965

0.978

0.983

β = 0.50

6.2 × 10−4 1.977

1.5 × 10−4 1.984

3.98 × 10−5 1.992

1.00 × 10−5

β = 0.75

1.4 × 10−5 2.975

1.8 × 10−6 2.899

2.4 × 10−7 2.994

3.996 × 10−8

β = 0.25

Table 4 Maximum absolute error and numerical convergence order are given respectively for the 2nd model (54)–(55) for τ = 0.2, r = 0.15, k = 1. n=8

n = 16

n = 32

n = 64

β = 0.25

3.3 × 10−4 0.988

1.6 × 10−4 0.996

3.3 × 10−5 0.999

4.2 × 10−5

β = 0.50

2.2 × 10−5 1.976

5.6 × 10−6 1.985

1.4 × 10−6 1.995

3.6 × 10−7

β = 0.75

7.2 × 10−5 1.997

1.8 × 10−5 1.999

4.5 × 10−6 2.01

1.1 × 10−6

Acknowledgment This work was supported by Act 211 Government of the Russian Federation program 02.A03.21.0006 on 27.08.2013. References [1] R.P. Agarwal, Yong Zhou, Yunyun He, Existence of fractional neutral functional differential equations, Comput. Math. Appl. 59 (3) (2010) 1095–1100, http://dx.doi.org/10.1016/j.camwa.2009.05.010. [2] M. Benchohra, J. Henderson, S.K. Ntouyas, A. Ouahab, Existence results for fractional order functional differential equations with infinite delay, J. Math. Anal. Appl. 338 (2008) 1340–1350. [3] Sachin Bhalekar, Varsha Daftardar-Gejji, Fractional ordered Liu system with time-delay, Commun. Nonlinear Sci. Numer. Simul. 15 (8) (2010) 2178–2191, http://dx.doi.org/10.1016/j.cnsns.2009.08.015. [4] Sachin Bhalekar, Varsha Daftardar-Gejji, Dumitru Baleanu, Richard Magin, Fractional Bloch equation with delay, Comput. Math. Appl. (ISSN 0898-1221) 61 (5) (2011) 1355–1365, http://dx.doi.org/10.1016/j.camwa.2010.12.079. [5] A.H. Bhrawy, M.A. Zaky, Numerical algorithm for the variable-order Caputo fractional functional differential equation, Nonlinear Dyn. (ISSN 1573-269X) 85 (3) (2016) 1815–1823, http://dx.doi.org/10.1007/s11071-016-2797-y. [6] J.C. Butcher, Numerical Methods for Ordinary Differential Equations, John Wiley & Sons, 2008. [7] Fulai Chen, Yong Zhou, Attractivity of fractional functional differential equations, Comput. Math. Appl. 62 (3) (2011) 1359–1369, http://dx.doi.org/ 10.1016/j.camwa.2011.03.062. [8] Charles W. Clenshaw, Alan R. Curtis, A method for numerical integration on an automatic computer, Numer. Math. 2 (1) (1960) 197–205, http:// dx.doi.org/10.1007/BF01386223. [9] V. Daftardar-Gejji, A. Babakhani, Analysis of a system of fractional differential equations, J. Math. Anal. Appl. 293 (2) (2004) 511–522, http://dx.doi.org/ 10.1016/j.jmaa.2004.01.013. [10] V. Daftardar-Gejji, Y. Sukale, S. Bhalekar, Solving fractional delay differential equations: a new approach, Fract. Calc. Appl. Anal. 18 (2) (2015) 400–418, http://dx.doi.org/10.1515/fca-2015-0026. [11] Jiqin Deng, Hailiang Qu, New uniqueness results of solutions for fractional differential equations with infinite delay, Comput. Math. Appl. 60 (8) (2010) 2253–2259, http://dx.doi.org/10.1016/j.camwa.2010.08.015. [12] K. Diethelm, N.J. Ford, Analysis of fractional differential equations, J. Math. Anal. Appl. 265 (2) (2002) 229–248, http://dx.doi.org/10.1006/ jmaa.2000.7194. [13] Y. Ding, H. Ye, A fractional-order differential equation model of HIV infection of CD4+ T-cells, Math. Comput. Model. 50 (3–4) (2009) 386–392, http://dx.doi.org/10.1016/j.mcm.2009.04.019. [14] E.H. Doha, A.H. Bhrawy, S.S. Ezz-Eldien, A Chebyshev spectral method based on operational matrix for initial and boundary value problems of fractional order, Comput. Math. Appl. (ISSN 0898-1221) 62 (5) (2011) 2364–2373, http://dx.doi.org/10.1016/j.camwa.2011.07.024. [15] E.P. Doolan, J.J.H. Miller, W.H.A. Schilders, Uniform Numerical Methods for Problems with Initial and Boundary Layers, Boole Press Ltd, 1980. [16] A.C. Fowler, Asymptotic methods for delay equations, J. Eng. Math. 53 (3–4) (2005) 271–290, http://dx.doi.org/10.1007/s10665-005-9016-z. [17] Ernst Hairer, S.P. Norset, Gerhard Wanner, Solving Ordinary Differential Equations I, 1993. [18] Ying Luo, YangQuan Chen, Fractional order [proportional derivative] controller for a class of fractional order systems, Automatica (ISSN 0005-1098) 45 (10) (2009) 2446–2450, http://dx.doi.org/10.1016/j.automatica.2009.06.022. [19] John C. Mason, David Cristopher Handscomb, Chebyshev Polynomials, Chapman and Hall/CRC, 2003. [20] Behrouz P. Moghaddam, Zeynab S. Mostaghim, A numerical method based on finite difference for solving fractional delay differential equations, J. Taibah Univ. Sci. 7 (3) (2013) 120–127, http://dx.doi.org/10.1016/j.jtusci.2013.07.002. [21] Concepción A. Monje, Blas M. Vinagre, Vicente Feliu, YangQuan Chen, Tuning and auto-tuning of fractional order controllers for industry applications, Control Eng. Pract. 16 (7) (2008) 798–812, http://dx.doi.org/10.1016/j.conengprac.2007.08.006.

276

V.G. Pimenov, A.S. Hendy / Applied Numerical Mathematics 118 (2017) 266–276

[22] M.L. Morgado, N.J. Ford, P.M. Lima, Analysis and numerical methods for fractional differential equations with delay, J. Comput. Appl. Math. 252 (2013) 159–168, http://dx.doi.org/10.1016/j.cam.2012.06.034. [23] Guojun Peng, Synchronization of fractional order chaotic systems, Phys. Lett. A 363 (5–6) (2007) 426–432, http://dx.doi.org/10.1016/ j.physleta.2006.11.053. [24] I. Podlubny, Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, Mathematics in Science and Engineering, vol. 198, Academic Press, San Diego, CA, USA, 1998. [25] Higinio Ramos, Jesús Vigo-Aguiar, A fourth-order Runge–Kutta method based on BDF-type Chebyshev approximations, J. Comput. Appl. Math. 204 (1) (2007) 124–136, http://dx.doi.org/10.1016/j.cam.2006.04.033. [26] H. Saeedi, M. Mohseni Moghadam, N. Mollahasani, G.N. Chuev, A CAS wavelet method for solving nonlinear Fredholm integro-differential equations of fractional order, Commun. Nonlinear Sci. Numer. Simul. 16 (3) (2011) 1154–1163. [27] Zhixin Tai, Xingcheng Wang, Controllability of fractional-order impulsive neutral functional infinite delay integrodifferential systems in Banach spaces, Appl. Math. Lett. 22 (11) (2009) 1760–1765, http://dx.doi.org/10.1016/j.aml.2009.06.017. [28] J.A. Tenreiro Machado, Time-delay and fractional derivatives, Adv. Differ. Equ. 2011 (2011) 934094, http://dx.doi.org/10.1155/2011/934094. [29] Z. Wang, X. Huang, G.D. Shi, Analysis of nonlinear dynamics and chaos in a fractional order financial system with time delay, Comput. Math. Appl. 62 (3) (2011) 1531–1539. [30] Zhen Wang, Xia Huang, Jianping Zhou, A numerical method for delayed fractional-order differential equations: based on G-L definition, Appl. Math. Inf. Sci. 7 (2L) (2013) 525–529, http://dx.doi.org/10.12785/amis/072L22. [31] Z. Yang, J. Cao, Initial value problems for arbitrary order fractional differential equations with delay, Commun. Nonlinear Sci. Numer. Simul. 18 (11) (2013) 2993–3005, http://dx.doi.org/10.1016/j.cnsns.2013.03.006. [32] Yong Zhou, Feng Jiao, Jing Li, Existence and uniqueness for fractional neutral differential equations with infinite delay, Nonlinear Anal., Theory Methods Appl. 71 (7–8) (2009) 3249–3256, http://dx.doi.org/10.1016/j.na.2009.01.202.