NUMERICAL COMPUTATIONS OF THE FRACTIONAL DERIVATIVE IN IVPS, EXAMPLES IN MATLAB AND MATHEMATICA

. The paper concerns a numerical method that deals with the computations of the fractional derivative in Caputo and Riemann-Liouville definitions. The method can be applied in time stepping processes of initial value problems. The name of the method is SubIval, which is an acronym of its previous name – the subinterval-based method. Its application in solving systems of fractional order state equations is presented. The method has been implemented into an ActiveX DLL. Exemplary MATLAB and Mathematica codes are given, which provide guidance on how the DLL can be used.


Introduction
The increasing popularity of fractional calculus in science is owed to its many applications e.g. in circuit theory [8,19,20], control theory [16], fractional order filter design [10], electromagnetic field analysis [6] and temperature field computations [3]. Fractional calculus introduces the concept of a fractional order derivative and integral (or, in general, an integro-derivative). Among the many definitions that can be found in literature [9], the paper concerns the application of the most commonly used fractional derivatives, which are the Riemann-Liouville definition [14] and that of Caputo [4]. The first, for   (0, 1), is: For that same range Caputo's definition of the fractional derivative is defined by: x t x (2) Moreover, in the most common occurrencethe fractional derivative appears as a derivative in time, expressing a component with memory.
The various successes in attempts to apply the fractional time derivative are not the only reason why its study has been very popular of late. The appearance of the fractional time derivative introduces a component, which adds some additional computational problems. Recently, however, many methods have been introduced, which allow to deal with the computations of the fractional derivative. Lead researchers mention:  the Adomian decomposition method [13],  the CAS wavelet method [18],  the differential transform method [2],  the Taylor expansion method [7],  the collocation method [17],  in general -Fractional Linear Multistep Methods [12],  methods that base on the fractional difference operator [5,11].
The development of clearly described methods dealing with fractional derivative computations allows to conveniently study the application of the fractional derivative.
The current paper concerns the application of a method called SubIval, which originates from it being called the subinterval-based method in its first introduction in [21]. It is a numerical method, generally dealing with the computation of the fractional derivative in time stepping processes of initial value problems.

Method overview
To better describe the fundamentals of SubIvalthe following notation is used for Caputo's fractional derivative: while for the Riemann-Liouville fractional derivative: where instead of the integration bounds t a and t b the integration where: with L s,i being Lagrange basis polynomials: defined on axes of local variables  s (introduced in order to avoid numerical errors), with c s being the normalizing coefficient: An assumption is made that only  M contains the rightmost time node t = t now (whose variables are treated implicitly). This When extracting the node values from the integrodifferentiations one obtains the relation: , and: )).
Finally, when knowing the Lagrange polynomial coefficients one can use the formula for monomial integrodifferentiation: for Caputo's fractional derivative, where the auxiliary variables: , and: while the incomplete beta function can be computed according to the formula: If the Riemann-Liouville definition is used then the above formulae can also be used with the difference that the following coefficient needs to be added: connecting both of the fractional derivative definitions [20]. In conclusion -when for Caputo's derivative one obtains: , then for the Riemann-Liouville fractional derivative: .
So far it has been explained how to obtain the approximation (12)  Typically the maximum polynomial order p s can be the same for all subinterval pairs. It's variability (through adaptivity) has not yet been discussed and is a subject for future research.  A time stepping solver does not need to contain only one subinterval dynamics processe.g. it is possible to perform two of them in cooperation, leading to a predictor-corrector method [21].

Relation to state equations
This chapter gives directions on how SubIval can be applied to solve fractional order state equations. The following system is considered: where the left-hand side contains a vector of the fractional derivatives of the state variables: When SubIval is applied, for a currently computed state vector x j , in the time instance t j , the left-hand side turns equation (22) into: where a diag is a diagonal matrix filled with the a coefficients obtained for each respective state variable; b on the other hand is a vector containing the b coefficients (i.e. b C for Caputo's definition, b C + b RL for the Riemann-Liouville fractional derivative). Therefore, in a computed time step the state vector values can be obtained by solving the following system of equations:

Example
The simple circuit problem given in figure 3 is considered.

Fig. 3. Fractional order circuit example -RLC circuit with fractional order capacitor and coil
The circuit contains both a fractional order coil and a fractional order capacitor. Both are of the order .
The circuit equations can be given in a matrix form: where the state vector is: is a vector of the remaining useful variables: and v(t) is generally a vector of source time functionsin this case, however, represented by a single function:

ActiveX SubIval DLL
The ActiveX DLL has been built in such a way as to support standard integer, double and string operations between the higher level program and the DLL. This allows it to be used in a larger variety of computational programs even in their older versions.
The DLL has been written in C# and even though nowadays many programs allow the usage of external .NET libraries older software did not support such functionality. For example, MATLAB has had this functionality introduced in version 7.8 (2009a) [15].
The code applying SubIval should have an order as suggested by the following guidelines (assuming an application in MATLAB): first a COM object (in here named hSubIval) needs to be obtained by creating a COM object (note that the SubIval DLL needs to be registered first): hSubIval = actxserver('SubIvalNS.SubIvalLauncher'); now from this moment on all functions are run as methods of the hSubIval object, one can define how many independent subinterval dynamics processes will be run independently by executing the first function, e.g. here one initializes one process: hSubIval.Initialize(1); subsequentlyeach process needs to be initialized by referring to it through its zero based index p_id; the number of variables under fractional derivatives n needs to be given, along with the polynomial order p, the initial time instance t0 and the type of subinterval dynamics process (0 stands for constant step, which is the only one discussed in this paper): hSubIval.
InitializeComputer(p_id, n, p, t0, 0); each of the fractional derivative orders needs to be given by entering its zero based index j and value alfaj: hSubIval.InitAlpha(p_id, alfaj, j); then the initial values need to be entered, also through a zero based index entry: hSubIval.Putxi(p_id, x0j, j); at this point the initialization is complete and the time stepping process can begin; in every new time instance tnext one needs to run: hSubIval.newt(p_id, tnext); to perform the core SubIval computations one needs to execute: hSubIval.ComputeSimple(p_id); in the case of a predictor-corrector schemethe functions are different, this, however, will not be discussed in this paper, after the computations had been performed the a and b coefficients are available; to obtain them one can run: hSubIval.getab(p_id, j); for the variable with the zero based index j, after a solution has been obtained then each state variable value can be entered also with Putxi.

Summary
A numerical method for the computation of the fractional derivative in initial value problems has been discussed. The method is now known as SubIval. Its brief explanation has been given in section 2. In section 3 a simple circuit example has been introduced, which requires dealing with the fractional derivative.
The state equations have been formulated through a matrix method. A short introduction has been given to the functions that need to be executed from the SubIval DLL through MATLAB. Exemplary scripts have been presented in both MATLAB and Mathematica. These scripts allow to obtain the time functions for the state variables of the discussed problem. In this paper only a constant step variant of SubIval has been used. An adaptive time step predictor-corrector scheme cooperating with the SubIval DLL will be discussed in a future paper.