An effective numerical method for solving fractional delay differential equations using fractional-order Chelyshkov functions

In this paper, the fractional-order Chelyshkov functions (FCHFs) and Riemann-Liouville fractional integrals are utilized to ﬁnd numerical solutions to fractional delay diﬀerential equations, by transforming the problem into a system of algebraic equations with unknown FCHFs coeﬃcients. An error bound of FCHFs approximation is estimated and its convergence is also demonstrated. The eﬀectiveness and accuracy of the presented method are established through several examples. The resulting solution is accurate and agrees with the exact solution, even if the exact solution is not a polynomial. Moreover, comparisons between the obtained numerical results and those recently reported in the literature are shown.

A delay differential equation is a type of differential equation where the derivative of the unknown function at the current time depends not only on the solution at the current time, but also on the solution at previous times.Fractional delay differential equations can be considered as a generalization of delay differential equations of integer orders.They have been used for modeling in different fields of science, such as control [20], biology [9], optics [41] and electrical networks [16].Many researchers have examined the existence and uniqueness of solutions to fractional delay differential equations [1,27].
Generally, It is a challenge to find exact solutions for most fractional delay differential equations, hence numerical or approximate methods should be developed to solve these equations.Various techniques have been presented to address this type of problem in recent decades.In [31], a numerical method based on finite differences has been presented to approximate solutions of fractional delay differential equations.The development of a predictor-corrector method with error analysis has been introduced in [14].The Bernoulli wavelet operational matrix of fractional integration together with the collocation method has been employed in [42] to provide an approximate solution of fractional delay differential equations.In [21], the collocation method, optimization techniques and generalized Laguerre polynomials are applied to approximate a solution of delay fractional differential equations.A numerical technique based on hybrid functions of block-pulse and fractional-order Fibonacci polynomials has been given in [43].They used Riemann-Liouville fractional integral operational matrix and delay operational matrix of the hybrid functions together with the collocation method to reduce the problem to a system of algebraic equations.
In [15], a computational technique based on the combination of fractional Bessel functions with block-pulse functions to produce a sparse operational matrix to reduce the problem to a system of algebraic equations.A spectral operational matrices-based algorithm has been presented in [48].These matrices have been generated by using shifted Gegenbauer polynomials.Instead of using the operational matrix of fractional integration, which requires some approximation, a numerical method based on the explicit formula of Riemann-Liouville fractional integral of Legendre wavelets functions has been given in [50].A numerical method based on applying the exact formula for the Riemann-Liouville fractional integral of Taylor wavelet functions together with the collocation method has been proposed in [47].Recently, Avci [10] has introduced a numerical spectral technique depending on the operational matrix of fractional integration for fractional-order Taylor basis.
Approximation by orthogonal functions has grown significantly as a useful tool for solving fractional differential equations.They convert the fractional problem into a set of algebraic problems.Some popular orthogonal families are Chebyshev polynomials, Legendre polynomials, Laguerre polynomials, Taylor polynomials, Chelyshkov polynomials and block-pulse functions.Chelyshkov polynomials were presented in [13] and have been used successfully for solving different kinds of fractional differential equations.Solving fractional differential equations with basis functions of integer-order may lead to some issues when the solutions contain terms with fractional powers [19].The major weaknesses are that the rate of convergence is weak and a large number of basis functions are required to provide accurate results.These drawbacks can be avoided by approximating the solutions using orthogonal functions of fractional order.
Based on the above considerations, the motivation of this paper is to extend the application of FCHFs to present a new numerical method for solving fractional delay differential equations.In order to accomplish this goal, the operational matrix of fractional integral for FCHFs is derived and used together with the collection method to convert the frac-tional delay problem into a set of algebraic equations.To the best of our knowledge, this is the first time that the operational matrix of fractional integration for FCHFs is applied to solve fractional delay differential equations The main advantages of the presented technique are: • FCHFs are fractional-order orthogonal functions.This feature leads to a good rate of convergence and reduces the number of basis functions required to get satisfactory results.• The solution procedure is simple since the fractional delay problem is reduced to a set of algebraic equations.• The computational errors are less than in the conventional methods where the obtained numerical results in some cases are the exact solutions and in others are in good agreement with the exact solutions.The structure of the paper is as follows.Section 2 presents some fundamental concepts of fractional calculus that will be used in our subsequent development.In Sect.3, the FCHFs are described, and the operational matrix of fractional integration is also constructed.The convergence of FCHFs is demonstrated in Sect. 4. Section 5 describes our new numerical method for solving fractional delay differential equations.Section 6 considers some numerical examples and reports our numerical results as well as some comparisons with other methods.The paper is finally summarized in Sect.7.

Preliminaries on fractional calculus
Here, we review some definitions and properties of fractional integration and differentiation that will be used in the next sections [36,39].

Definition 2.1
The Riemann-Liouville fractional integration of order ν ≥ 0 of a function g(t) is defined as where t ν-1 * g(t) is the convolution of t ν-1 and g(t).

Definition 2.2
The Caputo fractional derivative of order ν > 0 of a function g(t) is defined as where ν denotes the smallest integer greater than or equal to ν.
If β = 1, then the generalized Taylor's formula converts to the classical Taylor's formula.

Fractional-order Chelyshkov functions
The fractional-order Chelyshkov functions (FCHFs) can be defined by using Chelyshkov polynomials in the following form [5]: The set {ϕ β mr } m r=0 of FCHFs is a family of orthogonal functions over the interval ρ = [0, 1] with respect to the weight function w β (t) = t β-1 .The orthogonality property for these functions is It is clear from Equation (3.1) that every member of the set {ϕ β mr } m r=0 is of degree βm.This is the key distinction between the Chelyshkov polynomials and the other sets of orthogonal w β (ρ) be the set of all real-valued measurable functions u(t) on ρ that satisfy with inner product and norm given by ) w β (ρ), then for any function g(t) there exists a unique best approximation g m (t) ∈ A m that minimizing the distance to g(t), that is Since g m (t) ∈ A m , then there are unique coefficients 0 , 1 , . . ., m such that where L and β (t) are given by and

The operational matrix of fractional integration
This subsection is devoted to construct the fractional integration operational matrix of the FCHFs, which will be used to reduce the fractional delay differential equations into a set of algebraic equations.
Theorem 3.1 Let β (t) be a (m+1)×1 vector of FCHFs.The Riemann-Liouville fractional integral of order ν > 0 of the vector β (t) can be given by where (ν) = (λ rk ) m r,k=0 is the (m + 1) × (m + 1) operational matrix of fractional integral of order ν for FCHFs and its elements can be obtained by Proof Using Equation (2.1) to integrate Equation (3.1), yields By taking the Laplace transform of Equation (3.11), we get and Then Equation (3.12) becomes Now, by utilizing the inverse Laplace transform of Equation (3.15), we obtain Approximating t jβ+ν in terms of FCHFs, we get where jk can be determined from Equations (3.1) and (3.7) as follows From Equations (3.16), (3.17) and (3.18), we can write which leads to the desired result.

Convergence analysis
In this section, we discuss the convergence of the best approximation and show that by increasing the number of the FCHFs, the error in calculating the operational matrix of fractional integration tends to zero.
Proof Assume the generalized Taylor's formula of g(t) is denoted by ψ(t), then from Definition 2.3 the error bound is where is the best estimation of g(t) in A m and ψ(t) ∈ A m then by using Equation (4.2), we get The right-hand side of inequality (4.3) depends on both m and β, so for a fixed β, we have Theorem 4.2 Let the error vector E ν of the operational matrix of fractional integration (ν) be given by Proof Suppose that E ν = 0 , 1 , . . ., m T , then by using Equations (3.10), (3.16) and (3.18), we get Hence, Equation (4.8) becomes

Numerical method
The aim of this section is to present a new numerical method for solving fractional delay differential equations by using FCHFs and their operational matrix of fractional integration.
Step 6: Find the unknowns of the algebraic system obtained in Step 5. Output: An approximate solution given by equation (5.6) or (5.7).

Numerical results
To demonstrate the applicability and effectiveness of the introduced method, we use it to solve some examples, and we compare the numerical results it generates with those reported in the literature.All numerical results have been obtained using Mathematica 11 software and a laptop with Intel(R) Core(TM) i7-4600M 2.90 GHz CPU and 8.0 GB RAM.
For ν = 1 2 , we take m = 4, β = 1 2 , and again we can obtain the exact solution, where Similarly, by using Equation (6.10), subdivide the interval (0, 1] into four subintervals and consider the following cases of τ then collocate Equations ( 6.3) at the points t s = s 4 , s = 0, 1, . . ., 4, for each case of τ and solving the resulting systems for the unknowns 0 , 1 , 2 , 3 , 4 , we get the same solution for all cases of τ , thus which is the exact solution.Additionally, the total CPU time needed to find this solution using the proposed method for the four cases of τ in Equation (6.11) is 24.0608 seconds.
Problem (6.1) has been solved in [10,47,50].By comparing the obtained solution with those in [10,47,50], it can be observed that the introduced method is more accurate because we get the exact solution for ν = 1 and for any value of τ with three terms of FCHFs, whereas they used four terms of Legendre wavelets, six terms of Taylor wavelets and eight terms of fractional-order Taylor functions in [47,50] and [10] respectively to obtain only approximate solutions.
Example 2 Consider the following fractional delay differential equation [8,10]: The exact solution is given by y(t) = t 2 .By applying the method proposed in Sect.5, we obtain (6.15) Equations (6.15) transform Equation (6.14) to (2.5) (2.5) t  .17)which leads to the exact solution.Also, the total CPU time required to find the solution using the presented method in this example is 0.2969 seconds.This problem has been solved by using the fractional-order Taylor method [10], the Haar wavelet collocation technique [8] and the fractional backward difference method [32].Their solutions were approximate solutions.By comparing with these methods, it is clear that the proposed method is more effective, more accurate and less time-consuming than these methods.
Example 3 Consider the following fractional delay differential equation [50]: (2.8) (1.8) t 0.8 + t 2 -3t + 1, t ≥ 0, The exact solution to this problem is y(t) = t 2t -1. Figure 3 displays the approximate solution obtained for various values of β at m = 10 with the exact solution, while Fig. 4 plots the absolute error function at β = 0.2 and m = 10.It can be observed that our numerical solutions agree with the exact solution for all choices of β and high-accuracy numerical solutions can be obtained with some terms of FCHFs.Table 1 shows a comparison between the results obtained by the present method and those reported in [50] in terms of absolute errors, where M denotes the number of basis functions used in [50] for solving the problem.From Table 1, it can be seen that the absolute errors achieved by the introduced method are less than those in [50] for all corresponding cases which shows the advantage of using basis functions of fractional-order.In addition, we have used only eleven basis functions, while they used more than a hundred basis functions in [50].This demonstrates that the introduced method is in more agreement with the exact solution than [50] for this problem and can identify high-precision solutions at reasonable computational costs.Example 4 Consider the following fractional delay differential equation [15,50]: This problem has the exact solution y(t) = e -t when ν = 3. Figures 5 and 6 present graphically the obtained results for the approximate solution for various values of ν with the exact solution for ν = 3 and the absolute error at ν = 3, respectively.From Fig. 5, it can be seen that as ν approaches 3, the approximate solution converges to the exact solution for ν = 3 and from Fig. 6, we can see that the absolute error is at most of order 10 -14 , which means that the present method has good accuracy using a reasonable number of basis functions.Table 2 compares the results obtained by the present method with the collected results in [15] in terms of absolute errors.It can be observed that the results obtained by the present method are less than those in [15,22,38,46] for most of the listed values in that table by using a fewer number of basis functions for each case.Also, Problem (6.19) has been solved using the Legendre wavelet method [50] by selecting M = 14.The maximum absolute error achieved by this method is approximately 8.78 × 10 -12 while our maximum absolute error is approximately 1.74 × 10 -15 at m = 10.This means that the introduced method outperforms the method in [50] with respect to accuracy.These comparisons show that the presented method is more accurate, efficient, and coincidental with the exact solution than [15,22,38,46,50].2.20 ×10 -14 0 0 0 0 0 0.2 2.47 ×10 -10 5.11 ×10 -15 3.05 ×10 -7 1.86 ×10 -15 1.00 ×10 -10 8.54 ×10 -8 3.70 ×10 -7 0.4 4.12 ×10 -11 3.77 ×10 -15 9.81 ×10 -7 3.79 ×10 -13 0 5 .3 6 ×10 -6 2.38 ×10 -6 0.6 6.32 ×10 -11 2.66 ×10 -15 1.85 ×10 -6 8.11 ×10 -12 1.00 ×10 -10 5.95 ×10 -5 5.97 ×10 -6 0.8 2.35 ×10 -10 3.11 ×10 -15 2.64 ×10 -6 7.00 ×10 -11 1.00 ×10 -10 3.26 ×10 -4 3.48 ×10 -5 1.0 1.12 ×10 -9 6.52 ×10 -14 3.03 ×10 -6 3.71 ×10 -10 2.00 ×10 -10 1.21 ×10 -3 2.03 ×10 -4 CPU time 0.0625 s 0.1562 s ----- Example 5 Consider the following fractional delay differential equation [21]: We have y(t) = t q as the exact solution for this problem.By using the technique described in Sect.5, we get The fractional delay differential equation (6.20) can be written after using Equations (6.21) in the form Consider ν = 2 and q = 2 with the exact solution y(x) = t 2 .In order to solve this case, we take m = 2, and β = 1, hence By using Equation (6.23) and collocate Equation (6.22) at the points 0, 1 2 , 1, we get the following system 732 0 + 1 -3 2 -480 = 0, 1186 0 -1167 1 -419 2 + 2940 = 0, 964 0 -2113 1 + 1579 2 -1680 = 0. ( Solving Equations (6.24) gives the solution hence which is the exact solution.Moreover, the total CPU time required to find this solution using the introduced method in this example is 0.0156 seconds.Therefore, the presented method is accurate and time-efficient.Additionally, in case ν = 3 2 and q = 3 2 , we can again obtain the exact solution y(x) = t which leads to the exact solution for this case.Also, the total CPU time needed to find the solution using the proposed method for this case is 0.0937 seconds.The collocation Laguerre operational matrix technique [21] has been used to solve Problem (6.20).By comparing with our solution, it is clear that our method is more efficient than [21], because we utilized three terms of FCHFs to get the exact solution for ν = 2, while they used six terms of the generalized Laguerre polynomials to obtain only approximate solution.
Example 6 Consider the following fractional delay differential equation [33,42]: (6.28) The exact solution to this problem at ν = 1 is y(t) = e -t .Figures 7 and 8 display the approximate solution obtained for several values of ν and τ together with the exact solution for ν = 1, respectively.Figure 7 shows that as ν approaches 1, the approximate solution converges to that of the integer-order delay differential equation.From Fig. 8, it can be seen that the approximate solution obtained by our method coincides with the exact solution for all the used values of the delay τ .The absolute error of the presented method is plotted in Fig. 9.This figure shows that the proposed method is accurate because the error is at most of order 10 -13 .Table 3 presents the absolute errors for different values of m.It is   obvious that the absolute errors decrease when the value of m increases, which demonstrates the convergence of the introduced method to the exact solution.Table 4 compares the absolute errors obtained by the introduced method with the collected results in [42].
From Table 4, we can note that the presented technique is more accurate and efficient in comparison to the methods [33,42].
Example 7 Consider the following fractional delay differential equation [34,42,50]: The exact solution for this problem is y(t) = sin(t) when ν = 1.The approximate solution obtained for different values of ν together with the known exact solution at ν = 1 is plotted in Fig. 10.It can be noted that as ν tends to 1, the approximate solution converges to the exact solution for ν = 1.The absolute error function for ν = 1 is presented in Fig. 11.This figure shows that even if the exact solution is not a polynomial, an accurate approximation can be achieved using a sufficient number of FCHFs.In Table 5, we compare the absolute errors of the present method with those collected in [50] for the Bernoulli wavelet method [42], the spectral method based on a modification of hat functions [34] and the Legendre wavelet method [50].These results show that the proposed method is superior to the methods in [34,42,50] for most of the listed values of absolute errors by using a smaller number of basis functions than they have been used in [34,42,50].This example shows that the present method provides high efficiency regardless of the type of the exact solution for the problem.
The exact solution for this problem is y(t) = cos(t) when ν = 2. Figure 12 shows the obtained approximate solution for various choices of ν with the exact solution for ν = 2.It can be seen that as ν converges to 2, the approximate solution approaches the exact solution for ν = 2. Figure 13 displays the absolute error function for ν = 2.This graph shows that the presented method has a good convergence rate and its solutions are accurate.Table 6 presents a comparison between the obtained results and those reported in [42] in terms of absolute errors.It can be noted that our method outperforms the method in [42] since the absolute errors achieved by the introduced method are less than those in [42] for all corresponding cases by using a fewer number of basis functions for each case.This demonstrates the agreement of our method with the exact solution with reasonable computational costs.
Example 9 Consider the following fractional delay differential equation [44]: y(0) = 1.(6.31)This problem has the exact solution y(t) = e t when ν = 1.The approximate solution for different values of ν is plotted in Fig. 14 along with the exact solution for ν = 1.The absolute  error of the presented method is drawn in Fig. 15.These graphs demonstrate the convergence of the obtained approximate solutions to the exact solution.Table 7 compares the results obtained by the present method with the those in [44].It can be seen that our results are in agreement with the results in [44] for m = 8 and slightly better for m = 9.This means that the presented method is compatible with the method in [44] for this problem.
For all tested examples, the presented method has found the solutions with high accuracy within acceptable computation costs.We can observe that using fractional values of β leads to a good rate of convergence and a low number of basis functions is required to obtain satisfactory results when the solutions of fractional delay differential equations contain terms with fractional powers or when the fractional derivative ν is a proper fraction.

Conclusion
This paper has used FCHFs and their properties to provide an accurate numerical method for solving fractional delay differential equations.The fractional derivative was considered in the Caputo sense.The problem of finding a solution to a fractional delay differential equation was transformed into the problem of solving a system of algebraic equations.The convergence of the presented approximation has been discussed.The introduced method has been applied to solve many problems and the obtained results were compared with exact solutions and some other published numerical methods.It can be observed that the presented method has high accuracy by using a smaller number of basis functions than other methods used in the comparisons.This reveals that the presented method is efficient, accurate and its performance is quite satisfactory, which establishes the significance of the present method for solving fractional delay differential equations.

Figure 3 3 Figure 4
Figure 3 The graph of the approximate solution for various choices of β at m = 10 and the exact solution for Example 3

Figure 5 Figure 6
Figure 5 The graph of the approximate solution for various values of ν at m = 10, β = 1 and the exact solution for Example 4

Figure 7 Figure 8 6 Figure 9
Figure 7 The graph of the approximate solution for different values of ν at m = 10, β = ν, τ = 0.2 and the exact solution for Example 6

Figure 10 Figure 11
Figure 10 The graph of the approximate solution for several choices of ν at m = 10, β = ν and the exact solution for Example 7

Figure 12 Figure 13
Figure 12 The graph of the approximate solution for various values of ν at m = 10, β = 1 and the exact solution for Example 8

Figure 14 9 Figure 15
Figure 14 The graph of the approximate solution for various values of ν at m = 9, β = ν and the exact solution for Example 9

Table 2
Comparison of the absolute errors at β = 1 and ν = 3 for Example 4

Table 3
The absolute errors for different values of m at β = ν = 1 and τ = 0.2 for Example 6

Table 4
Comparison of the absolute errors at β = ν = 1 and τ = 0.2 for Example 6

Table 5
Comparison of the absolute errors at β = ν = 1 for Example 7

Table 6
Comparison of the absolute errors at β = 1 and ν = 2 for Example 8

Table 7
Comparison of the absolute errors at β = ν = 1 for Example 9