1. 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,
,
denote a constant, delay, and delay argument, respectively, with
.
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
and
where
,
,
represent the Caputo fractional derivatives of order
, respectively, with
and
,
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.
4. 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
and
We start by expanding
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
.
As a final step, by using the collocation points where , we get a system of linear algebraic equations with a vector of unknown constants . This system might be solved efficiently for the vector of coefficients .