© Hindawi Publishing Corp. RECENT APPLICATIONS OF FRACTIONAL CALCULUS TO SCIENCE AND ENGINEERING

This paper deals with recent applications of fractional calculus 
to dynamical systems in control theory, electrical circuits with 
fractance, generalized voltage divider, viscoelasticity, 
fractional-order multipoles in electromagnetism, electrochemistry, tracer in fluid flows, and model of neurons 
in biology. Special attention is given to numerical computation of 
fractional derivatives and integrals.


Introduction.
During the second half of the twentieth century, considerable amount of research in fractional calculus was published in engineering literature.Indeed, recent advances of fractional calculus are dominated by modern examples of applications in differential and integral equations, physics, signal processing, fluid mechanics, viscoelasticity, mathematical biology, and electrochemistry.There is no doubt that fractional calculus has become an exciting new mathematical method of solution of diverse problems in mathematics, science, and engineering.In order to stimulate more interest in the subject and to show its utility, this paper is devoted to new and recent applications of fractional calculus in science and engineering.
Historically, fractional calculus originated from the Riemann-Liouville definition of fractional integral of order α in the form (1.1) Often a D −α x = a J α x is called the Riemann-Liouville integral operator.When a = 0, (1.1) is the original Riemann definition of fractional integral, and if a = −∞, (1.1) represents the Liouville definition (see Debnath [12,13]).Integrals of this type were found to arise in the theory of linear ordinary differential equations where they are known as Euler transforms of the first kind and in Riesz's classic memoir [38] on the Cauchy problem for hyperbolic equations.
If a = 0 and x > 0, then the Laplace transform solution of the initial value problem D n y(x) = f (x), x > 0, y (k) where e −sx y(x)dx.(1.4)This inverse Laplace transform gives the solution of the initial value problem ( This is the Riemann-Liouville integral formula for an integer n.Replacing n by real α gives the Riemann-Liouville fractional integral (1.1) with a = 0. We consider a definite integral in the form and hence Thus, for a positive integer n, it follows that (1.9) Replacing n by α where Re α > 0 in (1.9) leads to the definition of the Riemann-Liouville fractional integral (1.10) Thus, this fractional-order integral formula is a natural extension of an iterated integral.The fractional integral formula (1.10) can also be obtained from the Euler integral formula x 0 (x − t) r t s dt = Γ (r + 1)Γ (s + 1) Γ (r + s + 2) x r +s+1 , r,s >−1. (1.11) Replacing r by n − 1 and s by n gives (1.12) Consequently, (1.9) follows from (1.12) when f (t) = t n and a = 0.In general, (1.12) gives (1.9) replacing t n by f (t).Hence, when n is replaced by α, we derive (1.10).
It may be interesting to point out that Euler's integral expression for the 2 F 1 (a,b,c; x) hypergeometric series can now be expressed as a fractional integral of order (c − b) where f (t) = t b−1 (1 − t) −a .One of the very useful results is the formula for the Laplace transform of the derivative (see Debnath [11]) of an integer-order n of a function f (t): (1.14) where f (n−k) (0) = c k represents the physically realistic given initial conditions.
Like the Laplace transform of integer-order derivative, it is easy to show that the Laplace transform of fractional-order derivative is given by where (1.17) represents the initial conditions which do not have obvious physical interpretation.Consequently, formula (1.16) has limited applicability for finding solutions of initial value problems in differential equations.We now replace α by an integer-order integral J n , and D n f (t) ≡ f (n) (t) is used to denote integer-order derivative of a function f (t).It turns out that This simply means that D n is the left inverse (not the right inverse) of J n .It also follows from (1.21) with α = n that Similarly, D α can also be defined as the left inverse of J α .With n−1 < α ≤ n, we define the fractional derivative of order α > 0 by where n is an integer and the identity operator I is defined by 16), Caputo and Mainardi [9] adopted as an alternative new definition of fractional derivative to solve initial value problems in viscoelasticity.This new definition was originally introduced by Caputo [7,8] in the form where n − 1 < α < n and n is an integer.According to this definition, That is, Caputo's fractional derivative of a constant is zero.It follows from (1.20) and (1.21) that unless f (t) and its first (n − 1) derivatives vanish at t = 0.
Furthermore, it follows from (1.21) and (1.19) that .24)This implies that This shows that Caputo's fractional derivative incorporates the initial values The Laplace transform of Caputo's fractional derivative (1.25) gives an interesting formula This is a natural generalization of the corresponding well-known formula for the Laplace transform of f n (t) when α = n, and can be used to solve initial value problems in differential equations with physically realistic initial conditions.The reader is referred to [12,13] for examples of applications of fractional calculus to ordinary and partial differential equations in applied mathematics and fluid mechanics.

The Weyl fractional integral and the Mellin transform.
In order to extend the idea of fractional integrals of ordinary functions to fractional integrals of generalized functions, it is convenient to adopt the convolution theory of distributions.In view of formula (1.1), 0 D −α t f (t) can be treated as the convolution of f (t) (assumed to vanish for t < 0) with the function φ α (t) defined by This function φ α can be extended to all complex values of α as a pseudofunction and becomes a distribution whose support is [0, ∞) except for the case α = 0, −1, −2,.... Similarly, the Weyl [46] fractional integral t W −α ∞ f (t) can also be regarded as the convolution of φ α (−t) with f (t) so that (2.2) Under suitable conditions by Fubini's theorem, a formal computation of the inner product gives This shows that D −α and W −α are adjoint operators (see [14]) in some sense.
If we can obtain a test function space Φ which is mapped by one of the operators D −α or W −α continuously into another test function space Ψ , then the adjoint operator would map generalized function on Ψ (the dual of Ψ ) into generalized functions on Φ, and this adjoint operator will define the other of D −α or W −α for generalized functions.
We next calculate the Mellin transforms of the Riemann-Liouville fractional integrals and derivatives.It is convenient to write where and H(t) is the Heaviside unit step function.
Using the definition of the Mellin transform of f (t) (see Debnath [11]), we obtain Using the properties of the Mellin transform, it turns out that where the Mellin transform f (p) of f (t) is defined by Similarly, using formula (1.20), we can obtain the Mellin transform of the Riemann-Liouville fractional derivative as (2.9) We next find the Mellin transform of the Weyl fractional integral where (2.11) Thus, the Mellin transform of (2.10) is Or, equivalently, Similarly, the Mellin transform of the Weyl fractional derivative is given by (2.14) Or, equivalently, We conclude that formulas (2.7) and (2.14) can be used to define fractional integrals of generalized functions as the Mellin transformation has been extended to generalized functions (see Zemanian [47]).

3.
Fractional-order dynamical systems in control theory.New and effective methods for the time-domain analysis of fractional-order dynamical systems are required for solving problems of control theory.As a new generalization of the classical P ID-controller, the idea of P I λ D µ -controller, involving fractional-order integrator and fractional-order differentiator, has been found to be a more efficient control of fractional-order dynamical systems.In his series of papers and books (see references of Podlubny's book [37]), Oustaloup [34,35,36] successfully used the fractional-order controller to develop the so-called CRONE-controller (Commande Robuste d'Ordre Non Entrier controller) which is an interesting example of application of fractional derivatives in control theory.He demonstrated the advantage of the CRONEcontroller compared to the classical P ID-controller and also showed that the P I λ D µ -controller has a better performance record when used for the control of fractional-order systems than the classical P ID-controller.
In the time domain, a dynamical system is described by the fractional-order differential equation (FDE) where The unit-impulse response y i (t) of the system is defined by the inverse Laplace transform of G n (s) so that and the unit-step response function is given by the integral of g n (t) so that We illustrate the above system response by simple examples.
Example 3.1.We consider a simple fractional-order transfer function The corresponding system responses are where Ᏹ 0 (t, x; α, β) is defined in terms of Mittag-Leffler's function (see Erdélyi [16]) E α,β (z) as and the Laplace transform is The function Ᏹ m satisfies the property (3.9) Example 3.2.We consider a more general fractional-order-controlled system with the transfer function The fractional differential equation in the time domain corresponding to the transfer function (3.10) is with the initial conditions The unit-impulse response y i (t) to the system is given by and the unit-impulse response to the system is Example 3.3 (the P I λ D µ -controller).The transfer function G c (s) of this controller is defined as the ratio of the controller output U(s) and error E(s) so that we can write The corresponding output u(t) = ᏸ −1 {U(s)} for the P I λ D µ -controller in the time domain is given by (3.16) When λ = µ = 1, the above results reduce to those of the classical P IDcontroller.When λ = 1 and µ = 0, we obtain the corresponding results for P I-controller, while λ = 0 and µ = 1 leads to the P D-controller.
All P ID-controllers are special cases of the fractional P I λ D µ -controller described by its transfer function (3.15).Numerous applications have demonstrated that P I λ D µ -controllers perform sufficiently better for the control of fractional-order dynamical systems than the classical P ID-controllers.

Electrical circuits with fractance.
Classical electrical circuits consist of resistors and capacitors and are described by integer-order models.However, circuits may have the so-called fractance which represents an electrical element with fractional-order impedance as suggested by Le Mehaute and Crepy [23].
We consider two kinds of fractances: (i) tree fractance and (ii) chain fractance.
Following Nakagawa and Sorimachi [30], we consider a tree fractance consisting of a finite self-similar circuit with resistors of resistance R and capacitors of capacitance C. The impedance of the fractance is given by The associated fractional-order transfer function of this tree fractance is For a chain fractance consisting of N pairs of resistor-capacitor connected in a chain, Oldham and Spanier [33] have shown that the transfer function is approximately given by It can be shown that this chain fractance behaves as a fractional-order integrator of order 1/2 in the time domain 6RC ≤ t < (1/6)N 2 RC.
Due to microscopic electrochemical processes at the electrode-electrolyte interface, electric batteries produce a limited amount of current.At metalelectrolyte interfaces the impedance function Z(ω) does not show the desired capacitive features for frequencies ω.Indeed, as ω → 0, Or, equivalently, in the Laplace transform space, the impedance function is This illustrates the fact that the electrode-electrolyte interface is an example of a fractional-order process.The value of η is closely associated with the smoothness of the interface as the surface is infinitely smooth as η → 1.
Kaplan et al. [21] proposed a physical model by the self-affine Cantor block with N-stage electrical circuit of fractance type.Under suitable assumptions, they found the importance of the fractance circuit in the form where η = 2 − log(N 2 )/ log a, K and a are constants, and N 2 > a implies 0 < η < 1.This shows that the model of Kaplan et al. is also an example of the fractional-order electric circuit.In a resistive-capacitive transmission line model, the inter conductor potential φ(x, t) or inter conductor current i(x, t) satisfies the classical diffusion equation where the diffusion constant κ is replaced by (RC) −1 , R and C denote the resistance and capacitance per unit length of the transmission line, and u(x, t) = φ(x, t) or i(x, t).
Using the initial and boundary conditions This confirms that the current field in the transmission line of infinite length is expressed in terms of the fractional derivative of order 1/2 of the potential φ(0,t).This is another example of the involvement of fractional-order derivative in the electric transmission line.[45] observed that both the tree fractance and chain fractance consist not only of resistors and capacitors properties, but also they exhibit electrical properties with noninteger-order impedance.He generalized the classical voltage divider in which the fractionalorder impedances F 1 and F 2 represent impedances not only on Westerlund's capacitors, classical resistors, and induction coils, but also impedances of tree fractance and chain fractance.The transfer function of Westerlund's voltage divider circuit is given by

Generalized voltage divider. Westerlund
where −2 < α < 2 and k is a constant that depends on the elements of the voltage divider.The negative values of α correspond to a highpass filter, while the positive values of α correspond to a lowpass filter.Westerlund considered some special cases of the transfer function (5.1) for voltage dividers that consist of different combinations of resistors (R), capacitors (C), and induction coils (L).
If U in (s) is the Laplace transform of the unit-step input signal u in (t), then the Laplace transform of the output signal U out (s) is given by The inverse Laplace transform of (5.2) is obtained from (3.8) to obtain the output signal (5.3) Although the exact solution for the output signal is obtained, this cannot describe physical properties of the signal.Some physically interesting properties of the output signal can be described for various values of α by evaluating the inverse Laplace transform in the complex s-plane.For 1 < |α| < 2, the output signal exhibits oscillations.

Fractional calculus in viscoelasticity.
Almost all deformed materials exhibit both elastic and viscous properties through simultaneous storage and dissipation of mechanical energy.So any viscoelastic material may be treated as a linear system with the stress (or strain) as excitation function (input) and the strain (or stress) as the response function (output).Several stress-strain relationships are all known in viscoelasticity, and they represent only mathematical models for an ideal solid or liquid.Neither of the classical laws (Hooke's law for elastic solids and Newton's law for viscous liquids) can adequately describe the real situation in the real world.On the other hand, Kelvin's model for viscoelasticity describing the connection of the Hooke elastic element and the Voigt viscoelastic element is given by where σ is the stress, ε is the strain, α, β, η, E 1 , and E 2 are constants related by Zener [48] formulated a mathematical model of viscoelasticity connecting the Hooke elastic element and the Maxwell viscoelastic element in the form Both Kelvin's and Zener's mathematical models provide a reasonable qualitative description; however, they are not satisfactory from a quantitative viewpoint (see Caputo and Mainardi [9]).Several authors including Scott Blair [42], Gerasimov [19], Slonimsky [43], Stiassnie [44], and Mainardi [28] suggested that integer-order models for viscoelastic materials seem to be inadequate from both qualitative and quantitative points of view.At the same time, they proposed fractional-order laws of deformation for modelling the viscoelastic behavior of real materials.Scott Blair [42] formulated a fractional-order model in the form where E and α are constants which depend on the nature of the material.On the other hand, Gerasimov [19] suggested a similar generalization of the fundamental law of deformation as where κ is referred to as the generalized viscosity of the material.Caputo and Mainardi [9] extended Zener's integer-order model of standard linear solid to a fractional-order model in the form It has been shown experimentally that law (6.6) is very useful for modelling of most viscoelastic materials (see Bagley and Torvik [3,5] and Rogers [40]).In addition to experimental findings, Bagley and Torvik proved that the fourparameter model (6.6) seems to be satisfactory for most real materials.
In a series of papers, Bagley and Torvik [3,5,6] and Koeller [22] proposed the fractional derivative approach to viscoelasticity in order to describe the properties of numerous viscoelastic materials.It has been demonstrated that this new approach leads to well-posed problems of motion of structures containing elastic and viscoelastic components even when incorporated into finiteelement formulations.They suggested the general form of the empirical model as where D α is defined by (1.20) with n = 1.Result (6.7) shows that the timedependent stress σ (t) is the superposition of an elastic term E 0 ε(t) and a viscoelastic term containing a fractional derivative of order α.Application of the Fourier transform, with respect to t in (−∞, ∞), to (6.7) gives where ω is the Fourier transform variable and E 0 , E 1 , and α are three parameters of the model.It is found that results of this model describe accurately the properties of viscoelastic materials as predicted by Bagley and Torvik [6].
According to the molecular theory for dilute polymer solutions due to Rouse [41], the shear modulus G(ω) is a complex function of real frequency ω so that where n is the number of molecules per unit volume of the polymer solution, k is the Boltzmann constant, T is the absolute temperature, τ p is the characteristic relaxation time of the solution, N is the number of submolecules in a polymer chain, and µ s is the steady-flow viscosity of the solvent in the solution.It is shown by Rouse [41] that where µ 0 is the steady-flow viscosity of the solution.
The spectrum of the stress time history can be calculated from where G(ω) is given in (6.11).Finally, in the time domain, the stress in a dilute polymer solution obtained by Rouse [41] consists of two terms as where the first term is the contribution from the solvent and the second fractional derivative term is the contribution from the solute chain molecules.Thus, the Rouse theory provides with a nonempirical basis for the presence of fractional derivatives along with the first derivative of classical viscoelasticity in the relation between stress and strain for some polymers.Subsequently, Ferry et al. [17] modified the Rouse theory in concentrated polymer solutions and polymer solids with no crosslinking, and obtained the following result: where M is the molecular weight, ρ is the density, and R is the universal gas constant.Thus, the fractional calculus approach to viscoelasticity for the study of viscoelastic material properties is justified, at least for polymer solutions and for polymer solids without crosslinking.
Bagley and Torvik [3] also used fractional calculus to construct stress-strain relationships for viscoelastic materials.They proposed five-parameter model in the form where b, E 0 , E 1 , α, and β are parameters.
Using the Fourier transform of (6.15) with respect to t in the five-parameter viscoelastic model gives .16)This shows that the frequency-dependent modulus Ẽ(ω) is a function of fractional powers of frequency.
At very low frequencies (the rubbery region), the modulus tends to the rubbery modulus E 0 .As the frequency increases into the transition region, the second term in the numerator of (6.16) produces the increase in the real and imaginary parts of the modulus associated with the transition region.Finally, as the frequency is advanced into the glassy region, the modulus tends to the constant value (E 1 /b) provided that α = β.
Based on this fractional calculus model, Bagley and Torvik [3] solved equations of motion for viscoelastically damped structures by the finite-element analysis.It is shown that very few empirical parameters are required to model viscoelastic materials and to compute the dynamic response of the structure under general loading conditions.Later on, Bagley and Torvik [6] also examined fractional calculus model of the viscoelastic phenomenon based on thermodynamic principles.In particular, they have shown that thermodynamic constraints on parameters of the model lead to the nonnegative rate of energy dissipation as well as the nonnegative internal work.Furthermore, these constraints ensure the model to predict realistic sinusoidal response as well as realistic relaxation and creep responses.Bagley and Torvik's analysis reveals that the fractional calculus models of viscoelastic behavior are much more satisfactory than the previously adopted classical models of viscoelasticity.
Subsequently, Bagley [2] showed that the frequency-dependent complex modulus for the four-parameter fractional model (α = β) is obtained from (6.16) in the form For low frequencies, the term [i + b(iω) α ] −1 in (6.17) can be replaced by 1 − b(iω) α for bω α  1. Retaining only those terms of order (iω) α in the resulting expression of (6.17) leads to Ẽ(ω) = E 0 + E 1 − bE 0 (iω) α , bω α  1. (6.18) On the other hand, the frequency-dependent complex modulus for the modified power law (see Bagley [2, equation ( 21 where E e is the rubbery modulus, E g is the glassy modulus, and τ 0 and n are parameters selected to fit the data.For low frequencies (ωτ 0 1), the approximate value of (6.19) is obtained by neglecting the integral in (6.19) so that At high frequencies, bω α 1 and (ωτ 0 ) n 1, the real parts of (6.18) and (6.20) are equal.When bω α  1, (6.18) approaches the value (E 1 /b).It is also shown by Bagley [2] that, for high frequencies, (ωτ 0 ) n  1, the incomplete Gamma function yields an asymptotic expansion for large (iωτ 0 ) so that the first term of the resulting expression (6.19) tends to the constant value E g .In summary, Bagley's analysis [2] suggests that the unmodified power law is a special case of the general fractional calculus model.The modified power law (6.19) is closely associated with the fractional calculus model of viscoelasticity.However, the two models have mathematically similar relaxation spectra.This similarity is found to generate an asymptotic equivalence of the models at long relaxation times and correspondingly low frequencies of motion in the lower transition and rubbery regions.On the other hand, the different behavior of the models at high frequencies in the transition and glassy regions is observed.
For particular numerical values of model's parameters, the difference is found to be small.

Fractional-order multipoles in electromagnetism.
It is well known that the axial multipole expansion of the electrostatic potential of electric charge distribution in three dimensions is where q is the so-called electric monopole moment, ε is constant permittivity of the homogeneous isotropic medium, r = (x 2 + y 2 + z 2 ) 1/2 , P n (cos θ) is the Legendre function of integer-order n.
In particular, the electrostatic potential functions for monopole (2 0 ), dipole (2 1 ), and quadrupole (2 2 ) are, respectively, given by Engheta [15] generalized the idea of the integer-order multipoles related to powers of 2 to the fractional-order multipoles that are called 2 α -poles.He obtained the potential function for 2 α -poles (0 < α < 1) along the z-axis, in terms of the Riemann-Liouville fractional derivatives in the form where l is a constant with dimension of length so that the usual dimension of the resulting volume charge density is Coulomb/m 3 .Evaluating the fractional derivative (7.3) yields the following result for the electrostatic potential: where P α (x) is the Legendre function of the first kind and of fractional degree α.

Fractional calculus in electrochemistry and tracer fluid flows.
Although the idea of a half-order fractional integral of the current field was known in electrochemistry, Oldham [31] has initiated a mathematical study of some semi-integral electroanalysis with some experimental support.Simultaneously, he and his associates (Oldham and Spanier [33]) have given considerable attention to their new approach to the solution of electrochemical problems that deal with diffusion phenomena.Subsequently, Goto and Ishii [20] developed the idea of semidifferential electroanalysis with the fractional-order diffusion equation that may occur in other fields including diffusion, heat conduction, and mass transfer.
Oldham and Spanier [32] also suggested the replacement of the classical integer-order Fick's law describing the diffusion of electroactive species toward the electrodes by a fractional-order integral law in the form where C 0 is the uniform concentration of electoactive species, κ is the diffusion coefficient, and K and R are constants.Incorporating a constant source S and a first-order removal term kC(x, t), Oldham and Spanier [33] considered the diffusion equation for the concentration C(x, t) in the form where the uniform steady state is described by Application of the Laplace transform to (8.2) and ( 8.3) gives The inverse Laplace transform gives the surface concentration C s (0,t) as This diffusion problem can be applied to modelling diffusion of atmospheric pollutants.We assume that C(z, t) represents the concentration of pollutant at height z at time t so that C(z, 0) = 0.With S ≡ 0 and the surface flux the ground-level concentration of the pollutant can be obtained from (8.5) in the form Or, equivalently, This enables us to calculate the pollutant generation rate J(0,t) from the ground-level pollution concentration C s (0,t).
In particular, if J(0,t) = J 0 = constant, then the ground-level pollution concentration C s (0,t) is This shows that the pollution concentration reaches a maximum constant value of J 0 / √ κk, and the time to reach one half of this maximum level is (0.23/k).
Dispersion is one of the most common observable phenomena in nature.When dyes or drugs are injected into blood vessels of animals, or any material (tracer) is added to fluid flows, the question is how the material is dispersed or spread in the fluid medium.Dispersion (or spreading) of tracers depends strongly on the scale of observation.In general, there are three different mechanisms of dispersion: molecular diffusion, variations in the permeability field (macrodispersion), and variations of the fluid velocity in a porous medium (microdispersion).These mechanisms take place at different scales.At large scale, dispersion is essentially controlled by permeability heterogenetics.We consider here how long-range correlations in the permeability field lead to transport equations involving fractional-order derivatives.
In a simple model of tracer dispersing into a straight cylindrical tube, usually the tracer is transported by the flowing fluid (pure convection or advection), but no molecular diffusion is involved in this dispersion process.In the case of steady laminar flow in a circular tube, it turns out that at a large-distance downstream from the injection point the concentration of material is approximately uniform over the tube cross section and the longitudinal spreading is governed by the diffusion equation.The relative importance of convection to diffusion transport is characterized by the Péclet number P e = UL/κ, where U is the fluid velocity, L is the characteristic length scale, and κ is molecular diffusion coefficient.
For a homogeneous fluid medium, the mean concentration of a tracer C(x, t) is governed by the convection-diffusion equation where u represents the mean fluid velocity.When dispersion is dominated by convection, that is, for high Péclet number, κ is proportional to u so that κ = γu, where γ is called the dispersivity constant which is the order of the pore size of the medium.Evidently, the transport equation is a diffusion equation with a constant fluid velocity.A spike (the Dirac delta function) of material injected at the entrance disperses and tends to the Gaussian form.The standard deviation of concentration which characterizes the dispersion is proportional to the square root of time so that σ x = √ 2κt.In terms of distance traveled l = ut, σ x = 2lγ.
We follow the theoretical model of convective fluid flow described in detail by Lenormand [25,26] to discuss the dispersion phenomenon due to the variation in stream tube cross section related to the permeability field K so that the local fluid velocity is directly proportional to K. In this fluid flow model, we assume that dispersion occurs due to the difference in filling time τ for each of the stream tube segments, and the probability distribution function of τ is related to the mean velocity u in the x-direction and statistical properties of the permeability field.In order to simplify notations, it is convenient to use the moments of the inverse permeability R = 1/K which is referred to as resistivity.
The correlation function of the function R between two segments at a distance h is defined by ρ(h) = R(x) • R(x + h) , where • denotes the mean value.The covariance is given by where distance h is considered as a continuous variable when h a or a discrete variable equal to ak with k denoting the number of segments between the two points of observation.
Since the displacement is governed by the time required to fill the different porous segments, it follows from Darcy's law that the variance is given by Introducing the flux through the porous sample by F(x,t), some material of mass m 0 is injected through the cross section A of the sample at time t = 0 so that the flux of injection is The spatial moment of order r is defined by weighted traveled distance with the resident mass dm = AC(x, t)dx normalized by the total mass m 0 as The arrival time moments of order r are defined by the integral over time t r weighted by the arrival mass dm = F(x,t)Adt, normalized by the total mass m 0 as The first moment of the arrival time t for a given stream tube is the sum of all the times required to fill the n segments (x = na).For any correlation, the average of a sum is the sum of the averages so that The variance of arrival times can be calculated by the formula where the dimensional covariance is To derive the fractional-order transport equation for the probability distribution F(x,t) of the arrival times, we transform the time variable t into t = t − t .The Fourier transform of F(x,t) with respect to time t is where ω is the Fourier transform variable.
Neglecting the effects of moments of arrival times higher than two, we use a Maclaurin series expansion of the Fourier transform about t = 0 so that and by invoking the time moments (the first centered moment is zero), Application of the Laplace transform to (8.22) with respect to x gives In this dispersion model, the spreading of the front is assumed to be small compared to displacement.This means that the second term in (8.23) is small compared to the first term.Consequently, we find or, equivalently, The use of the inverse Fourier transform to (8.25) yields This can be transformed into the equation Noting that the Dirac delta function represents the initial condition, application of the inverse Laplace transform gives the transport equation for F(x,t ) as We next use a permeability field with power law correlation with exponent β: where a is the length of a segment.We next compute the variance of arrival time for β ≠ 1 so that In the case of 0 < β < 1, the leading term is the power of the highest order (2 − β) and hence τ .(8.31)This is the same equation derived by Gelhar and Axness [18] based on local averaging and random distribution.However, the equation is different from the standard diffusion equation as the dispersive term is made up of the secondorder time derivatives instead of space derivatives.They used the statistical methods to explain the difference.
For long-range correlation (β < 1), the transport equation is nonlocal and has the form and in terms of fractional integral, or, equivalently, This is the fractional-order differential equation for the probability distribution function F(x,t ).
For short-range correlation (β > 1), the highest order in x is 1 and and hence, the equation for F(x,t) is where For β = 1, integration leads to a logarithm and hence this special case can be neglected.
For usual dispersion, the width of the dispersion front increases as the square root of time.A different scaling law exists for a displacement in a layered medium where front width is proportional to time t.A more general spreading law in t α has been derived for heterogeneous media.Maloy et al. [29] discovered that the dispersion front has a fractal structure at the scale of a large network.The fractal approach introduces a scale-dependent traveled distance which decreases the value of α either for diffusion or convection.Lenormand [24,27] suggested that modelling of fluid flow in a porous media presents three main properties including (i) very large range of length scales, (ii) strong heterogeneity, and (iii) in some case, instability of fluid flows.All these properties can be investigated by using a fractal approach and statistical physics.On the other hand, another new approach known as multifractal has also been suggested for describing heterogeneity and correlation of the permeability field.9. Fractional-order model of neurons in biology.Robinson [39] described the neurodynamics of the vestibulo-ocular reflex (VOR) model.The main function of VOR is to keep the retinal image stable by producing eye rotations which counterbalance head rotations.At lower frequencies (less than 0.3 Hz), the dynamics of canal afferents A(s) and vestibular and prepositus nuclei neurons V (s) reflect those of the canal receptors, and frequency response of neural discharge rate relative to angular velocity Ω(s) can be described approximately as a first-order highpass filter: where s(= iω = 2πiν) is the Laplace transform variable with ω in radians/s and ν in Hz, τ v is the vestibular time constant.Due to the action of the velocity storage integrator, the value of τ v can be larger for premotor neurons than for canal afferents.At high frequencies (greater than 0.3 Hz), motoneurons (M) dynamics could offset the mechanical lag of the eye, and the frequency response of neural discharge rate relative to eye angular position (E) can be approximately equal to a first-order leading function as where τ e is the eye time constant.Robinson's model is based on direct and integrated parallel pathways to the motoneurons.The first-order transfer functions given by (9.1) and (9.2) approximate time and frequency domain data from canal afferents, vestibular and prepositus nuclei neurons, and motoneurons.However, the fractional-order dynamics of vestibulo-oculomotor neurons suggests that fractional-order rather than integer-order forms of signal processing occur in the vestibulo-oculomotor system.
Anastasio [1] recognized some difficulties in classical integer-order models to describe the behavior of premotor neurons in the vestibulo-ocular reflex system.In order to overcome these difficulties, he proposed a fractional-order model in terms of the Laplace transform R(s) of the premotor neuron discharge rate r (t) and the Laplace transform Ω(s) of the angular velocity of the head in the form where τ 1 and τ 2 represent time constants of the neuron model, α d and α i are, respectively, fractional-order derivative and integral of the model.Assuming that α i > α d and γ = α i − α d > 0, we write so that the inverse Laplace transform of G(s) can be determined from (3.8) in terms of the Mittag-Leffler function in the form Since R(s) = G(s)Ω(s), the neuron discharge rate r (t) is given by the Laplace convolution Anastasio [1] suggested that the fractional-order derivatives can be used to describe the premotor neurons and motoneurons because the muscle and joint tissues throughout the musculoskeletal system behave as viscoelastic polymers and indicate s −k dynamics.This is likely to be compensated by s k dynamics of the associated motor and premotor neurons.In the Laplace transform domain, the fractional integrator (s −k ) is the inverse of the fractional differentiator (s k ).The transfer function assumes the form where H(s) represents the angular position of head.With regard to the possible physiological basis of fractional differentiation (s k ), it is apparent that the profuse branching of afferent terminals within the canal sensory neuroepithelium may constitute a sum of highpass filters and lead to s k dynamics.Similarly, the dendritic branching of vestibular and prepositus nuclei neurons and motoneurons may contribute to s k dynamics at these levels.On the other hand, fractional neural integration s −k could result from the joint activity of lowpass filters distributed over premotor neurons at the synaptic, cellular, and network levels.All these suggest that fractional differentiation and integration can effectively be used to describe various aspects of vestibulo-oculomotor dynamics.
In his work on the electric conductance of membranes of cells in biological systems, Cole [10] suggested the following form for the membrane reactance function X(ω) of frequency ω: where X 0 and α are constants.Thus, X(ω) corresponds to the Laplace transform function Y (s) so that where Y 0 is a constant.Cole also reported various experimental values for α that are found by other authors for different kinds of cells.This is another example of fractional-order models in biological systems.

Numerical computation of fractional derivatives and integrals.
It is often necessary to use numerical computation of fractional derivatives and integrals.For numerical treatment of fractional differential and integral equations, it is essential to have good approximation of the fractional differential operator D α and the fractional integral operator J α .
It is convenient to recall the classical formulas for the nth derivative of a function f (t) that can be approximated by the backward and central difference operators as follows: where the backward difference operator 0 ∆ n h and the central difference operator 0 δ n h are given by 0 Similarly, the Grünwald-Letnikov fractional derivative can be approximated by the backward difference operator as follows: where We next give some examples of computational solutions of fractional differential equations.
We consider the fractional-order initial value problem 0 D α t y(t) + ω 2 y(t) = f (t), t > 0, y(0) = y 0 , (10.5) where y 0 is constant.The first-order numerical approximation of this initial value problem is The results of computation from (10.4) and (10.5) are in excellent agreement with the analytical solution.The computational solution of this initial value problem has been given by Podlubny [37, Figure 8.4] and Bagley and Torvik's [4].

Call for Papers
As a multidisciplinary field, financial engineering is becoming increasingly important in today's economic and financial world, especially in areas such as portfolio management, asset valuation and prediction, fraud detection, and credit risk management.For example, in a credit risk context, the recently approved Basel II guidelines advise financial institutions to build comprehensible credit risk models in order to optimize their capital allocation policy.Computational methods are being intensively studied and applied to improve the quality of the financial decisions that need to be made.Until now, computational methods and models are central to the analysis of economic and financial decisions.
However, more and more researchers have found that the financial environment is not ruled by mathematical distributions or statistical models.In such situations, some attempts have also been made to develop financial engineering models using intelligent computing approaches.For example, an artificial neural network (ANN) is a nonparametric estimation technique which does not make any distributional assumptions regarding the underlying asset.Instead, ANN approach develops a model using sets of unknown parameters and lets the optimization routine seek the best fitting parameters to obtain the desired results.The main aim of this special issue is not to merely illustrate the superior performance of a new intelligent computational method, but also to demonstrate how it can be used effectively in a financial engineering environment to improve and facilitate financial decision making.In this sense, the submissions should especially address how the results of estimated computational models (e.g., ANN, support vector machines, evolutionary algorithm, and fuzzy models) can be used to develop intelligent, easy-to-use, and/or comprehensible computational systems (e.g., decision support systems, agent-based system, and web-based systems) This special issue will include (but not be limited to) the following topics: • Computational methods: artificial intelligence, neural networks, evolutionary algorithms, fuzzy inference, hybrid learning, ensemble learning, cooperative learning, multiagent learning ) where σ = τ σ is called the relaxation time and τ ε = b/m is referred to as the retardation time.This model includes the classical Stokes law when a = b = 0, the fractional-order Newton or Scott Blair model when a = m = 0 and b = E, the fractional Voigt model when a = 0, and the fractional Maxwell law when m = 0.

•
Application fields: asset valuation and prediction, asset allocation and portfolio selection, bankruptcy prediction, fraud detection, credit risk management • Implementation aspects: decision support systems, expert systems, information systems, intelligent agents, web service, monitoring, deployment, implementation