A numerical scheme for problems in fractional calculus

Abstract. In this paper, a new numerical method for solving the fractional differential equations with boundary value problems is presented. The method is based upon hybrid functions approximation. The properties of hybrid functions consisting of block-pulse functions and Bernoulli polynomials are presented. The Riemann-Liouville fractional integral operator for hybrid functions is given. This operator is then utilized to reduce the solution of the boundary value problems for fractional differential equations to a system of algebraic equations. Illustrative examples are included to demonstrate the validity and applicability of the technique.


Introduction
Fractional differential equations (FDEs) are generalized from integer order ones, which are obtained by replacing integer order derivatives with fractional ones. A history of the development of fractional differential operators can be found in [1,2].
FDEs have attracted considerable interest because of their ability to model complex phenomena such as continuum and statistical mechanics [3] visco-elastic materials [4], solid mechanics [5] and dynamics of interfaces between soft-nanoparticles and rough substrates [6]. In general, most of the FDEs do not have exact solutions. In particular, there is no known method for solving fractional boundary value problems exactly [7].
Most of the work published to date concerning exact and numerical solutions for FDEs is devoted to the initial value problems. The theory of boundary value problems for FDEs has received attention quite recently [7].
The available sets of orthogonal functions can be divided into three classes. The first class includes sets of piecewise constant basis functions (e.g., block-pulse, Haar, Walsh, etc.). The second class consists of sets of orthogonal polynomials (e.g., Chebyshev, Laguerre, Legendre, etc.). The third class is the set of sine-cosine functions in the Fourier series. Orthogonal functions have been used when dealing with various problems of the dynamical systems. The main advantage, of using orthogonal functions is that they reduce the dynamical system problems to those of solving a system of algebraic equations by using the operational matrix of integration.
In recent years, the hybrid functions consisting of the combination of block-pulse functions with Legendre polynomials, Chebyshev polynomials, Taylor series, Lagrange polynomials or Bernoulli polynomials [17][18][19][20][21][22] have been shown to be a mathematical power tool for discretization of selected problems. Among these hybrid functions, the hybrid functions of block-pulse and Bernoulli polynomials have been shown to be computationally more effective [21][22][23]. In 2015, to the best of our knowledge, the Riemann-Liouville fractional integral operator for hybrid of block-pulse functions and Bernoulli polynomials was derived directly for the first time [24,25].
In the present paper, a new numerical method for solving the initial and boundary value problems for fractional order differential equations is presented. The method is based upon Hybrid functions approximation. These hybrid functions, which consist of the hybrid of block-pulse and Bernoulli polynomials are given. We then obtain the Riemann-Liouville fractional integral operator for hybrid of block-pulse functions and Bernoulli polynomials. This operator is then utilized to reduce the solution of the fractional order differential equations to the solution of algebraic equations.

The fractional derivative and integral
There are various definitions of fractional derivative and integration. The widely used definition of a fractional derivative is the Caputo definition and a fractional integration is the Riemann-Liouville definition.
Definition 1. Caputo's fractional derivative of order q is defined as [10] ( where q > 0 is the order of the derivative and n is the smallest integer greater than q. The Caputo fractional derivative is a linear operation, namely where λ and µ are constants. Definition 2. The Riemann-Liouville fractional integral operator of order q is defined as [10] where t q−1 * f (t) is the convolution product of t q−1 and f (t). For the Riemann-Liouville fractional integral we have The Riemann-Liouville fractional integral is also a linear operation. The Caputo derivative and Riemann-Liouville integral satisfy the following properties [10]

Function Approximation
Let f ∈ L 2 [0, 1], the best approximation of f by using hybrid functions is given by where and

Riemann-Liouville fractional integral operator for hybrid of block-pulse functions and Bernoulli polynomials
We now derive the operator I α for B(t) in Eq. (10) given by To obtain I α b nm (t), we use the Laplace transform. By using Eq. (5), we have where µ c (t) is unit step function defined as By taking the Laplace transform from Eq. (13) and using Eq. (7), we get From the definition of Bernoulli polynomials in Eq. (6), we have From Eq. (1), Taking the inverse Laplace transform of Eq. (16) yields [24,25] where

Problem statement
We focus on the following problem: Caputo fractional differential equation with the boundary conditions For this problem we have: where G(t, s) is the Green function given by The existence and uniqueness results for this problem are given in [28].

Numerical method
In this section, we use the hybrid of block-pulse functions and Bernoulli polynomials for solving problem given in Eqs. (18) and (19). Here, we expand f (t) by the hybrid functions as by using Eq. (11), we get Substituting Eqs. (21) and (22) in Eq. (20), collocating the resulting equation at the Newtoncotes nodes t i , given by These equations give N(M + 1) algebraic equations, which can be solved for the unknown vector A T using Newton's iterative method.
We solve this problem by using hybrid functions with N = 1 and M = 5. Let Substituting Eqs. (32) and (33) in Eq. (31), we get  For these values, the exact solution is obtained by the present method.

Conclusion
A general formulation for the Riemann-Liouville fractional integral operator for hybrid of block-pulse functions and Bernoulli polynomials has been given. This operator is used to approximate numerical solution of FDEs. Our numerical findings are compared with exact solutions and with the solutions obtained by some other numerical methods. The solution obtained using the present method shows that this approach can solve the problem effectively.