A Novel Operational Matrix of Caputo Fractional Derivatives of Fibonacci Polynomials: Spectral Solutions of Fractional Differential Equations

Herein, two numerical algorithms for solving some linear and nonlinear fractional-order differential equations are presented and analyzed. For this purpose, a novel operational matrix of fractional-order derivatives of Fibonacci polynomials was constructed and employed along with the application of the tau and collocation spectral methods. The convergence and error analysis of the suggested Fibonacci expansion were carefully investigated. Some numerical examples with comparisons are presented to ensure the efficiency, applicability and high accuracy of the proposed algorithms. Two accurate semi-analytic polynomial solutions for linear and nonlinear fractional differential equations are the result.


Introduction
Fractional calculus is a vital branch of mathematical analysis which deals with derivatives and integrals to an arbitrary order (real or complex).There are intensive studies of fractional calculus from both theoretical and practical points of view.The applications of fractional calculus in many fields such as physics, economics, probability theory and statistics are numerous.For example, many physical and engineering models are described by fractional differential equations.For some applications of fractional calculus, one can refer to [1].
Due to the great importance of fractional differential equations (FDEs) and their applications in various disciplines, many researchers have shown considerable interest in developing numerical algorithms for solving these types of equations.Numerical solutions for various types of FDEs by applying several techniques were proposed, for instance, Adomian's decomposition method [2,3], the Taylor collocation method [4], the variational iteration method [5], the finite difference method [6,7] and the ultraspherical wavelets method [8,9].In addition, orthogonal polynomials have been widely used for obtaining numerical solutions for different types of FDEs.For example, Chebyshev polynomials of the second kind were employed by Sweilam et al. for solving space fractional-order diffusion equations [10].Moreover, Chebyshev polynomials of the first kind were used for solving some fractional differential equations in Khader and Sweilam [11].
The three well-known versions of spectral methods; namely collocation, tau and Galerkin methods play pivotal roles in the numerical solutions of various types of differential equations.The spectral solutions are expressed in terms of truncated series of various polynomials.Actually, this solution is written as ∑a k φ k , where {φ k } is a set of polynomials, which is often orthogonal.The expansion coefficients a k can be determined by applying an appropriate spectral method.The collocation approach is performed by enforcing the differential equation to be satisfied exactly at some selected points called nodes or collocation points.The tau method is a synonym for expanding the residual function as a series of orthogonal polynomials and then, applying the boundary conditions as constraints.The Galerkin method is basically based on selecting some combinations satisfying the underlying boundary (initial) conditions as basis functions, and then, enforcing the residual to be orthogonal with the basis functions.This approach is fruitfully applied in many studies to solve two point linear boundary value problems (BVPs), in one and two dimensions [12][13][14][15].
The utilization of operational matrices of derivatives of different orthogonal polynomials is very effective for treating almost all types of differential equations.To be more precise, several operational matrices are employed for obtaining numerical solutions of ordinary differential equations.In this regard, Abd-Elhameed in [16] has developed a novel harmonic numbers operational matrix of derivatives to solve linear and nonlinear sixth-order BVPs.Likewise, Abd-Elhameed [17] and Doha et al. [18] used tau and Galerkin operational matrices of derivatives, respectively, to solve the singular Lane-Emden differential equations.The approach of employing the operational matrices is not limited to application in ordinary differential equations, but it can be also followed to handle FDEs.In this respect, a number of articles were devoted to introducing numerical solutions for these equations.In this respect, a Legendre operational matrix of fractional derivatives were constructed and employed for solving some types of FDEs in [19].As well, Rostamy et al. [20] handled these equations with the aid of a Bernstein operational matrix of fractional derivatives.Other articles were interested in using the operational matrices of fractional integration to solve FDEs, i.e., Kazem [21].
The Fibonacci polynomials with their generalizations, and Fibonacci numbers and the golden ratio, are key elements in number theory.They have enormous applications in various disciplines such as computer science, statistics, physics, biology and graph theory [22].Several studies were peformed to theoretically discuss Fibonacci polynomials and generalized Fibonacci polynomials (see [23][24][25]).However, as far we know, articles interested in using these polynomials in different applications are scarce.For example, a collocation procedure for treating BVPs based on using the Fibonacci operational matrix of derivatives has been developed in [26,27].
The principal objective of the present article is to present and implement new numerical spectral solutions of some types of FDEs.For this purpose, the operational matrix of fractional derivatives of Fibonacci polynomials are constructed, and then employed along with the tau and collocation spectral methods for obtaining numerical solutions of FDEs.
The paper is organized as follows: first, in Section 2, some preliminaries including some fundamental definitions of the fractional calculus theory are presented.This section also includes an overview on Fibonacci polynomials and their relevant properties.Section 3 is interested in constructing Fibonacci operational matrix of the fractional derivatives in the Caputo sense.Section 4 is devoted to implementing and presenting two numerical algorithms for solving two-term FDEs.Section 5 is interested in investigating the convergence and error analysis of the suggested expansion of Fibonacci polynomials.Some numerical examples and discussions are presented in Section 6, hoping to illustrate the efficiency, simplicity and applicability of the two proposed algorithms.Finally, some conclusions are reported in Section 7.

Some Definitions and Properties of Fractional Calculus
In this subsection, some notations, definitions and preliminary facts about fractional calculus theory are presented.Definition 1.The Riemann-Liouville fractional integral operator I µ of order µ on the usual Lebesgue space L 1 [0,1] is defined as: where the following properties are satisfied: (i) where µ > −1 and v > −1.
Definition 2. The Riemann-Liouville fractional derivative of order µ > 0 is defined by: Definition 3. The fractional differential operator in Caputo sense is defined as: The operator D µ satisfies the following basic properties for n − 1 < µ < n: where the ceiling notation µ denotes the smallest integer greater than or equal to µ.For more properties of fractional derivatives and integrals, see for example [28,29].

Some Properties of Fibonacci Polynomials
It is well-known that the Fibonacci polynomials may be constructed by the following recurrence relation with the initial values: These polynomials have the following explicit form: Moreover, they can be represented by the following power form: where the notation z represents the largest integer less than or equal to z .
The first derivative of Fibonacci polynomials is related with the original Fibonacci polynomials by the following relation (see [30]): Furthermore, the polynomials x m can be written as a combination of Fibonacci polynomials as (see [30]): It is usefull to rewrite relations ( 2) and ( 3) as: respectively.

Construction of the Fibonacci Operational Matrix of the Caputo Fractional Derivative
Assume that u(x) is a square Lebesgue integrable function on (0,1).In addition, assume that u(x) can be expanded as: If we kept only the first (n + 1) terms of Fibonacci polynomials, then: where and Now, the first derivative of the vector Φ(x) can be written as: where is the Fibonacci operational matrix of derivatives of order (n + 1) × (n + 1).
Moreover, the elements of M (1) are given explicitly as: For example, for n = 6, the operational matrix M (1) is:

Introducing a Fibonacci Operational Matrix of Fractional Derivatives
The principal goal of this section is to introduce a generalization of the operational matrix of the integer case.First note that relation (7) leads to the relation: where s is any postive integer.The following theorem gives the Fibonacci operational matrix of fractional derivatives: is the Fibonacci polynomial vectors which defined in Equation ( 6), then for any α > 0, the following relation holds: where M (α) = m α ij is the (n + 1) × (n + 1) Fibonacci operational matrix of derivatives of order α in Caputo sense and it is given explicitly as: In fact, the elements m α ij are given explicitly by: where Proof.If we apply D α to Equation (4), then with the aid of relation (1), the following relation is obtained: The inversion formula (5) enables one-after performing some manipulations-to write: where ξ α (i,j) is given in (8).Now Equation ( 9) can be rewritten in the following vector form: and also we can write: Finally, Equation ( 10) together with Equation ( 11) yield the desired result.

Two New Matrix Algorithms for Solving Fractional-Order Differential Equations
This section is devoted to analyzing and presenting two numerical algorithms for solving fractional-order linear and nonlinear differential equations based on using the constructed Fibonacci operational matrix.The two algorithms, namely, the tau Fibonacci matrix method (TFMM) and the collocation Fibonacci matrix method (CFMM) are employed for solving linear and nonlinear multi-term orders fractional differential equations, respectively.

Use of TFMM for Handling Linear Fractional Differential Equations
Consider the linear fractional differential equation ( [31,32]): (12) subject to the initial conditions: or the boundary conditions: Now, assume that u(x) can be approximated as: In virtue of Theorem 1, D α u(x) and D β u(x) can be approximated as: and Making use of the approximations in ( 15)-( 17), the residual of ( 12) can be calculated by the formula: and accordingly the application of tau method (see [18]) yields: Moreover, the initial conditions (13) yield: while the boundary conditions ( 14) yield: From Equation ( 18) with ( 19) or ( 20), a linear system of equations in the unknown expansion coefficients c i of dimension (n + 1) can be generated.Thanks to the Gaussian elimination technique or any suitable technique, this system can be solved, and consequently, the approximate solution (15) can be found.

Use of CFMM for Handling Nonlinear Fractional Differential Equations
Consider the following nonlinear fractional-order differential equation: subject to the initial conditions (13) or the boundary conditions (14).We approximate u(x), D α u(x) and D β u(x) as in the previous section, hence the residual R (x) of Equation ( 21) is given by: The application of the spectral collocation method requires that R (x) must vanish at certain collocation points.We select the collocation points to be: i n+1 , i = 1, 2, . . ., n − q 1 , and therefore: Equations ( 22) with ( 19) or ( 20) generate a nonlinear system of equations in the unknown expansion coefficients c i of dimension (n + 1).Newton's iterative technique can be employed for solving this system, and consequently, the approximate solution is obtained from (15).

Remark 1.
It is worthy to mention here that the CFMM, which described above, can be also applied to linear fractional differential equations.

Convergence and Error Analysis
In this section, a comprehensive study of the convergence and error analysis of the suggested expansion of Fibonacci polynomials is performed.For this purpose, the following Lemmas are needed.

Lemma 1 (Byrd [33]
).If f(x) is an infinitely differentiable function at the origin, then f(x) can be expanded in terms of Fibonacci polynomials as: where

Lemma 2 (Abramowitz et al. [34]
).Let I v (x) denote the modified Bessel function of order v of the first kind.
The following identity holds: Lemma 3 (Luke [35]).The modified Bessel function of the first kind I v (x) satisfies the following inequality: denotes the golden ratio, then the following inequality is valid: Proof.We have: The following two theorems discuss the convergence and error analysis of the Fibonacci expansion.
, then we have: The series converges absolutely.
Proof.Starting from Lemma 1, we have: and therefore the application of Lemma 2 implies that: The inequality in (23) along with Lemma 3 leads to the inequality in the first part.Now, to prove that the series ∑ ∞ k=1 c k F k (x) converges absolutely, the comparison test is applied.Making use of the first part along with Lemma 4, yields = coshMe Mσ/2 and therefore the series converges absolutely.
Theorem 3. If f(x) satisfies the hypothesis of Theorem 2, and e n (x) = ∑ ∞ k=n+2 c k F k (x) then we have the following error estimate: Proof.By Theorem 2, we can write: which can be written as: where Γ(n + 1) and Γ(n + 1,µ) denote the well-known gamma function and the incomplete gamma function, respectively (see [36]).which completes the proof of the theorem.

Numerical Tests
This section presents some numerical experiments and comparisons to illustrate the efficiency and applicability of the two proposed algorithms in this paper.

Example 1 (Doha et al. [37]). Consider the following linear fractional initial value problem
We apply TFMM for the case n = 2.The residual of Equation ( 24) can be calculated by the formula: where the operational matrices M (2) and M ( 12 ) are given explicitly as follows: If we apply the method which described in Section 4.1 to Equation ( 24), then we get: Moreover, the initial conditions (25) yield: Equations ( 26) and ( 27) can be immediately solved to give c 1 = −1, c 2 =0, c 3 = 1, and consequently u(x) = x 2 which is the exact solution.

Example 2 (Doha et al. [38]
).Consider the following inhomogeneous Bagley-Trovik equation: where f(x) is chosen such that the exact solution of Equation ( 28) is u(x) = sin(αx).We apply the two methods, namely, TFMM, CFMM for different values of α and n.In Table 1, we display a comparison between the results obtained if the two methods TFMM, CFMM are applied with the results obtained by applying Chebyshev spectral method (CSM) which developed in [38].The displayed results in this table show that our algorithm gives a lesser error in almost all cases.
The exact solution of (29) in the case corresponds to q 1 = 2 and q 2 = 1 is u(x) = x(1 − e x-1 ).In Table 2, we compare our results with those obtained in [32,39].Figure 1 illustrates that the approximate solutions in the case corresponds to q 2 = 2 and for various values of q 1 near the value 1, have similar behavior.

Example 4 ([32,39]). Consider the following nonlinear boundary value problem:
(1,2), (0,1), ), ( sin subject to the homogenous boundary conditions: The exact solution of Equation ( 30) for the case q = 2 is u(x) = sin 2 (πx) We apply CFMM.Table 3 subject to the homogenous boundary conditions: The exact solution of Equation ( 30) for the case q = 2 is u(x) = sin 2 (πx) We apply CFMM.Table 3 lists the maximum pointwise error of Equation ( 30) for the case q = 2 and different values of n. Figure 2 illustrates that the behavior of the approximate solution for values of q near 2 is somehow close to the behavior of the exact solution.Example 4 ([32,39]).Consider the following nonlinear boundary value problem: (1,2), (0,1), ), ( sin subject to the homogenous boundary conditions: The exact solution of Equation ( 30) for the case q = 2 is u(x) = sin 2 (πx) We apply CFMM.Table 3 lists the maximum pointwise error of Equation ( 30) for the case q = 2 and different values of n. Figure 2 illustrates that the behavior of the approximate solution for values of q near 2 is somehow close to the behavior of the exact solution.

Example 5 ([40]
).Consider the fractional Van der Pol equation with fractional damping: subject to the homogenous initial conditions: We solve Equation (31) for the case corresponds to a = 1.31, ω = 0.5.In Table 4, we compare our results with the numerical solution obtained in [40] and the solution obtained from Mathematica 9 by applying the fourth-order Runge-Kutta method.In Table 5, we compare our results with the results obtained in [40] in case with no damping η = 0.  5 indicate that the numerical spectral solutions obtained by CFMM are in high agreement with the solutions of the fourth order Runge-Kutta method, also it is more accurate than the results obtained in Maleknejad [40] with the same number of retained modes.
We solve Equation (32) for the case corresponds to γ = 0.1, 0.5.In Table 6 we give a comparison between CFMM and the operational matrix method developed in Maleknejad [40] for the case corresponds to α = 1, γ = 0.1, 0.5.

Remark 3.
From the numerical results of Table 6, we conclude that CFMM with fewer number of terms is compared favorably with the fourth-order Runge-Kutta solution obtained from Mathematica.
subject to the initial conditions: u where, er f (x) = 2 √ π x 0 e −t 2 dt.
The exact solution of Equation ( 33) is u(x) = e x .To illustrate the advantages of our method if compared with methods which employ orthogonal polynomials, We solve Equation (33) by our proposed algorithm namely, TFMM, and also by using Legendre spectral tau method (LSTM).Table 7 presents a comparison between the maximum pointwise errors which resulted from the application of TFMM and LSTM for different values of n.In addition, a comparison between the running times (in seconds) if the two methods are applied is also displayed in Table 7.The results of this table indicate the advantages of our method if compared with LSTM.Moreover, they ensure that TFMM is efficient for higher order fractional differential equations.Remark 4: From the obtained numerical results in Tables 1-7 and the obtained estimates in Theorems 2 and 3, we list here the advantages of the current work: (1) The generation of Fibonacci Polynomials is very easy either from the recurrence relation or from the direct definition.(2) The implementation of Fibonacci polynomials is efficient compared to the orthogonal polynomials in other words.The running time of the "Mathematica" code is very short compared with using the orthogonal polynomials.(3) The speed of convergence of Fibonacci polynomials (inverse factorial rate of convergence) is very efficient.
Theorem 2 ensures this result.

Conclusions
In this work, two numerical algorithms for solving linear and nonlinear two-term fractional order differential equations are developed.The basic idea behind the proposed algorithms is based on constructing a new operational matrix of fractional derivatives of Fibonacci polynomials to reduce both linear and nonlinear fractional differential equations to systems of linear or nonlinear algebraic equations which can be solved efficiently.As far as we know, the derived Fibonacci operational matrix of fractional derivatives is novel and it is the first time to use such matrix in solving fractional-order differential equations.The tested numerical examples show the high efficiency and performance of the proposed algorithms.

Table 3 .
Maximum absolute error of Example 4.

Table 3 .
Maximum absolute error of Example 4.
Remark 2. The results of Table