Generalized Bessel Polynomial for Multi-Order Fractional Differential Equations

Abstract: The main goal of this paper is to define a simple but effective method for approximating solutions of multi-order fractional differential equations relying on Caputo fractional derivative and under supplementary conditions. Our basis functions are based on some original generalization of the Bessel polynomials, which satisfy many properties shared by the classical orthogonal polynomials as given by Hermit, Laguerre, and Jacobi. The main advantages of our polynomials are two-fold: All the coefficients are positive and any collocation matrix of Bessel polynomials at positive points is strictly totally positive. By expanding the unknowns in a (truncated) series of basis functions at the collocation points, the solution of governing differential equation can be easily converted into the solution of a system of algebraic equations, thus reducing the computational complexities considerably. Several practical test problems also with some symmetries are given to show the validity and utility of the proposed technique. Comparisons with available exact solutions as well as with several alternative algorithms are also carried out. The main feature of our approach is the good performance both in terms of accuracy and simplicity for obtaining an approximation to the solution of differential equations of fractional order.


Introduction
In recent years, fractional calculus has becoming an efficient and successful tool for the analysis of several physical-mathematical problems. The main reason for the increasing number of papers dealing with fractional problems is also explained by the intrinsic and natural possibility of the fractional calculus to take into account also some memory effects, which is quite impossible by using the ordinary differential operators [1]. In this work, we consider the nonlinear multi-order fractional differential equations (MOFDEs) of the form subjected by the following boundary or supplementary conditions H j pxpη j q, x p1q pη j q, . . . , x pm´1q pη j qq " d j j " 0, 1, . . . , m´1, where H j are linear functions, d j P R, and η j are some given points in r0, Ls. In (1), D γ ‹ denotes the Caputo fractional derivative operator such that m´1 ă γ ď m, m P N, and 0 ă β 1 ă β 2 ă . . . ă β ă γ are real constants. The function F can be either linear or nonlinear function of its arguments. In [2], some preliminary results both on the existence and uniqueness of the solution of MOFDEs (1) are obtained.
It is well-known that usually the exact solution of fractional differential equations cannot be obtained analytically. Therefore, many authors have recently developed some suitable numerical methods for such equations. Among the many approximation algorithms for (1) and (2), we mention the systems-based decomposition approach [3], the Adomian decomposition method [4], the spectral methods [5][6][7][8], the B-spline approach [9], and the generalized triangular function [10].
It is known that the traditional orthogonal polynomials such as Jacobi, Hermit, and Laguerre are solution of a second order differential equation. In addition, the derivatives of these polynomials constitute an orthogonal system. Moreover, there exist another set of polynomials with the two aforementioned properties. They satisfy the following second order differential equation where n is a positive integer. Krall and Frink [11] called these the Bessel polynomials thank to their close relation with the Bessel functions of half-integral order. In fact, they have shown that these polynomials occur naturally as the solutions of the classical wave equation in spherical coordinates. These polynomials also appear in the study of various mathematical topics including transcendental number theory [12,13] and student t-distributions [14]. These polynomials seem to have been considered first by Bochner [15] as pointed in Grosswald [16]. However, Krall and Frink considered them in more general setting so that to include also the polynomial solutions of the differential Equation (3). The properties of Bessel polynomials such as recurrence relations, generating functions, and orthogonality were investigated in [11]. The algebraic properties of these polynomials were considered by Grosswald [16]. Some more information about the theory and applications of Bessel polynomials can be found e.g., in [17].
In this research, we establish a new approximation algorithm based upon the Bessel polynomials to obtain a solution of a fractional model (1). In fact, one of our motivation comes from a recent paper [18], which proved the total positiveness of any collocation matrix of theses polynomials evaluated at positive (collocation) points. To the best of our knowledge, this is the first attempt to study these polynomials for approximating MOFDEs. In summary, the main idea behind the presented approximation algorithm based on using the Bessel polynomials with together the collocation points is that it transforms the differential operators in (1) to an equivalent algebraic form, thus greatly reducing the numerical efforts. It should be mentioned that our Bessel polynomials are different from those Bessel functions known as Bessel functions of the first kind that previously considered in the literature, see [19] for a recent review.
The content of the paper is structured as follows. In Section 2 some relevant properties of the Caputo fractional derivative and the generalized Taylor's formula for the Caputo derivative are presented. Section 3 is dedicated to the definitions of Bessel polynomials and their generalized fractional-order counterpart. Moreover, the results about the convergence and error bound of these polynomials are established. In Section 4, where a collocation method also shown to solve the MOFDEs. By using these Bessel basis functions along with collocation points, the proposed method converts the MOFDEs into a nonlinear matrix equation. Hence, the residual error function is introduced to assess the accuracy of Bessel-collocation scheme when the exact solutions are not available. In Section 5, some examples with various parameters together with error evaluation are given to show the utility and applicability of the method. The obtained results are interpreted through tables and figures. Finally, in Section 6, the report ends with a summary and conclusion.

Some Preliminaries
To continue, some definitions and theorems from fractional calculus theory are presented. Definition 1. Let f ptq be a m-times continuously differentiable function. The fractional derivative D q ‹ of f ptq of order q ą 0 in the Caputo's sense is defined as where The properties of the operator D q ‹ can be found in [1]. Besides the linearity, the following properties will be also used t β´q , for β P N 0 and β ě rqs, or β R N 0 and β ą tqu, 0, for β P N 0 and β ă rqs.
Now, we define a generalization of Taylor's formula which involves Caputo fractional derivatives (see a proof in [20]).

Fractional-Order Bessel Functions
In this section, definitions of Bessel polynomials as well as their generalized fractional-order version are introduced. Hence, some properties and convergence results for them are established.
Beside B 0 pxq and B 1 pxq, the next four of these polynomials are listed as follows The coefficients of these polynomials are positive with B n p0q " 1 and B 1 n p0q " npn`1q{2. The explicit expression for the Bessel polynomials as the unique solution of the given differential equation is defined by These polynomials form an orthogonal system with respect to the weight function wpxq " expp´2{xq on the unite circle C, i.e., 1 2πi where δ nm is the Kronecker delta function. Please note that the path of integration is not unique, and it can be replaced by an arbitrary curve surrounding x " 0. The same conclusion is true for the weight function wpxq. This implies that an arbitrary analytic function may be added to wpxq and wpxq may be multiplied by a nonzero constant. By means of the orthogonality relation (9), one may expand a function gpxq in terms of Bessel functions gpxq « 8 ÿ n"0 a n B n pxq, where the coefficients a n are a n " p´1q n`1 pn`1 2 q ż C B n pxq gpxq wpxqdx.

Fractional Bessel Polynomials
The fractional-order Bessel functions can be defined by introducing the change of variable x " t α {L, L, α ą 0 in (8). Let these polynomials will be denoted by B α n ptq " B n pxq. By generalizing these polynomials on the interval r0, Ls we obtain where η " 1 2L . It is not difficult to show that the set of fractional polynomial functions tB α 0 , B α 1 , . . .u is orthogonal on r0, Ls with respect to the weight function w α L ptq " t α´1 expp´2L{t α q. The fractional-order polynomials are useful in particular when the solutions of the underlying MOFDEs have fractional behavior.

Function Approximation and Convergence
Our goal is to obtain an approximate solution for the model problem (1) represented by the truncated Bessel series where the unknown coefficients a n , n " 0, 1, . . . , N must be sought. For this purpose, we express B α n ptq in the matrix representation as where T T T α ptq " and the lower triangular matrix D D D of size pN`1qˆpN`1q takes the form By expressing the relation (11) in a matrix form and exploiting (12), the approximate solution x N,α ptq in the matrix form can be rewritten as where the vector of unknown is A A A " ra 0 a 1 . . . a N s t . Our further aim is to establish the convergence results of the fractional Bessel polynomials. Roughly speaking, the next theorem shows that the approximate solution x N,α ptq converges to the solution xptq of differential Equation (1) as N Ñ 8, see e.g., [21] for a similar proof.

Theorem 2.
Let assume that D kα ‹ gptq P Cp0, Ls for k " 0, 1, . . . , N and let Suppose that g N,α ptq " B B B α ptq A A A is the best approximation out of S α N to g, then the following error bound holds: where M α ě |D Nα ‹ gptq|, t P p0, Ls.
Proof. According to Theorem 1, the generalized Taylor's formula for gptq can be represented as Using the fact that B B B α ptq A A A is the best approximation to g from S α N and G P S α N , we conclude that Employing the inequality´2 L α t α ď´2, which holds for all t P p0, Ls, one immediately find that expp´2 L t α q ď expp´2 L α´1 q. Thus, by inserting this inequality into (14) and then integrating we conclude that The proof is complete by taking the square roots of both sides.
Therefore, for obtaining an approximate solution of the form (11) for the solution of (1) the following collocation points are used on 0 ă t ď L,

The Collocation Scheme
To proceed, we approximate the solution xptq of MOFDEs (1) in terms of pN`1q-terms Bessel polynomials series denoted by x N,α ptq on the interval r0, Ls. In the matrix representation, we consider By placing the collocation points (15) into (16), we get to a system of matrix equations as Hence, we write the preceding equations compactly as where To handle the fractional derivative of order γ in (1), we differentiate both sides of (16), By means of the property (5) and (6), the calculation of D γ ‹ T T T α ptq can be easily obtained as follows To write the fractional derivative D γ ‹ involved in (1) in the matrix form, the collocation points (15) will be inserted into (18) which can be expressed equivalently as where Similarly, the fractional derivative operators D β j ‹ xptq in (1) for j " 1, . . . , can be approximated as where X X X pβ j q and T T T pβ j q are obtained as in (20) by replacing γ with β j . By inserting the collocation points into (1), we have the system Considering these equations in a matrix form and substituting the relations (17), (19), and (20) into the resulting system, a fundamental matrix equation is obtained to be solved. Let us assume that the function F in (21) is the linear form where c k ptq for k " 1, . . . , and c 0 ptq, hptq are given functions. In this case, the equations in (21) can be rewritten in the matrix representation as where the coefficient matrices C C C k , k " 0, 1, . . . with size pN`1qˆpN`1q and the vector H H H of size pN`1qˆ1 have the forms Substituting the relations (17), (19), and (20) where Obviously, (23) is a linear matrix equation and a n , n " 0, 1, . . . , N are the unknowns Bessel coefficients to be determined.
Next aim is to take into account the initial or boundary conditions (2). For the first condition xp0q " x 0 , we tend t Ñ 0 in (16) to get the following matrix representation For the remaining initial conditions, one needs to calculate the integer-order derivatives d k dt k T T T α ptq, k " 1, 2, . . . , n´1, which strictly depend on α as well as N. For example, by choosing α " 1{2 and N " 7 we get Differentiation twice with respect to t reveals that Now, by differentiating k times in (16), and defining with the limit t Ñ 0, we conclude for k " 1, 2, . . . , n´1 that Similarly, for the end conditions x pkq pLq " x Lk , k " 0, . . . , n´1, the following matrix expressions are obtained Now, we replace the first n rows of the augmented matrix rW W W; H H Hs in (23) by the row matrices r p X X X k ; x k s or r p X X X Lk ; x Lk s, k " 0, 1, . . . , n´1 to get the (nonlinear) algebraic system of equations Thus, the unknown Bessel coefficients in (16) will be known through solving this (nonlinear) system. This can be obtained by using the Newton's iterative algorithm.

Remark 1.
In numerical applications below, we frequently encounter the nonlinear terms like x s ptq for s " 2, 3 . . .. To approximate the nonlinear term x 2 ptq in terms of x 2 N,α ptq, the collocation points (15) will be substituted into x 2 N,α ptq. It can be easily seen that in the matrix representation we have " p X X X X X X.
Using (16), we further express the matrix p X X X as a product of three block diagonal matrices as follows Analogously, the higher-order nonlinear terms can be treated recursively X X X s " p p X X Xq s´1 X, s " 3, 4, . . ..

Error Estimation
In general, the exact solution of most MOFDEs cannot be explicitly obtained. Thus, we need some measurements to test the accuracy of the proposed scheme. Since the truncated Bessel series (11) as an approximate solution is satisfied in (1), our expectation is that the residual error function denoted by R N,α ptq becomes approximately small. Here, R N,α ptq : r0, Ls Ñ R obtained by inserting the computed approximated solution x N,α ptq into the differential equation (1). More precisely, for testing accuracy of some numerical models we calculate It should be noticed that the fractional derivatives of order γ, β j , j " 1, . . . , of the approximate solution x N,α ptq in (24) are calculated by using the properties (5) and (6). Obviously, the residual function is vanished at the collocation points (15), so our expectation is that R N,α ptq Ñ 0 as N tends to infinity. This implies that the smallness of the residual error function shows the closeness of the approximate solution to the true exact solution.

Illustrative Test Problems
Now, we show the benefits of the presented Bessel-collocation scheme by simulating some case examples including various linear and nonlinear initial and boundary value problems. The numerical models and calculations are verified through a comparison with existing computational schemes and experimental measurements. Our computations were carried out using MATLAB software version R2017a.
By employing N " 2 and L " 1, we are looking for an approximate solution in the form x N,α ptq " ř 2 n"0 a n B α n ptq. To this end, we calculate the unknown coefficients a 0 , a 1 , and a 2 . For this example we set α " 1 and the collocation points t 0 " 0, t 1 " 1 2 , t 3 " 0 are used. Using γ " 2 and β 1 " Afterwards, by inserting the obtained coefficients into x 2,1 ptq we get the approximate solution which is the desired exact solution.
The initial conditions are xp0q " 0, x p1q p0q " 0, and x p2q p0q " 2. It can be easily checked that the exact true solution is xptq " t 2 .
Therefore, we get which is obviously the exact solution.
In the next example, we show the advantage of using the fractional-order Bessel functions in the computations.
We first consider N " 5 and α " 1{2. The approximated solution x 5, 1 2 ptq for t P r0, 1s takes the form However, with a lower number of basis functions one can also obtain an accurate result. Using N " 2, α " 5{2 and N " 3, α " 5{6, the following approximations are obtained Moreover, to show the advantage of the presented approach and to validate our obtained approximated solutions, we make a comparison in terms of errors in the L 8 and L 2 norms in Table 1. We compare the Bessel-collocation approach and the Chelyshkov collocation spectral method [7]. In this comparison, we use different N " 1, 2, 3 and α " 1{2, 5{2, 5{6.
Let N " 8 and set α " 1. For γ " 2, β " 1, the approximate solution x 8,1 ptq of the model Problem 4 using Bessel functions in the interval 0 ď t ď 1 is In Table 2, we report the numerical results corresponding to these values of γ, β using different N " 8, 16 evaluated at some points t P r0, 1s. The corresponding absolute errors E N,α ptq :" |xptq´x N,α ptq| are also reported in this table. Moreover, the numerical results based on Haar wavelet operational matrices [22] are given in the last column of Table 2. As can see from Table 2, our approximate solutions agree with the results obtained in [22]. The next observation is that more accurate solutions are obtained if one increases the number of Bessel functions N.  In Figure 1, x 10,1 ptq is plotted when γ " 2 (β " 1) is fixed and different values of β " 0.25, 0.5, 0.75, 1 (γ " 1.25, 1.5, 1.75, 2) are examined. It is observed that as γ and β approached to 1 and 2 respectively, numerical solutions tend to the exact solutions.   [24,25] x p2q ptq`1 2 with initial conditions xp0q " 0, x p1q p0q " 0. The exact solution is xptq " t 3 .
Clearly, the exact solution is a third-degree polynomial. Therefore, we take N " 3 and α " 1, which are sufficient to get the desired approximations. Using the usual collocation points as in Problem 2 and similar to Problem 1, we get the final augmented matrix Solving the resulting linear system, we find Therefore, the approximated solution x 3,1 ptq is obtained as which is the exact solution.
Problem 6. Consider the fractional Riccati equation [23,26] D γ ‹ xptq`xptq´x 2 ptq " 0, 0 ă γ ď 1, on a long time interval with L " 5 and initial condition xp0q " 1{2. When γ " 1, the exact solution is xptq " We calculate the approximated solution x N,α ptq using N " 7 and γ equals to α " 1{4. Thus, we get To validate this solution, we also employ the old fractional-order Bessel polynomials as well as Chelyshkov and Legendre functions from the previous works [26,27] with the same parameters as above. The corresponding solutions take the forms respectively To further compare these collocation schemes based on various polynomials, we calculate the estimated residual errors obtained by the relation (24). The graphs of R N,α ptq on the interval r0, 5s correspond to γ, α " 1{4 and for N " 7 are shown in Figure 2. With respect to Figure 2, it is obviously seen that the residual error functions obtained by the presented Bessel-collocation method are smaller compared to the errors of other polynomial-based numerical collocation schemes. Problem 7. Consider the following nonlinear boundary value problem with variable coefficients [6] x p2q ptq`Γp ‹ xptq´rx p1q ptqs 2 " 2`1 10 t 2 , 0 ă t ă 1, with boundary conditions xp0q " 1 and xp1q " 2. The exact solution of this example is xptq " 1`t 2 .

1.0000000000012118324.
To show the gain of the proposed scheme, we compare our results with the collocation method based on Bernstein operational matrix (BOM) of fractional derivative from [6]. Table 3 reports the errors in L 8 and L 2 norms of the new Bessel-collocation procedure and the errors of the BOM algorithm. This comparison shows the thoroughness of the proposed method. Problem 8. We consider the following initial-value problem of multi-term nonlinear fractional differential equation [6] D γ ‹ xptq`D where 2 ă γ ă 3, 0 ă β 1 ă 1, and 1 ă β 2 ă 2 and the initial conditions are xp0q " 0, x p1q p0q " 0, and x p2q p0q " 0. An easy calculation shows that xptq " t 3 is the exact solution.
A comparison between our collocation scheme at C1 and the method of shifted Jacobi operational matrix (SJOM) [6] with N " 24 is made in Table 4. Besides the cases C1 and C2: pγ, β 1 , β 2 q " p2.000001, 0.000009, 1.000001q, the following values of pγ, β 1 , β 2 q are used in Table 5    Looking at Tables 4 and 5 reveals that our numerical solutions obtained via novel Bessel-collocation method are in excellent agreement with the corresponding exact solutions. Moreover, our proposed scheme is superior compared to the SJOM. Problem 9. We consider the fractional relaxation-oscillation equation [5,6] D γ ‹ xptq`xptq " 0, 0 ă γ ă 2, with the initial condition xp0q " 1. If γ ą 1 we also have x p1q p0q " 0. The exact solution in terms of Mittag-Leffler function is given by xptq " E γ p´t γ q. Here, E γ pzq " ř 8 First, we consider γ " 85{100 and set α equals to γ. We get the approximated solution x N,α ptq using N " 8 terms on r0, 1s as follows In Table 6, we calculate the maximum absolute errors using N " 8 and N " 10. In addition, a comparison is done in this table with the results obtained via SJOM [6]. Looking at Table 6 one can find that the achievement of good approximations to the exact solution is possible using only a few terms of fractional Bessel polynomials. In the next experiments, we investigate the impact of varying γ on the maximum absolute errors while N " 10 is fixed. Table 7 presents the L 8 errors for γ " 0.2, 0.4, 0.6, 0.8 as well as γ " 1.2, 1.4, 1.6, 1.8. In all cases, we exploit α " γ. Comparisons with existing approximation techniques based on operational matrix of fractional derivatives via B-spline functions [9] and shifted Jacobi functions [6] are also carried out in Table 7. Table 7. Comparison of L 8 error norms for N " 10 and various γ in test Problem 9.
The exact solution is xptq " t 3´t2 .
Obviously, these approximations are accurate up to machine epsilon. Table 8 reports the comparison of the absolute errors evaluated at some points t P r0, 1s obtained by the Bessel-collocation method. For comparison, the results obtained by the collocation method (CM) [29] and the reproducing kernel method (RKM) [28] are also shown in Table 8. The comparisons show that the proposed method is considerably more accurate than other methods.

Conclusions
A practical matrix approach based on novel (orthogonal) Bessel polynomials is presented to solve multi-order fractional-order differential equations (MOFDEs). Using the matrix representations of the generalized Bessel polynomials and their derivatives with the aid of collocation points, the scheme transforms MOFDEs to a fundamental matrix equation, which corresponds to a system of (non)linear algebraic equations. To assess the efficiency and accuracy of the presented technique, several numerical examples with initial and boundary conditions are investigated. Comparisons with the exact solutions and with various alternative numerical simulations and experimental measurements have also been made. Based on the experiments, it is found that the numerical approximations are in an excellent agreement, which demonstrate the reliability and the great potential of the presented technique.