Numerical Simulation of Fractional Delay Differential Equations Using the Operational Matrix of Fractional Integration for Fractional-Order Taylor Basis

: In this paper, we consider numerical solutions for a general form of fractional delay differential equations (FDDEs) with fractional derivatives deﬁned in the Caputo sense. A fractional integration operational matrix, created using a fractional Taylor basis, is applied to solve these FDDEs. The main characteristic of this approach is, by utilizing the operational matrix of fractional integration, to reduce the given differential equation to a set of algebraic equations with unknown coefﬁcients. This equation system can be solved efﬁciently using a computer algorithm. A bound on the error for the best approximation and fractional integration are also given. Several examples are given to illustrate the validity and applicability of the technique. The efﬁciency of the presented method is revealed by comparing results with some existing solutions, the ﬁndings of some other approaches from the literature and by plotting absolute error ﬁgures.


Introduction
Fractional calculus is a field of mathematical analysis which might be represented as the generalization of integration and differentiation to arbitrary orders.The history of this theory started with the speculations of G.W. Leibniz (1695) and L. Euler (1730).Despite their long history, fractional calculus and the corresponding fractional differential equations (FDEs) are only recently garnering attention and popularity.There are several different definitions of fractional derivatives in the literature, including Riemann-Liouville, Caputo, Grünwald-Letnikov, Weyl, Marchaud, and Prabhakar, to mention just a few of the well-known ones.The history of this topic can be found in [1][2][3][4].
In point of fact, fractional calculus provides mathematical modeling of various phenomena, sometimes in a stronger way than ordinary calculus.When compared with classical calculus, fractional calculus has memory effects which are useful for modeling long-time interactions.Due to this, it can illustrate many different kinds of dynamical and engineering models in a more accurate way.Over the last few decades, applications were investigated in numerous fields of science and engineering such as bioengineering [5], biology [6], chaotic systems [7], chemistry [8], control theory [9,10], economics [11], electrochemistry [12], finance [13], mechanics [14], optics [15], physics [16], rheology [17], social sciences [18], viscoelasticity [19], and so forth.
The importance and useful properties of FDEs leads them and their solution methods to attract widespread interest.Moreover, since analytical solutions such as those using matrix Mittag-Leffler functions [20] or Laplace transforms [21] are often impossible or impractical for more advanced FDEs, numerical methods are crucial for solving FDEs in practice.A variety of such methods have been established in the literature, such as the q-homotopy analysis transform method [22], the B-spline collocation method [23], the predictor-corrector method [24], the space-spectral time-fractional Adams-Bashforth-Moulton method [25], the fractional Taylor operational matrix method [26], the Bernoulli polynomials method [27], the differential transform method [28], and so on.Spectral methods are efficient numerical techniques used to solve differential equations (DEs) in applied mathematics.In 1938, spectral methods were introduced by Lanczos [29] by showing the powerful role of Fourier series and Chebyshev polynomials for solutions of some problems.Spectral methods have received considerable interest in recent years, solving many different types of integral and differential equations numerically, because of their easy applicability over both finite and infinite intervals.When the solution of a given problem is smooth, spectral methods have very good error properties, namely the so-called "exponential convergence".These highly accurate methods are based on expressing the approximate solution of a DE as a linear combination of a chosen set of orthogonal basis functions and choosing the coefficients in the sum in order to satisfy the solution of the DE [30].In general, there are three types of such methods; collocation, Tau, and Galerkin.In our work, for discretization of FDDEs, we will use a spectral collocation method with a fractional Taylor basis to approximate the functions.Collocation is based on interpolation; similarly to the finite difference method, the collocation spectral method uses collocation points, namely a set of grid points in the domain.
A delay differential equation (DDE) is a special kind of differential equation, which lies somewhere between ordinary differential equations (ODEs) and partial differential equations (PDEs); the derivative of the unknown function at a certain time is given in terms of the values of the function at previous times.Mathematical models that utilize DDEs turn out to be beneficial, particularly when the description of the investigated systems depends not only on the status of a system in that moment, but also in the past.There are different types of DDEs with different types of delays, such as constant delay, time-dependent delay, state-dependent delay, neutral delay, stochastic delay, and so on.DDEs play an important role in various applied sciences such as biology, electrical networks, environment science, life science, mechanics, physiology, and population dynamics; see, for instance, [31][32][33] and references therein.
Even though delays can be recognised everywhere in real world systems and there has been widespread interest in the study of DDEs for many years, FDDEs are a very recent topic.The FDDE is a generalization of the DDE to arbitrary non-integer order.Since a small delay can have a great impact, FDDEs have gained considerable interest for modeling some real world problems.FDDEs have many applications in different scientific disciplines such as biology, control theory, chemistry, economics, electrodynamics, finance, vibration theory, and so on [34][35][36].For pure mathematical results, we refer to [37] for existence and uniqueness, Ref. [38][39][40] for stability, and [41] for controllability.For most FDDEs, exact solutions are not known.Therefore, various numerical techniques such as the Legendre wavelet method [42], the Chebyshev operational matrix method [43], the Bernoulli wavelets method [44], the shifted Gegenbauer-Gauss collocation method [45], the Haar wavelet collocation method [46], the modified operational matrix method [47], the discontinuous Galerkin method [48], spectral methods [49], etc., have been developed and applied to provide approximate solutions.
The classical linear DDE is given as where ω is the solution to be determined, f is a known analytic function, and c, τ, ct − τ denote a constant, delay, and delay argument, respectively, with 0 < c ≤ 1.
In this paper, motivated by the powerful properties of fractional calculus, especially the non-local nature of fractional operators, we replace the integer order derivative of DDE (1) with a non-integer order Caputo derivative, and focus on numerical solutions for the following general form of FDDE which involves more than one fractional differential operator: with initial conditions ω r (0) = λ r , r = {0, 1, ..., µ − 1}, and ω(t) = g(t), t < 0.
where D µ , D β , D γ represent the Caputo fractional derivatives of order µ, β, γ, respectively, with µ ≥ β and γ > 0, ω is the solution to be determined, and f and g are known analytic functions.
The main purpose of this work is to find an approximate solution of the problem given in ( 2) and (3) using a spectral method with non-integer order Taylor basis.In order to realize this aim, we derive the operational matrix of fractional integration by using a fractional order Taylor basis.Then, by utilizing this matrix, we obtain a system of linear algebraic equations.The solution of this system for unknown coefficients will give an approximate solution for the FDDE.
Beside the advantages of using a Taylor basis, like finding fractional derivatives and approximating the functions easily, it has also led us to create an operational matrix without any approximation, that is, approximating the unknown function is the only approximation involved in this method.Obviously, in addition to the important properties of spectral methods mentioned above, using a Taylor basis is important in order to obtain better results.In general, if the solution of a given problem is smooth and analytic on the defined interval, spectral methods will obtain highly accurate results.However, these efficient methods become less accurate when applied to problems with discontinuous coefficients and complex geometries.Moreover, the accuracy of derivatives acquired by term-by-term differentiation of a truncated expansion naturally deteriorates [50].
The structure of this paper is as follows.Section 2 comprises a brief introduction to the key definitions and properties of fractional calculus.Section 3 introduces an operational matrix for fractional integration, constructed by the use of fractional Taylor vectors.In Section 4, we propose a numerical method to solve the given FDDE.In Section 5, we give an error bound for the best approximation and fractional integration.In Section 6, the proposed scheme is applied to six examples to indicate its applicability.Conclusions and suggestions for future work are presented in Section 7. A numerical algorithm for the present method is provided in Appendix A.

Fractional Calculus
Definition 1.Let ω(t) be an integrable function.The Riemann-Liouville (R-L) fractional integral to order µ of ω(t) is given as When employed to a power function, the R-L fractional integral gives the following result: The R-L fractional integral operator has a semigroup property, that is: and the operator is linear, namely for any two functions ω 1 , ω 2 and constants C 1 , C 2 .
Definition 2. The Caputo fractional derivative of ω(t) to order µ is given as The Newton-Leibniz identity yields a crucial relation among the R-L fractional integral and the Caputo fractional derivative as follows: The Taylor vector with non-integer order is represented as [51], for m ∈ N and σ > 0.

Function Approximation by Fractional Taylor Vector
Therefore, the function ω can be expanded in terms of the fractional Taylor vector as: Truncating the series in Equation ( 9) to a finite number of terms, we get where the vector of coefficients C is given by

Construction of Operational Matrix of Fractional Integration
In this part, we create the operational matrix of fractional integration by use of the fractional Taylor vector.
Let us define the R-L fractional integration to order µ of the Taylor vector given in (8) as where M (t,µ) is the operational matrix of fractional integration of order µ with (m + 1) × (m + 1) dimensions.By using feature (5) of the R-L fractional integral and the vector (8), we get Γ(2σ+µ+1) t 2σ+µ , ..., Therefore, we can write (13) as where .
If we define Θ (t,µ) = t µ P µ , then the fractional Taylor operational matrix of fractional integration can be expressed as

Method of Solution
In this section, we define the proposed technique for the construction of a numerical solution for a general form of FDDEs.The numerical algorithm, based on the fractional Taylor operational matrix of fractional integration, is given below.
Let us recall and focus on the fractional order delay differential equations given in ( 2) and ( 3): with initial conditions ω r (0 and We start by expanding D µ ω(t) by using a fractional Taylor basis vector as follows: As a next step, by applying the R-L fractional integral of order µ on (18) and using initial conditions (17), we get By using the fractional Taylor operational matrix of integration (15), we can rewrite (19) as Further, applying a differential operator of order β on both sides of ( 19), we have and therefore Similarly, we can approximate the delay term as Now, by taking a differential operator of order γ on both sides of ( 23), we get Hence, by substituting Equations ( 18), ( 20) and ( 22)-( 24) into Equation ( 16), we get an algebraic equation with unknown coefficient C T .
As a final step, by using the collocation points t r = r/m where r = 0, 1, . . ., m, we get a system of m + 1 linear algebraic equations with a vector of unknown constants C T = [c 0 , c 1 , c 2 , . . ., c m ].This system might be solved efficiently for the vector of coefficients C T .

Error Estimation
To obtain an upper bound for the best approximation, we use the generalized Taylor's formula.
In case µ = 1, the generalized Taylor's formula reduces to the classical one.

Error Bound for the Best Approximation
In this section, we obtain an error bound for the best approximation by the following theorem.
If ω m = C T T mσ is the best approximation of y from B, then where N µ ≥ D (m+1)µ ω(t) .
Proof of Theorem 2. Since ω N is the best approximation of ω from B, and by use of the generalized Taylor's formula, we get Now, by taking the square root of Equation ( 25), the desired result can be obtained.

Error Bound for Fractional Integration
The error bound for fractional integration can be obtained by the following theorem.
Theorem 3 ([51]).Let all conditions of Theorem 2 be satisfied and µ > 1.Then Proof of Theorem 3. By using the definition of the R-L fractional integral of order µ, we have With use of ( 25) and ( 26), the desired result can be obtained.
It is noted that the error investigation for the solution of the FDDE might be reduced to the error analysis of the function approximation and integration which are presented in Theorems 2 and 3, respectively.

Illustrative Examples
In this part, we present six examples to illustrate the applicability and accuracy of the introduced scheme.In each example, we apply the proposed method to solve given problems and compare the obtained results with analytical solutions and the results obtained by other methods from the literature.The acquired results reveal that the present method is very efficacious for solving FDDEs.

Example 1
In this example, the following generalized form of FDDE is considered [53]: with the initial conditions where t > 0, 0 < µ ≤ 3. The exact solution of ( 27) is revealed by ω(t) = e −t in the case when f (t) = e −t/2+0.3 and µ = 3.
In Table 1, we present the comparison of the absolute errors of the present method with that of the shifted Chebyshev polynomials method (SCPM) [53] for m = 6, 8.  1 indicates that the solutions derived by using the present approach are closer to the exact solution when compared with the method of [53].
In Figure 1, we present the numerical solution results for the problem ( 27) when µ = 2.65, 2.75, 2.85, 2.95, 3.0.By comparing these findings with the exact solution, we can state that as µ approaches 3, the approximate solutions converge to those of the integerorder DE.
In Table 2, we compare the absolute errors obtained by using our method for m = 7, σ = 1, with those results obtained using the Chebyshev polynomial method (CPM) [54], spectral collocation method (SCM) [43], and Bernoulli wavelet method (BWM) [44], where τ = 0.01.Table 2.The comparison of the absolute error results between our method and the CPM [54], SCM [43], and BWM [44] for τ = 0.01.In Figure 3, we graphically illustrate the comparison of absolute errors given in Table 2.In Figure 4, we illustrate a comparison between the exact solution and the approximate solution acquired by the proposed technique for the constant delay τ = 0.01.In Figure 5, we graphically illustrate the acquired absolute error by using the given technique when τ = 0.01.In Table 3, we present a comparison of the acquired absolute errors between the present method for m = 5 and the method of [44] when τ = 0.001.These results show that the proposed method is very efficient even for a small number of basis elements.The graphical representations of the exact solution and the approximate solution obtained using the given method when τ = 0.001 are demonstrated in Figure 6.In Figure 7, we graphically illustrate the absolute error acquired by the present scheme.The approximate solutions of the given problem (28) with constant delays τ = 0.001e −t , 0.002e −t , 0.01e −t , 0.02e −t , 0.03e −t are presented graphically in Figure 8.

Example 3
Consider an FDDE with constant delay as follows [46]: with The analytical solution is given by ω(t) = t 2 .We solve this problem by taking σ = 1/2.In Table 4, we present the obtained results for the absolute error and relative error when m = 4.A comparison of maximum absolute errors obtained by the proposed method for m = 4, the Haar wavelet collocation technique (HWCT) [46] for N = 2048, and the fractional backward difference method (FBDM) [55] for N = 640, are presented in Table 5.This results show the efficiency of the given technique even for a small m value.
Table 5.Comparison of the maximum absolute errors between the present method, the HWCT [46], and the FBDM [55].The graphical demonstration of the exact solution and the behaviour of the approximate solution obtained by the given numerical technique are presented in Figure 9.The absolute error results are also presented in Figure 10.

Example 4
In this example, the following form of fractional order differential equation with a proportional delay is considered [47]: with ω(0) = 0. where The exact solution is given by ω(t) = t + µt 3 2 .We solve the given problem (31) for µ = 1/2 and by choosing σ = 1/2, m = 4.A comparison between the acquired absolute errors of the present technique and the modified operational matrix method (MOMM) [47] are presented in Table 6.From these findings, we can conclude that the present method is very efficient even for small numbers m.In Figure 11, we illustrate the exact solution and the approximate solution acquired by using the given approach for the problem (31).The behaviour of the acquired absolute error for the problem (31) when we use the proposed scheme is presented in Figure 12.

Example 5
In this example, the following form of fractional order neutral functional-differential equation with proportional delay is considered [45]: with the initial conditions where t ∈ [0, R] and The exact solution of this problem is revealed by ω(t) = t 4 .
In order to solve this fractional differential equation with proportional delay, we use the given method with σ = 1 and m = 4.The obtained results for the absolute errors are presented and compared with the shifted Gegenbauer-Gauss collocation method (SGGC) [45] in Table 7.
The approximate solution obtained by using our technique and the exact solution of the given problem are graphically presented in Figure 13.Additionally, illustration of the absolute errors can be seen in Figure 14 where t ∈ [0, 1].In Figures 15 and 16, we illustrate the behaviours of the approximate solution and the absolute errors acquired by use of the present method, respectively, where t ∈ [0, 7].From these results, we can see that the given method provides highly accurate results also for wide intervals of t even for a small number of basis elements.

Example 6
As a last example, we consider the following form of FDDE with proportional delay [44]: The exact solution of this problem is given by ω(t) = sin(t) when µ = 1 for any a, b ∈ R and 0 < c < 1.
In Table 8, we give comparison results for the acquired absolute errors of the approximate solution between the proposed method, the Bernoulli wavelet method [44], the spectral method based on a modification of hat functions [49], and the Legendre wavelet method [42] for c = 1/2.From Table 8, we can see that the present method provides higher efficiency even for a small number of basis elements from the fractional Taylor vector.
The illustration of the approximate solution and the absolute error obtained by our method are presented in Figure 17 and Figure 18, respectively.From Figure 19, we can see the graphical representations of the numerical solutions acquired by using the proposed scheme when c = 1/10 and t ∈ [0, 5].Here, we can conclude that even for a wide range of t, we can obtain highly accurate results by use of the present method.
Figure 22 presents a graphical representation of the comparison between the exact and approximate solutions provided by use of the present method by taking m = 3, 4, 5, 6 when t ∈ [0, 8].From these illustrations and results, we can conclude that the given technique rapidly provides efficient results.

Conclusions and Future Work
In this study, a numerical technique based on an operational matrix of fractional integration is used to solve FDDEs.The main purpose of this approach is to reduce the given FDDE to algebraic equations, which simplifies the given problem.Hence, the numerical solution is obtained by solving this set of algebraic equations.Concerning the difficulty of solving FDDEs analytically, the proposed technique can have important contributions to the field of numerical methods by having very efficient results and applicability, even for a very low number of basis elements.The results provided by using the proposed technique have been shown and compared with some other techniques from the literature in Section 6.From these results, we can conclude that the given method has wide applicability for the solutions of many different kinds of FDDEs.We also observe that the choice of Taylor basis performs very well in approximating the unknown function.Another advantage is that the operational matrix of fractional integration M (t,µ) is a sparse matrix; that is, most of the elements are zero which shortens the required CPU time.All computational results were calculated using MATLAB version R2019A software on an Intel(R) Core(TM) i5-6500 CPU, 3.20 GHz computer with 6.00 GB RAM.For future work, this technique can be extended to apply to other DDEs such as PDEs with delays and variable-order DDEs.Furthermore, other orthogonal polynomials can be used for the approximation of unknown state functions.

Figure 1 .
Figure 1.The comparison between the exact solution and the approximate solutions obtained for the various values of µ.

Figure 2
Figure 2 illustrates the graphical representation of the absolute error acquired by using the present method for parameters m = 8 where µ = 3.

Figure 3 .
Figure 3.The graphical comparison of absolute errors acquired by our method, CPM, SCM and BWM, respectively.

Figure 4 .
Figure 4.The behaviour of approximate and exact solutions for Example 2 with τ = 0.01.

Figure 6 .
Figure 6.The behaviour of approximate and the exact solutions for Example 2 with τ = 0.001.

Figure 7 .
Figure 7.The absolute error of the approximate solution for Example 2 with τ = 0.001 .

Figure 8 .
Figure 8.The behaviour of the approximate solutions for Example 2 with different τ values.

Figure 9 .
Figure 9.The behaviour of approximate and the exact solutions for Example 3.

Figure 10 .
Figure 10.The behaviour of absolute error for Example 3.

Figure 11 .
Figure 11.The compared results of the approximate and exact solutions for Example 4.

Figure 12 .
Figure 12.The acquired absolute error of the approximate solution for Example 4.

Figure 13 .
Figure 13.The demonstration of the approximate and exact solutions for Example 5.

Figure 14 .
Figure 14.The demonstration of the absolute error for Example 5.

Figure 15 .
Figure 15.The compared results of the approximate and exact solutions for Example 5 where t ∈ [0, 7].

Figure 17 .Figure 18 .
Figure 17.The compared results of the approximate and exact solutions for Example 6 in the interval t ∈ [0, 1]

Figure 19 .
Figure 19.The compared results of the approximate and exact solutions for Example 6 when c = 1/10 in the interval t ∈ [0, 5].

Figures 20 and 21
Figures 20 and 21 graphically present the results for the approximate solutions and the absolute errors provided by use of our method, respectively, where t ∈ [0, 8].

Figure 20 .
Figure 20.The demonstration of the approximate and exact solutions for Example 6 in the interval t ∈ [0, 8].

Figure 22 .
Figure 22.The demonstration of the exact and approximate solutions for Example 6 with different m values where t ∈ [0, 8].

Table 3 .
[44]absolute error results obtained by the proposed method and the method of[44]where τ = 0.001.

Table 6 .
[47]arison of the absolute errors between the present method and the MOMM[47].

Table 7 .
[45]compared results for the acquired absolute errors between the present method and the SGGC[45].

Table 8 .
[42,44,49]s for the absolute errors of the approximate solutions provided by the present method and the methods of[42,44,49].