A Review on Fractional Differential Equations and a Numerical Method to Solve Some Boundary Value Problems

Fractional differential equations can describe the dynamics of several complex and nonlocal systems with memory. They arise in many scientific and engineering areas such as physics, chemistry, biology, biophysics, economics, control theory, signal and image processing, etc. Particularly, nonlinear systems describing different phenomena can be modeled with fractional derivatives. Chaotic behavior has also been reported in some fractional models. There exist theoretical results related to existence and uniqueness of solutions to initial and boundary value problems with fractional differential equations; for the nonlinear case, there are still few of them. In this work we will present a summary of the different definitions of fractional derivatives and show models where they appear, including simple nonlinear systems with chaos. Existing results on the solvability of classical fractional differential equations and numerical approaches are summarized. Finally, we propose a numerical scheme to approximate the solution to linear fractional initial value problems and boundary value problems.


Introduction
Fractional calculus is the theory of integrals and derivatives of arbitrary real (and even complex) order and was first suggested in works by mathematicians such as Leibniz, L'Hôpital, Abel, Liouville, Riemann, etc. The importance of fractional derivatives for modeling phenomena in different branches of science and engineering is due to their nonlocality nature, an intrinsic property of many complex systems. Unlike the derivative of integer order, fractional derivatives do not take into account only local characteristics of the dynamics but considers the global evolution of the system; for that reason, when dealing with certain phenomena, they provide more accurate models of real-world behavior than standard derivatives.
To illustrate this fact, we will retrieve an example from [1]. Recall the relationship between stress σ t ð Þ and strain ε t ð Þ in a material under the influence of external forces: is the Newton's law for a viscous liquid, with η the viscosity of the material, and is Hooke's law for an elastic solid, with E the modulus of elasticity. We can rewrite both Eqs. (1) and (2) as with α ¼ 0 for elastic solids and α ¼ 1 for a viscous liquid. But, in practice, there exist viscoelastic materials that have a behavior intermediate between an elastic solid and a viscous liquid, and it may be convenient to give sense to the operator d α dt α if 0 , α , 1.
There exist various definitions of fractional derivatives. All of them involve integral operators with different regularity properties, and some of them have singular kernels.
Next, we will briefly review the most frequently fractional derivatives cited in the bibliography (see [2] for a more complete review and [1,[3][4][5][6] for rigorous theoretical expositions and calculation methods).
The classical Cauchy formula for the n-fold iterated integral, with n ∈ N, is Recalling that gamma function verifies nΓ n ð Þ ¼ n!, an immediate generalization of this formula for a real order α is known as Riemann-Liouville fractional integral operator of order α (the term "fractional" is misleading but has a historical origin). From this, the Riemann-Liouville fractional derivative of order α, with n À 1 , α , n, is defined as while the Caputo fractional derivative of order α is defined as Both Eqs. (6) and (7) define left inverse operators for the integral operator of Riemann-Liouville of order α and are associated to the idea that "deriving α times may be thought as integrating n À α times and deriving n times." Of course these definitions aren't equivalent: clearly the domains of the operators RL 0 D α t : ½ and C 0 D α t : ½ are different; because of the different hypothesis about f , we need to impose to guarantee their existence. Besides that, with the appropriate conditions for f , More recently, the Caputo-Fabrizio fractional derivative of order α, with α ∈ 0; 1 ½ Þ, was defined as for M α ð Þ is a normalizing factor verifying M 0 ð Þ ¼ M 1 ð Þ ¼ 1. Let us point out that, both in Eqs. (7) and (9), the lower limit in the integral could be changed by any value a ∈ À∞; t ½ Þ, i.e., In [7] the authors prove that the operator defined in Eq. (9) verifies the following (convenient) properties: Caputo-Fabrizio definition was then generalized by Atangana and Baleanu, who gave the following definition of the Atangana-Baleanu fractional derivative in Riemann-Lioville sense: and the Atangana-Baleanu fractional derivative in Caputo sense replacing the exponential by E α z ð Þ ¼ ∑ ∞ k¼0 z k Γ αkþ1 ð Þ , the generalized Mittag-Leffler function.
Other types of fractional derivatives are Grünwald-Letnikov's, Hadamard's, Weyl's, etc. In every definition it is clear that fractional derivative operators are not local, since they need the information of f in a whole interval of integration. When defined with 0 as lower limit of integration, as we did, function f is usually assumed to be causal (i.e., f t ð Þ 0 for t , 0), but this limit can also be changed. Authors choose one definition or the other depending on the real-world phenomena they need to model; the scope of application of each operator is still unknown, and, in relation to some of them, there is neither an agreement about whether it is appropriate or not to call them derivatives (see, for a discussion on this topic, [8,9]) nor what are the criteria to decide on it ( [10,11]).
Caputo and Fabrizio ( [12]) proposed the following terms to recognize if an integral operator merits to be called a fractional derivative: 1. The fractional derivative must be a linear operator.
2. The fractional derivative of an analytic function must be analytic.
3. If the order of the derivative is a positive integer, the derivative must be the classical one.
4.If the order is null, the original function must be recovered.
Many applications of fractional calculus have been reported in areas as diverse as diffusion problems, hydraulics, potential theory, control theory, electrochemistry, viscoelasticity, and nanotechnology, among others (see, e.g., [13,14], for a profuse listing of application areas). In Section 2 we will briefly exemplify a few of these applications, in quite different fields, and in Section 3, we will even mention some cases of fractional nonlinear systems which exhibit chaotic behavior.
Theoretical results concerning existence and uniqueness of solutions to fractional differential equations have been also developed.
In [15][16][17] the authors state conditions to guarantee the existence and uniqueness of solution to problems like for 0 , α , 2. After rewriting the equation as an integral equation with a kernel whose norm is bounded in a proper Banach space, they use generalizations of the fixed-point theorem. The function F, besides being continuous, must satisfy certain conditions that substitute the classical Lipschitz's one.
Similar results are stated in [18] for a coupled system of fractional differential equations involving Riemann-Liouville derivatives.
Analytical calculus of fractional operators is, in general, difficult. In [19][20][21] a few examples of quite different analytical methods are presented.
In [22,23], existence and uniqueness, for the solution to a simple case, are proved, and explicit formulae are presented when g is continuous, causal, and null at the origin. The case of Caputo derivative is also considered. In all cases, the computation of the primitive of the data function is required.
Numerical methods have also been proposed to obtain approximate solution to fractional differential equations.
In [6] some numerical approximations to solutions to different fractional differential equations are presented and experimentally verified on various examples, and in [24,25] complete surveys on numerical methods are offered. But numerous articles appear continuously with new approximation methods: we finish this section commenting briefly some works on numerical methods of quite different nature.
In [26] a local fractional natural homotopy perturbation method is proposed to found the solution to partial fractional differential equations as a series.
A method based on a semi-discrete finite difference approximation in time, and Galerkin finite element method in space, is proposed in [27] to solve fractional partial differential equations arising in neuronal dynamics.
In [28] a new numerical approximation of Atangana-Baleanu integral, as the summation of the average of the given function and its fractional integral in Riemann-Liouville sense, is proposed.
Semi-discrete finite element methods are introduced in [29] to solve diffusion equations, and implicit numerical algorithms for the case of spatial and temporal fractional derivatives appeared in [30]. A high-speed numerical scheme for fractional differentiation and fractional integration is proposed in [31]. In [32], a new numerical method to solve partial differential equations involving Caputo derivatives of fractional variable order is obtained in terms of standard (integer order) derivatives.
In [33], a discrete form is proposed for solving time fractional convectiondiffusion equation. The Laplace transform is used to solve fractional differential equations in [34].
Finally, in Section 4, we will present a numerical method we have developed, based on wavelets, to solve initial and boundary value problems with linear fractional differential equations.

Some mathematical models with fractional derivatives
The purpose of this section is to highlight the role of fractional derivatives for modeling certain real evolution processes. We enumerate several mathematical models of different fields, found in the recent literature. In each of them, the fractional order of derivation is justified by the nature of the phenomenon that is described. Usually, in the papers, both the symbols ∂ α ∂t α and D α t are used to indistinctly represent any of the fractional derivatives, whose type is clarified in the text.
In [35] the authors review the evolution of the general fractional equation: for a . 0, 0 , β ≤ 2, where x ∈ S ⊂ R and t ∈ R . 0 denote the space and time variables. This equation is obtained from the classical D'Alembert wave equation by replacing the second-order time derivative with the Caputo fractional derivative of order β ∈ 0; 2 ð : The authors show that, for 1 , β , 2, the behavior of the fundamental solutions turns out to be intermediate between diffusion (for a viscous fluid) and wave propagation (for an elastic solid), thus justifying the attribute of fractional diffusive waves.
In [36] another approach for time fractional wave (Eq. (19)) is proposed. It is solved, for 0 , β , 1 and special initial conditions, by the method of separation of variables.
In [37] the authors study the particular linear fractional Klein-Gordon equation: considering the fractional Caputo derivative with 1 , α ≤ 2. They achieve a numerical solution by variational iteration method and multivariate Padé approximation.
In [38] the following mathematical model, using Fick's law of diffusion, is developed to study the effect of fractional advection diffusion equation (cross flow) for the calcium profile, considering the Caputo fractional derivative (0 , α ≤ 1): Here, the calcium concentration C x; t ð Þ varies in time and space, D is a diffusion constant, and C ∞ is calcium concentration at infinity and is assumed that, at initial state of time and at a long distance, calcium concentration vanishes or becomes zero. The authors note that the physical parameter α characterizes the cytosolic calcium ion in astrocytes.
In [39] the authors explain that arteries, like other soft tissues, exhibit viscoelastic behavior and part of the mechanical energy transferred to them is dissipative (viscosity) and the other part is stored in a reversible form (elasticity). They modified the standard model by a fractional-order one and test it in human arterial segments. They conclude that fractional derivatives, in Riemann-Liouville sense, are a good alternative to model arterial viscosity. The generalized Voigt model consists of a spring in parallel with two springpots (a neologism for a model that is between a spring-purely elastic-and a dashpot, purely viscous) of fractional orders α and β. The governing fractional-order differential equation is where E 0 is the elastic constant for a spring and η 1 and η 2 represent the viscosities of two springpots in parallel with the spring.
In Ref. [40] the authors developed an accurate and efficient numerical method for the fractional-order standard model described by Eq. (22) with α ¼ β, combining a fast convolution method with the spectral element discretization based on a general Jacobi polynomial basis that can be used to generate 3D polymorphic high-order elements. In that way they model complicated arterial geometries, such as patient-specific aneurysms, and apply it to 3D fluid-structure interaction simulations.
In [12] the authors use fractional derivatives to model the magnetic hysteresis, a phenomenon where the "memory" of the ferromagnetic material is crucial. They use a nonlinear model for the constitutive law of an isotropic ferromagnetic material: for λ . 0, H t ð Þ is the magnetic excitation field, M t ð Þ is the magnetization vector, θ t ð Þ is the temperature, θ c is the Curie temperature below which the hysteresis is observed, and C 0 is the tensor with the constitutive properties of the magnetic material. They compare the resulting behavior when D α t is the Caputo fractional derivative with the one that results when the derivative is the Caputo-Fabrizio one. By numerical simulations they obtain examples of the classical hysteresis cycles and conclude that Caputo derivative expresses a stronger memory than the Caputo-Fabrizio operator.
These are just a few examples of the huge variety of problems that can be modeled by means of fractional differential equations. The nonlocality of the associated operators is the key to the success in the description of these phenomena.

Some simple systems exhibiting chaos
Chaos theory is also an area where fractional derivatives play an important role. In this section we comment on some nonlinear systems modeled with fractional derivative, recently published, that exhibit chaos.
In [41] the authors studied a system based on the classical Lorenz one, but described by the Atangana-Baleanu fractional derivative (in the Caputo sense) with 0 , α , 1: Under certain assumptions on the physical problem, they proved existence of solutions, and, by means of an iterative algorithm, numerical evidence of chaos is shown when 0:25 , α , 0:3 and 0:4 , α , 0:5 for the usual set of parameters σ ¼ 10, ρ ¼ 8 3 , and β ¼ 28. In [42] a three-dimensional fractional-order dynamical system for cancer growth is proposed replacing the standard derivatives in the evolution equations: by the Caputo-Fabrizio and the Atangana-Baleanu (Caputo sense) derivatives. The system parameters are related to the rate of change in the population of the different cells: healthy and tumor ones. The authors prove that the system has a unique solution and show that the system exhibits chaos for a proper choice of the parameters values and initial conditions.
In [43] a fractional Lorenz system is studied considering the generalized Caputo derivative, defined, for 0 , α , 1, ρ . 0, as A detailed analysis of the stability of the system is performed. Adomian method is used to find semi-analytical solution to the fractional nonlinear equations. Chaotic behavior and strange attractors are numerically found for some values of α and ρ:

A numerical approximation scheme to solve linear fractional differential equations
After having exemplified several applications of fractional differential equations to different real-world problems, including chaotic ones, we will show a method that we have developed to obtain numerical solutions to linear fractional initial value and boundary value problems modeled with Caputo or Caputo-Fabrizio derivatives. The idea of the approximation scheme is to transform the derivatives into integral operators acting on the Fourier transform and to perform a wavelet decomposition of the data. The wavelet coefficients of the unknown are then recovered from a linear system of algebraic equations, and the solution is built up from its coefficients. The properties of the chosen wavelet basis guarantee numerical stability and efficiency of the approximation scheme. In the case of singular kernel, this procedure enables us to handle the singularity.
We note that choosing a ¼ À∞ in definition of Eqs. (4) or (5) and being f ∈ H 1 À∞; b ð Þ, the Sobolev space of functions with (weak) first derivative in L 2 À∞; b ð Þ, both derivatives can be expressed as a convolution. For the Caputo-Fabrizio fractional derivative of order 0 , α , 1, changing variables in Eq. (9), we have where k is a causal function, k t ð Þ ¼ e À αt 1Àα for t ≥ 0, and k t ð Þ ¼ 0 for t , 0.
Using the properties of the Fourier transform, we can rewrite the last equality:

An initial value problem
Let us consider the initial value problem (IVP): We look for f satisfying Eq. (31), where g is a causal and smooth function with g 0 ð Þ ¼ 0. Other situations where the initial condition is not null can also be faced adapting the following scheme. We will consider that the fractional derivative is the Caputo-Fabrizio one; the Caputo case can be solved similarly (see [44] for a detailed description). First we choose a wavelet basis with special properties: well localized in both time and frequency domain, smooth, band limited and infinitely oscillating with fast decay. The mother wavelet ψ ∈ S (the Schwartz space) and its Fourier transform satisfy supp is an orthonormal basis (BON) of L 2 R ð Þ associated to a multiresolution analysis (MRA). We denote, by the wavelet and scale subspaces, respectively, and decompose the space 46]). The sets Ω jÀ1 , Ω j , Ω jþ1 have little overlap, and W j is nearly a basis for the set of functions whose Fourier transform has support in Ω j . This property of the basis will be crucial in the procedure. Now we decompose the data g as where the first and second terms are the projections of g in V J and W j , respectively.
The properties of localization of the wavelets guarantee absolute convergence in each W j (see [47] for details). Now we choose the levels where the energy of g is concentrated, J min , J max ∈ Z so that for small ε, and truncate the component in each level so that the following approximation of g j arises: Afterwards we obtain the fractional Caputo-Fabrizio derivatives of the wavelet basis by means of Eq. (29): Let us consider for a moment that in Eq. (31) we have CF À∞ D α t , i.e., a ¼ À∞. Note that, since supp |v jk | ⊂ Ω j , v jk ∈ W jÀ1 ⋃W j ⋃W jþ1 , but, from the properties of the chosen basis, we can consider v jk ∈ W j . This fact enables us to work on each level j separately (details can be found in [44,48]).
Then, since the unknown f can be expressed as f t Finally, we replace this last expression in Eq. (31), where for simplicity we will first consider σ 1 ¼ 0 and look for the wavelet coefficients b jk that satisfy, for each j ∈ J min ; J max ½ : that, in matrix form, results in From the properties of the wavelet basis, and those of v jk , it results that M j is a diagonal dominant matrix and, consequently, the vector of coefficients b j can be computed in a stable and accurate way. The solution f can be obtained from Eq. (36). Moreover, it can be shown that f 0 ð Þ ¼ 0. To correct the effect of is the characteristic function of the interval 0; T ½ .
Finally, e f is an approximate solution to Eq. (31).
The error introduced in the approximation can be controlled and reduced: a more accurate truncated projection of the data into the wavelet subspaces can be considered, and the elements of the matrix can be computed with good precision since they can be expressed as integrals over compact subsets; finally, the matrix of the resulting linear system is a diagonal dominant matrix, and the solution can be computed accurately. In summary, the good properties of the basis and the operator guarantee that the resulting approximation scheme is efficient and numerically stable and no additional conditions need to be imposed.

Example 1
We illustrate the performance of the proposed approximation scheme by solving the IVP described by Eq. (31) for σ 0 ¼ 0:9, σ 1 ¼ 0, and α ¼ 0:5 and g a causal function defined as g t ð Þ ¼ v t ð Þ sin 2πt ð Þ cos 0:5πt ð Þ, where v is a smooth window in the interval 0; 4 ½ . Wavelet analysis indicates that the energy of the data g is concentrated in the subspaces W 0 and W 1 ; thus, we consider levels À1 ≤ j ≤ 2 for the reconstruction.
we obtain similar equations for the coefficients on each level: Once more the matrix of the resulting linear system is diagonal dominant, and the system can be solved efficiently (see [49] for details).
The following example shows the performance of the method for σ 1 ¼ 0:3.
The plots of the f (blue) and e f (green) appear in Figure 2. This scheme can be adapted to solve boundary value problems. We show the procedure finding the solution to the fractional heat equation.

Boundary value problem
We show how to adapt the scheme used for initial value problem for solving boundary value problems with fractional partial differential equations in an example.
We will consider a fractional heat problem where we have replaced the classical time derivative by the Caputo-Fabrizio fractional derivative of order α: This equation models the evolution of temperatures in a bar of length 1, constituted by a heterogeneous material which has "memory," due to the fluctuations introduced by elements at different dimension scales ( [7]).
The smooth and causal function g represents an external source. We look for smooth solutions u ∈ C 2 0; 1 ð ÞÂ 0; T ð Þby separating variables and pose where u k t ð Þ is the Fourier coefficients of u x; t ð Þ for each t ∈ 0; T ð Þ. For the temporal part of the function, after replacing in Eq. (36), we obtain an initial value problem like that described by Eq. (31) for each coefficient u k t ð Þ: f and e f of Example 2.
We show the approximate solution to Eq. (41) for α ¼ 0:5, T ¼ 3, and g x; t ð Þ ¼ v t ð Þe Àt=2 sin 5πt ð Þsin 2πx ð Þ, with v a smooth window in 0; 3 ½ .  In this case we only need to solve Eq. (43) for k ¼ 2, with B 2 ¼ v t ð Þe Àt=2 sin 5πt ð Þ: Wavelet analysis indicates that the 95% of the energy of B 2 is concentrated in subspaces W 1 , W 2 , and W 3 , and we obtained the following condition numbers for the band matrices of Eq. (39): cond ∞ M 0 À Á ¼ 1:1153, cond ∞ M 1 À Á ¼ 1:0663, cond ∞ M 2 À Á ¼ 1:0132, cond ∞ M 3 À Á ¼ 1:0098, and cond ∞ M 4 À Á ¼ 1:0076. We consider levels 0 ≤ j ≤ 4 for the reconstruction, and the mean square error in this case is u 2 À e u 2 k L 2 ¼ 3:5020 10 À4 . Figure 3 shows the approximate and the exact solution to Eq. (43). In Figure 4 we draw the approximate solution u of the heat problem described by Eq. (41), and  in Figure 5 the difference between the true solution u x; t ð Þ and its approximation is plotted. The mean square error obtained in this case is 4:0016 10 À4 : Finally, in Figure 6, we show the approximate solutions u 2 for different orders of derivation α, α ! 1, exhibiting the tendency to the solution of Eq. (43) with α ¼ 1, as expected (for α ¼ 1, Eq. (41) describes classical heat problem, for a bar made of a homogeneous material).

Conclusions
In this chapter we have presented a summary of some recent works showing the relevance and the intense research work in the area of fractional calculus and its applications. We have focused on nonlinear models describing different phenomena where fractional differentiation plays an important role. In the last section we have presented an approximation scheme that we have developed to solve linear initial and boundary value problems based on wavelet decomposition, and the performance of the method is illustrated by examples. Possible extensions and adaptation to nonlinear equations are still under study.