Open access peer-reviewed chapter

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

By María I. Troparevsky, Silvia A. Seminara and Marcela A. Fabio

Submitted: February 22nd 2019Reviewed: April 10th 2019Published: May 15th 2019

DOI: 10.5772/intechopen.86273

Downloaded: 1249


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.


  • fractional derivatives
  • fractional differential equations
  • wavelet decomposition
  • numerical approximation

1. 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 σtand strain εtin 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 Ethe modulus of elasticity. We can rewrite both Eqs. (1) and (2) as


with α=0for elastic solids and α=1for a viscous liquid. But, in practice, there exist viscoelasticmaterials 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 nN, is


Recalling that gamma function verifies 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 n1<α<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 ntimes.” Of course these definitions aren’t equivalent: clearly the domains of the operators Dtα0RL.and Dtα0C.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 α01, was defined as


for Mαis a normalizing factor verifying M0=M1=1.

Let us point out that, both in Eqs. (7) and (9), the lower limit in the integral could be changed by any value at, i.e.,




In [7] the authors prove that the operator defined in Eq. (9) verifies the following (convenient) properties:

Dtα0CFk=0, for any constant k.


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=0zkΓα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 fin a whole interval of integration. When defined with 0as lower limit of integration, as we did, function fis usually assumed to be causal (i.e., ft0for 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.

  5. For nN,0<α<1:DtαDtnft=DtnDtαft.

  6. Dtαftmust depend on the past history of f.

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

Dtα0Cft=Ftft0T,T<initial or boundary conditionsE16


Dtα0RLft=Ftft0T,T<initial or boundary conditionsE17

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 gis 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 convection-diffusion 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.


2. 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 Dtα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 xSRand tR>0denote 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 β02.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<β<1and 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 Cxtvaries in time and space, Dis a diffusion constant, and Cis 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 E0is the elastic constant for a spring and η1and η2represent 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,Htis the magnetic excitation field, Mtis the magnetization vector, θtis the temperature, θcis the Curie temperature below which the hysteresis is observed, and C0is the tensor with the constitutive properties of the magnetic material. They compare the resulting behavior when Dtα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.

3. 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.3and 0.4<α<0.5for the usual set of parameters σ=10,ρ=83,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, for0<α<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 ρ.

4. 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 fH1b, the Sobolev space of functions with (weak) first derivative in L2b, 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 kis a causal function, kt=eαt1αfor t0, and kt=0for t<0.

Consequently, since k̂ω=1αα+1α,


Using the properties of the Fourier transform, we can rewrite the last equality:


where hω=2πk̂ω.

Meanwhile, in the Caputo case, we have


where kt=1tαand k̂ω=Γ1α1α.

4.1 An initial value problem

Let us consider the initial value problem (IVP):


We look for fsatisfying Eq. (31), where gis a causal and smooth function with g0=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ψ̂jk=Ωjwhere Ωj=ω:2jπβω2j(π+β}, with 0<βπ3([45]).

The family ψjk/ψjk=2j2ψ2jtk,jkZis an orthonormal basis (BON) of L2Rassociated to a multiresolution analysis (MRA). We denote, by Wj=spanψjkjkZand VJ=j<JWj, the wavelet and scale subspaces, respectively, and decompose the space L2R=jZWj=j>nWj+Vn,nZ. We have also a scale function φV0so that φtkkZis a BON of V0([46]).

The sets Ωj1,Ωj,Ωj+1have little overlap, and Wjis 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 gas


where the first and second terms are the projections of gin VJand Wj, respectively.

The properties of localization of the wavelets guarantee absolute convergence in each Wj(see [47] for details).

Now we choose the levels where the energy of gis concentrated, Jmin,JmaxZso that


for small ε, and truncate the component in each level so that the following approximation of gjarises:


Afterwards we obtain the fractional Caputo-Fabrizio derivatives of the wavelet basis by means of Eq. (29):


(recall suppψ̂jk=Ωj).

Let us consider for a moment that in Eq. (31) we have −∞CFDtα, i.e., a=. Note that, since suppvjkΩj, vjkWj1WjWj+1, but, from the properties of the chosen basis, we can consider vjkWj. This fact enables us to work on each level jseparately (details can be found in [44, 48]).

Then, since the unknown fcan be expressed as ft=jZfjt,where fjt=kZbjkψjktand bjk=fψjk,we have


Finally, we replace this last expression in Eq. (31), where for simplicity we will first consider σ1=0and look for the wavelet coefficients bjkthat satisfy, for each jJminJmax:




that, in matrix form, results in


where bj=bjkkKj, cj=cjkkKj, MjRKj×KjandMjkl=vklψkl+σ0ψkj,ψkl.

From the properties of the wavelet basis, and those of vjk, it results that Mjis a diagonal dominant matrix and, consequently, the vector of coefficients bjcan be computed in a stable and accurate way. The solution fcan be obtained from Eq. (36). Moreover, it can be shown that f0=0. To correct the effect of having considered −∞CFDtαinstead of 0CFDtα, we set f=f.χ0T, where χ0Tis the characteristic function of the interval 0T. Finally, fis 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.

4.1.1 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.5and ga causal function defined as gt=vtsin2πtcos0.5πt,where vis a smooth window in the interval 04. Wavelet analysis indicates that the energy of the data gis concentrated in the subspaces W0and W1; thus, we consider levels 1j2for the reconstruction.

For this case, being σ1=0, σ011α, gC0, and g0=0, there exists a formula for the “exact” solution to Eq. (31) (see [22]). The approximate solution fto the IVP is plotted (in green) in Figure 1, together with the exact solution (in blue).

Figure 1.

The approximation and the exact solution for the Example 1.

If σ10, since ft=12πRf̂ωeiωt, 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.

4.1.2 Example 2

Now we consider IVP described by Eq. (31) for σ0=0.9,σ1=0.3,α=0.5, and ft=t2vt2sin2.5πt0.5cos6πt, with vas a smooth window over interval [0, 7]. Since σ10and we have no formula to test the performance of the approximation, for this example we will set f,calculate gfrom Eq. (31), and then apply the proposed method to recover f.

Choosing 1j2, it results in f=j=12fj. The plots of the f(blue) and f(green) appear in Figure 2.

Figure 2.

fandfof Example 2.

This scheme can be adapted to solve boundary value problems. We show the procedure finding the solution to the fractional heat equation.

4.2 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 grepresents an external source. We look for smooth solutions uC201×0Tby separating variables and pose


where uktis the Fourier coefficients of uxtfor each t0T.

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 ukt:


with Bkt=201gxtsinkπxdxand the Fourier coefficients of gxtfor each t0T. The uniqueness of solution is guaranteed because ux0=0, so Bk0is null.

We show the approximate solution to Eq. (41) for α=0.5,T=3,and gxt=vtet/2sin5πtsin2πx, with va smooth window in 03.

In this case we only need to solve Eq. (43) for k=2,with B2=vtet/2sin5πt.Wavelet analysis indicates that the 95% of the energy of B2is concentrated in subspaces W1,W2, and W3, and we obtained the following condition numbers for the band matrices of Eq. (39): condM0=1.1153, condM1=1.0663, condM2=1.0132, condM3=1.0098, and condM4=1.0076. We consider levels 0j4 for the reconstruction, and the mean square error in this case is u2u2L2=3.5020104.

Figure 3 shows the approximate and the exact solution to Eq. (43). In Figure 4 we draw the approximate solution uof the heat problem described by Eq. (41), and in Figure 5 the difference between the true solution uxtand its approximation is plotted. The mean square error obtained in this case is 4.0016104.

Figure 3.

The approximate and the exact solutions toEq. (41)withα=0.8,T=3,andB3=vtet2sin5πt.

Figure 4.

The approximate solutions toEq. (39)withα=0.8,T=3,gxt=vtet/2sin5πtsin2πx,andva smooth window in03.

Figure 5.

The difference between the true solutionuxttoEq. (39)and its approximation by means of the wavelet scheme.

Finally, in Figure 6, we show the approximate solutions u2for 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).

Figure 6.

The approximate solutions toEq. (41) withα1.

5. 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.



This work was partially supported by University of Buenos Aires through UBACyT (2018–2021) 20020170100350BA.

Conflict of interest

The authors declare no conflicts of interest regarding this chapter.

© 2019 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

María I. Troparevsky, Silvia A. Seminara and Marcela A. Fabio (May 15th 2019). A Review on Fractional Differential Equations and a Numerical Method to Solve Some Boundary Value Problems, Nonlinear Systems -Theoretical Aspects and Recent Applications, Walter Legnani and Terry E. Moschandreou, IntechOpen, DOI: 10.5772/intechopen.86273. Available from:

chapter statistics

1249total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Numerical Solutions to Some Families of Fractional Order Differential Equations by Laguerre Polynomials

By Adnan Khan, Kamal Shah and Danfeng Luo

Related Book

First chapter

Applications of 2D Padé Approximants in Nonlinear Shell Theory: Stability Calculation and Experimental Justification

By Igor Andrianov, Jan Awrejcewicz and Victor Olevs’kyy

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us