Chebyshev operational matrix for solving fractional order delay-differential equations using spectral collocation method

Abstract In this article, the solution for general form of fractional order delay-differential equations (GFDDEs) is introduced. The proposed GFDDEs have multi-term of integer and fractional order derivatives for delayed or non delayed terms. An operational matrix is presented for all terms. The spectral collocation method is used to solve the proposed GFDDEs as a matrix discretization method. The reliability and efficiency of the proposed scheme are demonstrated by some numerical experiments.


Introduction
Spectral methods are one of the principal methods of matrix discretization for the solutions of differential equations. The main feature of these methods lies in their accuracy for a given number of unknown (see, for instance Doha & Bhrawy, 2008;Doha, Bhrawy, & Saker, 2011;Ramos, Hussaini, Quarteroni, & Zang, 1988).
Few years ago, many authors presented numerical treatments for a kind of FDEs contain a delay term, called Fractional-Delay-differential equations (FDDEs). In Morgado, Ford, and Lima (2013) and Xu and Lin (2016), the existence and uniqueness of its solution presented. Legender pseudo-spectral method (Khader & Hendy, 2012), Bernoulli wavelets (Rahimkhani, Ordokhani, & Babolian, 2017), Hermite wavelet method (Saeed & ur Rehman, 2014), collocation method with shifted Jacobi polynomials (Muthukumar & Ganesh Priya, 2017) and collocation method with Chebyshev polynomials (Muthukumar & Ganesh Priya, 2017) were presented as a spectral numerical solution for FDDEs. All of the previous work considered as a single fractional term in the right hand side and a delay term in the left hand side free of integer or fractional order derivatives.
In this article, an operational matrix of fractional derivatives is presented for Chebyshev polynomials, and employs it to deal with a general form of GFDDEs. The proposed model of GFDDEs is chosen to be multi-term of fractional derivatives and also integer order derivatives, where the delayed terms are taken to be multi-term of fractional and integer order derivatives too. The presented operational matrix is used with the collocation method to solve the proposed GFDDEs as a matrix discretization method. The high accuracy of this method is verified through some numerical examples. The obtained numerical results are compared with other methods, where they show the proposed method gives good accuracy.
The article is organized as: Preliminary and notations for some of properties for Chebyshev polynomials and fractional calculus in Section 2. The introduced operational matrix of fractional derivatives is in Section 3, where the proposed model of GFDDEs and the fundamental matrix relations are in Section 4. In Section 5 the method of solution is presented. Section 6 contains the numerical results and the comparisons.

Preliminaries and notations
In this section, we present some definitions and some properties for the fractional derivative and the Chebyshev polynomial of the first kind.

The Caputo fractional derivative
Definition 1. The Caputo fractional derivative operator D of order is defined in the following form: where mÀ1 < m; m 2 N; x > 0: where k and l are constants.
for n 2 N 0 and n < de Cðnþ1Þ Cðnþ1ÀÞ x nÀ ; for n 2 N 0 and n ! de ; ( where de denote to the smallest integer greater than or equal to , and N ¼ f0; 1; 2; :::g:

Chebyshev polynomials of the first-kind
The Chebyshev polynomials T n ðxÞ of the first-kind are orthogonal polynomials of degree n in x defined on ½À1; 1 In this case, we are going to use the last row for odd values of N ¼ 2l þ 1; otherwise previous one will be the last row of matrix M (N ¼ 2l). Now, from Eq. (5) we can obtain the k th derivative of the matrix T(x) as: ÞM T ; k ¼ 0; 1; 2; ::::

Operational matrices
In this section, we introduce the generalize operational matrices for T(x), T ðkÞ ðxÞ; TðxÀsÞ; T ðsÞ ðxÀsÞ; D i TðxÞ and D a i TðxÀsÞ according to fractional calculus.
x nþ1À i ; :::; w Corollary 3.2. The a th i order fractional derivative of the delay vector TðxÀsÞ can be written in following form: Proof. By using Eq. (11) and by putting x ! ðxÀsÞ we get w where X a i x ð Þ ¼ 0; 0; :::0; x nÀa i ; ::::: if nÀ1 < a i < n; n 2 N; 0 0::: 0 0 0 0::: ; nÀ1 < a i < n; n 2 N; and if 0 < a i < 1; then 4. Application to fractional order delay-differential equation The general from of the linear fractional order Delaydifferential equation is presented in this section. In addition, the general operational matrices for all terms will obtained. Now, we consider the general fractional order Delay-differential equation: under the conditions where x 2 ½Às; 0 and Q i ðxÞ; Q Ã j ðxÞ; P j ðxÞ; P Ã s ðxÞ are known functions, i ¼ 0; 1; :::; n 1 ; j ¼ 0; 1; :::n 2 ; k ¼ 0; 1; :::; n 3 and s ¼ 0; 1; :::n 4 : Also m is the greatest integer order derivative exist, or the highest integer order greater than the fractional derivative. Let us write Eq. (21) in the form: where Therefore, we consider the approximate solution in the form: where a r are unknown Chebyshev coefficients and N is chosen any positive integer such that N ! m: From Eq. (25) we get, and where A ¼ 1 2 a 0 ; a 1 ; a 2 ; :::; a N ! T : By putting x ! xÀs in the relation Eq. (27) we obtain the matrix form Consequently, by substituting the matrix form Eq. (7) into (27) we have the matrix relation By substituting the matrix forms Eq. (11) into Eq. (28) we have the matrix relation By using Eqs. (17) and (29) we get Similar form Eqs. (10) and (30) we obtain To obtain the solution Eq. (25) for Eqs. (21) and (22), we can use the collocation points defined in this form: ; l ¼ 0; 1; 2; :::; N: By using Chebyshev collocation method with Eq. (35), we obtain the system Eq. (23) as: where Then, we can write the system Eq. (36) in the form where ;

Matrix representation of fractional differential-delay part
We can obtain the second term F(x) in Eq. (37) by using Eqs. (33) and (35) in this form:

Matrix representation of differential part
From Eqs. (31) and (35) the third term L(x) in Eq. (37) can be write as:

Matrix representation of differential-delay part
By using Eqs. (35) and (34)

Matrix relation for the conditions
Finally, we can obtain the matrix form for the conditions Eq. (22) by using Eq. (26) on the form: or where ; u i1 ; ::::; u iN ½ :

Error estimation
The collocation method is the most common spectral methods used to solve differential equations, it is easy to implement once the differentiation matrices are computed. The coefficient matrix of the collocation system is always full with a condition number behaving like OðN 2m Þ (m is the order of the differential equation) see Du (2016), Shen, Tang, andWang (2011), andWang, Zhang, andZhang (2014). In addition, if the exact solution exists then the error is estimated form: where y(x) the exact solution and y N ðxÞ approximate solution. We can easily check the accuracy of the suggested method. Since the truncated Chebyshev series Eq. (25) is an approximate solution of (21), when the solution y N ðxÞ and its derivatives are substituted in Eq. (21), the resulting relation must be satisfied approximately; that is, for x ¼ x p 2 ½0; 1; p ¼ 0; 1; 2; ::: and Eðx p Þ 10 £ p (£ p positive integer). If max 10 £ p ¼10 ÀL (£ positive integer) is prescribed, then the truncation limit N is increased until the difference Eðx p Þ at each of the points becomes smaller than the prescribed 10 L : On the other hand, the error can be estimated by the function If E N ðxÞ ! 0; when N is sufficiently large enough, then the error decreases.

Applications and numerical results
In this section, we introduce some numerical examples for fractional order delay-differential equation to illustrate the above results. All results are obtained by using Mathematica 7 programming.
which is the exact solution of the problem Eq. (52).
Example 2. Consider of the linear delay-differential equation (G€ ulsu & Sezer, 2006) 2y with the initial condition yð0Þ ¼ 1; y 0 ð0Þ ¼ 2 and the exact solution is yðxÞ ¼ 2e x þ x 2 À1: Then for N ¼ 5 with Eqs. (25) and (35), and the fundamental matrix equation of the problem is defined by After the augmented matrices of the system and condition are computed, we obtain the solution at N ¼ 5 as: A ¼ 2:03555 2:24761 1:0474 0:0846814 Â 0:0118535 0:00128608: Table 1 represents comparison between exact solution, approximate solution using the present method and solution presented by Taylor method (G€ ulsu & Sezer, 2006) at N ¼ 8. Figure 1 shows the behaviour of the exact solution, present method and Taylor method for at N ¼ 8. Example 3: Consider the linear fractional order delay-differential equation (Muthukumar & Ganesh Priya, 2017) The subjected initial condition is yð0Þ ¼ 0 and the exact solution is yðxÞ ¼ x 2 ; where: Q Ã 0 ðxÞ ¼ 1; P 0 ðxÞ ¼ 0; P Ã 0 ðxÞ ¼ 1 and gðxÞ ¼ 2x þ Cð3Þ Cð1:5Þ x 1:5 À1: Then the fundamental matrix equation of the problem Eq. (59) is defined by The approximate the solution y N ðxÞ introduced by present method with N ¼ 5 as: Then the solution is given as: This problem is mentioned in Muthukumar and Ganesh Priya (2017) using Jacobi collocation method, Table 2 shows the exact, present and Jacobi (Muthukumar & Ganesh Priya, 2017) solutions. Also, Figure 2 introduces the behaviour of the exact, approximate solutions at N ¼ 6.  Figure 2. The behaviour of the exact and approximate solutions at N ¼ 6.  Example 4. Consider the following linear initial value problem (Kazem, 2013) D y x ð Þ þ y x ð Þ ¼ 0: The subjected initial conditions are yð0Þ ¼ 1; y 0 ð0Þ ¼ 0; and the exact solution is given by: The fundamental matrix equation of the problem is defined by: By using the initial conditions we can write Eq. (64) as: where G ¼ 1; 0; 0; ::::; 0 ½ T : The approximate solution given with N ¼ 10. Table 3 gives the comparison between present method and (Kazem, 2013) with different of root mean square (RMS), where RMS is given as: where yðx k Þ and y N ðx k Þ are achieved by exact and numerical solution on x k and M is number of test points. In addition, Figure 3 shows the numerical results of y N ðxÞ for N ¼ 10 and ¼ 0:2; 0:5; 0:8 and 1.2.
Example 5. Consider the fractional-delay-differential equation (Rahimkhani et al., 2017) with the initial condition yð0Þ ¼ 0; y 0 ð0Þ ¼ 0 and the exact solution is yðxÞ ¼ x 2 Àx when a ¼ 1, s ¼ 0:01 by using Eqs. (25) and (35) with Q 0 ðxÞ ¼ 1; P Ã 0 ðxÞ ¼ 1; Q Ã 0 ðxÞ ¼ 1; gðxÞ ¼ :at N ¼ 6. Then the fundamental matrix equation of the problem is defined by Our numerical results are presented in Table 4. Figure 4 displays the approximate solutions obtained for values of a ¼ 0:5; 0:7; 0:9 and the exact solution with N ¼ 6 and s ¼ 0:01: From these results, it is seen that the approximate solutions converge to the exact solution. Figure 5 displays the approximate solutions obtained for different values of a with the exact solution.

Conclusion
In this work, the general form of fractional order Delay-differential equations (GFDDEs) is presented. The spectral collocation method is used for solving GFDDEs. All terms in the proposed equation reduced by an operational matrices based on Chebyshev polynomials to matrix form. The accuracy of this method is obtained by many numerical examples. Finally, we used the Mathematica 7 to calculate our results.

Disclosure Statement
No potential conflict of interest was reported by the authors.