Computation of numerical solutions to variable order fractional di ﬀ erential equations by using non-orthogonal basis

: In this work, we present some numerical results about variable order fractional di ﬀ erential equations (VOFDEs). For the said numerical analysis, we use Bernstein polynomials (BPs) with nonorthogonal basis. The method we use does not need discretization and neither collocation. Hence omitting the said two operations su ﬃ cient memory and time can be saved. We establish operational matrices for variable order integration and di ﬀ erentiation which convert the consider problem to some algebraic type matrix equations. The obtained matrix equations are then solved by Matlab 13 to get the required numerical solution for the considered problem. Pertinent examples are provided along with graphical illustration and error analysis to validate the results. Further some theoretical results for time complexity are also discussed.


Introduction
Calculus of arbitrary order derivatives and integrations has been found very applicable in mathematical modeling of various real world problems. Such differential and integral operators have greater degree of freedom. Therefore its dynamical behavior is global instead of local. Also various features related to represent memory and hereditary process can be comprehensively explained through the said area [1][2][3]. Therefore in last few decades Fractional order differential and integral equations (FODIEs) has got proper attention from researchers. Interesting applications regarding FODIEs have been investigated in variety of disciplines including engineering, technology and applied sciences [4][5][6]. Numerous features related to memory and hereditary characteristics corresponding to different materials and processes can be found [7]. Due to these applications FODIEs have appealed researchers for further analysis and progress.
Since arbitrary order differential equations are infact definite integrals which include the classical derivative as a special case. Therefore, arbitrary order derivative has not a unique definition. It is interesting that researchers have introduced various definitions for the said differential operators. Geometrically fractional order derivative provides a complete spectrum of a function which include the classical order curve as a special case. Among all the definitions of fractional derivative, the definitions of Caputo and Riemann-Liouville has gained more popularity among the researchers. On the other hand various aspects including qualitative analysis, stability theory and numerical treatment have been studied for FODIEs. For such investigations the authors have used fixed point approach, tools of nonlinear analysis to handle the mentioned aspects for various problems of FODIEs [8]. One of the most important area in the theory of FODIEs is known as numerical solutions of the mentioned area. In variety of situations it is difficult to obtain analytical solution of numerous FODIEs due to the complex behavior of fractional order. In such conditions approximating the solutions to that problem will be more suitable to be determined. For this purposes, large numbers of analytical and numerical methods have been established. In this regards, eigen function procedure [9], perturbation tools [10], iteration techniques [11], transform methods [12], decomposition schemes [13] have been established for analytical or semi-analytical results of FODIEs. Also for numerical solutions, difference methods [14], numerical method for multi-terms FODEs [15], Tau method [16], collocation techniques [17], wavelet analysis [18] have been introduced in literature. Also spectral methods based on operational matrices [19] have been constructed in large numbers. The aforementioned methods have been very well applied for classical differential and integral equations. Here we remark that the said tools have also applied in the area of fractional calculus very well. Among the mentioned methods, spectral techniques are the most powerful and significant for the numerical solutions to various problems. The said methods are based on some operational matrices of integration and differentiations. Based on these matrices, the proposed problem is converted to some algebraic equation known as Sylvester equation of the form where L, M, N are constant co-efficient matrices and X is an unknown matrix that has to be determined. Further a Sylvester equation is a linear matrix equation having a unique solution if the spectra of L and M has an empty intersection. If X has order m × n, then the order of matrix L will be m × n and that of M will be n × n, where the order of X and N will be m × n. Sylvester linear matrix equation plays very important role in control and stability theory of many applied problems. Keeping in mind the importance of numerical solutions to fractional order problems, various powerful techniques have been established in last two decades. For instance authors [20] have developed numerical scheme for two dimensional stochastic Volterra-Fredholm integral equations. In same line the authors [21] have developed a hybrid numerical scheme for partial-integro type problems. Further some more advanced numerical methods have been recently established for ordinary and partial fractional order problems. For the mentioned problems the authors have used various procedure including BPs and Bernoulli and some other complex tools. For the mentioned numerical methods we refer [22][23][24][25][26][27][28].
In present literature, the operational matrices have been developed by using different orthogonal polynomials like Legendre, Laguerre, Jacobi, etc. All the operational matrices have been obtained from the aforesaid polynomials by using descritization techniques to solve FODIEs subject to some initial or boundary conditions [29]. However, boundary value problems are rarely investigated which constitute a very applicable branch of applied analysis. Also some new operational matrices for BPs have been developed recently [30]. Keeping in mind that the said polynomials are non orthogonal, some researchers have used discretization together with collocation method to construct some operational matrices for arbitrary order differentiation and integration. But the mentioned approach is limited to only initial value problems. Recently some simple problems of fractional order differential equations (FODEs) have been investigated by using the said nonorthogonal polynomials. Since in most of the literature discretization and collocation techniques have jointly used which exploit extra memory and consume much time. For some frequent results in this regards see [31,32]. Therefore to omit discretization and collocation to save memory and time from wastage, we will directly construct the operational matrices for solving FODIEs numerically with some initial/boundary conditions by using BPs. Authors [33] have started the area of variable order integration and differentiation during 1993. Recently the mentioned area of variable order FODEs has attracted much attention. The reason is that such problems have more degree of freedom in choosing the most suitable order for the accurate description of a real world problems. Currently many valuable articles have been published in this regards (see [34][35][36][37][38]). Recently various articles related to fractional order dynamics of epidemiological disease, neural network and PD controller theory have been published (see [39][40][41][42][43]). In respect of numerical analysis BPs give more accurate results as compared to other polynomials. Because Bernstein polynomials have non orthogonal basis. For variable order the said polynomials are very rarely used. Further BPs have been developed by a Russian mathematician Sergie Natanovich Bernstein. According to the Weierstrass Approximation Theorem every continuous real valued function can be approximated uniformly with the help of polynomial function over R. In this regard, BPs play a very important role in function approximation. Though the BPs are non-orthogonal, but good approximation for real valued continues functions. BPs play significant roles in distribution functions theory. In recent time BPs estimators of density functions have got great popularity from researchers of statistics. Some authors have investigated various problems of estimating a multivariate distribution function by using BPs in multiple dimensions (see [44]).
Since it has been proved that spectral methods are stable and convergent. So far we know these type of methods have been used for traditional fractional order derivatives very well. But in case of variable order where the differential operators are more flexible and posses greater degree of freedom. The spectral methods have not so properly applied in past many years. Also the mentioned tools have been shamefully applied in many disciplines to perform simulations like heat conduction, quantum mechanics, fluid dynamics, weather prediction and so on (see [45]). Further the proposed method has some advantages like more accurate than finite difference method with the same number of degrees of freedom, can be used spatial filters of very high order easily, and obtain power spectra directly. Also the present method does not required any prior discretization of data like need by finite methods, Galerkin and wavelet methods. Since these methods are scaling techniques depend on scale level. So for larger scale level it need more time for compilation but produce more significant results (see [46]). Motivated by the aforesaid work, here we consider the following two cases of VOFDEs including initial and boundary value problems of the form where in both cases g : [0, 1] × R → R is linear continuous function. Here initially, we will construct operational matrices of fractional integration and differentiation by using BPs. Afterwards, with the help of developed matrices, we will transfer the proposed problems to algebraic equations of Sylvester type. The proposed problem (3.1) and (1.3) include some problems as special case like: ) becomes a decay model which has various applications in radioactivity process of various elements.
to two point boundary value problems been studied in [47].
We use the computational software Matlab 13 to find the unknown matrices for the required numerical solution. Also the results are demonstrated graphically. A comparison between integer order and numerical solutions at various fractional order is also given. Further some theoretical analysis in the time complexity is provided.

Elementary results
Here we recall some basic results from [48].
Definition 2.1. If s : [0, ∞) → R is a function whose integral converges overR + , then its fractional integral of order ϑ(t) > 0 is given by
BPs are formed by the linear combination of Bernstein basis [27]. The polynomial is defined as where k is the order of the polynomial. The compact form of BPs [29] is obtained by using binomial expansion as The set of BPs of degree m is expressed as is known as the Bernstein basis. Further some needful results are given bellow. Since, BPs are non-orthogonal, In this regard the inner product of two Bernstein basis is given by (2.6) Which upon representation in matrix form can be written as This results into For the given two functions φ, ψ ∈ L 2 [0, 1], we can approximate any function φ(t) in term of Bernstein basis as where j = 0, 1, 2, · · · . Hence we get G 0,0 G 0,1 · · · G 0,r · · · G 0,m G 1,0 G 1,1 · · · G 1,r · · · G 1,m G 2,0 G 2,1 · · · G 2,r · · · G 2,m . . . . . . · · · . . . · · · . . .
where S 1xm = s 0 s 1 · · · s m is the coefficient matrix and as (2.7) we provide the following two lemmas for our algorithm.
Lemma 2.4. [30] Let P i,m (t) be function vector defined in (2.4) the then the fractional order integration over the function is given as

Numerical scheme
Here we divide this section to two subsection. In one subsection we establish numerical scheme for initial values problems. In second subsection, we derive the scheme for boundary value problems under variable order.

Numerical algorithm for initial value problem
Initial value problems have their own importance in the field of mathematics. We discuss different cases of initial value problems to investigated for numerical analysis. Case I: . Applying ϕ(t) order integral on both sides, we get

By using function approximation consider
Hence using v(t) in the given (3.1) we can get required matrix equation. Case II: If we take g(t, v) = v(t), then one has and applying ϕ(t) order integral, we get the following function v(t) with the given initial condition

4)
from (3.6), we can write The given Eq (3.7) is a Sylvester type equation that can be solved by Matlab. Case III: If the given differential equation is non-homogeneous then the given linear function will be as then the given differential equation will be as In such case using the function v(t) in Eq (3.4), we get After rearrangement we have which implies that . Then taking ϕ(t) order integral of both sides yields Taking first order derivative of (3.13), we get applying initial condition we have, By using values of v 0 and v 1 in (3.13), we get Hence (3.14) The above equation is used to solve the (3.11) with the help of Matlab. Case V: Let the Eq (3.11) be considered as Let us assume that c 0 D ϕ(t) t v(t) = A M P T M (t) and applying ϕ(t) order integral, we get the following function v(t) , Taking derivative of (3.16), we get Applying initial conditions in (3.16) implies Then the given Eq (3.16) become, On using Lemma 2.4 and Lemma 2.5, one has from above equation we can have The given Eq (3.18) is the transformed Sylvester equation that can be solved by Matlab. Case VI: If the given differential equation is non-homogeneous then the given linear function will be as, and the given differential equation will be as (3.19) in such case using the function v(t) reference to Eq (3.17) we get it implies

Thus we have
Hence (3.21) is the required Sylvester equation.

General algorithm for boundary value problems
Here we establish algorithm for boundary value problems. Case I: We consider the following problem as: we will approximate the given VOFDEs in Bernstein Basis as per previous practice. Let us assume that applying ϕ(t) order integral on both sides of (3.23), we get applying initial conditions given in (3.22), we can find the values of d 0 and d 1 as Hence we get v 0 = d 0 . Similarly the value of d 1 is found in the same manner by applying the boundary putting the values of d 0 and d 1 in (3.24) we can find our function v(t) as where R M is the coefficient matrix. We proceed as Therefore, the given scalar problem in (3.22) is transformed to , which is the required matrix equation that can be solved for unknown matrix C M with Matlab. Case III: If in (3.22) the linear function g : [0, 1] → R is selected in a way that it become a nonhomogeneous differential equation, that is , let g(v(t)) = Kv(t) + f (t), then the given problem become a non-homogenous VOFDE as (3.34) applying fractional integral of order ϕ(t), we get with the given initial conditions We take some approximations (3.37) where by using Lemma 2.4 and Lemma 2.5, we get required approximations as Thus the Eq (3.36) becomes Also the given Eq (3.33) becomes From the Eq (3.39), we have The given equation in (3.40) is the required Sylvester equation, which can be solved for the unknown matrix C M by using Matlab.

Numerical examples
Here we enrich this part by providing some examples. In this section we provide some examples for the demonstration of our proposed single valued problems of fractional order via using the suggested scheme.
The exact solution at fixed fractional order ϕ ∈ (1, 2] of the problem is given by We approximate the solution at various scale level and different fractional order by using the proposed method. In Figures 1 and 2, we give the detail graphically. The exact solution of the problem at fixed order ϕ(t) = 2 is given by We approximate the solution at various scale level and different variable fractional order by using the proposed method in Figures 3 and 4 respectively as:  From Figure 3, we see as ϕ(t) → 2 the solution converges to exact solution at order 2. Also as we increase scale level the absolute error is decreasing and greater the scale level higher be the accuracy and vice versa in Figure 4.
The exact solution of the problem at ϕ(t) is given by We approximate the solution at various scale level and different variable fractional order by using the proposed method.  From Figure 5, we see as ϕ(t) → 2 the solution converges to exact solution at order 2. Also as we increase the scale level, the absolute error is decreasing and greater the scale level higher be the accuracy and vice versa. This process can be seen in Figure 6.  Table 1. The CPU time complexity has been recorded in seconds.

Conclusions and Discussion
We have established an algorithm for VOFDEs including both initial and BVPs. BPs have been applied to construct operational matrices for variable order derivative and integrations. Based on these matrices we have converted our considered problems to some algebraic type equations. By using Matlab we have solved the concerned algebraic equation to get the required numerical solution. VOFDEs have greater degree of freedom in choosing of order for more accurate results. Hence these can be used as a powerful tool to describe various real world problems more accurately as compared to constant fractional order. By graphical presentation we have shown the validity of our results. Further we see that our method is scaling based method greater the scale more will be the accuracy and vice versa. Further our method has been omitted discretization and collocation which help in saving memory from wasting. We have also computed some theoretical analysis of time complexity in Table 1 which indicates that as the scale level is enlarging, the corresponding time is also increase. This is the weakness of the method but not a major drawback as the precision and accuracy is increasing. In future the above technique will be used for numerical treatment of various dynamical systems under different kinds of fractional order operators.