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
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 . Recall the relationship between stress and strain 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
with for elastic solids and for a viscous liquid. But, in practice, there exist
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  for a more complete review and [1, 3, 4, 5, 6] for rigorous theoretical expositions and calculation methods).
The classical Cauchy formula for the -fold iterated integral, with , is
Recalling that gamma function verifies , 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 , 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 times and deriving times.” Of course these definitions aren’t equivalent: clearly the domains of the operators and are different; because of the different hypothesis about , we need to impose to guarantee their existence. Besides that, with the appropriate conditions for ,
More recently, the
for is a normalizing factor verifying .
, for any constant .
Caputo-Fabrizio definition was then generalized by Atangana and Baleanu, who gave the following definition of the
replacing the exponential by , 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 in a whole interval of integration. When defined with as lower limit of integration, as we did, function is usually assumed to be causal (i.e., for ), 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 () proposed the following terms to recognize if an integral operator merits to be called a fractional derivative:
The fractional derivative must be a linear operator.
The fractional derivative of an analytic function must be analytic.
If the order of the derivative is a positive integer, the derivative must be the classical one.
If the order is null, the original function must be recovered.
must depend on the past history of .
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.
for . 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 , besides being continuous, must satisfy certain conditions that substitute the classical Lipschitz’s one.
Similar results are stated in  for a coupled system of fractional differential equations involving Riemann-Liouville derivatives.
are proved, and explicit formulae are presented when 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  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  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  to solve fractional partial differential equations arising in neuronal dynamics.
In  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  to solve diffusion equations, and implicit numerical algorithms for the case of spatial and temporal fractional derivatives appeared in . A high-speed numerical scheme for fractional differentiation and fractional integration is proposed in . In , 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.
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 and are used to indistinctly represent any of the fractional derivatives, whose type is clarified in the text.
In  the authors review the evolution of the general fractional equation:
for , where and 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 The authors show that, for , 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  the authors study the particular linear fractional Klein-Gordon equation:
considering the fractional Caputo derivative with . They achieve a numerical solution by variational iteration method and multivariate Padé approximation.
In  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 ():
Here, the calcium concentration varies in time and space, is a diffusion constant, and 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  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
where is the elastic constant for a spring and and represent the viscosities of two springpots in parallel with the spring.
In Ref.  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  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 is the magnetic excitation field, is the magnetization vector, is the temperature, is the Curie temperature below which the hysteresis is observed, and is the tensor with the constitutive properties of the magnetic material. They compare the resulting behavior when 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  the authors studied a system based on the classical Lorenz one, but described by the Atangana-Baleanu fractional derivative (in the Caputo sense) with :
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 and for the usual set of parameters .
In  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  a fractional Lorenz system is studied considering the
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.
For the Caputo-Fabrizio fractional derivative of order , changing variables in Eq. (9), we have
where is a causal function, for , and for .
Consequently, since ,
Using the properties of the Fourier transform, we can rewrite the last equality:
Meanwhile, in the Caputo case, we have
where and .
4.1 An initial value problem
Let us consider the initial value problem (IVP):
We look for satisfying Eq. (31), where is a causal and smooth function with . 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  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 (the Schwartz space) and its Fourier transform satisfy where , with ().
The family is an orthonormal basis (BON) of associated to a multiresolution analysis (MRA). We denote, by and , the wavelet and scale subspaces, respectively, and decompose the space . We have also a scale function so that is a BON of ().
The sets have little overlap, and is nearly a basis for the set of functions whose Fourier transform has support in . This property of the basis will be crucial in the procedure.
Now we decompose the data as
where the first and second terms are the projections of in and , respectively.
The properties of localization of the wavelets guarantee absolute convergence in each (see  for details).
Now we choose the levels where the energy of is concentrated, so that
for small , and truncate the component in each level so that the following approximation of 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 , i.e., . Note that, since , , but, from the properties of the chosen basis, we can consider . This fact enables us to work on each level separately (details can be found in [44, 48]).
Then, since the unknown can be expressed as where and we have
Finally, we replace this last expression in Eq. (31), where for simplicity we will first consider and look for the wavelet coefficients that satisfy, for each :
that, in matrix form, results in
where , , .
From the properties of the wavelet basis, and those of , it results that is a diagonal dominant matrix and, consequently, the vector of coefficients can be computed in a stable and accurate way. The solution can be obtained from Eq. (36). Moreover, it can be shown that . To correct the effect of having considered instead of , we set , where is the characteristic function of the interval . Finally, 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.
4.1.1 Example 1
We illustrate the performance of the proposed approximation scheme by solving the IVP described by Eq. (31) for and a causal function defined as where is a smooth window in the interval . Wavelet analysis indicates that the energy of the data is concentrated in the subspaces and ; thus, we consider levels for the reconstruction.
For this case, being , , , and , there exists a formula for the “exact” solution to Eq. (31) (see ). The approximate solution to the IVP is plotted (in green) in Figure 1, together with the exact solution (in blue).
If , since , 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  for details).
The following example shows the performance of the method for .
4.1.2 Example 2
Now we consider IVP described by Eq. (31) for , and , with as a smooth window over interval [0, 7]. Since and we have no formula to test the performance of the approximation, for this example we will set calculate from Eq. (31), and then apply the proposed method to recover .
Choosing , it results in . The plots of the (blue) and (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.
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 ().
The smooth and causal function represents an external source. We look for smooth solutions by separating variables and pose
where is the Fourier coefficients of for each .
with and the Fourier coefficients of for each . The uniqueness of solution is guaranteed because , so is null.
We show the approximate solution to Eq. (41) for and , with a smooth window in .
In this case we only need to solve Eq. (43) for with Wavelet analysis indicates that the 95% of the energy of is concentrated in subspaces , and , and we obtained the following condition numbers for the band matrices of Eq. (39): , , , , and . We consider levels 4 for the reconstruction, and the mean square error in this case is .
Figure 3 shows the approximate and the exact solution to Eq. (43). In Figure 4 we draw the approximate solution of the heat problem described by Eq. (41), and in Figure 5 the difference between the true solution and its approximation is plotted. The mean square error obtained in this case is
Finally, in Figure 6, we show the approximate solutions for different orders of derivation , , exhibiting the tendency to the solution of Eq. (43) with , as expected (for , Eq. (41) describes classical heat problem, for a bar made of a homogeneous material).
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.