A New Operational Matrices-Based Spectral Method for Multi-Order Fractional Problems

: The operational matrices-based computational algorithms are the promising tools to tackle the problems of non-integer derivatives and gained a substantial devotion among the scientiﬁc community. Here, an accurate and e ﬃ cient computational scheme based on another new type of polynomial with the help of collocation method (CM) is presented for di ﬀ erent nonlinear multi-order fractional di ﬀ erentials (NMOFDEs) and Bagley–Torvik (BT) equations. The methods are proposed utilizing some new operational matrices of derivatives using Chelyshkov polynomials (CPs) through Caputo’s sense. Two di ﬀ erent ways are adopted to construct the approximated (AOM) and exact (EOM) operational matrices of derivatives for integer and non-integer orders and used to propose an algorithm. The understudy problems have been transformed to an equivalent nonlinear algebraic equations system and solved by means of collocation method. The proposed computational method is authenticated through convergence and error-bound analysis. The exactness and e ﬀ ectiveness of said method are shown on some fractional order physical problems. The attained outcomes are endorsing that the recommended method is really accurate, reliable and e ﬃ cient and could be used as suitable tool to attain the solutions for a variety of the non-integer order di ﬀ erential equations arising in applied sciences.


Introduction
The subject of fractional calculus (FC) and fractional or non-integer differential equations (FDEs) is an inspirational mathematical topic of research among scientific community. In 1695, this antique area is carried out from a famous dialogue between two conspicuous mathematicians Hospital and Leibniz. Indeed, the FC area is an old topic as traditional calculus proposed individually by Leibniz and Newton, but FC gained a vital importance among the researcher community. The devotion of scientific communities is pretty realistic due to its specific appearance in mathematical as well as other areas of science and engineering [1]. However, it converts the symbolization of derivatives for those sections where the order of the derivatives is non-integer or fractional. Recently, the subject of FC has been developed in different domains of sciences, from pure mathematical theories to modeling Symmetry 2020, 12, 1471 2 of 22 fractional order physical problems in different disciplines of engineering and applied sciences including ultra-capacitor, beam heating, transfer of heat in heterogeneous media, etc. [2]. Physical and natural mechanisms can be modeled using the differential equations and the solution of integer and non-integer order problems are worthy to compute because it saves both time and money [3][4][5]. The main impact of using FDEs states the non-local properties, which means that the present state and all its earlier states impact the subsequent state of the dynamical systems. The mentioned practical situations are normally due to the fact that a number of physical structures are associated to non-integer order dynamics and their deeds are directed by fractional differential equations (FDEs). The readers are referred to see the cited references [1,6] to review some fundamental concepts of fractional calculus and its potential applications as well as various techniques to compute fractional or non-integer order integrals/derivatives. The fractional derivative of n − 1 < δ < n in Caputo's sense [7] is used in this work and stated below in Equation (1).
if v = n ∈ N. (1) The differential operator C 0 D t satisfies the properties in (2) where, γ and λ are the constants.
Many useful processes in engineering, applied and physical sciences can be effectively modeled by means of fractional order differential equations. A noticeable knowledge and analysis related potential application and theory of FC is cited in [8]. Recently, FC is arising in various disciplines of sciences including modeling of different fluid flow, mathematical biology, modeling of water waves, control and optimization theory, dynamical, advection-diffusion, electrical circuits, networking, etc. [9][10][11]. It is an imperative fact that the results of physical models based on FDEs are additionally complex compared with integer-order differential equations (DEs) or physical and engineering problems. The devotion of mathematical communities is to accomplish the mathematical modeling for a given physical system in terms of differential or integral equations and fetch a competent solution to deliver a hypothetical study that reaches to the similar level of experimental outcomes which makes the FDEs an emerging and hot research area in recent times.
To achieve this high level of accuracy and better analysis, various kinds of techniques including, analytical, numerical, soliton and wavelet have been adopted in the past [12][13][14][15][16][17][18][19][20][21]. Literature survey proved that there are many lacks to obtain the better approximation and analysis of fractional ordered problems arising in applied sciences. Recently, approximation through orthogonal basis functions has acknowledged more suitable and gained a substantial devotion in different areas of engineering and mathematical science. Indeed, algorithms depending upon the orthogonal basis functions converts the nonlinear problems under study to a set of linear or nonlinear algebraic equations and the obtained system could be solved to approximate the given problem. Previously, various kinds of polynomials including Bernstein, Chebyshev, Chelyshkov, Gegenbauer and shifted Gegenbauer, Jacobi, Laguerre, Laurent, Legendre and a few others had been utilized by researchers to examine different complex nature physical problems [22][23][24][25][26][27][28][29][30][31]. Lately, alternative kinds of orthogonal polynomials have been familiarized and adopted to evaluate the dynamics of the several types of problems governed by differential/integral equations. An algorithm based on operational matrices of integration using Chelyshkov polynomials is proposed by Mohammadi [32]. A Chelyshkov wavelet-based technique is developed to compute the control theory problems by Moradi et al. [33]. Readers are suggested to find some associated research to such polynomials in [34][35][36][37][38].
The motivation of the present work is to develop and compare some operational matrices of derivatives and use to propose a methodology to investigate the fractional order problems. The novel operational matrices are developed with the aid of a new kind of polynomial. The matrices are obtained by using two different approaches and later utilized to propose a computational algorithm together with collocation technique named Chelyshkov polynomial method (CPM). The CPM with exact operational matrix of derivative is termed as CPMEOM while CPM with approximate operational matrix is CPMAOM. The detailed evaluation is given in upcoming sections. The nonlinear fractional problems (multi-order and BT equations) under study have been investigated to show the reliability of the method. The BT-equation arises in FDE-NNs architecture and the flow chart is available in Figure 1 while a detailed description about the problem is given in Section 5. The outcomes have been compared with existing data in literature and exact solutions, whereas a set of graphical plots and a tabular comparison is presented. It is apparent that the developed computational scheme is an effectual tool and could be prolonged to evaluate the physical, dynamical and complex natured nonlinear non-integer order problems. More precisely, the method can be extended by coupling with other methods as well as for problem governed by integral equations. Jacobi, Laguerre, Laurent, Legendre and a few others had been utilized by researchers to examine different complex nature physical problems [22][23][24][25][26][27][28][29][30][31]. Lately, alternative kinds of orthogonal polynomials have been familiarized and adopted to evaluate the dynamics of the several types of problems governed by differential/integral equations. An algorithm based on operational matrices of integration using Chelyshkov polynomials is proposed by Mohammadi [32]. A Chelyshkov waveletbased technique is developed to compute the control theory problems by Moradi et al. [33]. Readers are suggested to find some associated research to such polynomials in [34][35][36][37][38]. The motivation of the present work is to develop and compare some operational matrices of derivatives and use to propose a methodology to investigate the fractional order problems. The novel operational matrices are developed with the aid of a new kind of polynomial. The matrices are obtained by using two different approaches and later utilized to propose a computational algorithm together with collocation technique named Chelyshkov polynomial method (CPM). The CPM with exact operational matrix of derivative is termed as CPMEOM while CPM with approximate operational matrix is CPMAOM. The detailed evaluation is given in upcoming sections. The nonlinear fractional problems (multi-order and BT equations) under study have been investigated to show the reliability of the method. The BT-equation arises in FDE-NNs architecture and the flow chart is available in Figure 1 while a detailed description about the problem is given in Section 5. The outcomes have been compared with existing data in literature and exact solutions, whereas a set of graphical plots and a tabular comparison is presented. It is apparent that the developed computational scheme is an effectual tool and could be prolonged to evaluate the physical, dynamical and complex natured nonlinear non-integer order problems. More precisely, the method can be extended by coupling with other methods as well as for problem governed by integral equations.

Chelyshkov Polynomials
Herein, we will determine the basic definition of Chelyshkov polynomial as well as some important properties. The function approximation of a function through Chelyshkov polynomials is also being made.

Classical Concepts and Properties
The classical concepts related to Chelyshkov polynomials (a new class of orthogonal polynomials) are presented. In 2006, Chelyshkov [34] introduced this class of polynomials while Gokmen et al. [35] (in 2017) redefined said polynomials explicitly as stated in (3) and (4).

Chelyshkov Polynomials
Herein, we will determine the basic definition of Chelyshkov polynomial as well as some important properties. The function approximation of a function through Chelyshkov polynomials is also being made.

Classical Concepts and Properties
The classical concepts related to Chelyshkov polynomials (a new class of orthogonal polynomials) are presented. In 2006, Chelyshkov [34] introduced this class of polynomials while Gokmen et al. [35] (in 2017) redefined said polynomials explicitly as stated in (3) and (4).
Symmetry 2020, 12, 1471 4 of 22 In the above expression, the constant a j,i is defined as follows.
The orthogonality of Chelyshkov polynomials [34] over the interval [0, 1] can be proved using the expression below (see Equation (5)), wherein the weight function is taken as 1. In Equation (5), δ mn is Kronecker delta and defined as: 1, if m and n are same, 0, otherwise, .
Another way to compute the polynomials is use of Rodrigues formula [36] which is stated below: It is obvious from Chelyshkov polynomials definition that a specific digit value of N, gives exactly N degree polynomials of the given CPs or functions χ k,M (t), k ∈ N ≤ N. However, it makes a clear difference between CPs and other orthogonal function defined on [0, 1]. Therefore, estimation assisted by CPs brings more consistent results compared with further polynomials such as Legendre, shifted Chebyshev, second kind shifted Chebyshev polynomials, etc. where k is degree of kth polynomial. The Chelyshkov polynomials have correspondent properties to other orthogonal polynomials. However, CPs do not have hypergeometric sort of the solutions of the equations but can be identified in the form of Jacobi polynomials (JPs) which is stated in (6).
Therefore, the Chelyshkov polynomials are associated to an alternative class of Jacobi polynomials while keeping all distinctive features of regular orthogonal polynomials. Chelyshkov [34] has proved the properties of zero polynomial of the JPs [24] using the relation (6).

Corollary 1 ([34]
). The polynomials χ j,N (t) have j multiple zeros t = 0 and j − N distinct real zeros in the unit interval.
Hence, the polynomials χ 0,N (t) have exactly N simple roots in unit interval [0, 1], which has been used to construct the alternative Gaussian quadrature rule sated in Corollary 2.
It is exact for any polynomial of degree 2N if and only if t k are the zeros of the polynomial χ i,0 (t) and the weight factors are given in Equation (8) while the expression for omega_0 is stated in Equation (9).

Function Approximation
Let H = L 2 (Ω ) be the space of square integrable functions w.r.t Lebesgue measure on the interval Ω := [0, 1]. The inner product expression in the space can be stated as in (10) while the Norm is expressed in Equation (11).
Let H = L 2 (Ω ) be the space of square integrable functions w.r.t Lebesgue measure on the interval Ω := [0, 1]. The Equation (10) is the inner product expression in the space while the Norm is expressed in Equation (11).
Suppose that Since, H has a closed finite dimensional subspace H n , for every given function U ∈ H there exists a distinct best estimation U ∈ H n such that where, We can express a function U(t) on interval [0, 1) from L 2 (R) space by using Chelyshkov polynomials as in (12), c k we compute using relation (14).
In the above expression (12), C and ϕ(t) are N vectors defined in (13).

Operational Matrices of Derivatives
Herein, operational matrices of derivatives for both (positive integer and positive non-integer) order through a new family of functions. The matrices are reported via two non-similar methods. Theorem 1 concerns with an approximate operational-matrix (D A or AOM) of fractional-order derivative whilst Theorems 2 and 3 provide an exact operational matrix (D E or EOM) of non-integer order derivative. Theorem 1. Let ϕ(t)be the N Chelyshkov polynomial vector as stated in the Equation (13). Its fractional derivative of order ν can be expressed as follows: the D

(ν)
A is an N order matrix given as: n,m are (n, m)th components and computed via relation (16).
Applying the definition of fractional-order derivative to overhead expression, we have: We change the resulting expression by means of the Caputo's idea: Now we will expand the term t n+ j−ν through CPs and get the expression below, where β k is expanded through the relation given below: Inserting the above relation in Equation (17), we obtain the following expression: Now to construct the exact operational matrix for fractional-order derivative, first we define the following family of piece-wise functions defined on [0, 1] as: here, l = 0, 1, 2, . . . , N − 1. The nth piece-wise functions can be stated as a vector form: Theorem 2. Consider Θ(t) be a vector present in Equation (19) and where Λ is a matrix having order N × N which is stated as: Then the following relation must satisfied: Proof. We consider the following relation with the aid of Equation (12).

Lemma 2.
The fractional differentiation of order ν of Equation (19) in Caputo sense is defined as follows, where γ − 1 < ν < γ is a positive integer: here H ν is square having order N × N defined as: Proof. The relation can be proved with the aid of Lemma 1.

Theorem 3.
The differentiation of fractional order ν; γ − 1 < ν < γ, of Equation (13) can be described as: In Equation (25), M = Λ T and H ν are matrices presented above while D ν is the operational matrix for derivative of fractional or non-integer order ν for the Chelyshkov polynomials and needs to be identified.

Proof.
We have the subsequent form via Equation (20) In Equation (26), taking the fractional derivative w.r.t t of order v and we obtain the following form with the aid of Lemma 2.
The theorem stated below is of great significance in order to convert the differential equation into a system algebraic equation.
where, C = [ c ik ] N i,k=0 is the Chelyshkov operational matrix of product wherein.
The values of τ ijk are explained in the Lemma 2.5 of [26]. Proof. The readers are referred to see to cf. [26] to see the proof of above theorem.

Validation of Both Operational Matrices
Herein, we will validate the credibility of our proposed operational matrix of derivative by comparing both kinds of OMDs = 1/3.
Firstly, we calculate the operational matrix of derivative by considering ν = 1/3 and assuming polynomial vector of order 3. The detailed calculations are stated below while the Equation (33) is evident that the method to calculate the operational matrix of derivative is not accurate.
The particular value of t gives the following error in case of approximate operational matrix of fractional derivative.
On the other hand, the new kind of operational matrix of derivative is found and termed as exact operational matrix of fractional derivative. The detailed evaluations are given below.
It is noticed that the operational matrices of derivatives matrix obtained previously have issues about the accuracy and comprise some errors which means the method based on D A will also contain less accurateness. Hence, the approach to compute the new operational matrices of derivatives is found accurate and efficient whereas a method based on this kind of operational matrix i.e., D E will provide a better rate of accuracy as compared to traditional approximate matrices of derivative.

Convergence and Error-Bound Analysis
Some theorems related to convergence analysis and error-bound are presented in the current section [36]. Some classical ideas of norm, Kreyszig approximation property and Taylor series are utilized to complete the proofs. The norm L 2 I is denoted as . 2 and defined below:  1] be expended employing the CPs as: then, where, L = max t∈I U (N+1) (t) and I = [0, 1].
Proof. ([36]). Let the polynomial V(t) be defined as follows: The concept of Taylor expansion refers that there exists η ∈ (0, 1) such that: The Kreyszig approximation property gives the following expression: By means of the definition of L 2 − norm we get: .
Taking the square root completes the proof.
be expended by means of the Chelyshkov polynomials as: and have the following relation lim Proof. The proof of the theorem is presented by Hamid et al. [36].

Theorem 7.
Let the expansion of a continuous function say U(t) through Chelyshkov polynomials which converges uniformly. Therefore, the expansion will converge to U(t).

Proof. Let
where, c k = U, χ k,N (t) . Now, for fixed values of k, apply both sides of the above relation with χ k,N (t) and differentiating on unit interval [0, 1] which gives the following relation (37).
Hence, the specified functions U(t) and U(t) have equal Chelyshkov polynomial expansion and correspondingly U(t) = U(t) for t ∈ [0, 1].

Proposed Methodology and Numerical Experiments
Herein, a detailed description of the proposed methodology based on developed both approximate and exact operational matrices of derivatives is being made. The credibility of the algorithm is proved by talking some fractional (linear and nonlinear) problems. The solution of the Bagley-Torvik (BT) equation is presented to show the application of proposed algorithm in the neural network system. The solutions are asserted via a set of graphs and tabular form of the comparative study is also included in the same section.

Solution Procedure
Here we are presenting the methodology based on operational matrices of derivatives to solve the FDEs. The step by step procedure of our suggested method is given below. Consider the multi-order differential equation (MODE) as follows: where m ∈ Z + , 0 < ν ≤ 1, and 1 ≤ ν 1 ≤ 2 ν, v 1 ∈ Q are the order of the any differential equation and a, b, c, d, e are coefficients, while k is an integer showing the nonlinearity. The associated conditions with Equation (38) are: The proposed method has the following steps: Step I: Consider the following trial solution to investigate the solution of above problem. The λ and ϕ(t) are column vectors presented in above while λ needs to be computed.
Step II: The involved derivatives in the above-mentioned problem can be determined with the aid of theorems stated in Section 3 and related expressions are given below: By means of the suggested solution and Equations (41) and (42), we get the following form of the problem under study (see Equation (38)), along with the condition as follows: Step III: The system of algebraic equations can be constructed by means of collocation method. That is given below: where, O is the order of discussed problem.
Step IV: Putting the solution of obtained algebraic equation into the Equation (40) and we get the approximate solution of the problem given in Equations (38) and (39).

Algorithms Working Steps
The chart of the proposed algorithm and step by step procedure is presented in Figure 2.

Applications and Discussion
The section is devoted to a comprehensive discussion of obtained results by means of proposed method. The solutions are obtained via both approximate and exact operational matrices-based algorithms and efficiency of the method is verified by means of some numerical examples with various orders of fractional or integer order DEs. The simulations are performed by using Maple 2017. The graphical representation as well as a comparative study with exact solution and existing literature of the problems understudy is being made to show the appropriateness of the algorithms.
Subject to the initial conditions (ICs) and source term: Subject to the initial conditions (ICs) and source term: The solution of the Problem 1 is obtained for a = b = c = e = 1 while the derivative orders are taken as, v = 1.75, v 1 = 0.5 [39,40]. The D in the initial conditions is representing derivative w.r.t time. The solutions by means of using CPM via D E and D A are obtained for N = 3, and 100 digits are taken to the account. Firstly, we obtained operational matrices of derivatives as given below: where, D E and D A are respectively representing exact and approximate operational matrices of derivatives having a fractional order. The polynomial vector of order 4 for t is as follows: The proposed methodology in the previous section is adopted to examine the problem understudy and graphical illustration is with the help of set of graphs while a tabular form (see Tables 1-3) of comparison is presented with exact solution and previous literature [39,40]. It is noted that the CPM via D E brings more accurate and efficient results as compared to CPM via D A which is clear from the graphical plots. The Figure 3a is plotted for v = 1.75, v 1 = 0.5 and N = 3. One can observe that the CPM via D E is in excellent agreement with the exact solution while the CPM via D A converges to the exact solution as the value of N increases. The Figure 3b is plotted for integer values of derivatives and CPM by means of both kinds of OMs found in good agreement with exact solution which can also be noted form the last two columns of Table 2. The Tables 1 and 3      Problem 2. Consider the nonlinear multi-order fractional ordinary differential equation (NMOFDE) [39,40].

CPM via D E CPM via D
Subject to the initial conditions (ICs) and source term: The analysis of Problem 2 under study is reported for v = 2, v 1 = 0.333, v 2 = 1.234, while the values of coefficients a, b, c, and e are considered as 1. The simulations are being made for N = 3, and digits = 100. The operational matrices via both approaches has been computed for above mentioned values of fractional order derivatives and stated under:   In the above expressions, D A and D E are respectively representing approximate and exact operational matrices of derivatives having a fractional order while the Chelyshkov polynomial vector of order 3 for t is stated under: Again, exact solution is obtained using the algorithm based on D E while the algorithm based on D A still contains errors. The Figure 5a is presented to show the graphical behavior of the all three solutions. The tabular form of comparison is also included in the study. One can observe (see Table 4) for small values of N the error is decreased at a tangible level via the proposed method.
In the above expressions,  3   35  60  30 4  21  30  10  7 6 . solutions. The tabular form of comparison is also included in the study. One can observe (see Table  4) for small values of N the error is decreased at a tangible level via the proposed method.   The last column of the Table 5 witnesses that no error is found via CPMEOM. The readers can notice that for integer values of derivative involving both algorithms provide the same level of accuracy. The Figure 5b is an error curve for CPM via D A which is gradually decreasing. The Tables 4 and 6 are presented for various values of fractional derivatives and ensure the credibility of the newly developed OM (D E ) while, the Figure 5c is plotted to show the comparison with exact solution for arbitrary values of fractional derivatives. In account of set of graphs and tabular comparison, the suggested algorithm is found as an efficient tool to examine multi-order nonlinear problems of physical and complex nature.  The Bagley-Torvik (BT) equations are considered as an enormously essential part to examine the performance of various materials through the applications of fractional calculus and fractional differential equations [41,42]. It has increased its importance in various disciplines of industry, engineering and applied sciences. In the past, the said problem has been modeled for fractional orders 1 2 or 3 2 while the problem arises in frequency dependent damping materials. [41,42].

Problem 3. The mathematical form of Bagley-Torvik equation is stated in
Subject to the initial conditions (ICs), boundary conditions (BCs) and source term: In the above expressions (51)-(53), the nonlinear operator and unknown function are represented as n and U(t). The constant T is span input in the interval [0, T]. The constant coefficients are indicted as σ, µ, ω, while I m , m m are contents. The BT-equation also appear in the neural networking specifically in feed-forward artificial neural network (FFA-NN) and its graphical sketch is presented in the Figure 1. In addition, it is important to highlight that the said problem arises in fractional calculus as a fractional differential equation in neural networks (FDE-NNs). Readers are referred to see [41] for more details about said physical fractional problem of complex nature. It is noted that mathematical formulation of aforementioned problem is the linear combination of the discussed networks while the solution of unknown function could be appropriate unknown weights. In account of this presented model the numerical experiments have been performed when digits = 100, N = 2, σ = n = 1, ω = µ = 0.5 and the orders of derivatives are taken as: v = 2, v 1 = 1.5.
The operational matrices via both approaches have been computed for above mentioned values of fractional order derivatives and stated under: In the above expressions, D A and D E are respectively representing approximate and exact operational matrices of derivatives having a fractional order while the Chelyshkov polynomial vector of order 3 for t is stated under: The proposed scheme is adopted to examine the solution of physical problem understudy and noted that an exact solution is obtained using CPMEOM method for a very small value of N while the CPMAOM contains some error, which is asserted by means of graph (see Figure 6a) and tabular form of the comparative study is given in the Table 7. The graphical illustration of error analysis is presented in Figure 6b. The accuracy level of proposed technique is examined by calculating the L 2 , L ∞ and root mean square (RMS) errors while their descriptions are given below.
Symmetry 2020, 12, x FOR PEER REVIEW 24 of 28     The Tables 8-10 are constructed to present the analysis about L 2 , L ∞ and root mean square (RMS) errors while it is noticed that when increasing the number of nodes (10, 10 2 , 500 and 10 3 ) we are obtaining better rate of accuracy for CPM via D A but the CPM via D E is free of error due to the accuracy of the operational matrix. Hence, we found very efficient and accurate results by means of Chelyshkov polynomial method based on exact operational matrices of derivative (CPMEOM). The set of graphs and tabular forms of comparison are endorsing the credibility and efficiency of the novel matrix−based proposed method and could be extended for the problems arising in engineering and nonlinear sciences.

Conclusions
The current work is dedicated to the construction of a new computational algorithm based on operational matrices of derivatives. The OMDs are obtained through Chelyshkov polynomials and Symmetry 2020, 12, 1471 20 of 22 piecewise functions. A novel mathematical method is proposed depending upon developed matrices and successfully applied for some class of FDEs. The error bound and convergence analysis is reported to validate the consistency of the mathematical structure of CPM.
Hence, some concluding remarks are stated below: • The matrices are derived by means of two different approaches. The operational matrices obtained via Theorem 1 are approximate for non-integer orders of derivative but exact for any integer order derivative while the Theorem 3 provides the exact operational matrices for both inter and non-integer derivatives. The second tactic was extra effectual and consistent.

•
The solutions of the non-integer order problem have been attained via suggested method and very productive results have been attained. In account of obtained results, it is noticed that the CPM has an enhanced accuracy and efficiency rate compared with data in the literature.

•
It is also perceived that computational procedure is a well-organized mathematical tool to scrutinize the fractional order problems. Hence, we developed the mentioned method for a class of nonlinear fractional problems and applications in neural networking systems have been reported. The work is extendable to other applied and natural sciences problems.