Fractional Derivatives, Fractional Integrals, and Fractional Differential Equations in Matlab

differentiators/integrators as well as various fractional-order digital ﬁlters.


Fractional calculus fundamentals 2.1 Special functions
Here we should mention the most important function used in fractional calculus -Euler's gamma function, which is defined as This function is a generalization of the factorial in the following form: This gamma function is directly implemented in the Matlab with syntax []=gamma(). Another function, which plays a very important role in the fractional calculus, was in fact introduced in 1953. It is a two-parameter function of the Mittag-Leffler type defined as (Podlubny, 1999): The Mittag-Leffler function (3) can be expressed in the integral representation as ( 4 ) the contour C starts and ends at −∞ and circles around the singularities and branch points of the integrand. There are some relationships (given e.g. in (Djrbashian, 1993;Podlubny, 1999)): E 2,1 (z)=cosh( √ z), E 2,1 (−z 2 )=cos(z), E 1,2 (z)= e z − 1 z .

Fractional derivatives and integrals
Fractional calculus is a generalization of integration and differentiation to non-integer-order fundamental operator a D α t ,w h e r ea and t are the bounds of the operation and α ∈ R. The continuous integrodifferential operator is defined as : α > 0, 1: α = 0, t a (dτ) α : α < 0.
The three most frequently used definitions for the general fractional differintegral are: the Grünwald-Letnikov (GL) definition, the Riemann-Liouville (RL) and the Caputo definition (Miller and Ross, 1993;Oldham and Spanier, 1974;Podlubny, 1999). Other definitions are connected with well-known names as for instance Weyl, Fourier, Cauchy, Abel, etc.
The GL definition is given as where [.] means the integer part. The RL definition is given as where Γ(.) is the gamma function. The Caputo definition of fractional derivatives can be written as In this chapter we will consider mainly the GL, the RL, and the Caputo definitions. This consideration is based on the fact that, for a wide class of functions, the three best known definitions -GL, RL, and Caputo -are equivalent under some conditions (Podlubny, 1999).
As it was already mentioned, under the homogeneous initial conditions the Riemann-Liouville and the Caputo derivative are equivalent. Let us denote the Riemann-Liouville fractional derivative as RL a D α t t(t) and the Caputo definition as C a D α t f (t), then the relationship between them is: for f (k) (a)=0 (k = 0, 1, . . . , n − 1). The initial conditions for the fractional-order differential equations with the Caputo derivatives are in the same form as for the integer-order differential equations. It is an advantage because applied problems require definitions of fractional derivatives, where there are clear interpretations of initial conditions, which contain f (a), f ′ (a), f ′′ (a),etc.
In addition, for approximation purpose we consider the Laplace transform method of the fractional differintegral a D α t . The Laplace transform of the join operator of order α for the GL, RL and Caputo definitions, under the zero initial conditions is: Laplace transform technique is considered as an efficient way in solving differential equations with integer order and fractional order as well. For differential equations with fractional order, the Laplace transform technique works effectively only for relatively simple equations, because of the difficulties of calculating inversion of Laplace transforms. This problem was solved by applying a numerical inverse Laplace transform algorithms in fractional calculus (Sheng et al., 2011). Three numerical inverse Laplace transform algorithms in Matlab, named invlap(), gavsteh(),andnilt(), were described there and tested. The Laplace transforms for several known Mittag-Leffler type functions are summarized as follows (Magin, 2006;Podlubny, 1999):

Numerical methods for fractional calculus
There are many numerical methods appropriate for the fractional calculus computation. They were described in several works and a good survey can be found in (Vinagre et al., 2000;2003). In the above-mentioned papers are described the time domain and frequency methods in the form of IIR and FIR approximations together with illustrative examples. Some of the mentioned frequency methods in both forms of approximation have been realized as the Matlab routines in Duarte Valerio's toolbox called ninteger (see detail review in (Valério, 2005)). Also created in this toolbox was a Simulink block nid for fractional derivative and integral, where the order of derivative/integral and method of its approximation can be selected.

Grünwald-Letnikov method
For numerical calculation of fractional-order derivatives we can use the relation (13) derived from the GL definition (8). This approach is based on the fact that for a wide class of functions the three definitions -GL (8), RL (9), and Caputo's (10) -are equivalent if f (a)=0. The relation for the explicit numerical approximation of q-th derivative at the points kh, (k = 1,2,...) has the following form (Dorčák, 1994;Podlubny, 1999): where L m is the "memory length", t k = kh, h is the time step of calculation and c (q) j (j = 0, 1, . . . , k) are binomial coefficients. For their calculation we can use for instance the following expression (Dorčák, 1994): The binomial coefficients c (q) j (j = 0,1,...,k) can be also expressed by using a factorial. Writing the factorial as a gamma function (2), it allows the binomial coefficient to be generalized to non-integer arguments. We can write the relations: Obviously, for this simplification we pay a penalty in the form of some inaccuracy. If f (t) ≤ M, we can easily establish the following estimate for determining the memory length L m , providing the required accuracy ǫ: The above-described numerical method is called Power Series Expansion (PSE) of a generating function. It is important to note that PSE leads to approximation in the form of polynomials, that is, the discretized fractional operator is in the form of FIR filter, which has only zeros (Chen and Moore, 2002;Vinagre et al., 2003). The resulting discrete transfer function, approximating fractional-order operators, can be expressed in z-domain as follows: where T is the sample period, PSE{u} denotes the function resulting from applying the power series expansion to the function u, Y(z) is the Z transform of the output sequence y(kT), F(z) is the Z transform of the input sequence f (kT), n is the order of the approximation, and R is polynomial of degree n, respectively, in the variable z −1 ,a n dk = 1, 2, . . . . A Matlab routine dfod2() of this method (17) can be downloaded (see (Petráš, 2003b)).

Fig. 2. Fractional derivatives of function y=sin(t)
In Fig. 2 is depicted the fractional derivative of the sine function for fractional order of derivative 0 < α < 1and0< t < 2π.

Continuous and discrete-time approximation techniques
Other approach can be realized by Continued Fraction Expansion (CFE) of the generating function and then the approximated fractional operator is in the form of IIR filter, which has poles and zeros. Taking into account that our aim is to obtain equivalents to the fractional integrodifferential operators in the Laplace domain, s ±r , the result of such approximation for an irrational function, G(s), can be expressed in the form (Vinagre et al., 2003): where a i s and b i s are rational functions of the variable s, or are constants. The application of the method yields a rational function, which is an approximation of the irrational function G(s).
In other words, for evaluation purposes, the rational approximations obtained by CFE frequently converge much more rapidly than the PSE and have a wider domain of convergence in the complex plane. On the other hand, the approximation by PSE and the short memory principle is convenient for the dynamical properties consideration. For interpolation purposes, rational functions are sometimes superior to polynomials. This is, roughly speaking, due to their ability to model functions with poles. These techniques are based on the approximations of an irrational function, G(s), by a rational function defined by the quotient of two polynomials in the variable s in frequency s-domain: passing through the points (s i , G(s i )),..., (s i+m , G(s i+m )). The resulting discrete transfer function, approximating fractional-order operators, can be expressed as (Vinagre et al., 2003): where T is the sample period, CFE{u} denotes the function resulting from applying the continued fraction expansion to the function u, Y(z) is the Z transform of the output sequence y(kT), F(z) is the Z transform of the input sequence f (kT), p and n are the orders of the approximation, and P and Q are polynomials of degrees p and n, respectively, in the variable z −1 ,andk = 1, 2, . . . . In general, the discretization of fractional-order differentiator/integrator s ±r (r ∈ R) can be expressed by the generating function s ≈ ω(z −1 ). This generating function and its expansion determine both the form of the approximation and the coefficients (Lubich, 1986). In this section, for direct discretizing s r , (0 < r < 1), we will concentrate on the IIR form of discretization where as a generating function we will adopt an Al-Alaoui idea on mixed scheme of Euler and Tustin operators (Al-Alaoui, 1993;1997), but we will use a different ratio between both operators. The mentioned new operator, raised to power ±r, has the form (Petráš, 2003a): where a is the ratio term and r is the fractional order. The ratio term a is the amount of phase shift and this tuning knob is sufficient for most engineering problems being solved. In expanding the above in rational functions, we will use the CFE. It should be pointed out that, for control applications, the obtained approximate discrete-time rational transfer function should be stable. Furthermore, for a better fit to the continuous frequency response, it would be of high interest to obtain discrete approximations with poles an zeros interlaced along the line z ∈ (−1, 1) of the z plane. The direct discretization approximations proposed in this paper enjoy the desired properties.
The result of such approximation for an irrational function, G(z −1 ), can be expressed by G(z −1 ) in the CFE form (Vinagre et al., 2000): where a i and b i are either rational functions of the variable z −1 or constants. The application of the method yields a rational function, G(z −1 ), which is an approximation of the irrational function G(z −1 ). The resulting discrete transfer function, approximating fractional-order operators, can be expressed as: where CFE{u} denotes the continued fraction expansion of u; p and q are the orders of the approximation and P and Q are polynomials of degrees p and q. Normally, we can set p = q = n. A Matlab routine dfod1() for this described method (22) can be downloaded (see (Petráš, 2003a)).
EXAMPLE 2.2. Here we present some results for fractional order r = ±0.5 (half-order derivative/integral). The value of approximation order n is truncated to n = 3 and weighting factor a was chosen a = 1/3. Assume sampling period T = 0.001 sec. For r = 0.5 we have the following approximation of the fractional half-order derivative: A Matlab sequence for the the fractional half-order derivative approximation (23) follows: sys=dfod1(3, 0.001, 1/3, 0.5) bode(sys) step(sys) For r = −0.5 we have the following approximation of the fractional half-order integral: A Matlab sequence for the the fractional half-order integral approximation (24) is the following: sys=dfod1(3, 0.001, 1/3, -0.5) bode(sys) step(sys) The Bode plots and unit step response of the digital fractional-order integrator (24) and the analytical continuous solution of a fractional semi-derivative are depicted in Fig. 4. Poles and zeros of the transfer function (24) lie in a unit circle. For simulation purpose, here we also present the Oustaloup's Recursive Approximation algorithm (Oustaloup et al., 2000). The method is based on the approximation of a function of the form: for the frequency range selected as (ω b , ω h ) by a rational function: using the following set of synthesis formulas for zeros, poles and the gain: where ω h , ω b are the high and low transitional frequencies. An implementation of this algorithm in Matlab as a function ora_foc() is given in (Chen, 2003).

Podlubny's matrix approach
This approach is based on the fact that operation of the fractional calculus can be expressed by matrix (Podlubny et al., 2009). It follows from Podlubny (2000), that the left-sided Riemann-Liouville or Caputo fractional derivative v (α) (t)= 0 D α t v(t) can be approximated in all nodes of the equidistant discretization net t = jτ (j = 0,1,...,n) simultaneously with the help of the upper triangular strip matrix B (α) n as (Podlubny, 2000): Similarly, the right-sided Riemann-Liouville or Caputo fractional derivative v (α) (t)= t D α b v(t) can be approximated in all nodes of the equidistant discretization net t = jτ (j = 0,1,...,n) simultaneously with the help of the lower triangular strip matrix F (α) n (Podlubny, 2000): The symmetric Riesz derivative of order β can be approximated based on its definition as a combination of the approximations (30) and (33) for the left and right-sided Riemann-Liouville derivatives, or using the centered fractional differences approximation of the symmetric Riesz derivative suggested recently by Ortigueira (2006). The general formula is the same (Podlubny et al., 2009): In the first case, the approximation for the left-sided Caputo derivative is taken one step ahead, and the approximation for the right-sided Caputo derivative is taken one step back. This leads to the matrix In the second case (Ortigueira's definition) we have the following symmetric matrix: Both these approximations of symmetric Riesz derivatives give practically the same numerical results and in the case of numerical solution of partial fractional differential equations lead to a well-posed matrix of the resulting algebraic system. Similarly, if in addition to fractional-order time derivative we also consider symmetric fractional-order spatial derivatives, then we have to use all nodes at the considered time layer from the leftmost to the rightmost spatial discretization node. Let us consider the nodes (ih, jτ), j = 0, 1, 2, . . . , n, corresponding to all time layers at i-th spatial discretization node. It has been shown in Podlubny (2000) that all values of α-th order time derivative of u(x, t) at these nodes are approximated using the discrete analogue of differentiation of arbitrary order (Podlubny et al., 2009): In order to obtain a simultaneous approximation of α-th order time derivative of u(x, t) in all nodes, we need to arrange all function values u ij at the discretization nodes to the form of a column vector (Podlubny et al., 2009): In visual terms, we first take the nodes of n-th time layer, then the nodes of (n − 1)-th time layer, and so forth, and put them in this order in a vertical column stack.
The matrix that transforms the vector U nm to the vector U (α) t of the partial fractional derivative of order α with respect to time variable can be obtained as a Kronecker product of the matrix B (α) n , which corresponds to the fractional ordinary derivative of order α (recall that n is the number of time steps), and the unit matrix E m (recall that m is the number of spatial discretization nodes): T Similarly, the matrix that transforms the vector U to the vector U (β) x of the partial fractional derivative of order β with respect to spatial variable can be obtained as a Kronecker product of the unit matrix E n (recall that n is the number of spatial discretization nodes), and the matrix R (β) m , which corresponds to a symmetric Riesz ordinary derivative of order β (Ortigueira, 2006), (recall that m is the number of time steps): Having these approximations for partial fractional derivatives with respect to both variables, we can immediately discretize the general form of the fractional diffusion equation by simply replacing the derivatives with their discrete analogs. Namely, the equation (Podlubny et al., 2009): is discretized as B The initial and boundary conditions must be equal to zero. If it is not so, then an auxiliary unknown function must be introduced, which satisfies the zero initial and boundary conditions. In this way, the non-zero initial and boundary conditions moves to the right-hand side of the equation for the new unknown function.

Ordinary fractional differential equations
A general fractional-order system can be described by a fractional differential equation of the form a n D α n t y(t)+a n−1 D α n−1 t where D γ t ≡ 0 D γ t denotes the Grünwald-Letnikov, the Riemann-Liouville or Caputo fractional derivative (Podlubny, 1999). The corresponding transfer function of incommensurate real orders has the following form (Podlubny, 1999): where a k (k = 0,... n), b k (k = 0,... m) are constant, and α k (k = 0,... n), β k (k = 0,... m) are arbitrary real or rational numbers and without loss of generality they can be arranged as In the particular case of commensurate order systems, it holds that, α k = αk, β k = αk, (0 < α < 1), ∀k ∈ Z, and the transfer function has the following form: With N > M, the function G(s) becomes a proper rational function in the complex variable s α which can be expanded in partial fractions of the following form: where λ i (i = 1,2,...,N) are the roots of the pseudo-polynomial P(s α ) or the system poles which are assumed to be simple without loss of generality. The analytical solution of the system (47) can be expressed as where E μ,ν (z) is the Mittag-Leffler function defined as (3). A fractional-order plant to be controlled can be described by a typical n-term linear homogeneous fractional-order differential equation (FODE) in the time domain a n D α n t y(t)+···+ a 1 D α 1 t y(t)+a 0 D α 0 t y(t)=0 (49) where a k (k = 0, 1, ··· , n) are constant coefficients of the FODE; α k (k = 0, 1, 2, ··· , n) are real numbers. Without loss of generality, assume that α n > α n−1 > ... > α 0 ≥ 0. The analytical solution of the FODE (49) is given by general formula in the form (Podlubny, 1999): where (m; k 0 , k 1 ,... ,k n−2 ) are the multinomial coefficients.
Consider a control function which acts on the FODE system (49) as follows: a n D α n t y(t)+···+ a 1 D α 1 t y(t)+a 0 D α 0 t y(t)=u(t).
By the Laplace transform, we can get a fractional transfer function: The fractional-order linear time-invariant (LTI) system can also be represented by the following state-space model: where x ∈ R n , u ∈ R r and y ∈ R p are the state, input and output vectors of the system and A ∈ R n×n , B ∈ R n×r , C ∈ R p×n ,a n dq = [q 1 , q 2 ,...,q n ] T are the fractional orders. If q 1 = q 2 = ...q n ≡ α, system (53) is called a commensurate-order system, otherwise it is an incommensurate-order system. EXAMPLE 3.1. Let us consider a simple two-term fractional differential equation with zero initial condition in general form: Solution can be obtained by using the Laplace transform method. In the Laplace domain we can express the solution as and by using the useful relations (6) and (12), in the time domain we can write a general solution For obtaining the solution in Matlab we can use the following sequence of commands: clear all; close all; a=2; b=1; alpha=1.5; t=0:0.001:20; % for time step 0.001 and computational time 20 sec y=(1/a) * t.^(alpha). * mlf(alpha,alpha+1,((-b/a) * t.^(alpha))); plot(t,y); xlabel('Time [sec]'); ylabel('y(t)'); EXAMPLE 3.2. Let us consider a three-term fractional differential equation in the form a 2 D α 2 t y(t)+a 1 D α 1 t y(t)+a 0 y(t)=u(t).  (57), one can write where t k = kh (k = 1, 2, . . . , N) and q (α) j are binomial coefficients calculated according to (14). After some rearrangement of the terms in relation (58) we can obtain the numerical solution of the fractional differential equation (57) in the following form: where k = 1,2,...,N for N = T sim /h and where T sim is the total time of the calculation. The above approach is general and can be used for n-term fractional differential equation (44).
In this case we will use a matrix approach and the following syntax of the Matlab functions (Podlubny et al., 2008): and the fractional-order controller designed for desired value of phase margin Φ m = 45 o and gain margin equal to infinity, has the following form (Petráš, 2009): where G o (s) is the transfer function of the open control loop. Transfer function (61) corresponds to the fractional differential equation: The feedback control loop described above can be simulated in Matlab environment by using the approximation technique described in previous section, namely Oustaloup's recursive approximation function ora_foc() for desired frequency range ω b = 10 −2 , ω h = 10 2 ,a n d N = 6.

Partial fractional differential equations
For the solution of partial differential equations with derivatives of any order (integer and non-integer), with time and space fractional derivatives, we can use a matrix approach. EXAMPLE 4.1. Let us consider the following one-dimensional fractional diffusion equation with fractional derivatives with respect both to time and to the spatial variable: with the initial condition y(x,0)=h(x) and the boundary conditions y(0, t)=y(L, t)= K.N o t e t h a t∂ β y(x, t)/∂|x| β is symmetric Riesz derivative. Let us consider the following parameters: α = 1.2, β = 2.1, a 2 = 1, f (x, t)=5, L = 2, T sim = 3secandK = 0. The Matlab code is (for more details see also (Podlubny et al., 2008;2009) When we consider non-zero initial condition, e.g. in the form of auxiliary function y(x,0)= h(x)=x(x − 1) and non-zero boundary conditions in the form of the constant K,w h e r e y(0, t)=y(L, t)=K, we have to modify the following lines in above Matlab code: K=2; ... % initial conditions: k = 1:m; U0 = (k-1). * ((m-1) -k + 1) * h^2 + K; In Fig. 10(a) and Fig. 10(b) are depicted the solutions of the fractional partial differential equation (63) for zero and non-zero initial and boundary conditions, respectively.

Nonlinear fractional differential equations
A general numerical solution of the nonlinear fractional differential equation in the form can be expressed by using the relations (13) and (14) as: For the memory term expressed by the sum, a "short memory" principle can be used. Then the lower index of the sums in relations (64) will be v = 1fork < (L m /h) and v = k − (L m /h) for k > (L m /h), or without using the "short memory" principle, we put v = 1forallk.  (63) for parameters: α = 1.2, β = 2.1, a = 1, and f (x, t)=5 for L = 2andT sim = 3sec EXAMPLE 5.1. In this section, we will demonstrate the proposed numerical solution on the set of three nonlinear fractional differential equations, which are used for describing an economical system. The fractional-order financial system is described as follows (Chen, 2008): where the total order of the system is denoted byq =( q 1 , q 2 , q 3 ), a is the saving amount, b is the cost per investment, and c is the elasticity of demand of commercial market. The state variables are: x(t) is the interest rate, y(t) is the investment demand, and z(t) is the price index.

Conclusion
In this chapter are presented the methods for calculation of the fractional derivative and fractional integral together with methods for solution of the fractional differential equations in Matlab. Those methods are based on the numerical approximation of the fractional derivatives and integral in the continuos time, discrete time and frequency domains. For each mentioned methods are presented illustrative examples with a Matlab code. Moreover, all described problems can be solved also in Simulink environment but such approach was omitted in this chapter. Obviously the described problems could be solved via other methods and also by different Matlab functions but we shown the way how one can create its own Fractional Calculus Toolbox for Matlab from functions, which are freely downloadable from the MathWorks, Inc. Matlab Central File Exchange. It depends on the problem which should be solved in Matlab. Except above-mentioned open source routines, there are very useful and effective Matlab functions for solution and analysis of the fractional calculus problems, which are described in (Monje et al., 2010). Especially helpful is the linear fractional-order differential equations solver fode_sol() and also the set of functions already published in paper . It is worth mentioning that there are many additional useful functions at the Matlab Central File Exchange web site, which are not used in this chapter.
For instance variable-order, distributed-order, and random-order fractional equations, predictor corrector numerical method, impulse and step responses invariant discretization of fractional-order differentiators/integrators as well as various fractional-order digital filters. MATLAB is a software package used primarily in the field of engineering for signal processing, numerical data analysis, modeling, programming, simulation, and computer graphic visualization. In the last few years, it has become widely accepted as an efficient tool, and, therefore, its use has significantly increased in scientific communities and academic institutions. This book consists of 20 chapters presenting research works using MATLAB tools. Chapters include techniques for programming and developing Graphical User Interfaces (GUIs), dynamic systems, electric machines, signal and image processing, power electronics, mixed signal circuits, genetic programming, digital watermarking, control systems, time-series regression modeling, and artificial neural networks.