An Integral Operational Matrix of Fractional-Order Chelyshkov Functions and Its Applications

: Fractional differential equations have been applied to model physical and engineering processes in many fields of science and engineering. This paper adopts the fractional-order Chelyshkov functions (FCHFs) for solving the fractional differential equations. The operational matrices of fractional integral and product for FCHFs are derived. These matrices, together with the spectral collocation method, are used to reduce the fractional differential equation into a system of algebraic equations. The error estimation of the presented method is also studied. Furthermore, numerical examples and comparison with existing results are given to demonstrate the accuracy and applicability of the presented method.


Introduction
Fractional differential and integral operators are generalizations of differential and integral operators of integer order. In physical sciences, differential equations are used to model phenomena. However, differential equations of integer order cannot give acceptable results for complex systems. So, fractional differential equations are used to improve these models [1][2][3][4]. Several papers have been devoted for studying existence and uniqueness of solutions to the fractional differential equations [5,6]. Solving fractional differential equations is a desirable objective. There are some issues in solving these equations analytically except for a limited number of these equations. So, several studies have been reported for the development of new methods for finding numerical or approximate solutions, such as fractional differential transformation method [7], Adomian decomposition method [8], homotopy perturbation method [9], homotopy analysis method [10], variational iteration method [11], Galerkin method [12], collocation method [13], wavelet method [14], B-Spline operational matrix method [15] and Jacobi operational matrix method [16].
Spectral methods, finite differences and finite element methods are the main methods of discretization that provide a numerical solution of differential equations. The spectral methods have the advantage that they are accurate for a given number of unknowns [17,18]. These methods exhibit exponential rates of convergence/spectral accuracy for smooth problems, and this differs from finite difference and finite element methods, which offer only algebraic convergence rates [19]. In spectral methods, the solution of differential equations is approximated as an expansion in terms of orthogonal polynomials. The most commonly used spectral schemes are the collocation, tau and Galerkin approaches [20,21].
Orthogonal functions play an important role in finding numerical solution of differential equations. These functions simplify the treatment of differential equations via converting its solution into the solution of a system of algebraic equations. Some examples are shifted Legendre polynomials [22], Chebyshev wavelets [23], shifted Chebyshev polynomials [17], shifted Jacobi polynomials [24], block pulse operational matrix [25], etc. Orthogonal Chelyshkov polynomials were presented in [26]. Using these polynomials, a solution has been given of weakly singular integral equations in [27], Volterra-Fredholm integral equations in [28], linear functional integro-differential equations with variable coefficients in [29] and Volterra-Fredholm-Hammerstein integral equations in [30]. Furthermore, the operational matrix of fractional derivatives based on Chelyshkov polynomials for solving multi-order fractional differential equations has been presented in [31] and the operational matrix of fractional integration to solve a class of nonlinear fractional differential equations has been introduced in [32]. Recently, Talaei [33] has proposed a numerical algorithm based on FCHFs for solving linear weakly singular Volterra integral equation.
Solving fractional differential equations by using series expansion of the form such as Adomian's decomposition method, homotopy perturbation method and He's variational iteration method, is a common and efficient method [34]. So, building orthogonal functions of fractional-order may be useful in solving fractional differential equations more successfully.
The motivation of this paper is to construct the FCHFs based on Chelyshkov polynomials. Additionally, the operational matrices of the fractional integration and product for FCHFs are derived. In addition to that, an application of FCHFs for solving linear and nonlinear fractional differential equations is presented. The paper is organized as follows. Following this introduction, some definitions and preliminaries of fractional calculus, as well as the FCHFs and their properties are presented in Section 2. In Section 3, we investigate the error of the proposed method. The FCHFs operational matrices of fractional integration and product are given in Section 4. Section 5 presents a new technique to use the FCHFs operational matrices for solving multi-order fractional differential equation. The performance of the proposed algorithm is reported in Section 6 with satisfactory numerical results. Finally, the paper ends with conclusions in Section 7.

Fractional Calculus
Now, we will present some definitions and properties of fractional integration and differentiation which will be used throughout this article [35,36]. In the literature, one can find many definitions of fractional integration of order α ≥ 0, but these definitions may not equivalent to each other [37]. The definition of Riemann-Liouville is the more widely used, which can be stated as follows Definition 1. The Riemann-Liouville fractional integral operator of order α ≥ 0 of a function y(x) is given by For Riemann-Liouville fractional integral operator I α we have where λ 1 and λ 2 are constants.

Definition 2.
The fractional derivative of order α > 0 in the Caputo sense is defined in the form where α is the smallest integer greater than or equal to α.
In the case of η = 1, the generalized Taylor's formula reduces to the classical Taylor's formula.

Fractional Chelyshkov Functions
The set C n = {C ni } n i=0 of Chelyshkov polynomials is a set of orthogonal polynomials of degree n over the interval [0, 1] with the weight function ω(z) = 1 and has the form It can be seen that, all members of the set C n have the same degree n, which represents the main difference with other sets of orthogonal polynomials. The orthogonality condition of these polynomials is Chelyshkov polynomials C ni (z) are related to Jacobi polynomials P (α,β) m (z), where α, β > −1 and m ≥ 0 by the following relation The fractional-order Chelyshkov functions (FCHFs) can be built by substituting z = x η and η > 0 in the Chelyshkov polynomials. Suppose that C η ni (x) denotes the FCHFs C ni (x η ). Then, from Equation (12), C η ni (x) can be defined as These functions are orthogonal with respect to the weight function ω η (x) = x η−1 over the interval [0, 1] and satisfy the orthogonality property Let the FCHFs vector be which can take the form where the vector X η (x) and the upper triangular matrix E are given by

Function Approximation and Error Estimation
Assume that Λ = [0, 1] and define the weighted space L 2 ω η (Λ) by with inner product and norm Additionally, y n (x) ∈ S n , thus it can be expanded in terms of FCHFs as where the coefficient vector Q is given by The following theorem gives an upper bound for estimating the error.
Theorem 1. Assume that D jη y(x) ∈ C(0, 1] for j = 0, 1, . . . , n + 1 and y n (x) is the best approximation to y out of S n then the error bound is given by where Proof. Let p(x) be the generalized Taylor's formula of y(x), then from Definition 3 we get with error bound as However, p(x) ∈ S n and y n is the best approximation to y from S n then and by taking the square roots, we have the desired result.

The Fractional Integration Operational Matrix of FCHFs
In this section, the FCHFs operational matrices of fractional integration and product are obtained, which can be built easily and their numerical results are accurate.
where P (α) = [p ik ] n i,k=0 is the (n + 1) × (n + 1) operational matrix of fractional integration of order α > 0 in the Riemann-Liouville sense and its elements are given by Proof. By integrating Equation (18), we obtain and by using Equation (2), we get where The vector X η (x) can be approximated in terms of FCHFs as follows where a rk can be obtained from Equations (15) and (26) as follows Then, where A = [a rk ] n r,k=0 is a (n + 1) × (n + 1) matrix and its elements are given by Equation (39). From Equations (34), (35) and (40), we can write and by using Equations (20), (36) and (39) for multiplying the matrices E, B and A, we get Remark 1. It can be noted that the operational matrix of fractional integration in [32] is a special case of the operational matrix of fractional integration of FCHFs with η = 1.
In the following, we present the operational matrix of product of two FCHFs vectors, which will be useful to reduce the fractional differential equation into a set of algebraic equations. Theorem 3. Suppose U = [u 0 , u 1 , . . . , u n ] T is an arbitrary vector, then Proof. From Equations (18)- (20), we have Now, we approximate x (i+k)η for i, k = 0, 1, 2, . . . , n, in terms of FCHFs, as follows where Then, for every r = 0, 1, . . . , n, we get By using Equation (48) into Equation (45) we obtain which completes the proof.

Solution Method
In this section, the operational matrices of fractional integral and product for FCHFs together with the spectral collocation method are applied for solving the fractional differential equations.
By substituting as needed from Equations (52), (59) and (62) in Equation (50), the multi-order fractional differential Equation (50) is converted into the following algebraic equation To find the unknown vector Q, let Ψ n = {x i : x i = i n , i = 0, 1, .....n} be a set of equidistant nodes and collocate Equation (63) at the nodes x i , i = 0, 1, .....n, which gives Equation (64) represents a system of n + 1 nonlinear algebraic equations and can be solved to find the vector Q. So, an approximate solution of Equation (50) can be determine.

Illustrative Examples
In order to demonstrate the efficiency of the proposed method, we apply it to solve some linear and nonlinear fractional differential equations. A comparison with other methods reveals that the presented method is effective and accurate. Example 1. Consider the inhomogeneous Bagley-Torvik equation [16,22] The exact solution of this problem is y(x) = x + 1. By using the technique presented in Section 5, we have Substituting Equations (66) in (65), gives By taking n = 2, η = 1, we obtain Solving Equation (72) gives Q = 0, 0, 0 T , and therefore which is the exact solution. Additionally, for the same n = 2 with η = 1 2 we can obtain the exact solution, where Similarly, by collocating Equation (67) at the nodes x i = i 2 , i = 0, 1, 2 and solving the resulting system, we get the solution Q = 0, 0, 0 T , which yields Example 2. Consider the following linear fractional differential equation [38,39] The exact solution of this problem is y(x) = √ x. Using the proposed method yields Equation (79) can be written after using Equation (80) in the form With n = 1 and η = 1 2 , we get and the generated set of linear algebraic equations is The solution of Equation (84) Then which is the exact solution.
Example 3. Consider the nonlinear initial value problem [16,22,31] The exact solution of the problem is y(x) = x 2 . By applying the method described in Section 5, we obtain where Q 1 can be calculated from Equation (43). Equation (87) transforms Equation (86) to With n = 3 and η = 1, the resulting system of nonlinear algebraic equations can be written in the form The solution of system (89) is Q = 0, 0, 0, 0 T , which leads to the exact solution.

Example 4.
Consider the nonlinear multi-order fractional differential equation [40] D α y(x) The problem has the exact solution y(x) = x α+1 Γ(α + 2) . Applying the proposed technique in Section 5 for this problem, we get Then (90) takes the form Consider α = 1 with the exact solution y(x) = x 2 2 . For solving this case we take n = 2 and η = 1, which leads to the following system of nonlinear algebraic equations 1 36 Solving these equations yields Q = 0, 1 4 , 5 4 T , which leads to the exact solution.
Example 5. Consider the following fractional differential equation [31] Dy The exact solution of this problem is y(x) = x 2 √ x. Table 1 shows a comparison of the results obtained by the present method and those in [31] in terms of L ∞ ω η and L 2 ω η errors for different values of n and η. Note that the results for the method that we compare have been executed by the method in [31] and we use these results to make a direct comparison with the presented method. Symbol "-" means that the result for n is unavailable for that method. From Table 1, it can be seen that the errors achieved by the presented method are less than those in [31] for all values of n. In addition to that, when n increases, the errors are reduced until they become zero at n = 10, η = 0.25 and n = 13, η = 0.25. This means that the presented method is more accurate than that in [31] for this problem.
x 4−α , 0 < α < 1, The exact solution of this problem is y(x) = x 4 − 1 2 x 3 . Table 2 describes the results obtained by the present method and those in [31,41] in terms of L 2 ω η errors for different values of n with α = 1 4 .
It can be seen that the new method performs better than the method in [31] and much better than the technique in [41]. The error achieved by the presented method becomes smaller and smaller with the increment of n and converge to zero at n = 20. Which means that the introduced method is more coincidental with exact solutions than those in [31,41].

Conclusions
In this paper, the FCHFs have been adopted for solving fractional differential equations. The operational matrices of fractional integral and product for FCHFs have been derived. These matrices are used to approximate solutions of fractional differential equations where the fractional derivatives are considered in Caputo sense. The introduced method is a spectral collocation method which reduces the fractional differential equation to a system of algebraic equations. The efficiency and applicability of the presented method are tested on some problems. Comparison with some other methods has been performed. The numerical results show that the new method is more efficient, its performance is quite satisfactory, and only a small number of FCHFs is needed to obtain good results. Finally, we can say that employing fractional-order basis functions achieves more accurate results than the corresponding integer-order basis functions, which shows the applicability of this method for solving the fractional problems.