Numerical Study of Caputo Fractional-Order Differential Equations by Developing New Operational Matrices of Vieta–Lucas Polynomials

: In this article, we develop a numerical method based on the operational matrices of shifted Vieta–Lucas polynomials (VLPs) for solving Caputo fractional-order differential equations (FDEs). We derive a new operational matrix of the fractional-order derivatives in the Caputo sense, which is then used with spectral tau and spectral collocation methods to reduce the FDEs to a system of algebraic equations. Several numerical examples are given to show the accuracy of this method. These examples show that the obtained results have good agreement with the analytical solutions in both linear and non-linear FDEs. In addition to this, the numerical results obtained by using our method are compared with the numerical results obtained otherwise in the literature.


Introduction
Fractional calculus has been playing a very important role in scientific computations. Scientists are able to describe and model many physical phenomena with fractional-order differential equations. As a result, fractional-order differential operators are widely used to solve systems by developing more accurate models [1][2][3][4]. The nonlocal property of the fractional-order operators makes them more efficient for modeling the various problems of physics, fluid dynamics and their related disciplines [1,[5][6][7][8][9]. For example, consider a thin rigid plate of mass a 1 and area R immersed in a Newtonian fluid of infinite extent and connected by a massless spring of stiffness K to a fixed point. A force g(z) is applied to the plate. Assume that the spring has no effect on the fluid and that the area of the plate is large enough to produce the fluid adjacent to the plate, whereas stresses σ(z, x) can be defined by the following relation: where x is the distance of a point in the fluid from the spring to the submerged plate. By some assumptions discussed in [10], the dynamics of the system are given by where σ(z, 0) = Dυ(z). Equation (2), with some assumptions considered in [10], takes the following form of a Bagley-Torvik-type problem solved in (Section 6, Example 1): a 1 D 2 υ(z) + a 2 D 3/2 υ(z) + a 3 υ(z) = g(z), z ∈ [0, 1].
The existence and uniqueness results of fractional-order differential equations (FDEs) have been investigated extensively in the literature. Some of them are presented as follows: Fazli and Nieto [11] investigated the existence and uniqueness of the solution of FDEs of Bagley-Torvik type by considering the existence of coupled lower and upper solutions. Pang et al. [12] investigated the existence and uniqueness of the solution of the generalized FDEs with initial conditions by proposing a novel max-metric containing a Caputo derivative. Abbas [13] studied the existence and uniqueness of the solution of FDEs by using Banach's contraction principle together with Krasnoselskii fixed point theorem. For more works on existence and uniqueness results, we refer the reader to [14][15][16][17].
As most FDEs do not have closed-form solutions, different numerical techniques, including the finite difference method, variational iteration method and spectral methods, are preferably used. Among them, spectral methods have received considerable attention from the fractional community for solving FDEs, both ordinary and partial. Spectral methods are classified into three types, known as the collocation, tau and Galerkin methods. The basic idea of the spectral methods is to write the solution as a linear combination of basis vectors of global polynomials, typically Legendre, Jacobi and Chebyshev. The speed of convergence is considered the best advantage of the spectral methods, as the rate of convergence is exponential in these methods, which gives a high level of accuracy. Many efficient spectral techniques are obtained in the literature using the various global polynomials [18][19][20][21].
The construction of the operational matrices of fractional derivative operators defined with singular or nonsingular kernels has played a key role in the development of spectral methods. Many researchers have worked on the construction of the operational matrices of fractional derivatives using different types of global polynomials. For example, Benattia et al. [22] introduced the operational matrix of the fractional derivatives to develop a numerical method that is based on the Chebyshev wavelet for solving FDEs. Saadatmandi et al. [23] derived an operational matrix of derivatives of fractional order using the fractional-order Chebyshev functions. They also extended the results of [23] for solving the coupled system of FDEs with variable coefficients [24]. Additionally, Bharway et al. [25] introduced a new shifted Chebyshev operational matrix of fractional integration for solving linear FDEs. Moreover, Talib et al. [26] developed a new operational matrix based on the orthogonal shifted Legendre polynomials to numerically solve the fractional partial differential equations. Meanwhile, Rahimkhani et al. [27] introduced a Bernoulli wavelet operational matrix of fractional integration for obtaining the approximate solution of a fractional delay differential equation. Kazem et al. [18] derived an operational matrix that generalized the results presented in [19]. Recently, Dehastani et al. [28] calculated modified operational matrices of integration and pseudo-operational of fractional derivatives for the Lucas wavelet functions to compute the numerical solution of fractional Fredholm-Volterra integro-differential equations. Moreover, Dehastani et al. [29,30] also derived operational matrices of fractional-order derivatives and integration for fractional-order Bessel functions and fractional-order hybrid Bessel functions. In the derived numerical techniques, the operational matrices are applied to reduce the FDEs to a system of algebraic equations.
Dehestani et al. [31] also presented a novel collocation method based on the Genocchi wavelet for the numerical solution of FDEs and time-fractional partial differential equations with delay.
Motivated by the aforementioned works, we extend the study of the spectral methods by constructing a numerical algorithm that is based on the fractional-order derivative operational matrix of VLPs in Caputo sense, together with the spectral tau method and spectral collocation method. The basis vectors of VLPs are used to approximate the solution of the problems. The derivative terms are approximated by using the fractional-order derivative operational matrix of VLPs. It is important to mention that the proposed algorithm is computer-oriented and is capable of reducing the FDEs to a system of algebraic equations, which greatly simplifies the problems. Subsequently, we use the operational matrices approach together with the spectral tau method in the case of linear FDEs, and the operational matrices approach together with the spectral collocation method in the case of nonlinear FDEs. Our proposed numerical algorithm produces highly efficient numerical results as obtained otherwise in the literature [32][33][34].
The novel aspects of our proposed study are the development of the new fractionalorder derivative operational matrix of VLPs in Caputo sense and the construction of the numerical algorithm that is based on this newly developed operational matrix. To the best of our knowledge, this is the first result where the numerical algorithm is presented by using the operational matrix of VLPs. Moreover, the proposed numerical algorithm is fit to solve both linear and nonlinear FDEs with initial conditions. In addition, our proposed method has advantages over other methods, such as the Homotopy perturbation method, because, in our case, the perturbation, linearization or discretization are not necessary to be implemented.
The structure of this paper is set in the following way. In Section 2, we discuss the VLPs along with their properties. In Section 3 , the Vieta-Lucas operational matrix of fractional-order derivatives is derived. In Section 4, the numerical method is developed by using the operational matrices of VLPs. In Section 5, the error bound is determined. In Section 6, the accuracy and the stability of the proposed method are analyzed by taking some numerical examples. In Section 6, we conclude and give the summary of this paper.

Preliminaries
In this section, we summarize some definitions, properties and results of fractional calculus that are essential to construct the numerical algorithm to solve the linear and nonlinear FDEs. Definition 1. The Riemann-Liouville fractional integral operator of order α > 0, of a function υ, is defined as: Definition 2. The Caputo operator of the fractional derivative is defined as follows: where n − 1 < α ≤ n, n ∈ N, and υ ∈ C n [0, 1]. Hence, the Caputo operator follows:

Vieta-Lucas Polynomials
Vieta-Lucas polynomials belong to the class of orthogonal polynomials and can be created by using the recurrence relation [35]. Consider |z| ≤ 2; then, the Vieta-Lucas polynomials of degree n ∈ N 0 in the variable z can be defined as The VLPs can be created by using the following recurrence relation: VL n (z) = z VL n−1 (z) − VL n−2 (z), n = 2, 3, . . . , VL 0 (z) = 2, VL 1 (z) = z.
Moreover, VL n (z) can be expressed using the following power series formula: where n 2 is the ceiling function. Moreover, the orthogonality of VL n (z) can be expressed as: where 1 √ 4−z 2 is the weight function.

Shifted VLPs
As a new class of orthogonal polynomials, the shifted VLPs, VL * n (z) of degree n, defined on the closed interval [0, 1], can be obtained as follows: Moreover, VL * n (z) are created by the following formula: with the starting values Moreover, analytically, VL * n (z) can be expressed as: Let the function u(z) be Lebesgue-square-integrable on the interval [0, 1], which can be expressed in terms of VL n (z) as follows: where the undetermined coefficients, c j , j = 0, 1, 2, . . . , n, can be determined through the following expression: or where For approximation, we can take the first n + 1 terms of the series; therefore, u(z) can be expanded in the form where the shifted VLP coefficient vector C and the shifted VLP vector Ψ(z) are given by

Operational Matrices of Differentiation
Theorem 1. Let Ψ(z) be the shifted VLP vector defined in (18) and also suppose that α > 0; then, where P α is the (m + 1) × (m + 1) operational matrix of the fractional derivative of order α in the Caputo sense and is defined as follows: and ξ i,j,k is given by Proof. Applying the Caputo derivative to (12), we have Applying the linearity of the Caputo derivative, and using (5), we have and Now, approximate z i−k−α by (m + 1) terms of the series, as where and u(z) = z i−k−α . Now, inserting the value of u(z)and VL * j (z) into Equation (24), we obtain by inserting the value of z i−k−α into Equation (21), we obtain where Rewriting Equation (26) in vector form, we obtain For simplicity, we can write Equation (28) as: where Equations (22) and (29) prove the required result.

Application of Operational Matrices Method
In this section, we apply the Vieta-Lucas operational matrix method to find the analytical-approximate solution of linear and nonlinear FDEs.

Consider the linear FDE
with initial conditions where b l , for l = 1, . . . , k + 2, are real constant coefficients and also n < α ≤ n + 1, The unknown function υ(z) and the source term g(z) can be approximated as: where H = [h 0 , . . . , h M ] T is known, and C = [c 0 , . . . , c M ] T is an unknown to be determined. Now, using Equations (29) and (33), we have Using Equations (33)-(36), the residual R(x) for Equation (31) can be written as Using the spectral tau method [36], a system of linear equations is generated by applying Moreover, by substituting Equation (33) in the initial conditions given in Equation (32), we obtain Equations (38) and (39) generate the (M − n) and (n + 1) set of linear equations, respectively. This system of linear equations can then be solved easily for the unknown coefficients. Consequently, we can approximate υ(z) given in Equation (33).

Error Estimate
Lemma 1 ([37]). The following assumptions for the function g(z), such that g(k) = b k , must hold true: 1.
The function g(z) is positive, decreasing and continuous for z ≥ m.

2.
∑ b m is convergent, and P m = ∑ ∞ k=m+1 b k . Then, Proof. It is evident that the shifted VLPs are orthogonal on the interval [0, 1] with respect to the weight function, w(z) = 1 √ z−z 2 . Hence, these polynomials form a complete L 2 w (∆) orthogonal set, where L 2 w (∆) represents the space of functions defined as v : ∆ → R. Thus, the error in space Using Equations (13) and (17) in Equation (47) Now, by applying the orthogonality property of shifted VLPs to Equation (48) Now, by using the substitution 4x − 2 = 2cos(θ) in Equation (15), the coefficients, b k , k = 0, 1, · · · , m, can be determined as Equation (50) can also be expressed as Using v (z) ≤ N, and the properties of trigonometric functions, we may express Equation (51) as Now, using Equation (52) in Equation (49) Now, by using the Lemma 1, we have

Illustrative Examples
In this section, we give some numerical examples to show the accuracy of our proposed method.

Example 1. Consider the following fractional Bagley-Torvik equation
subject to the initial conditions with integer order The source term g(z) is as follows: g(z) = 1 + z.
The exact solution of the problem in Example 1 is: Now, we apply the technique that is described in Section 4.1 by choosing the first three terms of VLPs. We may write the approximation solution as Now, we have Vieta-Lucas polynomials and The Vieta-Lucas operational matrices can be expressed as The residuals can be evaluated as where C = [c 0 , c 1 , c 2 ]. Now, using initial conditions, we have: Moreover, using the inner product of the residual with the Vieta-Lucas polynomials, we obtain a system of equations. If we take one equation from this system and two equations from Equation (65), then, by simultaneously solving these equations, we obtain c 0 = 0.75, c 1 = 0.25, c 2 = 0, hence which is the exact solution.

Remark 1.
The numerical results computed using our method are compared with the method of [33] by choosing various n. We observe that our proposed method produces efficient numerical results as compared to the numerical results obtained by using the method of [33] (see Tables 1 and 2 and Figure 1). Moreover, for a small value of n = 2, the exact solution and the approximate solution computed by using our method coincide (see Table 1 and Figure 1).

Remark 2.
The numerical results of Example 1 computed at n = 2 by using our method are compared with the results obtained by using the methods of [32,34] at n = 6. We observe that, for a small value of n = 2, the approximate solution obtained using our method coincides with the exact solution of Example 1 (see Tables 3 and 4). However, the exact solution and the approximate solution computed by using the methods of [32,34] coincide at n = 6. This shows that our proposed method is numerically more efficient. Example 2. Consider the following linear initial value problem [38]: The exact solution of the problem is [39]. To solve the problem, we use the technique described in Section 4.1. The absolute error for α = 0.85 and n = 2 , 5 and 8 is shown in Table 5. An error plot is also shown in Figure 2 for these values.
We can see in Table 5 that a good approximation has been achieved by using some initial terms of VLPs. Moreover, the numerical results for υ(z) when n = 10 and α = 0.5, 0.65, 0.8, 0.95 and 1 are plotted in Figure 3. The exact solution for α = 1 is υ(z) = exp(−z). It can be noted that the numerical solution converges to the analytical solution when α approaches 1. We also analyze the nonlocal behavior of the fractional derivative by computing the results at various non-integer values of α, which highlights the advantage of using the fractional derivatives, as the next state of the system depends not only upon its current state by also upon all of its historical states.  Example 3. Consider the following initial value problem [10]: subject to the initial condition with integer order The source term is given as 8Γ( 13 4 ) .
The exact solution at a 1 = −1 is given below: We test the behavior of our proposed method by solving Example (3) at various values of n. In Table 6, we list the L ∞ and L 2 errors for different values of n. We compare the numerical results obtained by using our proposed method with the numerical results obtained in [10]. It can be observed that the errors computed by using our method are much smaller than those computed by using the method presented in [10]; see Table 6. This highlights the efficiency of our method for this problem. Note that the symbol " − " means that the result for n is unavailable for the method [10]. Example 4. Consider the following nonlinear initial value problem [40]: The exact solution of the problem is υ(z) = z 8 − 3z (4+α/2) + 9 4 z α [39]. We have solved the problem using the technique described in Section 4.2. The absolute error for α = 0.85 and n = 2, 5 and 8 is shown in Table 5. An error plot is also shown in Figures 4 and  5 for these values. We can see in Table 7 that a good approximation has been achieved. Numerical results for υ(z) when n = 6 and α = 0.6, 0.7, 0.8, 0.9 and 1 are plotted in Figure 6, along with the exact solutions at the given values of α. It can be noted that, as α approaches 1, the solution of the FDEs approaches that of the integer-order differential equations. We also analyze the nonlocal behavior of the fractional derivative by computing the results at various non-integer values of α, which highlights the advantage of using the fractional derivatives, as the next state of the system depends not only upon its current state but also upon all of its historical states (Table 8).    Exact Sol at =0.7 Approx. Sol at =0.7 Exact Sol at =0.8 Approx. Sol at =0.8 Exact Sol at =0.9 Approx. Sol at =0.9 Exact Sol at =1 Approx. Sol at =1 We solved this problem by using the same technique as described in Section 4.2 with n = 3.
The exact solution of the problem is υ(z) = z 2 , whereas the source term is g(z) = z 4 . The operational matrices can be expressed as Now, using the initial condition, we obtain three equations: Meanwhile, using the technique in Section 4.2, we obtain the following equation: Now, we collocate Equation (74) at the first root of P 4 (z), and we obtain z 0 = 0.06698. Now, solving Equations (73) and (74), we obtain which is the exact solution.
Example 6. Consider the following non-homogenous multi-order fractional problem: subject to the following initial conditions The source term is as below: The exact solution corresponding to α = 2 ,a = c = −1, b = 2, d = 0, β 0 = 0, β 1 = 1, β 2 = 1 2 is given below: We can observe in Example 6 that a good approximation of the function has been achieved while using n = 7 as a scale level. The absolute error Table 9 at different scale levels is given below. (see Figures 7)  Moreover, a comparison between the exact solution and approximate solution at different scale levels for the values of z is given in Table 10.

Conclusions
In the present study, we introduce a new fractional-order derivative operational matrix of VLPs in Caputo sense. The newly derived operational matrix is used to develop a computer-oriented numerical algorithm to solve the linear and nonlinear FDEs that include the Caputo fractional-order derivative. The proposed numerical algorithm has the advantage of transforming the problems into a system of algebraic equations that are easy to solve using any computational software. To the best of our knowledge, this is the first result where the numerical algorithm is presented using the operational matrix of VLPs and the solution of the problems is approximated using its basis vectors.
We tested the accuracy and efficiency of the algorithm by solving various linear and nonlinear FDEs with initial conditions. We found that with an increase in the values of n, the approximate solutions were in good agreement with the exact solutions. We also demonstrated the high efficiency of the method by determining the amount of absolute error and observed that as we increased n, this amount was decreased significantly. In addition to this, the numerical efficiency was also demonstrated by comparing the results obtained by using our method with results obtained otherwise in the literature [10,[32][33][34]. We observed that our method produced more efficient results.