FRACALMO PRE-PRINT www.fracalmo.org
Probability distributions generated by
fractional diffusion equations1
Francesco MAINARDI(1), Paolo PARADISI(2) and Rudolf GORENFLO(3)
(1) Department of Physics, University of Bologna, and INFN,
Via Irnerio 46, I-40126 Bologna, Italy.
francesco.mainardi@unibo.it francesco.mainardi@bo.infn.it
(2) ISAC: Istituto per le Scienze dell’Atmosfera e del Clima del CNR,
Strada Provinciale Lecce-Monteroni Km 1.200, I-73100 Lecce, Italy.
p.paradisi@isac.cnr.ir
(3) Department of Mathematics and Computer Science,
Freie Universität Berlin, Arnimallee 3, D-14195 Berlin, Germany.
gorenflo@mi.fu-berlin.de
Contents
Abstract . . . . . . . . . . . . . . . . . . . . . . . . p. 2
1. Introduction . . . . . . . . . . . . . . . . . . . . . . . p. 2
2. The Standard Diffusion Equation . . . . . . . . . . . . . . p. 4
3. The Time-Fractional Diffusion Equation . . . . . . . . . . . p. 8
4. The Cauchy Problem for the Time-Fractional Diffusion Equation p.10
5. The Signalling Problem for the Time-Fractional Diffusion Equation p.13
6. The Cauchy Problem for the Symmetric Space-Fractional
Diffusion Equation . . . . . . . . . . . . . . . . . . . . p.15
7. Conclusions . . . . . . . . . . . . . . . . . . . . . . . p.21
A. The Riemann-Liouville Fractional Calculus . . . . . . . . . p.22
B. The Stable Probability Distributions . . . . . . . . . . . . p.31
References . . . . . . . . . . . . . . . . . . . . . . . p.41
1This paper is based on an invited talk given by Francesco Mainardi at the International
Workshop on Econophysics held at Bolyai College, Eötvös University, Budapest, on
July 21-27, 1997. The paper was originally edited as a contribution for the book J.
Kertesz and I. Kondor (Editors), Econophysics: an Emerging Science, Kluwer
Academic Publishers, Dordrecht (NL) that should contain selected papers presented at
that Workshop and should have appeared in 1998 or 1999. Unfortunately the book was
not published. The present e-print is a revised version (with up-date annotations and
references) of that unpublished contribution, but essentially represents our knowledge of
that early time.
http://arXiv.org/abs/0704.0320v1
Abstract
Fractional calculus allows one to generalize the linear, one-dimensional,
diffusion equation by replacing either the first time derivative or the second
space derivative by a derivative of fractional order. The fundamental
solutions of these generalized diffusion equations are shown to provide
probability density functions, evolving on time or variable in space, which
are related to the peculiar class of stable distributions. This property is
a noteworthy generalization of what happens for the standard diffusion
equation and can be relevant in treating financial and economical problems
where the stable probability distributions are known to play a key role.
1 Introduction
Non-Gaussian probability distributions are becoming more common as data
models, especially in economics where large fluctuations are expected. In
fact, probability distributions with heavy tails are often met in economics
and finance, which suggests to enlarge the arsenal of possible stochastic
models by non-Gaussian processes. This conviction started in the early
sixties after the appearance of a series of papers by Mandelbrot and
his associates, who point out the importance of non-Gaussian probability
distributions, formerly introduced by Pareto and Lévy, and related scaling
properties, to analyse economical and financial variables, as reported in
the recent book by Mandelbrot (1997). Some examples of such variables
are common stock prices changes, changes in other speculative prices, and
interest rate changes. In this respect many works by different authors have
recently appeared, see e.g. the recent books by Bouchaud & Potter (1997),
Mantegna & Stanley (1998) and the references therein quoted.
It is well known that the fundamental solution (or Green function) of
the Cauchy problem for the standard linear diffusion equation provides at
any time the probability density function (pdf) in space of the Gauss (or
normal) law. This law exhibits all moments finite thanks to its exponential
decay at infinity. In particular, the space variance of the Green function
is proportional to the first power of time, a noteworthy property that
can be understood by means of an unbiased random walk model for the
Brownian motion, see e.g. Feller (1957). Less known is the property for
which the fundamental solution of the Signalling problem for the same
diffusion equation, provides at any position a unilateral pdf in time, known
as Lévy law, using the terminology of Feller (1966-1973). Because of its
algebraic decay at infinity as t−3/2 , this law has all moments of integer
order divergent, and consequently its expectation value and variance are
infinite.
Both the Gauss and Lévy laws belong to the general class of stable probability
distributions, which are characterized by an index α (0 < α ≤ 2), called
index of stability or characteristic exponent. In particular, the index of the
Gauss law is 2 , whereas that of the Lévy law is 1/2 .
In this paper we consider two different generalizations of the diffusion
equation by means of fractional calculus, which allows us to replace either the
first time derivative or the second space derivative by a suitable fractional
derivative. Correspondingly, the generalized equation will be referred to
as the time-fractional diffusion equation or the symmetric, space-fractional
diffusion equation. Here we show how the fundamental solutions of this
equation for the Cauchy and Signalling problems provide probability density
functions related to certain stable distributions, so providing a natural
generalization of what occurs for the standard diffusion equation.
The plan of the paper is as follows. First of all, for the sake of convenience
and completeness, we provide the essential notions of Riemann-Liouville
Fractional Calculus and Lévy Stable Probability Distributions in Appendix
A and B, respectively.
In Section 2, we recall the basic results for the standard diffusion
equation concerning the fundamental solutions of the Cauchy and Signalling
problems. In particular we provide the derivation of these solutions by the
Fourier and Laplace transforms and the interpretation in terms of Gauss
and Lévy stable pdf , respectively.
In Section 3, we consider the time-fractional diffusion equation and we
formulate for it the basic Cauchy and Signalling problems to be treated in the
subsequent two sections. Here we adopt the Riemann-Liouville approach to
Fractional Calculus, and the related definition for the Caputo time-fractional
derivative of a causal function of time.
In Section 4, we solve the Cauchy problem for the time-fractional diffusion
equation by using the technique of Fourier transform and we derive the
corresponding fundamental solution in terms of a special function of Wright
type in the similarity variable. In this case the solution can be interpreted
as a noteworthy symmetric pdf in space with all moments finite, evolving
in time. In particular, its space variance turns out to be proportional to a
power of time equal to the order of the time-fractional derivative.
In Section 5, we derive the fundamental solution for the Signalling problem
of the time-fractional diffusion equation by using the technique of Laplace
transform. In this case the solution, still expressed in terms of a special
function of Wright type, can be interpreted as a unilateral stable pdf in
time, depending on position, with index of stability given by half of the
order of the time-fractional derivative.
In Section 6, we consider the symmetric, space-fractional diffusion equation.
Here we adopt the Riesz approach to Fractional Calculus, and the related
definition for the symmetric space-fractional derivative of a function of a
single space variable. Here we treat the Cauchy problem by technique
of Fourier transform and we derive the series representation of the
corresponding Green function. In this case the fundamental solution is
interpreted in terms of a symmetric stable pdf in space, evolving in time,
with index of stability given by the order of the space-fractional derivative.
To approximate such evolution we propose a random walk model, discrete
in space and time, which is based on the Grünwald-Letnikov approximation
of the fractional derivative.
Finally, Section 7 is devoted to conclusions and remarks on related work.
2 The standard diffusion equation
For the standard diffusion equation we mean the linear partial differential
equation
u(x, t) = D
u(x, t) , u = u(x, t) , (2.1)
where D denotes a positive constant with the dimensions L2 T−1 , x and t
are the space-time variables, and u = u(x, t) is the field variable, which is
assumed to be a causal function of time, i.e. vanishing for t < 0 .
The typical physical phenomenon related to such an equation is the heat
conduction in a thin solid rod extended along x , so the field variable u is
the temperature.
In order to guarantee the existence and the uniqueness of the solution,
we must equip (1.1) with suitable data on the boundary of the space-time
domain. The basic boundary-value problems for diffusion are the so-called
Cauchy and Signalling problems. In the Cauchy problem, which concerns
the space-time domain −∞ < x < +∞ , t ≥ 0 , the data are assigned at
t = 0+ on the whole space axis (initial data). In the Signalling problem,
which concerns the space-time domain x ≥ 0 , t ≥ 0 , the data are assigned
both at t = 0+ on the semi-infinite space axis x > 0 (initial data) and at
x = 0+ on the semi-infinite time axis t > 0 (boundary data); here, as mostly
usual, the initial data are assumed to be vanishing.
Denoting by g(x) and h(t) two given, sufficiently well-behaved functions, the
basic problems are thus formulated as following:
a) Cauchy problem
u(x, 0+) = g(x) , −∞ < x < +∞ ; u(∓∞, t) = 0 , t > 0 ; (2.2a)
b) Signalling problem
u(x, 0+) = 0 , x > 0 ; u(0+, t) = h(t) , u(+∞, t) = 0 , t > 0 . (2.2b)
Hereafter, for both the problems, we derive the classical results which will be
properly generalized for the fractional diffusion equation in the subsequent
sections.
Let us begin with the Cauchy problem. It is well known that this initial value
problem can be easily solved making use of the Fourier transform and its
fundamental solution can be interpreted as a Gaussian pdf in x. Adopting
the notation g(x) ÷ ĝ(κ) with κ ∈ R and
ĝ(κ) = F [g(x)] =
e+iκx g(x) dx ,
g(x) = F−1 [ĝ(κ)] = 1
e−iκx ĝ(κ) dκ ,
the transformed solution satisfies the ordinary differential equation of the
first order (
+ κ2 D
û(κ, t) = 0 , û(κ, 0+) = ĝ(κ) , (2.3)
and consequently it turns out to be
û(κ, t) = ĝ(κ) e−κ
2 D t . (2.4)
Then, introducing
Gdc (x, t) ÷ Ĝdc (κ, t) = e−κ
2 D t , (2.5)
where the upper index d refers to (standard) diffusion, the required solution,
obtained by inversion of (2.4), can be expressed in terms of the space
convolution u(x, t) =
−∞ Gdc (ξ, t) g(x − ξ) dξ , where
Gdc (x, t) =
t−1/2 e−x
2/(4D t) . (2.6)
Here Gdc (x, t) represents the fundamental solution (or Green function) of
the Cauchy problem, since it corresponds to g(x) = δ(x) . It turns out
to be a function in x , even and normalized, i.e. Gdc (x, t) = Gdc (|x|, t) and∫ +∞
−∞ Gdc (x, t) dx = 1 . We also note the identity
|x| Gdc (|x|, t) =
Md(ζ) , (2.7)
where ζ = |x|/(
D t1/2) is the well-known similarity variable and
Md(ζ) =
2/4 . (2.8)
We note that Md(ζ) satisfies the normalization condition
d(ζ) dζ = 1 .
The interpretation of the Green function (2.6) in probability theory is
straightforward since we easily recognize
Gdc (x, t) = pG(x;σ) :=
2/(2σ2) , σ2 = 2D t , (2.9)
where pG(x;σ) denotes the well-known Gauss or normal pdf spread out
over all real x (the space variable), whose moment of the second order, the
variance, is σ2 . The associated cumulative distribution function (cdf) is
known to be
PG(x;σ) :=
′;σ) dx′ =
1 + erf
, (2.10)
where erf (z) := (2/
0 exp (−u2) du denotes the error function.
Furthermore, the moments of even order of the Gauss pdf turn out to be∫ +∞
2n pG(x;σ) dx = (2n − 1)!!σ2n , so
x2n Gdc (x, t) dx = (2n − 1)!! (2D t)n , n = 1, 2, . . . . (2.11)
Let us now consider the Signalling problem. This initial-boundary value
problem can be easily solved by making use of the Laplace transform.
Adopting the notation h(t) ÷ h̃(s) with s ∈ C and
h̃(s) = L [h(t)] =
e−st h(t) dt ,
h(t) = L−1
h̃(t)
est h̃(s) ds ,
where Br denotes the Bromwich path, the transformed solution of the
diffusion equation satisfies the ordinary differential equation of the second
order
ũ(x, s) = 0 , ũ(0+, s) = h̃(s) , ũ(+∞, s) = 0 . (2.12)
and consequently it turns out to be
ũ(x, s) = h̃(s) e−(x/
D) s1/2 . (2.13)
Then introducing
Gds (x, t) ÷ G̃ds (x, s) = e−(x/
D) s1/2 , (2.14)
the required solution, obtained by inversion of (2.13), can be expressed in
terms of the time convolution, u(x, t) =
0 Gds (x, τ)h(t − τ) dτ , where
Gds (x, t) =
t−3/2 e−x
2/(4D t) . (2.15)
Here Gds (x, t) represents the fundamental solution (or Green function) of the
Signalling problem, since it corresponds to h(t) = δ(t) . We note that
Gds (x, t) = pLS(t;µ) :=
2π t3/2
e−µ/(2t) , t ≥ 0 , µ = x
, (2.16)
where pLS(t;µ) denotes the one-sided Lévy-Smirnov pdf spread out over all
non negative t (the time variable). The associated cdf is, see e.g. Feller
(1966-1971) and Prüss (1993),
PL(t;µ) :=
′;µ) dt′ = erfc
= erfc
, (2.17)
where erfc (z) := 1 − erf (z) denotes the complenatary error function.
The Lévy-Smirnov pdf has all moments of integer order infinite, since it
decays at infinity as t−3/2 . However, we note that the absolute moments of
real order ν are finite only if 0 ≤ ν < 1/2 . In particular, for this pdf the mean
is infinite, for which we can take the median as expectation value. From
PLs(tmed;µ) = 1/2 , it turns out that tmed ≈ 2µ , since the complementary
error function gets the value 1/2 as its argument is approximatively 1/2.
We note that in the common domain x > 0 , t > 0 the Green functions of
the two basic problems satisfy the identity
xGdc (x, t) = tGds (x, t) , (2.18)
that we refer to as the reciprocity relation between the two fundamental
solutions of the diffusion equation. Furthermore, in view of (2.7) and (2.18)
we recognize the role of the function of the similarity variable, Md(ζ) ,
in providing the two fundamental solutions; we shall refer to it as to the
normalized auxiliary function of the diffusion equation for both the Cauchy
and Signalling problems.
3 The time-fractional diffusion equation
By the time-fractional diffusion equation we mean the linear evolution
equation obtained from the classical diffusion equation by replacing the first-
order time derivative by a fractional derivative (in the Caputo sense) of order
α with 0 < α ≤ 2. In our notation it reads
, u = u(x, t) , 0 < α ≤ 2 , (3.1)
where D denotes a positive constant with the dimensions L2 T−α . From
Appendix A we recall the definition of the Caputo fractional derivative of
order α > 0 for a (sufficiently well-behaved) causal function f(t) , see (A.9),
Dα∗ f(t) :=
Γ(m − α)
(t − τ)m−α f (m)(τ) dτ , (3.2)
where m = 1, 2, . . . , and 0 ≤ m − 1 < α ≤ m . According to (3.2) we thus
need to distinguish the cases 0 < α ≤ 1 and 1 < α ≤ 2 . In the the latter case
(3.1) may be seen as a sort of interpolation between the standard diffusion
equation and the standard wave equation. Introducing
Φλ(t) :=
tλ−1+
, λ > 0 , (3.3)
where the suffix + is just denoting that the function is vanishing for t < 0 ,
we easily recognize that the equation (3.1) assumes the explicit forms :
if 0 < α ≤ 1 ,
Φ1−α(t) ∗
Γ(1 − α)
(t − τ)−α
dτ = D
; (3.4)
if 1 < α ≤ 2 ,
Φ2−α(t) ∗
Γ(2 − α)
(t − τ)1−α
dτ = D ∂
. (3.5)
Extending the classical analysis for the standard diffusion equation (2.1) to
the above integro-differential equations (3.4-5), the Cauchy and Signalling
problems are thus formulated as in equations (2.2), i.e.
a) Cauchy problem
u(x, 0+) = g(x) , −∞ < x < +∞ ; u(∓∞, t) = 0 , t > 0 ; (3.6a)
b) Signalling problem
u(x, 0+) = 0 , x > 0 ; u(0+, t) = h(t) , u(+∞, t) = 0 , t > 0 . (3.6b)
However, if 1 < α ≤ 2 , the presence in (3.5) of the second order time
derivative of the field variable requires to specify the initial value of the first
order time derivative ut(x, 0
+) , since in this case two linearly independent
solutions are to be determined. To ensure the continuous dependence of our
solution on the parameter α also in the transition from α = 1− to α = 1+ ,
we agree to assume ut(x, 0
+) = 0 .
We recognize that our fractional diffusion equation (3.1), when subject to
the conditions (3.6), is equivalent to the integro-differential equation
u(x, t) = g(x) +
(t − τ)α−1
dτ , (3.7)
where 0 < α ≤ 2 . Such integro-differential equation has been investigated
by several authors, including Schneider & Wyss (1989), Fujita (1990), Prüss
(1993) and Engler (1997).
In view of our subsequent analysis we find it convenient to put
, 0 < ν < 1 . (3.8)
In fact the analysis of the time-fractional diffusion equation turns out to
be easier if we adopt as a key parameter the half of the order of the
time-fractional derivative. In future we shall provide the symbol α with
other relevant meanings, as the index of stability of a stable probability
distribution or the order of the space derivative in the space-fractional
diffusion equation.
Henceforth, we agree to insert the parameter ν in the field variable, i.e.
u = u(x, t; ν) . By denoting the Green functions of the Cauchy and Signalling
problems by Gc(x, t; ν) and Gs(x, t; ν) , respectively, the solutions of the two
basic problems are obtained by a space or time convolution, u(x, t; ν) =∫ +∞
−∞ Gc(ξ, t; ν) g(x−ξ) dξ , u(x, t; ν) =
0 Gs(x, τ ; ν)h(t−τ) dτ , respectively.
It should be noted that Gc(x, t; ν) = Gc(|x|, t; ν) , since the Green function
turns out to be an even function of x .
In the following two sections we shall compute the two fundamental solutions
with the same techniques (based on Fourier and Laplace transforms) used
for the standard diffusion equation and we shall provide their interpretation
in terms of probability distributions. Most of the presented results are based
on the papers by Mainardi (1994), (1995), (1996), (1997) and by Mainardi
& Tomirotti (1995), (1997).
4 The Cauchy problem for the time-fractional
diffusion equation
For the fractional diffusion equation (3.1) subject to (3.6a) the application
of the Fourier transform leads to the ordinary differential equation of order
α = 2ν ,
+ κ2 D
û(κ, t; ν) = 0 , û(κ, 0+; ν) = ĝ(κ) , (4.1)
Using the results of Appendix A, see (A.22-30), the transformed solution is
û(κ, t; ν) = ĝ(κ)E2ν
−κ2 D t2ν
, (4.2)
where E2ν(·) denotes the Mittag-Leffler function of order 2ν , and conse-
quently for the Green function we have
Gc(x, t; ν) = Gc(|x|, t; ν) ÷ Ĝc(k, t; ν) = E2ν
−κ2D t2ν
. (4.3)
Since the Green function is a real and even function of x, its (exponential)
Fourier transform can be expressed in terms of the cosine Fourier transform
and thus is related to its spatial Laplace transform as follows
Ĝc(k, t; ν) = 2
Gc(x, t; ν) cos κx dx =
G̃c(s, t; ν)
s=+ik
+ G̃c(s, t; ν)
s=−ik
(4.4)
Indeed, a split occurs also in (4.3) according to the duplication formula for
the Mittag-Leffler function, see (A.26),
Ĝc(k, t; ν) = E2ν(−κ2 D t2ν) =
[Eν(+iκ
D tν) + Eν(−iκ
D tν)]/2 .
(4.5)
When ν 6= 1/2 the inversion of the Fourier transform in (4.5) cannot be
obtained by using a standard table of Fourier transform pairs; however, for
any ν ∈ (0, 1) such inversion can be achieved by appealing to the Laplace
transform pair (A.37) with r = |x| , and s = ±iκ . In fact, taking into
account the scaling property of the Laplace transform, we obtain from (4.5)
and (A.37)
Gc(|x|, t; ν) =
( |x|√
, (4.6)
where M(ζ; ν) is the special function of Wright type, defined by (A.31-33),
, (4.7)
the similarity variable. We note the identity
|x| Gc(|x|, t; ν) =
M(ζ; ν) , (4.8)
which generalizes to the time-fractional diffusion equation the identity (2.7)
of the standard diffusion equation. Since
0 M(ζ; ν) dζ = 1 , see (A.40),
the function M(ζ; ν) is the normalized auxiliary function of the fractional
diffusion equation.
We note that for the time-fractional diffusion equation the fundamental
solution of the Cauchy problem is still a bilateral symmetric pdf in x (with
two branches, for x > 0 and x < 0 , obtained one from the other by
reflection), but is no longer of Gaussian type if ν 6= 1/2 . In fact, for large
|x| each branch exhibits an exponential decay in the ”stretched” variable
|x|1/(1−ν) as can be derived from the asymptotic representation (A.36) of the
auxiliary function M(·; ν) . In fact, by using (4.7-8) and (A.36), we obtain
Gc(x, t; ν) ∼ a∗(t) |x|(ν−1/2)/(1−ν) exp
−b∗(t)|x|1/(1−ν)
, (4.9)
as |x| → ∞ , where a∗(t) and b∗(t) are certain positive functions of time.
Furthermore, the exponential decay in x provided by (4.9) ensures that all
the absolute moments of positive order of Gc(x, t; ν) are finite. In particular,
using (4.8) and (A.39) it turns out that the moments (of even order) are
x2n Gc(x, t; ν) dx =
Γ(2n + 1)
Γ(2νn + 1)
(Dt2ν)n , n = 0 , 1 , 2 , . . . (4.10)
The formula (4.10) provides a generalization of the corresponding formula
(2.11) valid for the standard diffusion equation, ν = 1/2 . Furthermore, we
recognize that the variance associated to the pdf is now proportional to Dt2ν ,
which for ν 6= 1/2 implies a phenomenon of anomalous diffusion. According
to a usual terminology in statistical mechanics, the anomalous diffusion is
said to be slow if 0 < ν < 1/2 and fast if 1/2 < ν < 1 .
In Figure 1, as an example, we compare versus |x| , at fixed t , the
fundamental solutions of the Cauchy problem with different ν (ν =
1/4 , 1/2 , 3/4 ). We consider the range 0 ≤ |x| ≤ 4 and assume D = t = 1 .
0 1 2 3 4
Figure 1: The Cauchy problem for the time-fractional diffusion equation.
The fundamental solutions versus |x| with a) ν = 1/4 , b) ν = 1/2 , c)
ν = 3/4 .
We note the different behaviour of the pdf in the cases of slow diffusion (ν =
1/4 ) and fast diffusion (ν = 3/4 ) with respect to the Gaussian behaviour
of the standard diffusion (ν = 1/2). In the limiting cases ν = 0 and ν = 1
we have
Gc(x, t; 0) =
e−|x|
, Gc(x, t; 1) =
δ(x −
D t) + δ(x +
. (4.11)
We also recognize from the appendix B that for 1/2 ≤ ν < 1 any branch
of the fundamental solution is proportional to the corresponding positive
branch of an extremal stable pdf with index of stability α = 1/ν , which
exhibits an exponential decay at infinity. In fact, applying (B.29) with
α = 1/ν and y = ζ = |x|/(
Dtν) , from (4.7-8) we obtain
Gc(|x|, t; ν) =
|x|/(
D tν) ; − (2 − 1/ν)
· p1/ν (|x|; +1, 1, 0) , 1 < 1/ν ≤ 2 .
(4.12)
We also note that the stable distribution in (4.12) satisfies the condition
p1/ν (x; +1, 1, 0) dx = ν , 1 < 1/ν ≤ 2 . (4.13)
5 The Signalling problem for the time-fractional
diffusion equation
For the fractional diffusion equation (3.1) subject to (3.6b) the application
of the Laplace transform leads to the ordinary differential equation of order
ũ(x, s; ν) , ũ(0+, s; ν) = h̃(s) , ũ(+∞, s; ν) = 0 . (5.1)
Thus the transformed solution reads
ũ(x, s; ν) = h̃(s) e−(x/
D) sν , (5.2)
so for the Green function we have
Gs(x, t; ν) ÷ G̃s(x, s; ν) = e−(x/
D) sν . (5.3)
When ν 6= 1/2 the inversion of this Laplace transform cannot be obtained by
looking in a standard table of Laplace transform pairs. Also here we appeal
to a Laplace transform pair related to the Wright-type function M(ζ; ν). In
fact, using (A.40) with r = t , and taking into account the scaling property
of the Laplace transform, we obtain
Gs(x, t; ν) = ν
D t1+ν
. (5.4)
Introducing the similarity variable ζ = x/(
Dtν) , we recognize the identity
tGs(x, t; ν) = ν ζ M(ζ; ν) , (5.5)
which is the counterpart for the Signalling problem of the identity (4.8) valid
for the Cauchy problem.
Comparing (5.5) with (4.8) we obtain the reciprocity relation between the
two fundamental solutions of the time-fractional diffusion equation, in the
common domain x > 0 , t > 0 ,
2ν xGc(x, t; ν) = tGs(x, t; ν) . (5.6)
The interpretation of Gs(x, t; ν) as a one-sided stable pdf in time is
straightforward: in this respect we need to apply (B.28), with index of
stability α = ν and variable y = ζ−1/ν = t (
D/x)1/ν , in (5.5). We obtain
Gs(x, t; ν) =
; − ν
 = pν (t; −1, 1, 0) . (5.7)
In Figure 2, as an example, we compare versus t , at fixed x , the fundamental
solutions of the Signalling problem with different ν (ν = 1/4 , 1/2 , 3/4 ). We
consider the range 0 ≤ t ≤ 3 and assume D = x = 1 .
We note the different behaviour of the pdf in the cases of slow diffusion
(ν = 1/4 ) and fast diffusion (ν = 3/4 ) with respect to the Lévy pdf for the
standard diffusion (ν = 1/2). In the limiting cases ν = 0 , 1 , we have
Gs(x, t; 0) = δ(t) , Gs(x, t; 1) = δ(t − x/
D) . (5.8)
0 1 2 3
Figure 2: The Signalling problem for the time-fractional diffusion equation.
The fundamental solutions versus t with a) ν = 1/4 , b) ν = 1/2 , c)
ν = 3/4 .
6 The Cauchy problem for the symmetric space-
fractional diffusion equation
The symmetric space-fractional diffusion equation is obtained from the
classical diffusion equation by replacing the second-order space derivative by
a symmetric space-fractional derivative (explained below) of order α with
0 < α ≤ 2 . In our notation we write this equation as
∂|x|α
, u = u(x, t;α) , x ∈ R , t ∈ R+0 , 0 < α ≤ 2 , (6.1)
where D is a positive coefficient with the dimensions Lα T−1 . The
fundamental solution for the Cauchy problem, Gc(x, t;α) is the solution of
(6.1), subject to the initial condition u(x, 0+;α) = δ(x) .
The symmetric space-fractional derivative of any order α > 0 of a sufficiently
well-behaved function φ(x) , x ∈ R , may be defined as the pseudo-
differential operator characterized in its Fourier representation by
d|x|α
φ(x) ÷ −|κ|α φ̂(κ) , x , k ∈ R , α > 0 . (6.2)
According to a usual terminology, −|κ|α is referred to as the symbol of our
pseudo-differential operator, the symmetric space-fractional derivative, of
order α . Here, we have adopted the notation introduced by Zaslavski, see
e.g. Saichev & Zaslavski (1997).
In order to properly introduce this kind of fractional derivative we need
to consider a peculiar approach to fractional calculus different from the
Riemann-Liouville one, already treated in Appendix A. This approach is
indeed based on the so-called Riesz potentials (or integrals), that we prefer
to consider later.
At first, let us see how things become highly transparent by using an
heuristic argument, originally due to Feller (1952). The idea is to start
from the positive definite differential operator
A := −
÷ κ2 = |κ|2 , (6.3)
whose symbol is |κ|2 , and form positive powers of this operator as pseudo-
differential operators by their action in the Fourier-image space, i.e.
Aα/2 :=
= |κ|α α > 0 . (6.4)
Thus the operator −Aα/2 can be interpreted as the required fractional
derivative, i.e.
Aα/2 ≡ − d
d|x|α
, α > 0 . (6.5)
We note that the operator just defined must not be confused with a power
of the first order differential operator d
for which the symbol is −iκ .
After the above considerations it is straightforward to obtain the Fourier
image of the Green function of the Cauchy problem for the space-fractional
diffusion equation. In fact, applying the Fourier transform to the equation
(6.1), subject to the initial condition u(x, 0+;α) = δ(x) , and accounting for
(6.2), we obtain
Gc(x, t;α) = Gc(|x|, t;α) ÷ Ĝc(k, t;α) = e−D t |κ|
, 0 < α ≤ 2 . (6.6)
We easily recognize that the Fourier transform of the Green function
corresponds to the canonic form of a symmetric stable distribution of index
of stability α and scaling factor γ = (Dt)1/α , see (B.8). Therefore we have
Gc(x, t;α) = pα(x; 0, γ, 0) , γ = (Dt)1/α . (6.7)
For α = 1 and α = 2 we easily obtain the explicit expressions of the
corresponding Green functions since in these cases they correspond to the
Cauchy and Gauss distributions,
Gc(x, t; 1) =
x2 + (D t)2
, (6.8)
see (B.5), and
Gc(x, t; 2)) =
2/(4D t) , (6.9)
in agreement with (2.6).
We easily recognize that
(D t)1/α
(6.10)
is the similarity variable for the space-fractional diffusion equation, in terms
of which we can express the Green function for any α ∈ (0, 2] . Indeed, we
recognize that
Gc(x, t;α) =
(D t)1/α
qα(η; 0) , (6.11)
where qα(η; 0) denotes the symmetric stable distribution of order α with
Feller-type characteristic function, see (B.14-15). Now we can express the
Green function using the Feller series expansions (B.21-22) with θ = 0 . We
obtain:
for 0 < α < 1 ,
qα(η; 0) = −
Γ(nα + 1)
, (6.12a)
for 1 < α ≤ 2 ,
qα(η; 0) =
(−1)m
Γ[(2m + 1)/α]
(2m)!
η2m . (6.12b)
In the limiting case α = 1 the above series reduce to geometrical series and
therefore are no longer convergent in all of C . In particular, they represent
the expansions of the function q1(η; 0) = 1/[π(1+η
2)] , convergent for η > 1
and 0 < η < 1 , respectively.
We also note that for any α ∈ (0, 2] the functions qα(η; 0) exhibit at the
origin the value qα(0; 0) = Γ(1/α)/(π α) , and at the queues, excluding the
Gaussian case α = 2 , the algebraic asymptotic behaviour, as η → ∞ ,
qα(η; 0) ∼
Γ(α + 1) sin
η−(α+1) , 0 < α < 2 . (6.13)
In Figure 3, as an example, we compare versus x , at fixed t , the fundamental
solutions of the Cauchy problem with different α (α = 1/2 , 1 , 3/2 , 2 ). We
consider the range −6 ≤ x ≤ +6 and assume D = t = 1 .
-6 -4 -2 0 2 4 6
-6 -4 -2 0 2 4 6
Figure 3: The Cauchy problem for the simmetric space-fractional diffusion
equation. The fundamental solutions versus x : plate a) α = 1/2
(continuous line), α = 1 (dashed line); plate b) α = 3/4 (continuous line),
α = 2 (dashed line).
Let us now express more properly our operator (6.4) (with symbol |κ|α)
as inverse of a suitable integral operator Iα whose symbol is |κ|−α . This
operator can be found in the approach by Marcel Riesz to Fractional
Calculus, see e.g. Samko, Kilbas & Marichev (1987-1993) and Rubin (1996).
We recall that for any α > 0 , α 6= 1 , 3 , 5 , . . . and for a sufficiently well-
behaved function φ(x) , x ∈ R , the Riesz integral or Riesz potential Iα and
its image in the Fourier domain read
Iα φ(x) :=
2Γ(α) cos(πα/2)
|x − ξ|α−1 φ(ξ) dξ ÷ φ̂(κ)
. (6.14)
On its turn, the Riesz potential can be written in terms of two Weyl integrals
Iα± according to
Iα φ(x) =
2 cos(πα/2)
Iα+φ(x) + I
−φ(x)
, (6.15)
where 
Iα+ φ(x) :=
(x − ξ)α−1 φ(ξ) dξ ,
Iα− φ(x) :=
(ξ − x)α−1 φ(ξ) dξ .
(6.16)
Then, at least in a formal way, the space-fractional derivative (6.2) turns
out to be defined as the opposite of the (left) inverse of the Riesz fractional
integral, i.e.
d|x|α
φ(x) := −I−α φ(x) = −
2 cos(πα/2)
I−α+ φ(x) + I
− φ(x)
. (6.17)
Notice that (6.14) and (6.17) become meaningless when α is an integer odd
number. However, for our range of interest 0 < α ≤ 2 , the particular case
α = 1 can be singled out since the corresponding Green function is already
known, see (6.8). Thus, excluding the case α = 1 , our space-fractional
diffusion equation (6.1) can be re-written, x ∈ R , t ∈ R+0 , as
= −D I−α u , u = u(x, t;α) , 0 < α ≤ 2 , α 6= 1 , (6.18)
where the operator I−α is defined by (6.16-17).
Here, in order to evaluate the fundamental solution of the Cauchy problem,
interpreted as a probability density, we propose a numerical approach,
original as far as we know, based on a (symmetric) random walk model,
discrete in space and time, see also Gorenflo & Mainardi (1998a), Gorenflo
& Mainardi (1998b) and Gorenflo, De Fabritiis & Mainardi (1999). We shall
see how things become highly transparent, in that we properly generalize
the classical random-walk argument of the standard diffusion equation
to our space-fractional diffusion equation (6.18). So doing we are in
position to provide a numerical simulation of the related (symmetric) stable
distributions in a way analogous to the standard one for the Gaussian law.
The essential idea is to approximate the left inverse operators I−α± by the
Grünwald-Letnikov scheme, on which the reader can inform himself in the
treatises on fractional calculus, see e.g. Oldham & Spanier (1974), Samko,
Kilbas & Marichev (1987-1993), Miller & Ross (1993), or in the recent review
article by Gorenflo (1997). If h denotes a ”small” positive step-length, these
approximating operators read
± φ(x) :=
(−1)k
φ(x ∓ kh) . (6.19)
Assume, for simplicity, D = 1 , and introduce grid points xj = j h with
h > 0 , j ∈ Z , and time instances tn = n τ with τ > 0 , n ∈ N0 . Let there
be given probabilities pj,k ≥ 0 of jumping from point xj at instant tn to
point xk at instant tn+1 and define probabilities yj(tn) of the walker being
at point xj at instant tn. Then, by
yk(tn+1) =
pj,k uj(tn) ,
pj,k =
pj,k = 1 , (6.20)
with pj,k = pk,j , a symmetric random walk (more precisely a symmetric
random jump) model is described. With the approximation
yj(tn) ≈
∫ (xj+h/2)
(xj−h/2)
u(x, tn) dx ≈ hu(xj , tn) , (6.21)
and introducing the ”scaling parameter”
2 | cos(απ/2)|
, (6.22)
we have solved
yj(tn+1) − yj(tn)
= − hI−α yj(tn) , (6.23)
for yj(tn+1) . So we have proved to have a consistent (for h → 0) symmetric
random walk approximation to (6.18) by taking
i) for 0 < α < 1 , 0 < µ ≤ 1/2 ,
−α yj(tn) = µ
+ yj(tn) + hI
− yj(tn)
pj,j = 1 − 2µ , pj,j±k = µ
)∣∣ , k ≥ 1 ;
(6.24)
ii) for 1 < α ≤ 2 , 0 < µ ≤ 1/(2α) ,


−α yj(tn) = µ
+ yj+1(tn) + hI
− yj−1(tn)
pj,j = 1 − 2µ α , pj,j±1 = µ
pj,j±k = µ
)∣∣∣ , k ≥ 2 .
(6.25)
We note that our random walk model is not only symmetric, but also
homogeneous, the transition probabilities pj,j±k not depending on the index
In the special case α = 2 we recover from (6.25) the well-known three-point
approximation of the heat equation, because pj,j±k = 0 for k ≥ 2 . This
means that for approximation of common diffusion only jumps of one step
to the right or one to the left or jumps of width zero occur, whereas for
0 < α < 2 (α 6= 1) arbitrary large jumps occur with power-like decaying
probability, as it turns out from the asymptotic analysis for the transition
probabilities given in (6.24-25). In fact, as k → ∞ , one finds
pj,j+k ∼
(τ/hα)
Γ(α + 1) sin
k−(α+1) , 0 < α < 2 . (6.26)
This result thus provides the discrete counterpart of the asymptotic
behaviour of the long power-law tails of the symmetric stable distributions,
as foreseen by (6.13) when 0 < α < 2 .
7 Conclusions
We have treated two generalizations of the standard, one-dimensional,
diffusion equation, namely, the time-fractional diffusion equation and the
symmetric space-fractional diffusion equation. For these equations we have
derived the fundamental solutions using the transform methods of Fourier
and Laplace, and exhibited their connections to extremal and symmetric
stable probability densities, evolving on time or variable in space. For the
symmetric space-fractional diffusion equation we have presented a stationary
(in time), homogeneous (in space) symmetric random walk model, discrete
in space and time, the step-lengths of the spatial grid and the time lapses
between transitions properly scaled. In the limit of infinitesimally fine
discretization this model (based on the Grünwald-Letnikov approximation
to fractional derivatives) is consistent with the continuous diffusion process,
i.e. convergent if interpreted as a difference scheme in the sense of numerical
analysis2.
From the mathematical viewpoint the field of such ”fractional” general-
izations is fascinating as there several mathematical disciplines meet and
come to a fruitful interplay: e.g. probability theory and stochastic processes,
2Further generalizations have been considered by us and our collaborators in other
papers, in which we have given a derivation of discrete random walk models related to
more general space-time fractional diffusion equations. For a comprehensive analysis, see
Gorenflo et al. (2002). Readers interested to the fundamental solutions of these fractional
diffusion equations are referred to the paper by Mainardi et al. (2001) where analytical
expressions and numerical plots are found.
integro-differential equations, transform theory, special functions, numerical
analysis. As one may take from our References, one can observe that since
some decades there is an ever growing interest in using the concepts of
fractional calculus among physicists and economists. Among economists we
like to refer the reader to a collection of papers on the topic of ”Fractional
Differencing and Long Memory Processes”, edited by Baillie & King (1996).
Appendix A: The Riemann-Liouville Fractional
Calculus
Fractional calculus is the field of mathematical analysis which deals with the
investigation and applications of integrals and derivatives of arbitrary order.
The term fractional is a misnomer, but it is retained following the prevailing
use. This appendix is mostly based on the recent review by Gorenflo &
Mainardi (1997). For more details on the classical treatment of fractional
calculus the reader is referred to Erdélyi (1954), Oldham & Spanier (1974),
Samko et al. (1987-1993) and Miller & Ross (1993).
According to the Riemann-Liouville approach to fractional calculus, the
notion of fractional Integral of order α (α > 0) is a natural consequence
of the well known formula (usually attributed to Cauchy), that reduces the
calculation of the n−fold primitive of a function f(t) to a single integral of
convolution type. In our notation the Cauchy formula reads
Jnf(t) := fn(t) =
(n − 1)!
(t − τ)n−1 f(τ) dτ , t > 0 , n ∈ N , (A.1)
where N is the set of positive integers. From this definition we note that
fn(t) vanishes at t = 0 with its derivatives of order 1, 2, . . . , n − 1 . For
convention we require that f(t) and henceforth fn(t) be a causal function,
i.e. identically vanishing for t < 0. In a natural way one is led to extend
the above formula from positive integer values of the index to any positive
real values by using the Gamma function. Indeed, noting that (n − 1)! =
Γ(n) , and introducing the arbitrary positive real number α , one defines the
Fractional Integral of order α > 0 :
Jα f(t) :=
(t − τ)α−1 f(τ) dτ , t > 0 , α ∈ R+ , (A.2)
where R+ is the set of positive real numbers. For complementation we define
J0 := I (Identity operator), i.e. we mean J0 f(t) = f(t) . Furthermore, by
Jαf(0+) we mean the limit (if it exists) of Jαf(t) for t → 0+ ; this limit
may be infinite.
We note the semigroup property JαJβ = Jα+β , α , β ≥ 0 , which implies
the commutative property JβJα = JαJβ , and the effect of our operators Jα
on the power functions
Jαtγ =
Γ(γ + 1)
Γ(γ + 1 + α)
tγ+α , α ≥ 0 , γ > −1 , t > 0 . (A.3)
These properties are of course a natural generalization of those known when
the order is a positive integer.
Introducing the Laplace transform by the notation L {f(t)} :=∫∞
−st f(t) dt = f̃(s) , s ∈ C , and using the sign ÷ to denote a Laplace
transform pair, i.e. f(t) ÷ f̃(s) , we note the following rule for the Laplace
transform of the fractional integral,
Jα f(t) ÷ f̃(s)
, α ≥ 0 , (A.4)
which is the generalization of the case with an n-fold repeated integral.
After the notion of fractional integral, that of fractional derivative of order
α (α > 0) becomes a natural requirement and one is attempted to substitute
α with −α in the above formulas. However, this generalization needs some
care in order to guarantee the convergence of the integrals and preserve the
well known properties of the ordinary derivative of integer order.
Denoting by Dn with n ∈ N , the operator of the derivative of order n ,
we first note that Dn Jn = I , Jn Dn 6= I , n ∈ N , i.e. Dn is left-inverse
(and not right-inverse) to the corresponding integral operator Jn . In fact
we easily recognize from (A.1) that
Jn Dn f(t) = f(t) −
f (k)(0+)
, t > 0 . (A.5)
As a consequence we expect that Dα is defined as left-inverse to Jα. For
this purpose, introducing the positive integer m such that m − 1 < α ≤ m ,
one defines the Fractional Derivative of order α > 0 :
Dα f(t) := Dm Jm−α f(t) , m − 1 < α ≤ m , m ∈ N , (A.6)
namely
Dα f(t)=
Γ(m − α)
(t − τ)α+1−m
, m − 1 < α < m,
f(t) , α = m.
(A.6′)
Defining for complementation D0 = J0 = I , then we easily recognize that
Dα Jα = I , α ≥ 0 , and
Dα tγ =
Γ(γ + 1)
Γ(γ + 1 − α)
tγ−α , α ≥ 0 , γ > −1 , t > 0 . (A.7)
Of course, these properties are a natural generalization of those known when
the order is a positive integer.
Note the remarkable fact that the fractional derivative Dα f is not zero for
the constant function f(t) ≡ 1 if α 6∈ N . In fact, (A.7) with γ = 0 teaches
us that
Dα1 =
Γ(1 − α)
, α ≥ 0 , t > 0 . (A.8)
This, of course, is ≡ 0 for α ∈ N, due to the poles of the gamma function in
the points 0,−1,−2, . . .. We now observe that an alternative definition of
fractional derivative, originally introduced by Caputo (1967) (1969) in the
late sixties and adopted by Caputo and Mainardi (1971) in the framework
of the theory of Linear Viscoelasticity, is
Dα∗ f(t) := J
m−α Dm f(t) m − 1 < α ≤ m , m ∈ N , (A.9)
namely
D ∗α f(t) =
Γ(m − α)
f (m)(τ)
(t − τ)α+1−m
dτ , m − 1 < α < m,
f(t) , α = m.
(A.9′)
This definition is of course more restrictive than (A.6), in that requires
the absolute integrability of the derivative of order m. Whenever we use
the operator Dα∗ we (tacitly) assume that this condition is met. We easily
recognize that in general
Dα f(t) := Dm Jm−α f(t) 6= Jm−α Dm f(t) := Dα∗ f(t) , (A.10)
unless the function f(t) along with its first m − 1 derivatives vanishes at
t = 0+. In fact, assuming that the passage of the m-derivative under the
integral is legitimate, one recognizes that, for m − 1 < α < m and t > 0 ,
Dα f(t) = Dα∗ f(t) +
Γ(k − α + 1)
f (k)(0+) , (A.11)
and therefore, recalling the fractional derivative of the power functions (A.7),
f(t) −
f (k)(0+)
= Dα∗ f(t) . (A.12)
The alternative definition (A.9) for the fractional derivative thus incorpo-
rates the initial values of the function and of its integer derivatives of lower
order. The subtraction of the Taylor polynomial of degree m − 1 at t = 0+
from f(t) means a sort of regularization of the fractional derivative. In
particular, according to this definition, the relevant property for which the
fractional derivative of a constant is still zero can be easily recognized, i.e.
Dα∗ 1 ≡ 0 , α > 0 . (A.13)
We now explore the most relevant differences between the two fractional
derivatives (A.6) and (A.9). We agree to denote (A.9) as the Caputo
fractional derivative to distinguish it from the standard Riemann-Liouville
fractional derivative (A.6). We observe, again by looking at (A.7), that
Dαtα−1 ≡ 0 , α > 0 , t > 0 .
From above we thus recognize the following statements about functions
which for t > 0 admit the same fractional derivative of order α , with
m − 1 < α ≤ m , m ∈ N ,
Dα f(t) = Dα g(t) ⇐⇒ f(t) = g(t) +
α−j , (A.14)
Dα∗ f(t) = D
∗ g(t) ⇐⇒ f(t) = g(t) +
m−j . (A.15)
In these formulas the coefficients cj are arbitrary constants.
For the two definitions we also note a difference with respect to the formal
limit as α → (m − 1)+ ; from (A.6) and (A.9) we obtain respectively,
Dα f(t) → Dm J f(t) = Dm−1 f(t) ; (A.16)
Dα∗ f(t) → J Dm f(t) = Dm−1 f(t) − f (m−1)(0+) . (A.17)
We now consider the Laplace transform of the two fractional derivatives.
For the standard fractional derivative Dα the Laplace transform, assumed to
exist, requires the knowledge of the (bounded) initial values of the fractional
integral Jm−α and of its integer derivatives of order k = 1, 2, . . . ,m−1 . The
corresponding rule reads, in our notation,
Dα f(t) ÷ sα f̃(s) −
Dk J (m−α) f(0+) sm−1−k , (A.18)
where m − 1 < α ≤ m .
The Caputo fractional derivative appears more suitable to be treated by
the Laplace transform technique in that it requires the knowledge of the
(bounded) initial values of the function and of its integer derivatives of
order k = 1, 2, . . . ,m− 1 , in analogy with the case when α = m . In fact, by
using (A.4) and noting that
Jα Dα∗ f(t) = J
α Jm−α Dm f(t) = Jm Dm f(t) = f(t) −
f (k)(0+)
(A.19)
we easily prove the following rule for the Laplace transform,
Dα∗ f(t) ÷ sα f̃(s) −
f (k)(0+) sα−1−k , m − 1 < α ≤ m . (A.20)
Indeed, the result (A.20), first stated by Caputo (1969) by using the
Fubini-Tonelli theorem, appears as the most ”natural” generalization of the
corresponding result well known for α = m .
Gorenflo and Mainardi (1997) have pointed out the major utility of the
Caputo fractional derivative in the treatment of differential equations of
fractional order for physical applications. In fact, in physical problems,
the initial conditions are usually expressed in terms of a given number
of bounded values assumed by the field variable and its derivatives of
integer order, no matter if the governing evolution equation may be a
generic integro-differential equation and therefore, in particular, a fractional
differential equation3.
We now analyze the most simple differential equations of fractional order,
including those which, by means of fractional derivatives, generalize the well-
known ordinary differential equations related to relaxation and oscillation
3We note that the Caputo fractional derivative was so named after the book by
Podlubny (1999). It coincides with that introduced, independently and a few later,
by Dzherbashyan and Nersesyan (1968) as a regularization of the Riemann-Liouville
fractional derivative. Nowadays, some Authors refer to it as the Caputo-Dzherbashyan
fractional derivative. The prominent role of this fractional derivative in treating initial
value problems was recognized in interesting papers by Kochubei (1989), (1990).
phenomena. Generally speaking, we consider the following differential
equation of fractional order α > 0 ,
Dα∗ u(t) = D
u(t) −
u(k)(0+)
= −u(t) + q(t) , t > 0 , (A.21)
where u = u(t) is the field variable and q(t) is a given function. Here m is
a positive integer uniquely defined by m − 1 < α ≤ m , which provides the
number of the prescribed initial values u(k)(0+) = ck , k = 0, 1, 2, . . . ,m−1 .
Implicit in the form of (A.21) is our desire to obtain solutions u(t) for which
the u(k)(t) are continuous. In particular, the cases of fractional relaxation
and fractional oscillation are obtained for 0 < α < 1 and 1 < α < 2 ,
respectively
The application of the Laplace transform through the Caputo formula (A.20)
yields
ũ(s) =
sα−k−1
sα + 1
sα + 1
q̃(s) . (A.22)
Now, in order to obtain the Laplace inversion of (A.22), we need to recall
the Mittag-Leffler function of order α > 0 , Eα(z) . This function, so named
from the great Swedish mathematician who introduced it at the beginning
of this century, is defined by the following series and integral representation,
valid in the whole complex plane,
Eα(z) =
Γ(αn + 1)
σα−1 e σ
σα − z
dσ , α > 0 . (A.23)
Here Ha denotes the Hankel path, i.e. a loop which starts and ends at −∞
and encircles the circular disk |σ| ≤ |z|1/α in the positive sense. It turns out
that Eα(z) is an entire function of order ρ = 1/α and type 1 .
The Mittag-Leffler function provides a simple generalization of the expo-
nential function, to which it reduces for α = 1 . Particular cases from which
elementary functions are recovered, are
= cosh z , E2
= cos z , z ∈ C , (A.24)
E1/2(±z1/2) = ez
1 + erf (±z1/2)
= ez erfc (∓z1/2) , z ∈ C , (A.25)
where erf (erfc) denotes the (complementary) error function. defined as
erf (z) :=
du , erfc (z) := 1 − erf (z) , z ∈ C .
A noteworthy property of the Mittag-Leffler function is based on the
following duplication formula
Eα(z) =
Eα/2(+z
1/2) + Eα/2(−z1/2)
. (A.26)
In (A.25-26) we agree to denote by z1/2 the main branch of the complex
root of z .
The Mittag-Leffler function is connected to the Laplace integral through the
equation ∫ ∞
e−u Eα (u
α z) du =
1 − z
α > 0 . (A.27)
The integral at the L.H.S. was evaluated by Mittag-Leffler who showed that
the region of its convergence contains the unit circle and is bounded by the
line Re z1/α = 1 . The above integral is fundamental in the evaluation of the
Laplace transform of Eα (−λ tα) with α > 0 and λ ∈ C . In fact, putting in
(A.27) u = st and uα z = −λ tα with t ≥ 0 and λ ∈ C , we get the Laplace
transform pair
Eα (−λ tα) ÷
sα + λ
, Re s > |λ|1/α . (A.28)
Then, using (A.28), we put for k = 0, 1, . . . ,m − 1 ,
uk(t) := J
keα(t) ÷
sα−k−1
sα + 1
, eα(t) := Eα(−tα) , (A.29)
and, from inversion of the Laplace transforms in (A.22), we find
u(t) =
ck uk(t) −
q(t − τ)u′0(τ) dτ . (A.30)
In particular, the formula (A.30) encompasses the solutions for α = 1 , 2 ,
since e1(t) = exp(−t) , e2(t) = cos t . When α is not integer, namely for
m − 1 < α < m , we note that m − 1 represents the integer part of α
(usually denoted by [α]) and m the number of initial conditions necessary
and sufficient to ensure the uniqueness of the solution u(t). Thus the m
functions uk(t) = J
keα(t) with k = 0, 1, . . . ,m−1 represent those particular
solutions of the homogeneous equation which satisfy the initial conditions
+) = δk h , h, k = 0, 1, . . . ,m − 1 , and therefore they represent the
fundamental solutions of the fractional equation (A.21), in analogy with the
case α = m . Furthermore, the function uδ(t) = −u′0(t) = −e′α(t) represents
the impulse-response solution.
The Mittag-Leffler function of order less than one turns out to be related
through the Laplace integral to another special function of Wright type,
denoted by M(z, ν) with 0 < ν < 1 , following the notation introduced
by Mainardi (1994, 1995). Since this function turns out to be relevant in
the general framework of fractional calculus with special regard to stable
probability distributions, we are going to summarize its basing properties.
For more details on this function, see Mainardi (1997), Appendix A.
Let us first recall the more general Wright function Wλ,µ(z) , z ∈ C , with
λ > −1 and µ > 0 . This function, so named from the British mathematician
who introduced it between 1933 and 1941, is defined by the following series
and integral representation, valid in the whole complex plane,
Wλ,µ(z) =
n! Γ(λn + µ)
eσ + zσ
−λ dσ
, (A.31)
where Ha denotes the Hankel path. It is possible to prove that the Wright
function is entire of order 1/(1+λ) , hence of exponential type if λ ≥ 0 . The
case λ = 0 is trivial since W0,µ(z) = e
z/Γ(µ) . The case λ = −ν , µ = 1 − ν
with 0 < ν < 1 provides the function M(z, ν) of special interest for us.
Specifically, we have
M(z; ν) := W−ν,1−ν(−z) =
W−ν,0(−z) , 0 < ν < 1 , (A.32)
and therefore from (A.31-32)
M(z; ν) =
(−z)n−1
(n − 1)!
Γ(ν n) sin (ν n π)
eσ − zσ
, 0 < ν < 1 .
(A.33)
In the series representation we have used the reflection formula for the
Gamma function, Γ(x) Γ(1−x) = π/ sin πx . Explicit expressions of M(z; ν)
in terms of simpler known functions are expected in particular cases when
ν is a rational number. Relevant cases are ν = 1/2 , 1/3 for which
M(z; 1/2) =
− z2/4
, (A.34)
M(z; 1/3) = 32/3 Ai
z/31/3
, (A.35)
where Ai denotes the Airy function.
When the argument is real and positive, i.e. z = r > 0 , the existence of
the Laplace transform of M(r; ν) is ensured by the asymptotic behaviour,
as derived by Mainardi & Tomirotti (1995), as r → +∞ ,
M(r/ν; ν) ∼ a(ν) r(ν − 1/2)/(1 − ν) exp
−b(ν) r1/(1 − ν)
, (A.36)
where a(ν) = 1/
2π (1 − ν) , b(ν) = (1 − ν)/ν .
It is an instructive exercise to derive the Laplace transform by interchanging
the Laplace integral with the Hankel integral in (A.33) and recalling the
integral representation (A.23) of the Mittag-Leffler function. We obtain the
Laplace transform pair
M(r; ν) ÷ Eν(−s) , 0 < ν < 1 . (A.37)
For ν = 1/2 , (A.37) with (A.25) and (A.34) provides the result, see e.g.
Doetsch (1974),
M(r; 1/2) :=
− r2/4
÷ E1/2(−s) := exp
erfc (s) . (A.38)
It would be noted that, since M(r, ν) is not of exponential order,
transforming term-by-term the Taylor series of M(r; ν) yields a series of
negative powers of s , which represents the asymptotic expansion of Eν(−s)
as s → ∞ in a certain sector around the real axis.
We also note that (A.37) with (A.23) allows us to compute the moments of
any real order δ ≥ 0 of M(r; ν) in the positive real axis. We obtain
r δ M(r; ν) dr =
Γ(δ + 1)
Γ(νδ + 1)
, δ ≥ 0 . (A.39)
When δ is integer we note that the moments are provided by the derivatives
of the Mittag-Leffler function in the origin, i.e.
rn M(r; ν) dr = lim
(−1)n
Eν(−s) =
Γ(n + 1)
Γ(νn + 1)
, (A.40)
where n = 0, 1, 2, . . . . The normalization condition
0 M(r; ν) dr =
Eν(0) = 1 is recovered for n = 0 . The relation with the Mittag-Leffler
function stated in (A.40) can be extended to the moments of non integer
order if we replace the ordinary derivative, of order n, with the corresponding
fractional derivative, of order δ 6= n, in the Caputo sense.
Another exercise on the function M concerns the inversion of the Laplace
transform exp(−sν) , either by the complex integral formula or by the formal
series method. We obtain the Laplace transform pair
M (1/rν ; ν) ÷ exp (−sν) , 0 < ν < 1 . (A.41)
For ν = 1/2 , (A.41) with (A.34) provides the known result, see e.g. Doetsch
(1974),
2 r3/2
M(1/r1/2; 1/2) :=
π r3/2
exp [− 1/(4r)] ÷ exp
−s1/2
. (A.42)
We recall that a rigorous proof of (A.41) was formerly given by Pollard
(1946), based on a formal result by Humbert (1945). The Laplace transform
pair was also obtained by Mikusiński (1959) and, albeit unaware of the
previous results, by Buchen & Mainardi (1975) in a formal way.
Appendix B: The Stable Probability Distributions
The stable distributions are a fascinating and fruitful area of research in
probability theory; furthermore, nowadays, they provide valuable models in
physics, astronomy, economics, and communication theory.
The general class of stable distributions was introduced and given this name
by the French mathematician Paul Lévy in the early 1920’s, see Lévy (1924,
1925). The inspiration for Lévy was the desire to generalize the celebrated
Central Limit Theorem, according to which any probability distribution
with finite variance belongs to the domain of attraction of the Gaussian
distribution.
Formerly, the topic attracted only moderate attention from the leading
experts, though there were also enthusiasts, of whom the Russian
mathematician Alexander Yakovlevich Khintchine should be mentioned first
of all. The concept of stable distributions took full shape in 1937 with the
appearance of Lévy’s monograph, see Lévy (1937-1954), soon followed by
Khintchine’s monograph, see Khintchine (1938).
The theory and properties of stable distributions are discussed in some
classical books on probability theory including Gnedenko & Kolmogorov
(1949-1954), Lukacs (1960-1970), Feller (1966-1971), Breiman (1968-1992),
Chung (1968-1974) and Laha & Rohatgi (1979). Also treatises on fractals
devote particular attention to stable distributions in view of their properties
of scale invariance, see e.g. Mandelbrot (1982) and Takayasu (1990). Sets of
tables and graphs have been provided by Mandelbrot & Zarnfaller (1959),
Fama & Roll (1968), Bo’lshev & Al. (1968) and Holt & Crow (1973).
Only recently, monographs devoted solely to stable distributions and related
stochastic processes have been appeared, i.e. Zolotarev (1983-1986), Janicki
& Weron (1994), Samorodnitsky & Taqqu (1994), Uchaikin & Zolotarev
(1999). We now can cite the paper by Mainardi, Luchko & Pagnini (2001)
where the reader can find (convergent and asymptotic) representations and
plots of the symmetric and non-symmetric stable densities generated by
fractional diffusion equations.
Stable distributions have three exclusive properties, which can be briefly
summarized stating that they 1) are invariant under addition, 2) possess
their own domain of attraction, and 3) admit a canonic characteristic
function.
Let us now illustrate the above properties which, providing necessary and
sufficient conditions, can be assumed as equivalent definitions for a stable
distribution. We recall the basic results without proof.
A random variable X is said to have a stable distribution P (x) = Prob {X ≤
x} if for any n ≥ 2 , there is a positive number cn and a real number dn such
X1 + X2 + . . . + Xn
= cn X + dn , (B.1)
where X1,X2, . . . Xn denote mutually independent random variables with
common distribution P (x) with X . Here the notation
= denotes equality
in distribution, i.e. means that the random variables on both sides have the
same probability distribution.
When mutually independent random variables have a common distribution
[shared with a given random variable X], we also refer to them as
independent, identically distributed (i.i.d) random variables [independent
copies of X]. In general, the sum of i.i.d. random variables becomes
a random variable with a distribution of different form. However, for
independent random variables with a common stable distribution, the sum
obeys to a distribution of the same type, which differs from the original
one only for a scaling (cn) and possibly for a shift (dn). When in (B.1) the
dn = 0 the distribution is called strictly stable.
It is known, see Feller (1966-1971), that the norming constants in (B.1) are
of the form
cn = n
1/α with 0 < α ≤ 2 . (B.2)
The parameter α is called the characteristic exponent or the index of stability
of the stable distribution.
We agree to use the notation X ∼ Pα(x) to denote that the random variable
X has a stable probability distribution with characteristic exponent α . We
simply refer to P (x) , p(x) := dP/dx (probability density function = pdf)
and X as α-stable distribution, density, random variable, respectively.
The definition (B.1) with the theorem (B.2) can be stated in an alternative
version that needs only two i.i.d. random variables. see also Lukacs (1960-
1970). A random variable X is said to have a stable distribution if for any
positive numbers A and B, there is a positive number C and a real number
D such that
AX1 + B X2
= C X + D , (B.3)
where X1 and X2 are independent copies of X . Then there is a number
α ∈ (0, 2] such that the number C in (B.3) satisfies Cα = Aα + Bα .
For a strictly stable distribution (B.3) holds with D = 0 . This implies that
all linear combinations of i.i.d. random variables obeying to a strictly stable
distribution is a random variable with the same type of distribution.
A stable distribution is called symmetric if the random variable −X has the
same distribution. Of course, a symmetric stable distribution is necessarily
strictly stable.
Noteworthy examples of stable distributions are provided by the Gaussian
(or normal) law (with α = 2) and by the Cauchy-Lorentz law (α = 1). The
corresponding pdf are known to be
pG(x;σ, µ) :=
e−(x − µ)
2/(2σ2) , x ∈ R , (B.4)
where σ2 denotes the variance and µ the mean, and
pC(x; γ, δ) :=
(x − δ)2 + γ2
, x ∈ R , (B.5)
where γ denotes the semi-interquartile range and δ the ”shift”.
Another (equivalent) definition states that stable distributions are the only
distributions that can be obtained as limits of normalized sums of i.i.d.
random variables. A random variable X is said to have a domain of
attraction,i.e. if there is a sequence of i.i.d. random variables Y1, Y2, . . .
and sequences of positive numbers {γn} and real numbers {δn}, such that
Y1 + Y2 + . . . Yn
d⇒X . (B.6)
The notation
d⇒ denotes convergence in distribution.
It is clear that the previous definition (B.1) yields (B.6), e.g. , by taking the
Yis to be independent and distributed like X . The converse is easy to show,
see Gnedenko & Kolmogorov (1949-1954). Therefore we can alternatively
state that a random variable X is said to have a stable distribution if it has
a domain of attraction.
When X is Gaussian and the Yis are i.i.d. with finite variance, then (B.6)
is the statement of the ordinary Central Limit Theorem. The domain
of attraction of X is said normal when γn = n
1/α ; in general, γn =
n1/α h(n) where h(x) , x > 0 , is a slow varying function at infinity, that
is, lim
h(ux)/h(x) = 1 for all u > 0 , see Feller (1971). The function
h(x) = log x , for example, is slowly varying at infinity.
Another definition specifies the canonic form that the characteristic function
(cf) of a stable distribution of index α must have. Recalling that the cf is
the Fourier transform of the pdf , we use the notation p̂α(κ) := 〈exp (iκX)〉 ÷
pα(x) . We first note that a stable distribution is also infinitely divisible, i.e.
for every positive integer n its cf can be expressed as the nth power of
some cf . In fact, using the characteristic function, the relation (B.1) is
transformed into
[p̂α(κ)]
n = p̂α(cn κ) e
idnκ . (B.7)
The functional equation (B.7) can be solved completely and the solution is
known to be
p̂α(κ;β, γ, δ) = exp {iδκ − γα |κ|α [1 + i (sign κ)β ω(|κ|, α)]} , (B.8)
where
ω(|κ|, α) =
tan (α π/2) , if α 6= 1 ,
−(2/π) log |κ| , if α = 1 . (B.9)
Consequently a random variable X is said to have a stable distribution if
there are four real parameters α, β, γ, δ with 0 < α ≤ 2 , −1 ≤ β ≤ +1 ,
γ > 0 , such that its characteristic function has the canonic form (B.8-9).
Then we write pα(x;β, γ, δ)÷ p̂α(κ;β, γ, δ) and X ∼ Pα(x;β, γ, δ) , so partly
following the notation of Holt & Crow (1973) and Samorodnitsky & Taqqu
(1994).
We note in (B.8-9) that β appears with different signs for α 6= 1 and α = 1 .
This minor point has been the source of great confusion in the literature, see
Hall (1980) for a discussion. The presence of the logarithm for α = 1 is the
source of many difficulties, so this case has often to be treated separately.
The cf (B.8-9) turns out to be a useful tool for studying α-stable distri-
butions and for providing an interpretation of the additional parameters,
β (skewness parameter), γ (scale parameter) and δ (shift parameter), see
Samorodnitsky & Taqqu (1994). When α = 2 the cf refers to the Gaussian
distribution with variance σ2 = 2 γ2 and mean µ = δ ; in this case the value
of the skewness parameter β is not specified because tan π = 0 , and one
conventionally takes β = 0 .
One easily recognizes that a stable distribution is symmetric if and only if
β = δ = 0 and is symmetric about δ if and only if β = 0 . Stable distributions
with extremal values of the skewness parameter are called extremal. One
can prove that all the extremal stable distributions with 0 < α < 1 are
one-sided, the support being R+0 if β = −1 , and R
0 if β = +1 .
For the stable distributions Pα(x;β, γ, δ) we now consider the asymptotic
behaviour of the tail probabilities, T+(λ) := Prob {X > λ} and T−(λ) :=
Prob {X < −λ} , as λ → ∞ . For the Gaussian case α = 2 the result is well
known, see e.g. Feller (1957),
α = 2 : T±(λ) ∼ 1
2/(4γ2)
, λ → ∞ . (B.10)
Because of the above exponential decay all the moments of the corresponding
pdf turn out to be finite, which is an exclusive property of this stable
distribution. For all the other stable distributions the singularity of the
characteristic function in the origin is responsible for the algebraic decay of
the tail probabilities as indicated below, see e.g. Samorodnitsky & Taqqu
(1994),
0 < α < 2 : lim
λα T±(λ) = Cα γ
α (1 ∓ β)/2 , (B.11)
where
x−α sin x dx
1 − α
Γ(2 − a) cos (απ/2)
, if α 6= 1 ,
2/π , if α = 1 .
(B.12)
We note that for extremal distributions (β = ±1) the above algebraic decay
holds true only for one tail, the left one if β = +1 , the right one if β = −1 .
The other tail is either identically zero if 0 < α < 1 (the distribution is
one-sided !), or exhibits an exponential decay if 1 ≤ α < 2 . Because of the
algebraic decay we recognize that
0 < α < 2 :
|x|>λ
pα(x;β, γ, δ) dx = O(λ
−α) , (B.13)
so the absolute moments of a stable non-Gaussian pdf turn out to be finite
if their order ν is 0 ≤ ν < α and infinite if ν ≥ α . We are now convinced
that the Gaussian distribution is the unique stable distribution with finite
variance. Furthermore, when α ≤ 1 , the first absolute moment 〈|X|〉 is
infinite as well, so we need to use the median to characterize the expected
value.
There is however a fundamental property shared by all the stable
distributions that we like to point out: for any α the stable pdf are unimodal
and indeed bell-shaped, i.e. their n-th derivative has exactly n zeros, see
Gawronski (1964).
We now come back to the cf of a stable distribution, in order to provide for
α 6= 1 and δ = 0 a simpler canonic form which allow us to derive convergent
and asymptotic power series for the corresponding pdf . We first note that
the two parameters γ and δ in (B.8), being related to a scale transformation
and a translation, are not so essential since they do not change the shape
of distributions. If we take γ = 1 and δ = 0 , we obtain the so-called
standardized form of the stable distribution and X ∼ Pα(x;β, 1, 0) is referred
to as the α-stable standardized random variable. Furthermore, we can choose
the scale parameter γ in such a way to get from (B.8-9) the simplified canonic
form used by Feller (1952, 1966-1971) and Takayasu (1990) for strictly stable
distributions (δ = 0) with α 6= 1 , which reads in an ad hoc notation,
q̂α(κ; θ) :=
eiκ y pα(y; θ) dy = exp
−|κ|α e±i θ π/2
, (B.14)
where the symbol ± takes the sign of κ . This canonic form, that we refer to
as the Feller canonic form, is derived from (B.8-9) if in addition to α 6= 1
and δ = 0 we require
γα = cos
, tan
= β tan
. (B.15)
Here θ is the skewness parameter instead of β and its domain is restricted
in the following region (depending on α)
|θ| ≤
α , if 0 < α < 1 ,
2 − α , if 1 < α < 2 . (B.16)
Thus, when we use the Feller canonic form for strictly stable distributions
with index α 6= 1 and skewness θ , we implicitly select the scale parameter
γ (0 < γ ≤ 1), which is related to α , β and θ by (B.15). Specifically, the
random variable Y ∼ Qα(y; θ) turns out to be related to the standardized
random variable X ∼ Pα(x;β, 1, 0) by the following relations
Y = X/γ , pα(x;β, 1, 0) = γ qα(y = γx; θ) , (B.17)
with 
γ = [cos (θπ/2)]1/α ,
θ = (2/π) arctan [β tan (απ/2)] ,
tan (θπ/2)
tan (απ/2)
(B.18)
We recognize that qα(y, θ) = qα(−y,−θ) , so the symmetric stable
distributions are obtained if and only if θ = 0 . We note that for the
symmetric stable distributions we get the identity between the standardized
and the Lévy canonic forms, since in (B.18) β = θ = 0 implies γ = 1 .
A particular but noteworthy case is provided by p2(x; 0, 1, 0) = q2(y; 0) ,
corresponding to the Gaussian distribution with variance σ2 = 2 .
The extremal stable distributions, corresponding to β = ±1 , are now
obtained for θ = ±α if 0 < α < 1 , and for θ = ∓(2 − α) if 1 < α < 2 ; for
them the scaling parameter turns out to be γ = [cos (|α|π/2)]1/α . It may be
an instructive exercise to carry out the inversion of the Fourier transform
when α = 1/2 and θ = −1/2 . In this case we obtain the analytical expression
for the corresponding extremal stable pdf , known as the (one-sided) Lévy-
Smirnov density,
q1/2(y;−1/2) =
y−3/2 e−1/(4y) , y ≥ 0 . (B.19)
The standardized form for this distribution can be easily obtained from
(B.19) using (B.17-18) with α = 1/2 and θ = −1/2 . We get γ =
[cos (−π/4)]2 = 1/2 , β = −1 , so
p1/2(x;−1, 1, 0) =
q1/2(x/2;−1/2) =
x−3/2 e−1/(2x) , (B.20)
where x ≥ 0 , in agreement with Holt & Crow (1973) [§2.13, p. 147].
Feller (1952) has obtained from (B.14) the following representations by
convergent power series for the stable distributions valid for y > 0 , with
0 < α < 1 (negative powers),
qα(y; θ) =
(−y−α)n Γ(nα + 1)
(θ − α)
, (B.21)
1 < α ≤ 2 (positive powers),
qα(y; θ) =
(−y)n
Γ(n/α + 1)
(θ − α)
. (B.22)
The values for y < 0 can be obtained from (B.21-22) using the identity
qα(−y; θ) = qα(y;−θ) , y > 0 . As a consequence of the convergence in all of
C of the series in (B.21-22) we recognize that the restrictions of the functions
y qα(y; θ) on the two real semi-axis turn out to be equal to certain entire
functions of argument 1/|y|α for 0 < α < 1 and argument |y| for 1 < α ≤ 2 .
It has be shown, see e.g. Bergström (1952), Chao Chung-Jeh (1953), that the
two series in (B.21-22) provide also the asymptotic (divergent) expansions to
the stable pdf with the ranges of α interchanged from those of convergence.
From (B.21-22) a relation between stable pdf with index α and 1/α can be
derived as noted in Feller (1966-1971). Assuming 1/2 < α < 1 and y > 0 ,
we obtain
q1/α(y
−α; θ) = qα(y; θ
∗) , θ∗ = α(θ + 1) − 1 . (B.23)
A quick check shows that θ∗ falls within the prescribed range, |θ∗| ≤ α ,
provided that |θ| ≤ 2 − 1/α .
We now consider two particular cases of the Feller series (B.21-22), of
particular interest for us, which turn out to be related to the entire function
of Wright type, M(z; ν) with 0 < ν < 1 , reported in Appendix A. These
cases correspond to the following extremal distributions
Φ1(y) := qα(y;−α) , y > 0 , 0 < α < 1 , (B.24)
Φ2(y) := qα(y;α − 2) , y > 0 , 1 < α ≤ 2 , (B.25)
for which the Feller series (B.21-22) reduce to
Φ1(y) =
(−1)n−1 y−αn−1 Γ(nα + 1)
sin (nπα) , y > 0 , (B.26)
Φ2(y) =
(−1)n−1 yn−1 Γ(n/α + 1)
, y > 0 . (B.27)
In fact, recalling the series representation of the general Wright function,
Wλ,µ(z) with λ > −1 , µ > 0 , see (A.31), and the definition of the function
M(z; ν) with 0 < ν < 1 , see (A.32-33), we recognize that
Φ1(y) =
W−α,0(−y−α) =
M(y−α;α) , y > 0 , (B.28)
Φ2(y) =
W−1/α,0(−y) =
M(y; 1/α) , y > 0 . (B.29)
We would like to remark that the above relations with the Wright functions
have been noted also by Engler (1997).
It is worth to point out that, whereas Φ1(y) totally represents the one-
sided stable pdf qα(y;−α) , 0 < α < 1 , with support in R+0 , Φ2(y) is the
restriction on the positive axis of qα(y;α− 2) , 1 < α ≤ 2 , whose support is
all of R . Since the function M(z; ν) turns out to be normalized in R+0 , see
(A.39-40), we also note
Φ1(y) dy = 1 ;
Φ2(y) dy = 1/α . (B.30)
Using the results (A.41) and (A.37) we can easily evaluate the Laplace
transforms of Φ1(y) and Φ2(y) , respectively. We obtain
L[Φ1(y)] = Φ̃1(s) = exp (−sα) , 0 < α < 1 , (B.31)
L[Φ2(y)] = Φ̃2(s) =
E1/α (−s) , 1 < α ≤ 2 , (B.32)
where E1/α(·) denotes the Mittag-Leffler function of order 1/α , see (A.23).
It is an instructive exercise to derive the asymptotic behaviours of Φ1(y) and
Φ2(y) as y → 0+ and y → +∞ . By using the expressions (B.28−29) in terms
of the function M and recalling the series and asymptotic representations of
this function, see (A.33) and (A.36), we obtain
Φ1(y) =
y−(2−α)/[2(1−α)] e−c1 y
−α/(1−α)
, as y → 0+ ,
Γ(1 − α)
y−α−1 [1 + O (y−α)] , as y → +∞ ,
(B.33)
Φ2(y) =
Γ(1 − 1/α)
[1 + O (y)] , as y → 0+ ,
y(2−α)/[2(α−1)] e−c2 y
α/(α−1)
, as y → +∞ ,
(B.34)
where c1 , c2 are positive constants depending on α . We note that the
exponential decay is found for Φ1(y) as y → 0+ but as y → +∞ for Φ2(y) .
Explicit expressions for stable pdf can be derived form those for the function
M(z; ν) when ν = 1/2 and ν = 1/3 , given in Appendix A, see (A.34-
35). Of course the ν = 1/2 expression can be used to recover the well-
known (symmetric) Gaussian distribution q2(y; 0) accounting for (B.29), and
the (one-sided) Lévy distribution q1/2(y;−1/2), see (B.19), accounting for
(B.28). The ν = 1/3 expression provides, accounting for (B.28),
q1/3(y;−1/3) = 3−1/3 y−4/3 Ai
(3y)−1/3
y−3/2 K1/3
(B.35)
where Ai denotes the Airy function and K1/3 the modified Bessel function of
the second kind of order 1/3 . The equivalence between the two expressions
in (B.35) can be proved in view of the relation, see Abramowitz & Stegun
(1965-1972) [(10.4.14)],
Ai (z) =
. (B.36)
The case α = 1/3 has also been discussed by Zolotarev (1983-1986), who
has quoted the corresponding expression of the pdf in terms of K1/3 .
A general representation of all stable distributions (thus including the
extremal distributions above considered) in terms of special functions has
been only recently achieved by Schneider (1986). In his remarkable
(but almost ignored) article, Schneider has established that all the stable
distributions can be characterized in terms of a general class of special
functions, the so-called Fox H functions, so named after Charles Fox (1961).
For details on Fox H functions, see e.g. the books Mathai & Saxena (1978),
Srivastava & Al. (1982) and the most recent paper by Kilbas and Saigo
(1999). These functions are expressed in terms of special integrals in the
complex-plane, the Mellin-Barnes integrals4.
4The names refer to the two authors, who in the first 1910’s developed the theory of
these integrals using them for a complete integration of the hypergeometric differential
equation. However, as pointed out in the the Bateman Project Handbook on High
Transcendental Functions, see Erdelyi (1953), these integrals were first used by S. Pincherle
in 1888. For a revisited analysis of the pioneering work of Pincherle (1853-1936, Professor
of Mathematics at the University of Bologna from 1880 to 1928) we refer to the paper by
Mainardi and Pagnini (2003).
References
[1] Abramowitz, M. and Stegun, I.A. (Editors) : Handbook of Mathematical
Functions, Dover, New York 1965. [reprint, 1972]
[2] Baillie, R.T. and King, M.L. (Editors) : Fractional Differencing and
Long Memory Processes, Journal of Econometrics, 73 1-324 (1996).
[3] Bergström, H. : On some expansions of stable distribution functions,
Ark. Mat., 2 375-378 (1952).
[4] Bol’shev, L.N., Zolotarev, V.M., Kedrova, E.S. and Rybinskaya, M.A.
: Tables of cumulative functions of on-sided stable distributions, Theor.
Probability Appl., 15 299-309 (1968).
[5] Bouchaud, J.-P and Potters, M. : Theory of Financial Risk: from Sta-
tistical Physics to Risk Management, Cambridge Univ. Press, Cambridge
2000. [English enlarged version of Théories des Risques Financiers, CEA,
Aléa, Saclay 1997.]
[6] Breiman, L. : Probability, SIAM, Philadelphia 1992. [The SIAM edition
is an unabridged, corrected republication of the 1st edn, Addison-Wesley,
Reading Mass. 1968]
[7] P.W. Buchen, and Mainardi, F. : Asymptotic expansions for transient
viscoelastic waves, J. de Mécanique, 14 597-608 (1975).
[8] Caputo, M. : Linear models of dissipation whose Q is almost frequency
independent, Part II., Geophys. J. R. Astr. Soc., 13 529-539 (1967).
[9] Caputo, M.: Elasticità e Dissipazione, Zanichelli, Bologna 1969.
[10] Caputo, M. and Mainardi, F. : Linear models of dissipation in anelastic
solids, Rivista del Nuovo Cimento (Ser. II), 1 161-198 (1971).
[11] Chung-Jeh, Chao : Explicit formula for the stable law of distribution,
Acta Math. Sinica, 3 177-185 (1953). [in Chinese with English summary]
[12] Chung, K.L. : A Course in Probability Theory, 2nd edn, Academic
Press, New York 1974. [1st edn, Harcourt Brace Jowvanovich, 1968]
[13] Dzherbashyan, M.M. and Nersesyan, A.B. : Fractional derivatives and
the Cauchy problem for differential equations of fractional order. Izv.
Acad. Nauk Armjanskvy SSR, Matematika, 3 3–29 (1968). [In Russian]
[14] Engler, H. : Similarity solutions for a class of hyperbolic integrodiffer-
ential equations, Differential Integral Eqns, 10 815-840 (1997).
[15] Erdélyi, A. (Editor) : Higher Transcendental Functions, Bateman
Project, McGraw-Hill, New York 1953; Vol. 1, Ch. 1, §1.19, p. 49.
[16] Erdélyi, A. (Editor) : Higher Transcendental Functions, Bateman
Project, McGraw-Hill, New York 1955; Vol. 3, Ch. 18, pp. 206-227.
[17] Erdélyi, A. (Editor) : Tables of Integral Transforms, Bateman Project
McGraw-Hill, New York 1954; Vol. 2, Ch. 13, pp. 181-212.
[18] Fama, E. and Roll, R. : Some properties of symmetric stable
distributions, J. Amer. Statist. Assoc., 63 817-836 (1968).
[19] Feller, W. : On a generalization of Marcel Riesz’ potentials and
the semi-groups generated by them, Meddelanden Lunds Universitets
Matematiska Seminarium (Comm. Sém. Mathém. Université de Lund),
Tome suppl. dédié a M. Riesz, Lund (1952), pp. 73-81.
[20] Feller, W. : An Introduction to Probability Theory and its Applications,
Vol. 1, 3rd edn, Wiley, New York 1968. [1st edn, 1957]
[21] Feller, W. : An Introduction to Probability Theory and its Applications,
Vol. 2, 2nd edn, Wiley, New York 1971. [1st edn. 1966]
[22] Fox, C. : The G and H functions as symmetrical Fourier kernels, Trans.
Amer. Math. Soc., 98 395-429 (1961).
[23] Fujita, Y. : Integrodifferential equation which interpolates the heat
equation and the wave equation, Osaka J. Math., 27 309-321, 797-804
(1990). [2 papers]
[24] Gawronski, W. : On the bell-shape of stable distributions, Annals of
Probability, 12 230-242 (1984).
[25] Gnedenko, B.V. and Kolmogorov, A.N. : Limit Distributions for Sums
of Independent Random Variables, Addison-Wesley, Cambridge, Mass.
1954. [English Transl. from the Russian edition 1949, with notes by K.L.
Chung, revised 1968]
[26] Gorenflo, R. : Fractional calculus: some numerical methods, in
Carpinteri, A. and Mainardi, F. (Eds.), Fractals and Fractional Calculus
in Continuum Mechanics, CISM Courses and Lectures # 378, Springer
Verlag, Wien 1997, pp. 277-290. [Reprinted in www.fracalmo.org]
[27] Gorenflo, R. and Mainardi, F. : Fractional calculus: integral and
differential equations of fractional order, in Carpinteri, A. and Mainardi,
F. (Eds.), Fractals and Fractional Calculus in Continuum Mechanics,
CISM Courses and Lectures # 378, Springer Verlag, Wien 1997, pp. 223-
276. [Reprinted in www.fracalmo.org]
[28] Gorenflo, R. and Mainardi, F. : Fractional calculus and stable
probability distributions, Archives of Mechanics, 50 377-388 (1998a).
[29] Gorenflo, R. and Mainardi, F. : Random walk models for space-
fractional diffusion processes, Fractional Calculus and Applied Analysis,
1 167-190 (1998b).
[30] Gorenflo, R. De Fabritiis, G. and Mainardi, F. : Discrete random walk
models for symmetric Lévy-Feller diffusion processes, Physica A, 269 79–
89 (1999).
[31] Gorenflo, R., Mainardi, F., Moretti, D., Pagnini, G. and Paradisi,
P. : Discrete random walk models for space-time fractional diffusion,
Chemical Physics, 284 521-541 (2002). Special issue on Strange Kinetics
Guest Editors: R. Hilfer, R. Metzler, A. Blumen, J. Klafter. [E-print
arXiv:cond-mat/0702072]
[32] Hall, P. : A comedy of errors: the canonical form for a stable
characteristic function, Bull. London Math. Soc., 13 23-27 (1980).
[33] Holt, D.R. and Crow, E.L. : Tables and graphs of the stable probability
density functions, J. Res. Nat. Bureau Standards, 77B 143-198 (1973).
[34] Humbert, P. : Nouvelles correspondances symboliques, Bull. Sci.
Mathém. (Paris, II ser.), 69 121-129 (1945).
[35] Janicki, A. and Weron, A. : Simulation and Chaotic Behavior of α-
Stable Stochastic Processes, Marcel Dekker, New York 1994.
[36] Khintchine, A.Y. : Limit Laws for Sums of Independent Variables,
ONTI, Moscow 1938 [in Russian]
[37] Kilbas, A.A. and Saigo, M. : On the H functions, Journal of Applied
Mathematics and Stochastic Analysis, 12 191-204 (1999).
[38] Kochubei, A.N. : A Cauchy problem for evolution equations of
fractional order, Differential Equations, 25 967–974 (1989). [English
translation from the Russian Journal Differentsial’nye Uravneniya]
[39] Kochubei, A.N. : Fractional order diffusion, Differential Equations,
26 485–492 (1990). [English translation from the Russian Journal
Differentsial’nye Uravneniya]
[40] Laha, R.G. and Rohatgi, V.K. : Probability Theory, Wiley, New York
1979.
[41] Lévy, P. : Théorie des erreurs. La Loi de Gauss et les lois exceptionelles,
Bull. Soc. Math. France, 52 49-85 (1924).
[42] Lévy, P. : Calcul des probabilités, Gauthier-Villars, Paris 1925: Part II,
Chap. 6.
[43] Lévy, P. : Théorie de l’addition des variables aléatoires, 2nd edn,
Gauthier-Villars, Paris 1954. [1st edn, 1937]
[44] Lukacs, E. : Characteristic Functions, 2nd edn, Griffin, London 1970.
[1st edn, 1960]
[45] Mainardi, F. : On the initial value problem for the fractional diffusion-
wave equation, in S. Rionero and T. Ruggeri (Eds), Waves and Stability
in Continuous Media, World Scientific, Singapore 1994, pp. 246-251.
[46] Mainardi, F. : The time fractional diffusion-wave equation, Radiofisika,
38 20-36 (1995). [English Translation: Radiophysics & Quantum
Electronics]
[47] Mainardi, F. : Fractional relaxation-oscillation and fractional diffusion-
wave phenomena, Chaos, Solitons & Fractals, 7 1461-1477 (1996).
[48] Mainardi, F. : Fractional calculus: some basic problems in continuum
and statistical mechanics, in Carpinteri, A. and Mainardi, F. (Eds),
Fractals and Fractional Calculus in Continuum Mechanics, CISM Courses
and Lectures # 378, Springer-Verlag, Wien 1997, pp. 291-348. [Reprinted
in www.fracalmo.org]
[49] Mainardi, F., Luchko, Yu. and Pagnini, G. : The fundamental solution
of the space-time fractional diffusion equation, Fractional Calculus and
Applied Analysis, 4 153-192 (2001). [E-print arXiv:cond-mat/0702419]
[50] Mainardi, F. and Pagnini, G. : Salvatore Pincherle, the pioneer of the
Mellin-Barnes integrals, Jour. Computational and Applied Mathematics,
153 331-342 (2003).
[51] Mainardi, F. and Tomirotti, M. : On a special function arising in the
time fractional diffusion-wave equation, in P. Rusev, I. Dimovski and V.
Kiryakova (Eds), Transform Methods and Special Functions, Sofia 1994,
Science Culture Technology, Singapore 1995, pp. 171-183.
[52] Mainardi, F. and Tomirotti, M. : Seismic pulse propagation with
constant Q and stable probability distributions, Annali di Geofisica, 40
1311-1328 (1997).
[53] Mandelbrot, B.B. and Zarnfaller, F. : Five place tables of certain stable
distributions, Technical Report RC-421, IBM Thomas J. Watson Research
Center, Yorktown Heights, New York, Dec 31, 1959.
[54] Mandelbrot, B.B. : The Fractal Geometry of Nature, Freeman, San
Francisco 1982.
[55] Mandelbrot, B.B. : Fractals and Scaling in Finance, Springer-Verlag,
New York 1997.
[56] Mantegna, R.N. and Stanley, H.E. : An Introduction to Econophysics,
Cambridge University Press, Cambridge 2000.
[57] Mathai, A.M. and Saxena, R.K. : The H-function with Applications in
Statistics and Other Disciplines, Wiley Eastern Ltd, New Delhi 1978.
[58] Mikusiński, J. : On the function whose Laplace transform is
exp (−sαλ) , Studia Math., 18 191-198 (1959).
[59] Miller, K.S. and Ross, B. : An Introduction to the Fractional Calculus
and Fractional Differential Equations, Wiley, New York 1993.
[60] Oldham, K.B. and Spanier, J. : The Fractional Calculus, Academic
Press, New York 1974.
[61] Pincherle, S. : Sulle funzioni ipergeometriche generalizzate, Atti R.
Accademia Lincei, Rend. Cl. Sci. Fis. Mat. Nat. (4), 4, 694-700, 792-799
(1888).
[62] Podlubny, I. : Fractional Differential Equations, Academic Press, San
Diego 1999. [Mathematics in Science and Engineering, Vol. 198]
[63] Pollard, H. : The representation of exp (−xλ) as a Laplace integral,
Bull. Amer. Math. Soc., 52, 908-910 (1946).
[64] Prüss, J. : Evolutionary Integral Equations and Applications,
Birkhäuser, Basel 1993.
[65] Rubin, B. : Fractional Integrals and Potentials, Pitman Monographs
and Surveys in Pure and Appl. Mathematics # 82, Longman, London
1996.
[66] Saichev, A.I and Zaslavsky, G.M. : Fractional kinetic equations:
solutions and applications, Chaos, 7 753-764 (1997).
[67] Samko, S.G., Kilbas, A.A. and Marichev, O.I. : Fractional Integrals and
Derivatives, Theory and Applications, Gordon and Breach, Amsterdam
1993. [Engl. Transl. from the Russian edition, 1987]
[68] Samorodnitsky, G. and Taqqu, M.S. : Stable non-Gaussian Random
Processes, Chapman & Hall, New York 1994.
[69] Schneider, W.R. : Stable distributions: Fox function representation and
generalization, in S. Albeverio, G. Casati and D. Merlini (Eds), Stochastic
Processes in Classical and Quantum Systems, Lecture Notes in Physics #
262, Springer Verlag, Berlin 1986, 497-511.
[70] Schneider, W.R. and Wyss, W. : Fractional diffusion and wave
equations, J. Math. Phys., 30 134-144 (1989).
[71] Srivastava, H.M., Gupta, K.C. and Goyal, S.P. : The H-Functions of
One and Two Variables with Applications, South Asian Publ., New Delhi
1982.
[72] H. Takayasu, H. : Fractals in the Physical Sciences, Manchester Univ.
Press, Manchester and New York 1990.
[73] Uchaikin, V.V. and Zolotarev, V.M. : Chance and Stability. Stable
Distributions and their Applications, VSP, Utrecht 1999. [Series ”Modern
Probability and Statistics”, No 3]
[74] Zolotarev, V.M. : One-dimensional stable distributions, Amer. Math.
Soc., Providence, R.I. 1986. [English Transl. from the Russian edition,
1982]
	Introduction
	The standard diffusion equation
	The time-fractional diffusion equation
	The Cauchy problem for the time-fractional diffusion equation
	The Signalling problem for the time-fractional diffusion equation
	The Cauchy problem for the symmetric space-fractional diffusion equation
	Conclusions

  Fractional calculus allows one to generalize the linear, one-dimensional,
diffusion equation by replacing either the first time derivative or the second
space derivative by a derivative of fractional order. The fundamental solutions
of these equations provide probability density functions, evolving on time or
variable in space, which are related to the class of stable distributions. This
property is a noteworthy generalization of what happens for the standard
diffusion equation and can be relevant in treating financial and economical
problems where the stable probability distributions play a key role.

Introduction . . . . . . . . . . . . . . . . . . . . . . . p. 2
2. The Standard Diffusion Equation . . . . . . . . . . . . . . p. 4
3. The Time-Fractional Diffusion Equation . . . . . . . . . . . p. 8
4. The Cauchy Problem for the Time-Fractional Diffusion Equation p.10
5. The Signalling Problem for the Time-Fractional Diffusion Equation p.13
6. The Cauchy Problem for the Symmetric Space-Fractional
Diffusion Equation . . . . . . . . . . . . . . . . . . . . p.15
7. Conclusions . . . . . . . . . . . . . . . . . . . . . . . p.21
A. The Riemann-Liouville Fractional Calculus . . . . . . . . . p.22
B. The Stable Probability Distributions . . . . . . . . . . . . p.31
References . . . . . . . . . . . . . . . . . . . . . . . p.41
1This paper is based on an invited talk given by Francesco Mainardi at the International
Workshop on Econophysics held at Bolyai College, Eötvös University, Budapest, on
July 21-27, 1997. The paper was originally edited as a contribution for the book J.
Kertesz and I. Kondor (Editors), Econophysics: an Emerging Science, Kluwer
Academic Publishers, Dordrecht (NL) that should contain selected papers presented at
that Workshop and should have appeared in 1998 or 1999. Unfortunately the book was
not published. The present e-print is a revised version (with up-date annotations and
references) of that unpublished contribution, but essentially represents our knowledge of
that early time.
http://arXiv.org/abs/0704.0320v1
Abstract
Fractional calculus allows one to generalize the linear, one-dimensional,
diffusion equation by replacing either the first time derivative or the second
space derivative by a derivative of fractional order. The fundamental
solutions of these generalized diffusion equations are shown to provide
probability density functions, evolving on time or variable in space, which
are related to the peculiar class of stable distributions. This property is
a noteworthy generalization of what happens for the standard diffusion
equation and can be relevant in treating financial and economical problems
where the stable probability distributions are known to play a key role.
1 Introduction
Non-Gaussian probability distributions are becoming more common as data
models, especially in economics where large fluctuations are expected. In
fact, probability distributions with heavy tails are often met in economics
and finance, which suggests to enlarge the arsenal of possible stochastic
models by non-Gaussian processes. This conviction started in the early
sixties after the appearance of a series of papers by Mandelbrot and
his associates, who point out the importance of non-Gaussian probability
distributions, formerly introduced by Pareto and Lévy, and related scaling
properties, to analyse economical and financial variables, as reported in
the recent book by Mandelbrot (1997). Some examples of such variables
are common stock prices changes, changes in other speculative prices, and
interest rate changes. In this respect many works by different authors have
recently appeared, see e.g. the recent books by Bouchaud & Potter (1997),
Mantegna & Stanley (1998) and the references therein quoted.
It is well known that the fundamental solution (or Green function) of
the Cauchy problem for the standard linear diffusion equation provides at
any time the probability density function (pdf) in space of the Gauss (or
normal) law. This law exhibits all moments finite thanks to its exponential
decay at infinity. In particular, the space variance of the Green function
is proportional to the first power of time, a noteworthy property that
can be understood by means of an unbiased random walk model for the
Brownian motion, see e.g. Feller (1957). Less known is the property for
which the fundamental solution of the Signalling problem for the same
diffusion equation, provides at any position a unilateral pdf in time, known
as Lévy law, using the terminology of Feller (1966-1973). Because of its
algebraic decay at infinity as t−3/2 , this law has all moments of integer
order divergent, and consequently its expectation value and variance are
infinite.
Both the Gauss and Lévy laws belong to the general class of stable probability
distributions, which are characterized by an index α (0 < α ≤ 2), called
index of stability or characteristic exponent. In particular, the index of the
Gauss law is 2 , whereas that of the Lévy law is 1/2 .
In this paper we consider two different generalizations of the diffusion
equation by means of fractional calculus, which allows us to replace either the
first time derivative or the second space derivative by a suitable fractional
derivative. Correspondingly, the generalized equation will be referred to
as the time-fractional diffusion equation or the symmetric, space-fractional
diffusion equation. Here we show how the fundamental solutions of this
equation for the Cauchy and Signalling problems provide probability density
functions related to certain stable distributions, so providing a natural
generalization of what occurs for the standard diffusion equation.
The plan of the paper is as follows. First of all, for the sake of convenience
and completeness, we provide the essential notions of Riemann-Liouville
Fractional Calculus and Lévy Stable Probability Distributions in Appendix
A and B, respectively.
In Section 2, we recall the basic results for the standard diffusion
equation concerning the fundamental solutions of the Cauchy and Signalling
problems. In particular we provide the derivation of these solutions by the
Fourier and Laplace transforms and the interpretation in terms of Gauss
and Lévy stable pdf , respectively.
In Section 3, we consider the time-fractional diffusion equation and we
formulate for it the basic Cauchy and Signalling problems to be treated in the
subsequent two sections. Here we adopt the Riemann-Liouville approach to
Fractional Calculus, and the related definition for the Caputo time-fractional
derivative of a causal function of time.
In Section 4, we solve the Cauchy problem for the time-fractional diffusion
equation by using the technique of Fourier transform and we derive the
corresponding fundamental solution in terms of a special function of Wright
type in the similarity variable. In this case the solution can be interpreted
as a noteworthy symmetric pdf in space with all moments finite, evolving
in time. In particular, its space variance turns out to be proportional to a
power of time equal to the order of the time-fractional derivative.
In Section 5, we derive the fundamental solution for the Signalling problem
of the time-fractional diffusion equation by using the technique of Laplace
transform. In this case the solution, still expressed in terms of a special
function of Wright type, can be interpreted as a unilateral stable pdf in
time, depending on position, with index of stability given by half of the
order of the time-fractional derivative.
In Section 6, we consider the symmetric, space-fractional diffusion equation.
Here we adopt the Riesz approach to Fractional Calculus, and the related
definition for the symmetric space-fractional derivative of a function of a
single space variable. Here we treat the Cauchy problem by technique
of Fourier transform and we derive the series representation of the
corresponding Green function. In this case the fundamental solution is
interpreted in terms of a symmetric stable pdf in space, evolving in time,
with index of stability given by the order of the space-fractional derivative.
To approximate such evolution we propose a random walk model, discrete
in space and time, which is based on the Grünwald-Letnikov approximation
of the fractional derivative.
Finally, Section 7 is devoted to conclusions and remarks on related work.
2 The standard diffusion equation
For the standard diffusion equation we mean the linear partial differential
equation
u(x, t) = D
u(x, t) , u = u(x, t) , (2.1)
where D denotes a positive constant with the dimensions L2 T−1 , x and t
are the space-time variables, and u = u(x, t) is the field variable, which is
assumed to be a causal function of time, i.e. vanishing for t < 0 .
The typical physical phenomenon related to such an equation is the heat
conduction in a thin solid rod extended along x , so the field variable u is
the temperature.
In order to guarantee the existence and the uniqueness of the solution,
we must equip (1.1) with suitable data on the boundary of the space-time
domain. The basic boundary-value problems for diffusion are the so-called
Cauchy and Signalling problems. In the Cauchy problem, which concerns
the space-time domain −∞ < x < +∞ , t ≥ 0 , the data are assigned at
t = 0+ on the whole space axis (initial data). In the Signalling problem,
which concerns the space-time domain x ≥ 0 , t ≥ 0 , the data are assigned
both at t = 0+ on the semi-infinite space axis x > 0 (initial data) and at
x = 0+ on the semi-infinite time axis t > 0 (boundary data); here, as mostly
usual, the initial data are assumed to be vanishing.
Denoting by g(x) and h(t) two given, sufficiently well-behaved functions, the
basic problems are thus formulated as following:
a) Cauchy problem
u(x, 0+) = g(x) , −∞ < x < +∞ ; u(∓∞, t) = 0 , t > 0 ; (2.2a)
b) Signalling problem
u(x, 0+) = 0 , x > 0 ; u(0+, t) = h(t) , u(+∞, t) = 0 , t > 0 . (2.2b)
Hereafter, for both the problems, we derive the classical results which will be
properly generalized for the fractional diffusion equation in the subsequent
sections.
Let us begin with the Cauchy problem. It is well known that this initial value
problem can be easily solved making use of the Fourier transform and its
fundamental solution can be interpreted as a Gaussian pdf in x. Adopting
the notation g(x) ÷ ĝ(κ) with κ ∈ R and
ĝ(κ) = F [g(x)] =
e+iκx g(x) dx ,
g(x) = F−1 [ĝ(κ)] = 1
e−iκx ĝ(κ) dκ ,
the transformed solution satisfies the ordinary differential equation of the
first order (
+ κ2 D
û(κ, t) = 0 , û(κ, 0+) = ĝ(κ) , (2.3)
and consequently it turns out to be
û(κ, t) = ĝ(κ) e−κ
2 D t . (2.4)
Then, introducing
Gdc (x, t) ÷ Ĝdc (κ, t) = e−κ
2 D t , (2.5)
where the upper index d refers to (standard) diffusion, the required solution,
obtained by inversion of (2.4), can be expressed in terms of the space
convolution u(x, t) =
−∞ Gdc (ξ, t) g(x − ξ) dξ , where
Gdc (x, t) =
t−1/2 e−x
2/(4D t) . (2.6)
Here Gdc (x, t) represents the fundamental solution (or Green function) of
the Cauchy problem, since it corresponds to g(x) = δ(x) . It turns out
to be a function in x , even and normalized, i.e. Gdc (x, t) = Gdc (|x|, t) and∫ +∞
−∞ Gdc (x, t) dx = 1 . We also note the identity
|x| Gdc (|x|, t) =
Md(ζ) , (2.7)
where ζ = |x|/(
D t1/2) is the well-known similarity variable and
Md(ζ) =
2/4 . (2.8)
We note that Md(ζ) satisfies the normalization condition
d(ζ) dζ = 1 .
The interpretation of the Green function (2.6) in probability theory is
straightforward since we easily recognize
Gdc (x, t) = pG(x;σ) :=
2/(2σ2) , σ2 = 2D t , (2.9)
where pG(x;σ) denotes the well-known Gauss or normal pdf spread out
over all real x (the space variable), whose moment of the second order, the
variance, is σ2 . The associated cumulative distribution function (cdf) is
known to be
PG(x;σ) :=
′;σ) dx′ =
1 + erf
, (2.10)
where erf (z) := (2/
0 exp (−u2) du denotes the error function.
Furthermore, the moments of even order of the Gauss pdf turn out to be∫ +∞
2n pG(x;σ) dx = (2n − 1)!!σ2n , so
x2n Gdc (x, t) dx = (2n − 1)!! (2D t)n , n = 1, 2, . . . . (2.11)
Let us now consider the Signalling problem. This initial-boundary value
problem can be easily solved by making use of the Laplace transform.
Adopting the notation h(t) ÷ h̃(s) with s ∈ C and
h̃(s) = L [h(t)] =
e−st h(t) dt ,
h(t) = L−1
h̃(t)
est h̃(s) ds ,
where Br denotes the Bromwich path, the transformed solution of the
diffusion equation satisfies the ordinary differential equation of the second
order
ũ(x, s) = 0 , ũ(0+, s) = h̃(s) , ũ(+∞, s) = 0 . (2.12)
and consequently it turns out to be
ũ(x, s) = h̃(s) e−(x/
D) s1/2 . (2.13)
Then introducing
Gds (x, t) ÷ G̃ds (x, s) = e−(x/
D) s1/2 , (2.14)
the required solution, obtained by inversion of (2.13), can be expressed in
terms of the time convolution, u(x, t) =
0 Gds (x, τ)h(t − τ) dτ , where
Gds (x, t) =
t−3/2 e−x
2/(4D t) . (2.15)
Here Gds (x, t) represents the fundamental solution (or Green function) of the
Signalling problem, since it corresponds to h(t) = δ(t) . We note that
Gds (x, t) = pLS(t;µ) :=
2π t3/2
e−µ/(2t) , t ≥ 0 , µ = x
, (2.16)
where pLS(t;µ) denotes the one-sided Lévy-Smirnov pdf spread out over all
non negative t (the time variable). The associated cdf is, see e.g. Feller
(1966-1971) and Prüss (1993),
PL(t;µ) :=
′;µ) dt′ = erfc
= erfc
, (2.17)
where erfc (z) := 1 − erf (z) denotes the complenatary error function.
The Lévy-Smirnov pdf has all moments of integer order infinite, since it
decays at infinity as t−3/2 . However, we note that the absolute moments of
real order ν are finite only if 0 ≤ ν < 1/2 . In particular, for this pdf the mean
is infinite, for which we can take the median as expectation value. From
PLs(tmed;µ) = 1/2 , it turns out that tmed ≈ 2µ , since the complementary
error function gets the value 1/2 as its argument is approximatively 1/2.
We note that in the common domain x > 0 , t > 0 the Green functions of
the two basic problems satisfy the identity
xGdc (x, t) = tGds (x, t) , (2.18)
that we refer to as the reciprocity relation between the two fundamental
solutions of the diffusion equation. Furthermore, in view of (2.7) and (2.18)
we recognize the role of the function of the similarity variable, Md(ζ) ,
in providing the two fundamental solutions; we shall refer to it as to the
normalized auxiliary function of the diffusion equation for both the Cauchy
and Signalling problems.
3 The time-fractional diffusion equation
By the time-fractional diffusion equation we mean the linear evolution
equation obtained from the classical diffusion equation by replacing the first-
order time derivative by a fractional derivative (in the Caputo sense) of order
α with 0 < α ≤ 2. In our notation it reads
, u = u(x, t) , 0 < α ≤ 2 , (3.1)
where D denotes a positive constant with the dimensions L2 T−α . From
Appendix A we recall the definition of the Caputo fractional derivative of
order α > 0 for a (sufficiently well-behaved) causal function f(t) , see (A.9),
Dα∗ f(t) :=
Γ(m − α)
(t − τ)m−α f (m)(τ) dτ , (3.2)
where m = 1, 2, . . . , and 0 ≤ m − 1 < α ≤ m . According to (3.2) we thus
need to distinguish the cases 0 < α ≤ 1 and 1 < α ≤ 2 . In the the latter case
(3.1) may be seen as a sort of interpolation between the standard diffusion
equation and the standard wave equation. Introducing
Φλ(t) :=
tλ−1+
, λ > 0 , (3.3)
where the suffix + is just denoting that the function is vanishing for t < 0 ,
we easily recognize that the equation (3.1) assumes the explicit forms :
if 0 < α ≤ 1 ,
Φ1−α(t) ∗
Γ(1 − α)
(t − τ)−α
dτ = D
; (3.4)
if 1 < α ≤ 2 ,
Φ2−α(t) ∗
Γ(2 − α)
(t − τ)1−α
dτ = D ∂
. (3.5)
Extending the classical analysis for the standard diffusion equation (2.1) to
the above integro-differential equations (3.4-5), the Cauchy and Signalling
problems are thus formulated as in equations (2.2), i.e.
a) Cauchy problem
u(x, 0+) = g(x) , −∞ < x < +∞ ; u(∓∞, t) = 0 , t > 0 ; (3.6a)
b) Signalling problem
u(x, 0+) = 0 , x > 0 ; u(0+, t) = h(t) , u(+∞, t) = 0 , t > 0 . (3.6b)
However, if 1 < α ≤ 2 , the presence in (3.5) of the second order time
derivative of the field variable requires to specify the initial value of the first
order time derivative ut(x, 0
+) , since in this case two linearly independent
solutions are to be determined. To ensure the continuous dependence of our
solution on the parameter α also in the transition from α = 1− to α = 1+ ,
we agree to assume ut(x, 0
+) = 0 .
We recognize that our fractional diffusion equation (3.1), when subject to
the conditions (3.6), is equivalent to the integro-differential equation
u(x, t) = g(x) +
(t − τ)α−1
dτ , (3.7)
where 0 < α ≤ 2 . Such integro-differential equation has been investigated
by several authors, including Schneider & Wyss (1989), Fujita (1990), Prüss
(1993) and Engler (1997).
In view of our subsequent analysis we find it convenient to put
, 0 < ν < 1 . (3.8)
In fact the analysis of the time-fractional diffusion equation turns out to
be easier if we adopt as a key parameter the half of the order of the
time-fractional derivative. In future we shall provide the symbol α with
other relevant meanings, as the index of stability of a stable probability
distribution or the order of the space derivative in the space-fractional
diffusion equation.
Henceforth, we agree to insert the parameter ν in the field variable, i.e.
u = u(x, t; ν) . By denoting the Green functions of the Cauchy and Signalling
problems by Gc(x, t; ν) and Gs(x, t; ν) , respectively, the solutions of the two
basic problems are obtained by a space or time convolution, u(x, t; ν) =∫ +∞
−∞ Gc(ξ, t; ν) g(x−ξ) dξ , u(x, t; ν) =
0 Gs(x, τ ; ν)h(t−τ) dτ , respectively.
It should be noted that Gc(x, t; ν) = Gc(|x|, t; ν) , since the Green function
turns out to be an even function of x .
In the following two sections we shall compute the two fundamental solutions
with the same techniques (based on Fourier and Laplace transforms) used
for the standard diffusion equation and we shall provide their interpretation
in terms of probability distributions. Most of the presented results are based
on the papers by Mainardi (1994), (1995), (1996), (1997) and by Mainardi
& Tomirotti (1995), (1997).
4 The Cauchy problem for the time-fractional
diffusion equation
For the fractional diffusion equation (3.1) subject to (3.6a) the application
of the Fourier transform leads to the ordinary differential equation of order
α = 2ν ,
+ κ2 D
û(κ, t; ν) = 0 , û(κ, 0+; ν) = ĝ(κ) , (4.1)
Using the results of Appendix A, see (A.22-30), the transformed solution is
û(κ, t; ν) = ĝ(κ)E2ν
−κ2 D t2ν
, (4.2)
where E2ν(·) denotes the Mittag-Leffler function of order 2ν , and conse-
quently for the Green function we have
Gc(x, t; ν) = Gc(|x|, t; ν) ÷ Ĝc(k, t; ν) = E2ν
−κ2D t2ν
. (4.3)
Since the Green function is a real and even function of x, its (exponential)
Fourier transform can be expressed in terms of the cosine Fourier transform
and thus is related to its spatial Laplace transform as follows
Ĝc(k, t; ν) = 2
Gc(x, t; ν) cos κx dx =
G̃c(s, t; ν)
s=+ik
+ G̃c(s, t; ν)
s=−ik
(4.4)
Indeed, a split occurs also in (4.3) according to the duplication formula for
the Mittag-Leffler function, see (A.26),
Ĝc(k, t; ν) = E2ν(−κ2 D t2ν) =
[Eν(+iκ
D tν) + Eν(−iκ
D tν)]/2 .
(4.5)
When ν 6= 1/2 the inversion of the Fourier transform in (4.5) cannot be
obtained by using a standard table of Fourier transform pairs; however, for
any ν ∈ (0, 1) such inversion can be achieved by appealing to the Laplace
transform pair (A.37) with r = |x| , and s = ±iκ . In fact, taking into
account the scaling property of the Laplace transform, we obtain from (4.5)
and (A.37)
Gc(|x|, t; ν) =
( |x|√
, (4.6)
where M(ζ; ν) is the special function of Wright type, defined by (A.31-33),
, (4.7)
the similarity variable. We note the identity
|x| Gc(|x|, t; ν) =
M(ζ; ν) , (4.8)
which generalizes to the time-fractional diffusion equation the identity (2.7)
of the standard diffusion equation. Since
0 M(ζ; ν) dζ = 1 , see (A.40),
the function M(ζ; ν) is the normalized auxiliary function of the fractional
diffusion equation.
We note that for the time-fractional diffusion equation the fundamental
solution of the Cauchy problem is still a bilateral symmetric pdf in x (with
two branches, for x > 0 and x < 0 , obtained one from the other by
reflection), but is no longer of Gaussian type if ν 6= 1/2 . In fact, for large
|x| each branch exhibits an exponential decay in the ”stretched” variable
|x|1/(1−ν) as can be derived from the asymptotic representation (A.36) of the
auxiliary function M(·; ν) . In fact, by using (4.7-8) and (A.36), we obtain
Gc(x, t; ν) ∼ a∗(t) |x|(ν−1/2)/(1−ν) exp
−b∗(t)|x|1/(1−ν)
, (4.9)
as |x| → ∞ , where a∗(t) and b∗(t) are certain positive functions of time.
Furthermore, the exponential decay in x provided by (4.9) ensures that all
the absolute moments of positive order of Gc(x, t; ν) are finite. In particular,
using (4.8) and (A.39) it turns out that the moments (of even order) are
x2n Gc(x, t; ν) dx =
Γ(2n + 1)
Γ(2νn + 1)
(Dt2ν)n , n = 0 , 1 , 2 , . . . (4.10)
The formula (4.10) provides a generalization of the corresponding formula
(2.11) valid for the standard diffusion equation, ν = 1/2 . Furthermore, we
recognize that the variance associated to the pdf is now proportional to Dt2ν ,
which for ν 6= 1/2 implies a phenomenon of anomalous diffusion. According
to a usual terminology in statistical mechanics, the anomalous diffusion is
said to be slow if 0 < ν < 1/2 and fast if 1/2 < ν < 1 .
In Figure 1, as an example, we compare versus |x| , at fixed t , the
fundamental solutions of the Cauchy problem with different ν (ν =
1/4 , 1/2 , 3/4 ). We consider the range 0 ≤ |x| ≤ 4 and assume D = t = 1 .
0 1 2 3 4
Figure 1: The Cauchy problem for the time-fractional diffusion equation.
The fundamental solutions versus |x| with a) ν = 1/4 , b) ν = 1/2 , c)
ν = 3/4 .
We note the different behaviour of the pdf in the cases of slow diffusion (ν =
1/4 ) and fast diffusion (ν = 3/4 ) with respect to the Gaussian behaviour
of the standard diffusion (ν = 1/2). In the limiting cases ν = 0 and ν = 1
we have
Gc(x, t; 0) =
e−|x|
, Gc(x, t; 1) =
δ(x −
D t) + δ(x +
. (4.11)
We also recognize from the appendix B that for 1/2 ≤ ν < 1 any branch
of the fundamental solution is proportional to the corresponding positive
branch of an extremal stable pdf with index of stability α = 1/ν , which
exhibits an exponential decay at infinity. In fact, applying (B.29) with
α = 1/ν and y = ζ = |x|/(
Dtν) , from (4.7-8) we obtain
Gc(|x|, t; ν) =
|x|/(
D tν) ; − (2 − 1/ν)
· p1/ν (|x|; +1, 1, 0) , 1 < 1/ν ≤ 2 .
(4.12)
We also note that the stable distribution in (4.12) satisfies the condition
p1/ν (x; +1, 1, 0) dx = ν , 1 < 1/ν ≤ 2 . (4.13)
5 The Signalling problem for the time-fractional
diffusion equation
For the fractional diffusion equation (3.1) subject to (3.6b) the application
of the Laplace transform leads to the ordinary differential equation of order
ũ(x, s; ν) , ũ(0+, s; ν) = h̃(s) , ũ(+∞, s; ν) = 0 . (5.1)
Thus the transformed solution reads
ũ(x, s; ν) = h̃(s) e−(x/
D) sν , (5.2)
so for the Green function we have
Gs(x, t; ν) ÷ G̃s(x, s; ν) = e−(x/
D) sν . (5.3)
When ν 6= 1/2 the inversion of this Laplace transform cannot be obtained by
looking in a standard table of Laplace transform pairs. Also here we appeal
to a Laplace transform pair related to the Wright-type function M(ζ; ν). In
fact, using (A.40) with r = t , and taking into account the scaling property
of the Laplace transform, we obtain
Gs(x, t; ν) = ν
D t1+ν
. (5.4)
Introducing the similarity variable ζ = x/(
Dtν) , we recognize the identity
tGs(x, t; ν) = ν ζ M(ζ; ν) , (5.5)
which is the counterpart for the Signalling problem of the identity (4.8) valid
for the Cauchy problem.
Comparing (5.5) with (4.8) we obtain the reciprocity relation between the
two fundamental solutions of the time-fractional diffusion equation, in the
common domain x > 0 , t > 0 ,
2ν xGc(x, t; ν) = tGs(x, t; ν) . (5.6)
The interpretation of Gs(x, t; ν) as a one-sided stable pdf in time is
straightforward: in this respect we need to apply (B.28), with index of
stability α = ν and variable y = ζ−1/ν = t (
D/x)1/ν , in (5.5). We obtain
Gs(x, t; ν) =
; − ν
 = pν (t; −1, 1, 0) . (5.7)
In Figure 2, as an example, we compare versus t , at fixed x , the fundamental
solutions of the Signalling problem with different ν (ν = 1/4 , 1/2 , 3/4 ). We
consider the range 0 ≤ t ≤ 3 and assume D = x = 1 .
We note the different behaviour of the pdf in the cases of slow diffusion
(ν = 1/4 ) and fast diffusion (ν = 3/4 ) with respect to the Lévy pdf for the
standard diffusion (ν = 1/2). In the limiting cases ν = 0 , 1 , we have
Gs(x, t; 0) = δ(t) , Gs(x, t; 1) = δ(t − x/
D) . (5.8)
0 1 2 3
Figure 2: The Signalling problem for the time-fractional diffusion equation.
The fundamental solutions versus t with a) ν = 1/4 , b) ν = 1/2 , c)
ν = 3/4 .
6 The Cauchy problem for the symmetric space-
fractional diffusion equation
The symmetric space-fractional diffusion equation is obtained from the
classical diffusion equation by replacing the second-order space derivative by
a symmetric space-fractional derivative (explained below) of order α with
0 < α ≤ 2 . In our notation we write this equation as
∂|x|α
, u = u(x, t;α) , x ∈ R , t ∈ R+0 , 0 < α ≤ 2 , (6.1)
where D is a positive coefficient with the dimensions Lα T−1 . The
fundamental solution for the Cauchy problem, Gc(x, t;α) is the solution of
(6.1), subject to the initial condition u(x, 0+;α) = δ(x) .
The symmetric space-fractional derivative of any order α > 0 of a sufficiently
well-behaved function φ(x) , x ∈ R , may be defined as the pseudo-
differential operator characterized in its Fourier representation by
d|x|α
φ(x) ÷ −|κ|α φ̂(κ) , x , k ∈ R , α > 0 . (6.2)
According to a usual terminology, −|κ|α is referred to as the symbol of our
pseudo-differential operator, the symmetric space-fractional derivative, of
order α . Here, we have adopted the notation introduced by Zaslavski, see
e.g. Saichev & Zaslavski (1997).
In order to properly introduce this kind of fractional derivative we need
to consider a peculiar approach to fractional calculus different from the
Riemann-Liouville one, already treated in Appendix A. This approach is
indeed based on the so-called Riesz potentials (or integrals), that we prefer
to consider later.
At first, let us see how things become highly transparent by using an
heuristic argument, originally due to Feller (1952). The idea is to start
from the positive definite differential operator
A := −
÷ κ2 = |κ|2 , (6.3)
whose symbol is |κ|2 , and form positive powers of this operator as pseudo-
differential operators by their action in the Fourier-image space, i.e.
Aα/2 :=
= |κ|α α > 0 . (6.4)
Thus the operator −Aα/2 can be interpreted as the required fractional
derivative, i.e.
Aα/2 ≡ − d
d|x|α
, α > 0 . (6.5)
We note that the operator just defined must not be confused with a power
of the first order differential operator d
for which the symbol is −iκ .
After the above considerations it is straightforward to obtain the Fourier
image of the Green function of the Cauchy problem for the space-fractional
diffusion equation. In fact, applying the Fourier transform to the equation
(6.1), subject to the initial condition u(x, 0+;α) = δ(x) , and accounting for
(6.2), we obtain
Gc(x, t;α) = Gc(|x|, t;α) ÷ Ĝc(k, t;α) = e−D t |κ|
, 0 < α ≤ 2 . (6.6)
We easily recognize that the Fourier transform of the Green function
corresponds to the canonic form of a symmetric stable distribution of index
of stability α and scaling factor γ = (Dt)1/α , see (B.8). Therefore we have
Gc(x, t;α) = pα(x; 0, γ, 0) , γ = (Dt)1/α . (6.7)
For α = 1 and α = 2 we easily obtain the explicit expressions of the
corresponding Green functions since in these cases they correspond to the
Cauchy and Gauss distributions,
Gc(x, t; 1) =
x2 + (D t)2
, (6.8)
see (B.5), and
Gc(x, t; 2)) =
2/(4D t) , (6.9)
in agreement with (2.6).
We easily recognize that
(D t)1/α
(6.10)
is the similarity variable for the space-fractional diffusion equation, in terms
of which we can express the Green function for any α ∈ (0, 2] . Indeed, we
recognize that
Gc(x, t;α) =
(D t)1/α
qα(η; 0) , (6.11)
where qα(η; 0) denotes the symmetric stable distribution of order α with
Feller-type characteristic function, see (B.14-15). Now we can express the
Green function using the Feller series expansions (B.21-22) with θ = 0 . We
obtain:
for 0 < α < 1 ,
qα(η; 0) = −
Γ(nα + 1)
, (6.12a)
for 1 < α ≤ 2 ,
qα(η; 0) =
(−1)m
Γ[(2m + 1)/α]
(2m)!
η2m . (6.12b)
In the limiting case α = 1 the above series reduce to geometrical series and
therefore are no longer convergent in all of C . In particular, they represent
the expansions of the function q1(η; 0) = 1/[π(1+η
2)] , convergent for η > 1
and 0 < η < 1 , respectively.
We also note that for any α ∈ (0, 2] the functions qα(η; 0) exhibit at the
origin the value qα(0; 0) = Γ(1/α)/(π α) , and at the queues, excluding the
Gaussian case α = 2 , the algebraic asymptotic behaviour, as η → ∞ ,
qα(η; 0) ∼
Γ(α + 1) sin
η−(α+1) , 0 < α < 2 . (6.13)
In Figure 3, as an example, we compare versus x , at fixed t , the fundamental
solutions of the Cauchy problem with different α (α = 1/2 , 1 , 3/2 , 2 ). We
consider the range −6 ≤ x ≤ +6 and assume D = t = 1 .
-6 -4 -2 0 2 4 6
-6 -4 -2 0 2 4 6
Figure 3: The Cauchy problem for the simmetric space-fractional diffusion
equation. The fundamental solutions versus x : plate a) α = 1/2
(continuous line), α = 1 (dashed line); plate b) α = 3/4 (continuous line),
α = 2 (dashed line).
Let us now express more properly our operator (6.4) (with symbol |κ|α)
as inverse of a suitable integral operator Iα whose symbol is |κ|−α . This
operator can be found in the approach by Marcel Riesz to Fractional
Calculus, see e.g. Samko, Kilbas & Marichev (1987-1993) and Rubin (1996).
We recall that for any α > 0 , α 6= 1 , 3 , 5 , . . . and for a sufficiently well-
behaved function φ(x) , x ∈ R , the Riesz integral or Riesz potential Iα and
its image in the Fourier domain read
Iα φ(x) :=
2Γ(α) cos(πα/2)
|x − ξ|α−1 φ(ξ) dξ ÷ φ̂(κ)
. (6.14)
On its turn, the Riesz potential can be written in terms of two Weyl integrals
Iα± according to
Iα φ(x) =
2 cos(πα/2)
Iα+φ(x) + I
−φ(x)
, (6.15)
where 
Iα+ φ(x) :=
(x − ξ)α−1 φ(ξ) dξ ,
Iα− φ(x) :=
(ξ − x)α−1 φ(ξ) dξ .
(6.16)
Then, at least in a formal way, the space-fractional derivative (6.2) turns
out to be defined as the opposite of the (left) inverse of the Riesz fractional
integral, i.e.
d|x|α
φ(x) := −I−α φ(x) = −
2 cos(πα/2)
I−α+ φ(x) + I
− φ(x)
. (6.17)
Notice that (6.14) and (6.17) become meaningless when α is an integer odd
number. However, for our range of interest 0 < α ≤ 2 , the particular case
α = 1 can be singled out since the corresponding Green function is already
known, see (6.8). Thus, excluding the case α = 1 , our space-fractional
diffusion equation (6.1) can be re-written, x ∈ R , t ∈ R+0 , as
= −D I−α u , u = u(x, t;α) , 0 < α ≤ 2 , α 6= 1 , (6.18)
where the operator I−α is defined by (6.16-17).
Here, in order to evaluate the fundamental solution of the Cauchy problem,
interpreted as a probability density, we propose a numerical approach,
original as far as we know, based on a (symmetric) random walk model,
discrete in space and time, see also Gorenflo & Mainardi (1998a), Gorenflo
& Mainardi (1998b) and Gorenflo, De Fabritiis & Mainardi (1999). We shall
see how things become highly transparent, in that we properly generalize
the classical random-walk argument of the standard diffusion equation
to our space-fractional diffusion equation (6.18). So doing we are in
position to provide a numerical simulation of the related (symmetric) stable
distributions in a way analogous to the standard one for the Gaussian law.
The essential idea is to approximate the left inverse operators I−α± by the
Grünwald-Letnikov scheme, on which the reader can inform himself in the
treatises on fractional calculus, see e.g. Oldham & Spanier (1974), Samko,
Kilbas & Marichev (1987-1993), Miller & Ross (1993), or in the recent review
article by Gorenflo (1997). If h denotes a ”small” positive step-length, these
approximating operators read
± φ(x) :=
(−1)k
φ(x ∓ kh) . (6.19)
Assume, for simplicity, D = 1 , and introduce grid points xj = j h with
h > 0 , j ∈ Z , and time instances tn = n τ with τ > 0 , n ∈ N0 . Let there
be given probabilities pj,k ≥ 0 of jumping from point xj at instant tn to
point xk at instant tn+1 and define probabilities yj(tn) of the walker being
at point xj at instant tn. Then, by
yk(tn+1) =
pj,k uj(tn) ,
pj,k =
pj,k = 1 , (6.20)
with pj,k = pk,j , a symmetric random walk (more precisely a symmetric
random jump) model is described. With the approximation
yj(tn) ≈
∫ (xj+h/2)
(xj−h/2)
u(x, tn) dx ≈ hu(xj , tn) , (6.21)
and introducing the ”scaling parameter”
2 | cos(απ/2)|
, (6.22)
we have solved
yj(tn+1) − yj(tn)
= − hI−α yj(tn) , (6.23)
for yj(tn+1) . So we have proved to have a consistent (for h → 0) symmetric
random walk approximation to (6.18) by taking
i) for 0 < α < 1 , 0 < µ ≤ 1/2 ,
−α yj(tn) = µ
+ yj(tn) + hI
− yj(tn)
pj,j = 1 − 2µ , pj,j±k = µ
)∣∣ , k ≥ 1 ;
(6.24)
ii) for 1 < α ≤ 2 , 0 < µ ≤ 1/(2α) ,


−α yj(tn) = µ
+ yj+1(tn) + hI
− yj−1(tn)
pj,j = 1 − 2µ α , pj,j±1 = µ
pj,j±k = µ
)∣∣∣ , k ≥ 2 .
(6.25)
We note that our random walk model is not only symmetric, but also
homogeneous, the transition probabilities pj,j±k not depending on the index
In the special case α = 2 we recover from (6.25) the well-known three-point
approximation of the heat equation, because pj,j±k = 0 for k ≥ 2 . This
means that for approximation of common diffusion only jumps of one step
to the right or one to the left or jumps of width zero occur, whereas for
0 < α < 2 (α 6= 1) arbitrary large jumps occur with power-like decaying
probability, as it turns out from the asymptotic analysis for the transition
probabilities given in (6.24-25). In fact, as k → ∞ , one finds
pj,j+k ∼
(τ/hα)
Γ(α + 1) sin
k−(α+1) , 0 < α < 2 . (6.26)
This result thus provides the discrete counterpart of the asymptotic
behaviour of the long power-law tails of the symmetric stable distributions,
as foreseen by (6.13) when 0 < α < 2 .
7 Conclusions
We have treated two generalizations of the standard, one-dimensional,
diffusion equation, namely, the time-fractional diffusion equation and the
symmetric space-fractional diffusion equation. For these equations we have
derived the fundamental solutions using the transform methods of Fourier
and Laplace, and exhibited their connections to extremal and symmetric
stable probability densities, evolving on time or variable in space. For the
symmetric space-fractional diffusion equation we have presented a stationary
(in time), homogeneous (in space) symmetric random walk model, discrete
in space and time, the step-lengths of the spatial grid and the time lapses
between transitions properly scaled. In the limit of infinitesimally fine
discretization this model (based on the Grünwald-Letnikov approximation
to fractional derivatives) is consistent with the continuous diffusion process,
i.e. convergent if interpreted as a difference scheme in the sense of numerical
analysis2.
From the mathematical viewpoint the field of such ”fractional” general-
izations is fascinating as there several mathematical disciplines meet and
come to a fruitful interplay: e.g. probability theory and stochastic processes,
2Further generalizations have been considered by us and our collaborators in other
papers, in which we have given a derivation of discrete random walk models related to
more general space-time fractional diffusion equations. For a comprehensive analysis, see
Gorenflo et al. (2002). Readers interested to the fundamental solutions of these fractional
diffusion equations are referred to the paper by Mainardi et al. (2001) where analytical
expressions and numerical plots are found.
integro-differential equations, transform theory, special functions, numerical
analysis. As one may take from our References, one can observe that since
some decades there is an ever growing interest in using the concepts of
fractional calculus among physicists and economists. Among economists we
like to refer the reader to a collection of papers on the topic of ”Fractional
Differencing and Long Memory Processes”, edited by Baillie & King (1996).
Appendix A: The Riemann-Liouville Fractional
Calculus
Fractional calculus is the field of mathematical analysis which deals with the
investigation and applications of integrals and derivatives of arbitrary order.
The term fractional is a misnomer, but it is retained following the prevailing
use. This appendix is mostly based on the recent review by Gorenflo &
Mainardi (1997). For more details on the classical treatment of fractional
calculus the reader is referred to Erdélyi (1954), Oldham & Spanier (1974),
Samko et al. (1987-1993) and Miller & Ross (1993).
According to the Riemann-Liouville approach to fractional calculus, the
notion of fractional Integral of order α (α > 0) is a natural consequence
of the well known formula (usually attributed to Cauchy), that reduces the
calculation of the n−fold primitive of a function f(t) to a single integral of
convolution type. In our notation the Cauchy formula reads
Jnf(t) := fn(t) =
(n − 1)!
(t − τ)n−1 f(τ) dτ , t > 0 , n ∈ N , (A.1)
where N is the set of positive integers. From this definition we note that
fn(t) vanishes at t = 0 with its derivatives of order 1, 2, . . . , n − 1 . For
convention we require that f(t) and henceforth fn(t) be a causal function,
i.e. identically vanishing for t < 0. In a natural way one is led to extend
the above formula from positive integer values of the index to any positive
real values by using the Gamma function. Indeed, noting that (n − 1)! =
Γ(n) , and introducing the arbitrary positive real number α , one defines the
Fractional Integral of order α > 0 :
Jα f(t) :=
(t − τ)α−1 f(τ) dτ , t > 0 , α ∈ R+ , (A.2)
where R+ is the set of positive real numbers. For complementation we define
J0 := I (Identity operator), i.e. we mean J0 f(t) = f(t) . Furthermore, by
Jαf(0+) we mean the limit (if it exists) of Jαf(t) for t → 0+ ; this limit
may be infinite.
We note the semigroup property JαJβ = Jα+β , α , β ≥ 0 , which implies
the commutative property JβJα = JαJβ , and the effect of our operators Jα
on the power functions
Jαtγ =
Γ(γ + 1)
Γ(γ + 1 + α)
tγ+α , α ≥ 0 , γ > −1 , t > 0 . (A.3)
These properties are of course a natural generalization of those known when
the order is a positive integer.
Introducing the Laplace transform by the notation L {f(t)} :=∫∞
−st f(t) dt = f̃(s) , s ∈ C , and using the sign ÷ to denote a Laplace
transform pair, i.e. f(t) ÷ f̃(s) , we note the following rule for the Laplace
transform of the fractional integral,
Jα f(t) ÷ f̃(s)
, α ≥ 0 , (A.4)
which is the generalization of the case with an n-fold repeated integral.
After the notion of fractional integral, that of fractional derivative of order
α (α > 0) becomes a natural requirement and one is attempted to substitute
α with −α in the above formulas. However, this generalization needs some
care in order to guarantee the convergence of the integrals and preserve the
well known properties of the ordinary derivative of integer order.
Denoting by Dn with n ∈ N , the operator of the derivative of order n ,
we first note that Dn Jn = I , Jn Dn 6= I , n ∈ N , i.e. Dn is left-inverse
(and not right-inverse) to the corresponding integral operator Jn . In fact
we easily recognize from (A.1) that
Jn Dn f(t) = f(t) −
f (k)(0+)
, t > 0 . (A.5)
As a consequence we expect that Dα is defined as left-inverse to Jα. For
this purpose, introducing the positive integer m such that m − 1 < α ≤ m ,
one defines the Fractional Derivative of order α > 0 :
Dα f(t) := Dm Jm−α f(t) , m − 1 < α ≤ m , m ∈ N , (A.6)
namely
Dα f(t)=
Γ(m − α)
(t − τ)α+1−m
, m − 1 < α < m,
f(t) , α = m.
(A.6′)
Defining for complementation D0 = J0 = I , then we easily recognize that
Dα Jα = I , α ≥ 0 , and
Dα tγ =
Γ(γ + 1)
Γ(γ + 1 − α)
tγ−α , α ≥ 0 , γ > −1 , t > 0 . (A.7)
Of course, these properties are a natural generalization of those known when
the order is a positive integer.
Note the remarkable fact that the fractional derivative Dα f is not zero for
the constant function f(t) ≡ 1 if α 6∈ N . In fact, (A.7) with γ = 0 teaches
us that
Dα1 =
Γ(1 − α)
, α ≥ 0 , t > 0 . (A.8)
This, of course, is ≡ 0 for α ∈ N, due to the poles of the gamma function in
the points 0,−1,−2, . . .. We now observe that an alternative definition of
fractional derivative, originally introduced by Caputo (1967) (1969) in the
late sixties and adopted by Caputo and Mainardi (1971) in the framework
of the theory of Linear Viscoelasticity, is
Dα∗ f(t) := J
m−α Dm f(t) m − 1 < α ≤ m , m ∈ N , (A.9)
namely
D ∗α f(t) =
Γ(m − α)
f (m)(τ)
(t − τ)α+1−m
dτ , m − 1 < α < m,
f(t) , α = m.
(A.9′)
This definition is of course more restrictive than (A.6), in that requires
the absolute integrability of the derivative of order m. Whenever we use
the operator Dα∗ we (tacitly) assume that this condition is met. We easily
recognize that in general
Dα f(t) := Dm Jm−α f(t) 6= Jm−α Dm f(t) := Dα∗ f(t) , (A.10)
unless the function f(t) along with its first m − 1 derivatives vanishes at
t = 0+. In fact, assuming that the passage of the m-derivative under the
integral is legitimate, one recognizes that, for m − 1 < α < m and t > 0 ,
Dα f(t) = Dα∗ f(t) +
Γ(k − α + 1)
f (k)(0+) , (A.11)
and therefore, recalling the fractional derivative of the power functions (A.7),
f(t) −
f (k)(0+)
= Dα∗ f(t) . (A.12)
The alternative definition (A.9) for the fractional derivative thus incorpo-
rates the initial values of the function and of its integer derivatives of lower
order. The subtraction of the Taylor polynomial of degree m − 1 at t = 0+
from f(t) means a sort of regularization of the fractional derivative. In
particular, according to this definition, the relevant property for which the
fractional derivative of a constant is still zero can be easily recognized, i.e.
Dα∗ 1 ≡ 0 , α > 0 . (A.13)
We now explore the most relevant differences between the two fractional
derivatives (A.6) and (A.9). We agree to denote (A.9) as the Caputo
fractional derivative to distinguish it from the standard Riemann-Liouville
fractional derivative (A.6). We observe, again by looking at (A.7), that
Dαtα−1 ≡ 0 , α > 0 , t > 0 .
From above we thus recognize the following statements about functions
which for t > 0 admit the same fractional derivative of order α , with
m − 1 < α ≤ m , m ∈ N ,
Dα f(t) = Dα g(t) ⇐⇒ f(t) = g(t) +
α−j , (A.14)
Dα∗ f(t) = D
∗ g(t) ⇐⇒ f(t) = g(t) +
m−j . (A.15)
In these formulas the coefficients cj are arbitrary constants.
For the two definitions we also note a difference with respect to the formal
limit as α → (m − 1)+ ; from (A.6) and (A.9) we obtain respectively,
Dα f(t) → Dm J f(t) = Dm−1 f(t) ; (A.16)
Dα∗ f(t) → J Dm f(t) = Dm−1 f(t) − f (m−1)(0+) . (A.17)
We now consider the Laplace transform of the two fractional derivatives.
For the standard fractional derivative Dα the Laplace transform, assumed to
exist, requires the knowledge of the (bounded) initial values of the fractional
integral Jm−α and of its integer derivatives of order k = 1, 2, . . . ,m−1 . The
corresponding rule reads, in our notation,
Dα f(t) ÷ sα f̃(s) −
Dk J (m−α) f(0+) sm−1−k , (A.18)
where m − 1 < α ≤ m .
The Caputo fractional derivative appears more suitable to be treated by
the Laplace transform technique in that it requires the knowledge of the
(bounded) initial values of the function and of its integer derivatives of
order k = 1, 2, . . . ,m− 1 , in analogy with the case when α = m . In fact, by
using (A.4) and noting that
Jα Dα∗ f(t) = J
α Jm−α Dm f(t) = Jm Dm f(t) = f(t) −
f (k)(0+)
(A.19)
we easily prove the following rule for the Laplace transform,
Dα∗ f(t) ÷ sα f̃(s) −
f (k)(0+) sα−1−k , m − 1 < α ≤ m . (A.20)
Indeed, the result (A.20), first stated by Caputo (1969) by using the
Fubini-Tonelli theorem, appears as the most ”natural” generalization of the
corresponding result well known for α = m .
Gorenflo and Mainardi (1997) have pointed out the major utility of the
Caputo fractional derivative in the treatment of differential equations of
fractional order for physical applications. In fact, in physical problems,
the initial conditions are usually expressed in terms of a given number
of bounded values assumed by the field variable and its derivatives of
integer order, no matter if the governing evolution equation may be a
generic integro-differential equation and therefore, in particular, a fractional
differential equation3.
We now analyze the most simple differential equations of fractional order,
including those which, by means of fractional derivatives, generalize the well-
known ordinary differential equations related to relaxation and oscillation
3We note that the Caputo fractional derivative was so named after the book by
Podlubny (1999). It coincides with that introduced, independently and a few later,
by Dzherbashyan and Nersesyan (1968) as a regularization of the Riemann-Liouville
fractional derivative. Nowadays, some Authors refer to it as the Caputo-Dzherbashyan
fractional derivative. The prominent role of this fractional derivative in treating initial
value problems was recognized in interesting papers by Kochubei (1989), (1990).
phenomena. Generally speaking, we consider the following differential
equation of fractional order α > 0 ,
Dα∗ u(t) = D
u(t) −
u(k)(0+)
= −u(t) + q(t) , t > 0 , (A.21)
where u = u(t) is the field variable and q(t) is a given function. Here m is
a positive integer uniquely defined by m − 1 < α ≤ m , which provides the
number of the prescribed initial values u(k)(0+) = ck , k = 0, 1, 2, . . . ,m−1 .
Implicit in the form of (A.21) is our desire to obtain solutions u(t) for which
the u(k)(t) are continuous. In particular, the cases of fractional relaxation
and fractional oscillation are obtained for 0 < α < 1 and 1 < α < 2 ,
respectively
The application of the Laplace transform through the Caputo formula (A.20)
yields
ũ(s) =
sα−k−1
sα + 1
sα + 1
q̃(s) . (A.22)
Now, in order to obtain the Laplace inversion of (A.22), we need to recall
the Mittag-Leffler function of order α > 0 , Eα(z) . This function, so named
from the great Swedish mathematician who introduced it at the beginning
of this century, is defined by the following series and integral representation,
valid in the whole complex plane,
Eα(z) =
Γ(αn + 1)
σα−1 e σ
σα − z
dσ , α > 0 . (A.23)
Here Ha denotes the Hankel path, i.e. a loop which starts and ends at −∞
and encircles the circular disk |σ| ≤ |z|1/α in the positive sense. It turns out
that Eα(z) is an entire function of order ρ = 1/α and type 1 .
The Mittag-Leffler function provides a simple generalization of the expo-
nential function, to which it reduces for α = 1 . Particular cases from which
elementary functions are recovered, are
= cosh z , E2
= cos z , z ∈ C , (A.24)
E1/2(±z1/2) = ez
1 + erf (±z1/2)
= ez erfc (∓z1/2) , z ∈ C , (A.25)
where erf (erfc) denotes the (complementary) error function. defined as
erf (z) :=
du , erfc (z) := 1 − erf (z) , z ∈ C .
A noteworthy property of the Mittag-Leffler function is based on the
following duplication formula
Eα(z) =
Eα/2(+z
1/2) + Eα/2(−z1/2)
. (A.26)
In (A.25-26) we agree to denote by z1/2 the main branch of the complex
root of z .
The Mittag-Leffler function is connected to the Laplace integral through the
equation ∫ ∞
e−u Eα (u
α z) du =
1 − z
α > 0 . (A.27)
The integral at the L.H.S. was evaluated by Mittag-Leffler who showed that
the region of its convergence contains the unit circle and is bounded by the
line Re z1/α = 1 . The above integral is fundamental in the evaluation of the
Laplace transform of Eα (−λ tα) with α > 0 and λ ∈ C . In fact, putting in
(A.27) u = st and uα z = −λ tα with t ≥ 0 and λ ∈ C , we get the Laplace
transform pair
Eα (−λ tα) ÷
sα + λ
, Re s > |λ|1/α . (A.28)
Then, using (A.28), we put for k = 0, 1, . . . ,m − 1 ,
uk(t) := J
keα(t) ÷
sα−k−1
sα + 1
, eα(t) := Eα(−tα) , (A.29)
and, from inversion of the Laplace transforms in (A.22), we find
u(t) =
ck uk(t) −
q(t − τ)u′0(τ) dτ . (A.30)
In particular, the formula (A.30) encompasses the solutions for α = 1 , 2 ,
since e1(t) = exp(−t) , e2(t) = cos t . When α is not integer, namely for
m − 1 < α < m , we note that m − 1 represents the integer part of α
(usually denoted by [α]) and m the number of initial conditions necessary
and sufficient to ensure the uniqueness of the solution u(t). Thus the m
functions uk(t) = J
keα(t) with k = 0, 1, . . . ,m−1 represent those particular
solutions of the homogeneous equation which satisfy the initial conditions
+) = δk h , h, k = 0, 1, . . . ,m − 1 , and therefore they represent the
fundamental solutions of the fractional equation (A.21), in analogy with the
case α = m . Furthermore, the function uδ(t) = −u′0(t) = −e′α(t) represents
the impulse-response solution.
The Mittag-Leffler function of order less than one turns out to be related
through the Laplace integral to another special function of Wright type,
denoted by M(z, ν) with 0 < ν < 1 , following the notation introduced
by Mainardi (1994, 1995). Since this function turns out to be relevant in
the general framework of fractional calculus with special regard to stable
probability distributions, we are going to summarize its basing properties.
For more details on this function, see Mainardi (1997), Appendix A.
Let us first recall the more general Wright function Wλ,µ(z) , z ∈ C , with
λ > −1 and µ > 0 . This function, so named from the British mathematician
who introduced it between 1933 and 1941, is defined by the following series
and integral representation, valid in the whole complex plane,
Wλ,µ(z) =
n! Γ(λn + µ)
eσ + zσ
−λ dσ
, (A.31)
where Ha denotes the Hankel path. It is possible to prove that the Wright
function is entire of order 1/(1+λ) , hence of exponential type if λ ≥ 0 . The
case λ = 0 is trivial since W0,µ(z) = e
z/Γ(µ) . The case λ = −ν , µ = 1 − ν
with 0 < ν < 1 provides the function M(z, ν) of special interest for us.
Specifically, we have
M(z; ν) := W−ν,1−ν(−z) =
W−ν,0(−z) , 0 < ν < 1 , (A.32)
and therefore from (A.31-32)
M(z; ν) =
(−z)n−1
(n − 1)!
Γ(ν n) sin (ν n π)
eσ − zσ
, 0 < ν < 1 .
(A.33)
In the series representation we have used the reflection formula for the
Gamma function, Γ(x) Γ(1−x) = π/ sin πx . Explicit expressions of M(z; ν)
in terms of simpler known functions are expected in particular cases when
ν is a rational number. Relevant cases are ν = 1/2 , 1/3 for which
M(z; 1/2) =
− z2/4
, (A.34)
M(z; 1/3) = 32/3 Ai
z/31/3
, (A.35)
where Ai denotes the Airy function.
When the argument is real and positive, i.e. z = r > 0 , the existence of
the Laplace transform of M(r; ν) is ensured by the asymptotic behaviour,
as derived by Mainardi & Tomirotti (1995), as r → +∞ ,
M(r/ν; ν) ∼ a(ν) r(ν − 1/2)/(1 − ν) exp
−b(ν) r1/(1 − ν)
, (A.36)
where a(ν) = 1/
2π (1 − ν) , b(ν) = (1 − ν)/ν .
It is an instructive exercise to derive the Laplace transform by interchanging
the Laplace integral with the Hankel integral in (A.33) and recalling the
integral representation (A.23) of the Mittag-Leffler function. We obtain the
Laplace transform pair
M(r; ν) ÷ Eν(−s) , 0 < ν < 1 . (A.37)
For ν = 1/2 , (A.37) with (A.25) and (A.34) provides the result, see e.g.
Doetsch (1974),
M(r; 1/2) :=
− r2/4
÷ E1/2(−s) := exp
erfc (s) . (A.38)
It would be noted that, since M(r, ν) is not of exponential order,
transforming term-by-term the Taylor series of M(r; ν) yields a series of
negative powers of s , which represents the asymptotic expansion of Eν(−s)
as s → ∞ in a certain sector around the real axis.
We also note that (A.37) with (A.23) allows us to compute the moments of
any real order δ ≥ 0 of M(r; ν) in the positive real axis. We obtain
r δ M(r; ν) dr =
Γ(δ + 1)
Γ(νδ + 1)
, δ ≥ 0 . (A.39)
When δ is integer we note that the moments are provided by the derivatives
of the Mittag-Leffler function in the origin, i.e.
rn M(r; ν) dr = lim
(−1)n
Eν(−s) =
Γ(n + 1)
Γ(νn + 1)
, (A.40)
where n = 0, 1, 2, . . . . The normalization condition
0 M(r; ν) dr =
Eν(0) = 1 is recovered for n = 0 . The relation with the Mittag-Leffler
function stated in (A.40) can be extended to the moments of non integer
order if we replace the ordinary derivative, of order n, with the corresponding
fractional derivative, of order δ 6= n, in the Caputo sense.
Another exercise on the function M concerns the inversion of the Laplace
transform exp(−sν) , either by the complex integral formula or by the formal
series method. We obtain the Laplace transform pair
M (1/rν ; ν) ÷ exp (−sν) , 0 < ν < 1 . (A.41)
For ν = 1/2 , (A.41) with (A.34) provides the known result, see e.g. Doetsch
(1974),
2 r3/2
M(1/r1/2; 1/2) :=
π r3/2
exp [− 1/(4r)] ÷ exp
−s1/2
. (A.42)
We recall that a rigorous proof of (A.41) was formerly given by Pollard
(1946), based on a formal result by Humbert (1945). The Laplace transform
pair was also obtained by Mikusiński (1959) and, albeit unaware of the
previous results, by Buchen & Mainardi (1975) in a formal way.
Appendix B: The Stable Probability Distributions
The stable distributions are a fascinating and fruitful area of research in
probability theory; furthermore, nowadays, they provide valuable models in
physics, astronomy, economics, and communication theory.
The general class of stable distributions was introduced and given this name
by the French mathematician Paul Lévy in the early 1920’s, see Lévy (1924,
1925). The inspiration for Lévy was the desire to generalize the celebrated
Central Limit Theorem, according to which any probability distribution
with finite variance belongs to the domain of attraction of the Gaussian
distribution.
Formerly, the topic attracted only moderate attention from the leading
experts, though there were also enthusiasts, of whom the Russian
mathematician Alexander Yakovlevich Khintchine should be mentioned first
of all. The concept of stable distributions took full shape in 1937 with the
appearance of Lévy’s monograph, see Lévy (1937-1954), soon followed by
Khintchine’s monograph, see Khintchine (1938).
The theory and properties of stable distributions are discussed in some
classical books on probability theory including Gnedenko & Kolmogorov
(1949-1954), Lukacs (1960-1970), Feller (1966-1971), Breiman (1968-1992),
Chung (1968-1974) and Laha & Rohatgi (1979). Also treatises on fractals
devote particular attention to stable distributions in view of their properties
of scale invariance, see e.g. Mandelbrot (1982) and Takayasu (1990). Sets of
tables and graphs have been provided by Mandelbrot & Zarnfaller (1959),
Fama & Roll (1968), Bo’lshev & Al. (1968) and Holt & Crow (1973).
Only recently, monographs devoted solely to stable distributions and related
stochastic processes have been appeared, i.e. Zolotarev (1983-1986), Janicki
& Weron (1994), Samorodnitsky & Taqqu (1994), Uchaikin & Zolotarev
(1999). We now can cite the paper by Mainardi, Luchko & Pagnini (2001)
where the reader can find (convergent and asymptotic) representations and
plots of the symmetric and non-symmetric stable densities generated by
fractional diffusion equations.
Stable distributions have three exclusive properties, which can be briefly
summarized stating that they 1) are invariant under addition, 2) possess
their own domain of attraction, and 3) admit a canonic characteristic
function.
Let us now illustrate the above properties which, providing necessary and
sufficient conditions, can be assumed as equivalent definitions for a stable
distribution. We recall the basic results without proof.
A random variable X is said to have a stable distribution P (x) = Prob {X ≤
x} if for any n ≥ 2 , there is a positive number cn and a real number dn such
X1 + X2 + . . . + Xn
= cn X + dn , (B.1)
where X1,X2, . . . Xn denote mutually independent random variables with
common distribution P (x) with X . Here the notation
= denotes equality
in distribution, i.e. means that the random variables on both sides have the
same probability distribution.
When mutually independent random variables have a common distribution
[shared with a given random variable X], we also refer to them as
independent, identically distributed (i.i.d) random variables [independent
copies of X]. In general, the sum of i.i.d. random variables becomes
a random variable with a distribution of different form. However, for
independent random variables with a common stable distribution, the sum
obeys to a distribution of the same type, which differs from the original
one only for a scaling (cn) and possibly for a shift (dn). When in (B.1) the
dn = 0 the distribution is called strictly stable.
It is known, see Feller (1966-1971), that the norming constants in (B.1) are
of the form
cn = n
1/α with 0 < α ≤ 2 . (B.2)
The parameter α is called the characteristic exponent or the index of stability
of the stable distribution.
We agree to use the notation X ∼ Pα(x) to denote that the random variable
X has a stable probability distribution with characteristic exponent α . We
simply refer to P (x) , p(x) := dP/dx (probability density function = pdf)
and X as α-stable distribution, density, random variable, respectively.
The definition (B.1) with the theorem (B.2) can be stated in an alternative
version that needs only two i.i.d. random variables. see also Lukacs (1960-
1970). A random variable X is said to have a stable distribution if for any
positive numbers A and B, there is a positive number C and a real number
D such that
AX1 + B X2
= C X + D , (B.3)
where X1 and X2 are independent copies of X . Then there is a number
α ∈ (0, 2] such that the number C in (B.3) satisfies Cα = Aα + Bα .
For a strictly stable distribution (B.3) holds with D = 0 . This implies that
all linear combinations of i.i.d. random variables obeying to a strictly stable
distribution is a random variable with the same type of distribution.
A stable distribution is called symmetric if the random variable −X has the
same distribution. Of course, a symmetric stable distribution is necessarily
strictly stable.
Noteworthy examples of stable distributions are provided by the Gaussian
(or normal) law (with α = 2) and by the Cauchy-Lorentz law (α = 1). The
corresponding pdf are known to be
pG(x;σ, µ) :=
e−(x − µ)
2/(2σ2) , x ∈ R , (B.4)
where σ2 denotes the variance and µ the mean, and
pC(x; γ, δ) :=
(x − δ)2 + γ2
, x ∈ R , (B.5)
where γ denotes the semi-interquartile range and δ the ”shift”.
Another (equivalent) definition states that stable distributions are the only
distributions that can be obtained as limits of normalized sums of i.i.d.
random variables. A random variable X is said to have a domain of
attraction,i.e. if there is a sequence of i.i.d. random variables Y1, Y2, . . .
and sequences of positive numbers {γn} and real numbers {δn}, such that
Y1 + Y2 + . . . Yn
d⇒X . (B.6)
The notation
d⇒ denotes convergence in distribution.
It is clear that the previous definition (B.1) yields (B.6), e.g. , by taking the
Yis to be independent and distributed like X . The converse is easy to show,
see Gnedenko & Kolmogorov (1949-1954). Therefore we can alternatively
state that a random variable X is said to have a stable distribution if it has
a domain of attraction.
When X is Gaussian and the Yis are i.i.d. with finite variance, then (B.6)
is the statement of the ordinary Central Limit Theorem. The domain
of attraction of X is said normal when γn = n
1/α ; in general, γn =
n1/α h(n) where h(x) , x > 0 , is a slow varying function at infinity, that
is, lim
h(ux)/h(x) = 1 for all u > 0 , see Feller (1971). The function
h(x) = log x , for example, is slowly varying at infinity.
Another definition specifies the canonic form that the characteristic function
(cf) of a stable distribution of index α must have. Recalling that the cf is
the Fourier transform of the pdf , we use the notation p̂α(κ) := 〈exp (iκX)〉 ÷
pα(x) . We first note that a stable distribution is also infinitely divisible, i.e.
for every positive integer n its cf can be expressed as the nth power of
some cf . In fact, using the characteristic function, the relation (B.1) is
transformed into
[p̂α(κ)]
n = p̂α(cn κ) e
idnκ . (B.7)
The functional equation (B.7) can be solved completely and the solution is
known to be
p̂α(κ;β, γ, δ) = exp {iδκ − γα |κ|α [1 + i (sign κ)β ω(|κ|, α)]} , (B.8)
where
ω(|κ|, α) =
tan (α π/2) , if α 6= 1 ,
−(2/π) log |κ| , if α = 1 . (B.9)
Consequently a random variable X is said to have a stable distribution if
there are four real parameters α, β, γ, δ with 0 < α ≤ 2 , −1 ≤ β ≤ +1 ,
γ > 0 , such that its characteristic function has the canonic form (B.8-9).
Then we write pα(x;β, γ, δ)÷ p̂α(κ;β, γ, δ) and X ∼ Pα(x;β, γ, δ) , so partly
following the notation of Holt & Crow (1973) and Samorodnitsky & Taqqu
(1994).
We note in (B.8-9) that β appears with different signs for α 6= 1 and α = 1 .
This minor point has been the source of great confusion in the literature, see
Hall (1980) for a discussion. The presence of the logarithm for α = 1 is the
source of many difficulties, so this case has often to be treated separately.
The cf (B.8-9) turns out to be a useful tool for studying α-stable distri-
butions and for providing an interpretation of the additional parameters,
β (skewness parameter), γ (scale parameter) and δ (shift parameter), see
Samorodnitsky & Taqqu (1994). When α = 2 the cf refers to the Gaussian
distribution with variance σ2 = 2 γ2 and mean µ = δ ; in this case the value
of the skewness parameter β is not specified because tan π = 0 , and one
conventionally takes β = 0 .
One easily recognizes that a stable distribution is symmetric if and only if
β = δ = 0 and is symmetric about δ if and only if β = 0 . Stable distributions
with extremal values of the skewness parameter are called extremal. One
can prove that all the extremal stable distributions with 0 < α < 1 are
one-sided, the support being R+0 if β = −1 , and R
0 if β = +1 .
For the stable distributions Pα(x;β, γ, δ) we now consider the asymptotic
behaviour of the tail probabilities, T+(λ) := Prob {X > λ} and T−(λ) :=
Prob {X < −λ} , as λ → ∞ . For the Gaussian case α = 2 the result is well
known, see e.g. Feller (1957),
α = 2 : T±(λ) ∼ 1
2/(4γ2)
, λ → ∞ . (B.10)
Because of the above exponential decay all the moments of the corresponding
pdf turn out to be finite, which is an exclusive property of this stable
distribution. For all the other stable distributions the singularity of the
characteristic function in the origin is responsible for the algebraic decay of
the tail probabilities as indicated below, see e.g. Samorodnitsky & Taqqu
(1994),
0 < α < 2 : lim
λα T±(λ) = Cα γ
α (1 ∓ β)/2 , (B.11)
where
x−α sin x dx
1 − α
Γ(2 − a) cos (απ/2)
, if α 6= 1 ,
2/π , if α = 1 .
(B.12)
We note that for extremal distributions (β = ±1) the above algebraic decay
holds true only for one tail, the left one if β = +1 , the right one if β = −1 .
The other tail is either identically zero if 0 < α < 1 (the distribution is
one-sided !), or exhibits an exponential decay if 1 ≤ α < 2 . Because of the
algebraic decay we recognize that
0 < α < 2 :
|x|>λ
pα(x;β, γ, δ) dx = O(λ
−α) , (B.13)
so the absolute moments of a stable non-Gaussian pdf turn out to be finite
if their order ν is 0 ≤ ν < α and infinite if ν ≥ α . We are now convinced
that the Gaussian distribution is the unique stable distribution with finite
variance. Furthermore, when α ≤ 1 , the first absolute moment 〈|X|〉 is
infinite as well, so we need to use the median to characterize the expected
value.
There is however a fundamental property shared by all the stable
distributions that we like to point out: for any α the stable pdf are unimodal
and indeed bell-shaped, i.e. their n-th derivative has exactly n zeros, see
Gawronski (1964).
We now come back to the cf of a stable distribution, in order to provide for
α 6= 1 and δ = 0 a simpler canonic form which allow us to derive convergent
and asymptotic power series for the corresponding pdf . We first note that
the two parameters γ and δ in (B.8), being related to a scale transformation
and a translation, are not so essential since they do not change the shape
of distributions. If we take γ = 1 and δ = 0 , we obtain the so-called
standardized form of the stable distribution and X ∼ Pα(x;β, 1, 0) is referred
to as the α-stable standardized random variable. Furthermore, we can choose
the scale parameter γ in such a way to get from (B.8-9) the simplified canonic
form used by Feller (1952, 1966-1971) and Takayasu (1990) for strictly stable
distributions (δ = 0) with α 6= 1 , which reads in an ad hoc notation,
q̂α(κ; θ) :=
eiκ y pα(y; θ) dy = exp
−|κ|α e±i θ π/2
, (B.14)
where the symbol ± takes the sign of κ . This canonic form, that we refer to
as the Feller canonic form, is derived from (B.8-9) if in addition to α 6= 1
and δ = 0 we require
γα = cos
, tan
= β tan
. (B.15)
Here θ is the skewness parameter instead of β and its domain is restricted
in the following region (depending on α)
|θ| ≤
α , if 0 < α < 1 ,
2 − α , if 1 < α < 2 . (B.16)
Thus, when we use the Feller canonic form for strictly stable distributions
with index α 6= 1 and skewness θ , we implicitly select the scale parameter
γ (0 < γ ≤ 1), which is related to α , β and θ by (B.15). Specifically, the
random variable Y ∼ Qα(y; θ) turns out to be related to the standardized
random variable X ∼ Pα(x;β, 1, 0) by the following relations
Y = X/γ , pα(x;β, 1, 0) = γ qα(y = γx; θ) , (B.17)
with 
γ = [cos (θπ/2)]1/α ,
θ = (2/π) arctan [β tan (απ/2)] ,
tan (θπ/2)
tan (απ/2)
(B.18)
We recognize that qα(y, θ) = qα(−y,−θ) , so the symmetric stable
distributions are obtained if and only if θ = 0 . We note that for the
symmetric stable distributions we get the identity between the standardized
and the Lévy canonic forms, since in (B.18) β = θ = 0 implies γ = 1 .
A particular but noteworthy case is provided by p2(x; 0, 1, 0) = q2(y; 0) ,
corresponding to the Gaussian distribution with variance σ2 = 2 .
The extremal stable distributions, corresponding to β = ±1 , are now
obtained for θ = ±α if 0 < α < 1 , and for θ = ∓(2 − α) if 1 < α < 2 ; for
them the scaling parameter turns out to be γ = [cos (|α|π/2)]1/α . It may be
an instructive exercise to carry out the inversion of the Fourier transform
when α = 1/2 and θ = −1/2 . In this case we obtain the analytical expression
for the corresponding extremal stable pdf , known as the (one-sided) Lévy-
Smirnov density,
q1/2(y;−1/2) =
y−3/2 e−1/(4y) , y ≥ 0 . (B.19)
The standardized form for this distribution can be easily obtained from
(B.19) using (B.17-18) with α = 1/2 and θ = −1/2 . We get γ =
[cos (−π/4)]2 = 1/2 , β = −1 , so
p1/2(x;−1, 1, 0) =
q1/2(x/2;−1/2) =
x−3/2 e−1/(2x) , (B.20)
where x ≥ 0 , in agreement with Holt & Crow (1973) [§2.13, p. 147].
Feller (1952) has obtained from (B.14) the following representations by
convergent power series for the stable distributions valid for y > 0 , with
0 < α < 1 (negative powers),
qα(y; θ) =
(−y−α)n Γ(nα + 1)
(θ − α)
, (B.21)
1 < α ≤ 2 (positive powers),
qα(y; θ) =
(−y)n
Γ(n/α + 1)
(θ − α)
. (B.22)
The values for y < 0 can be obtained from (B.21-22) using the identity
qα(−y; θ) = qα(y;−θ) , y > 0 . As a consequence of the convergence in all of
C of the series in (B.21-22) we recognize that the restrictions of the functions
y qα(y; θ) on the two real semi-axis turn out to be equal to certain entire
functions of argument 1/|y|α for 0 < α < 1 and argument |y| for 1 < α ≤ 2 .
It has be shown, see e.g. Bergström (1952), Chao Chung-Jeh (1953), that the
two series in (B.21-22) provide also the asymptotic (divergent) expansions to
the stable pdf with the ranges of α interchanged from those of convergence.
From (B.21-22) a relation between stable pdf with index α and 1/α can be
derived as noted in Feller (1966-1971). Assuming 1/2 < α < 1 and y > 0 ,
we obtain
q1/α(y
−α; θ) = qα(y; θ
∗) , θ∗ = α(θ + 1) − 1 . (B.23)
A quick check shows that θ∗ falls within the prescribed range, |θ∗| ≤ α ,
provided that |θ| ≤ 2 − 1/α .
We now consider two particular cases of the Feller series (B.21-22), of
particular interest for us, which turn out to be related to the entire function
of Wright type, M(z; ν) with 0 < ν < 1 , reported in Appendix A. These
cases correspond to the following extremal distributions
Φ1(y) := qα(y;−α) , y > 0 , 0 < α < 1 , (B.24)
Φ2(y) := qα(y;α − 2) , y > 0 , 1 < α ≤ 2 , (B.25)
for which the Feller series (B.21-22) reduce to
Φ1(y) =
(−1)n−1 y−αn−1 Γ(nα + 1)
sin (nπα) , y > 0 , (B.26)
Φ2(y) =
(−1)n−1 yn−1 Γ(n/α + 1)
, y > 0 . (B.27)
In fact, recalling the series representation of the general Wright function,
Wλ,µ(z) with λ > −1 , µ > 0 , see (A.31), and the definition of the function
M(z; ν) with 0 < ν < 1 , see (A.32-33), we recognize that
Φ1(y) =
W−α,0(−y−α) =
M(y−α;α) , y > 0 , (B.28)
Φ2(y) =
W−1/α,0(−y) =
M(y; 1/α) , y > 0 . (B.29)
We would like to remark that the above relations with the Wright functions
have been noted also by Engler (1997).
It is worth to point out that, whereas Φ1(y) totally represents the one-
sided stable pdf qα(y;−α) , 0 < α < 1 , with support in R+0 , Φ2(y) is the
restriction on the positive axis of qα(y;α− 2) , 1 < α ≤ 2 , whose support is
all of R . Since the function M(z; ν) turns out to be normalized in R+0 , see
(A.39-40), we also note
Φ1(y) dy = 1 ;
Φ2(y) dy = 1/α . (B.30)
Using the results (A.41) and (A.37) we can easily evaluate the Laplace
transforms of Φ1(y) and Φ2(y) , respectively. We obtain
L[Φ1(y)] = Φ̃1(s) = exp (−sα) , 0 < α < 1 , (B.31)
L[Φ2(y)] = Φ̃2(s) =
E1/α (−s) , 1 < α ≤ 2 , (B.32)
where E1/α(·) denotes the Mittag-Leffler function of order 1/α , see (A.23).
It is an instructive exercise to derive the asymptotic behaviours of Φ1(y) and
Φ2(y) as y → 0+ and y → +∞ . By using the expressions (B.28−29) in terms
of the function M and recalling the series and asymptotic representations of
this function, see (A.33) and (A.36), we obtain
Φ1(y) =
y−(2−α)/[2(1−α)] e−c1 y
−α/(1−α)
, as y → 0+ ,
Γ(1 − α)
y−α−1 [1 + O (y−α)] , as y → +∞ ,
(B.33)
Φ2(y) =
Γ(1 − 1/α)
[1 + O (y)] , as y → 0+ ,
y(2−α)/[2(α−1)] e−c2 y
α/(α−1)
, as y → +∞ ,
(B.34)
where c1 , c2 are positive constants depending on α . We note that the
exponential decay is found for Φ1(y) as y → 0+ but as y → +∞ for Φ2(y) .
Explicit expressions for stable pdf can be derived form those for the function
M(z; ν) when ν = 1/2 and ν = 1/3 , given in Appendix A, see (A.34-
35). Of course the ν = 1/2 expression can be used to recover the well-
known (symmetric) Gaussian distribution q2(y; 0) accounting for (B.29), and
the (one-sided) Lévy distribution q1/2(y;−1/2), see (B.19), accounting for
(B.28). The ν = 1/3 expression provides, accounting for (B.28),
q1/3(y;−1/3) = 3−1/3 y−4/3 Ai
(3y)−1/3
y−3/2 K1/3
(B.35)
where Ai denotes the Airy function and K1/3 the modified Bessel function of
the second kind of order 1/3 . The equivalence between the two expressions
in (B.35) can be proved in view of the relation, see Abramowitz & Stegun
(1965-1972) [(10.4.14)],
Ai (z) =
. (B.36)
The case α = 1/3 has also been discussed by Zolotarev (1983-1986), who
has quoted the corresponding expression of the pdf in terms of K1/3 .
A general representation of all stable distributions (thus including the
extremal distributions above considered) in terms of special functions has
been only recently achieved by Schneider (1986). In his remarkable
(but almost ignored) article, Schneider has established that all the stable
distributions can be characterized in terms of a general class of special
functions, the so-called Fox H functions, so named after Charles Fox (1961).
For details on Fox H functions, see e.g. the books Mathai & Saxena (1978),
Srivastava & Al. (1982) and the most recent paper by Kilbas and Saigo
(1999). These functions are expressed in terms of special integrals in the
complex-plane, the Mellin-Barnes integrals4.
4The names refer to the two authors, who in the first 1910’s developed the theory of
these integrals using them for a complete integration of the hypergeometric differential
equation. However, as pointed out in the the Bateman Project Handbook on High
Transcendental Functions, see Erdelyi (1953), these integrals were first used by S. Pincherle
in 1888. For a revisited analysis of the pioneering work of Pincherle (1853-1936, Professor
of Mathematics at the University of Bologna from 1880 to 1928) we refer to the paper by
Mainardi and Pagnini (2003).
References
[1] Abramowitz, M. and Stegun, I.A. (Editors) : Handbook of Mathematical
Functions, Dover, New York 1965. [reprint, 1972]
[2] Baillie, R.T. and King, M.L. (Editors) : Fractional Differencing and
Long Memory Processes, Journal of Econometrics, 73 1-324 (1996).
[3] Bergström, H. : On some expansions of stable distribution functions,
Ark. Mat., 2 375-378 (1952).
[4] Bol’shev, L.N., Zolotarev, V.M., Kedrova, E.S. and Rybinskaya, M.A.
: Tables of cumulative functions of on-sided stable distributions, Theor.
Probability Appl., 15 299-309 (1968).
[5] Bouchaud, J.-P and Potters, M. : Theory of Financial Risk: from Sta-
tistical Physics to Risk Management, Cambridge Univ. Press, Cambridge
2000. [English enlarged version of Théories des Risques Financiers, CEA,
Aléa, Saclay 1997.]
[6] Breiman, L. : Probability, SIAM, Philadelphia 1992. [The SIAM edition
is an unabridged, corrected republication of the 1st edn, Addison-Wesley,
Reading Mass. 1968]
[7] P.W. Buchen, and Mainardi, F. : Asymptotic expansions for transient
viscoelastic waves, J. de Mécanique, 14 597-608 (1975).
[8] Caputo, M. : Linear models of dissipation whose Q is almost frequency
independent, Part II., Geophys. J. R. Astr. Soc., 13 529-539 (1967).
[9] Caputo, M.: Elasticità e Dissipazione, Zanichelli, Bologna 1969.
[10] Caputo, M. and Mainardi, F. : Linear models of dissipation in anelastic
solids, Rivista del Nuovo Cimento (Ser. II), 1 161-198 (1971).
[11] Chung-Jeh, Chao : Explicit formula for the stable law of distribution,
Acta Math. Sinica, 3 177-185 (1953). [in Chinese with English summary]
[12] Chung, K.L. : A Course in Probability Theory, 2nd edn, Academic
Press, New York 1974. [1st edn, Harcourt Brace Jowvanovich, 1968]
[13] Dzherbashyan, M.M. and Nersesyan, A.B. : Fractional derivatives and
the Cauchy problem for differential equations of fractional order. Izv.
Acad. Nauk Armjanskvy SSR, Matematika, 3 3–29 (1968). [In Russian]
[14] Engler, H. : Similarity solutions for a class of hyperbolic integrodiffer-
ential equations, Differential Integral Eqns, 10 815-840 (1997).
[15] Erdélyi, A. (Editor) : Higher Transcendental Functions, Bateman
Project, McGraw-Hill, New York 1953; Vol. 1, Ch. 1, §1.19, p. 49.
[16] Erdélyi, A. (Editor) : Higher Transcendental Functions, Bateman
Project, McGraw-Hill, New York 1955; Vol. 3, Ch. 18, pp. 206-227.
[17] Erdélyi, A. (Editor) : Tables of Integral Transforms, Bateman Project
McGraw-Hill, New York 1954; Vol. 2, Ch. 13, pp. 181-212.
[18] Fama, E. and Roll, R. : Some properties of symmetric stable
distributions, J. Amer. Statist. Assoc., 63 817-836 (1968).
[19] Feller, W. : On a generalization of Marcel Riesz’ potentials and
the semi-groups generated by them, Meddelanden Lunds Universitets
Matematiska Seminarium (Comm. Sém. Mathém. Université de Lund),
Tome suppl. dédié a M. Riesz, Lund (1952), pp. 73-81.
[20] Feller, W. : An Introduction to Probability Theory and its Applications,
Vol. 1, 3rd edn, Wiley, New York 1968. [1st edn, 1957]
[21] Feller, W. : An Introduction to Probability Theory and its Applications,
Vol. 2, 2nd edn, Wiley, New York 1971. [1st edn. 1966]
[22] Fox, C. : The G and H functions as symmetrical Fourier kernels, Trans.
Amer. Math. Soc., 98 395-429 (1961).
[23] Fujita, Y. : Integrodifferential equation which interpolates the heat
equation and the wave equation, Osaka J. Math., 27 309-321, 797-804
(1990). [2 papers]
[24] Gawronski, W. : On the bell-shape of stable distributions, Annals of
Probability, 12 230-242 (1984).
[25] Gnedenko, B.V. and Kolmogorov, A.N. : Limit Distributions for Sums
of Independent Random Variables, Addison-Wesley, Cambridge, Mass.
1954. [English Transl. from the Russian edition 1949, with notes by K.L.
Chung, revised 1968]
[26] Gorenflo, R. : Fractional calculus: some numerical methods, in
Carpinteri, A. and Mainardi, F. (Eds.), Fractals and Fractional Calculus
in Continuum Mechanics, CISM Courses and Lectures # 378, Springer
Verlag, Wien 1997, pp. 277-290. [Reprinted in www.fracalmo.org]
[27] Gorenflo, R. and Mainardi, F. : Fractional calculus: integral and
differential equations of fractional order, in Carpinteri, A. and Mainardi,
F. (Eds.), Fractals and Fractional Calculus in Continuum Mechanics,
CISM Courses and Lectures # 378, Springer Verlag, Wien 1997, pp. 223-
276. [Reprinted in www.fracalmo.org]
[28] Gorenflo, R. and Mainardi, F. : Fractional calculus and stable
probability distributions, Archives of Mechanics, 50 377-388 (1998a).
[29] Gorenflo, R. and Mainardi, F. : Random walk models for space-
fractional diffusion processes, Fractional Calculus and Applied Analysis,
1 167-190 (1998b).
[30] Gorenflo, R. De Fabritiis, G. and Mainardi, F. : Discrete random walk
models for symmetric Lévy-Feller diffusion processes, Physica A, 269 79–
89 (1999).
[31] Gorenflo, R., Mainardi, F., Moretti, D., Pagnini, G. and Paradisi,
P. : Discrete random walk models for space-time fractional diffusion,
Chemical Physics, 284 521-541 (2002). Special issue on Strange Kinetics
Guest Editors: R. Hilfer, R. Metzler, A. Blumen, J. Klafter. [E-print
arXiv:cond-mat/0702072]
[32] Hall, P. : A comedy of errors: the canonical form for a stable
characteristic function, Bull. London Math. Soc., 13 23-27 (1980).
[33] Holt, D.R. and Crow, E.L. : Tables and graphs of the stable probability
density functions, J. Res. Nat. Bureau Standards, 77B 143-198 (1973).
[34] Humbert, P. : Nouvelles correspondances symboliques, Bull. Sci.
Mathém. (Paris, II ser.), 69 121-129 (1945).
[35] Janicki, A. and Weron, A. : Simulation and Chaotic Behavior of α-
Stable Stochastic Processes, Marcel Dekker, New York 1994.
[36] Khintchine, A.Y. : Limit Laws for Sums of Independent Variables,
ONTI, Moscow 1938 [in Russian]
[37] Kilbas, A.A. and Saigo, M. : On the H functions, Journal of Applied
Mathematics and Stochastic Analysis, 12 191-204 (1999).
[38] Kochubei, A.N. : A Cauchy problem for evolution equations of
fractional order, Differential Equations, 25 967–974 (1989). [English
translation from the Russian Journal Differentsial’nye Uravneniya]
[39] Kochubei, A.N. : Fractional order diffusion, Differential Equations,
26 485–492 (1990). [English translation from the Russian Journal
Differentsial’nye Uravneniya]
[40] Laha, R.G. and Rohatgi, V.K. : Probability Theory, Wiley, New York
1979.
[41] Lévy, P. : Théorie des erreurs. La Loi de Gauss et les lois exceptionelles,
Bull. Soc. Math. France, 52 49-85 (1924).
[42] Lévy, P. : Calcul des probabilités, Gauthier-Villars, Paris 1925: Part II,
Chap. 6.
[43] Lévy, P. : Théorie de l’addition des variables aléatoires, 2nd edn,
Gauthier-Villars, Paris 1954. [1st edn, 1937]
[44] Lukacs, E. : Characteristic Functions, 2nd edn, Griffin, London 1970.
[1st edn, 1960]
[45] Mainardi, F. : On the initial value problem for the fractional diffusion-
wave equation, in S. Rionero and T. Ruggeri (Eds), Waves and Stability
in Continuous Media, World Scientific, Singapore 1994, pp. 246-251.
[46] Mainardi, F. : The time fractional diffusion-wave equation, Radiofisika,
38 20-36 (1995). [English Translation: Radiophysics & Quantum
Electronics]
[47] Mainardi, F. : Fractional relaxation-oscillation and fractional diffusion-
wave phenomena, Chaos, Solitons & Fractals, 7 1461-1477 (1996).
[48] Mainardi, F. : Fractional calculus: some basic problems in continuum
and statistical mechanics, in Carpinteri, A. and Mainardi, F. (Eds),
Fractals and Fractional Calculus in Continuum Mechanics, CISM Courses
and Lectures # 378, Springer-Verlag, Wien 1997, pp. 291-348. [Reprinted
in www.fracalmo.org]
[49] Mainardi, F., Luchko, Yu. and Pagnini, G. : The fundamental solution
of the space-time fractional diffusion equation, Fractional Calculus and
Applied Analysis, 4 153-192 (2001). [E-print arXiv:cond-mat/0702419]
[50] Mainardi, F. and Pagnini, G. : Salvatore Pincherle, the pioneer of the
Mellin-Barnes integrals, Jour. Computational and Applied Mathematics,
153 331-342 (2003).
[51] Mainardi, F. and Tomirotti, M. : On a special function arising in the
time fractional diffusion-wave equation, in P. Rusev, I. Dimovski and V.
Kiryakova (Eds), Transform Methods and Special Functions, Sofia 1994,
Science Culture Technology, Singapore 1995, pp. 171-183.
[52] Mainardi, F. and Tomirotti, M. : Seismic pulse propagation with
constant Q and stable probability distributions, Annali di Geofisica, 40
1311-1328 (1997).
[53] Mandelbrot, B.B. and Zarnfaller, F. : Five place tables of certain stable
distributions, Technical Report RC-421, IBM Thomas J. Watson Research
Center, Yorktown Heights, New York, Dec 31, 1959.
[54] Mandelbrot, B.B. : The Fractal Geometry of Nature, Freeman, San
Francisco 1982.
[55] Mandelbrot, B.B. : Fractals and Scaling in Finance, Springer-Verlag,
New York 1997.
[56] Mantegna, R.N. and Stanley, H.E. : An Introduction to Econophysics,
Cambridge University Press, Cambridge 2000.
[57] Mathai, A.M. and Saxena, R.K. : The H-function with Applications in
Statistics and Other Disciplines, Wiley Eastern Ltd, New Delhi 1978.
[58] Mikusiński, J. : On the function whose Laplace transform is
exp (−sαλ) , Studia Math., 18 191-198 (1959).
[59] Miller, K.S. and Ross, B. : An Introduction to the Fractional Calculus
and Fractional Differential Equations, Wiley, New York 1993.
[60] Oldham, K.B. and Spanier, J. : The Fractional Calculus, Academic
Press, New York 1974.
[61] Pincherle, S. : Sulle funzioni ipergeometriche generalizzate, Atti R.
Accademia Lincei, Rend. Cl. Sci. Fis. Mat. Nat. (4), 4, 694-700, 792-799
(1888).
[62] Podlubny, I. : Fractional Differential Equations, Academic Press, San
Diego 1999. [Mathematics in Science and Engineering, Vol. 198]
[63] Pollard, H. : The representation of exp (−xλ) as a Laplace integral,
Bull. Amer. Math. Soc., 52, 908-910 (1946).
[64] Prüss, J. : Evolutionary Integral Equations and Applications,
Birkhäuser, Basel 1993.
[65] Rubin, B. : Fractional Integrals and Potentials, Pitman Monographs
and Surveys in Pure and Appl. Mathematics # 82, Longman, London
1996.
[66] Saichev, A.I and Zaslavsky, G.M. : Fractional kinetic equations:
solutions and applications, Chaos, 7 753-764 (1997).
[67] Samko, S.G., Kilbas, A.A. and Marichev, O.I. : Fractional Integrals and
Derivatives, Theory and Applications, Gordon and Breach, Amsterdam
1993. [Engl. Transl. from the Russian edition, 1987]
[68] Samorodnitsky, G. and Taqqu, M.S. : Stable non-Gaussian Random
Processes, Chapman & Hall, New York 1994.
[69] Schneider, W.R. : Stable distributions: Fox function representation and
generalization, in S. Albeverio, G. Casati and D. Merlini (Eds), Stochastic
Processes in Classical and Quantum Systems, Lecture Notes in Physics #
262, Springer Verlag, Berlin 1986, 497-511.
[70] Schneider, W.R. and Wyss, W. : Fractional diffusion and wave
equations, J. Math. Phys., 30 134-144 (1989).
[71] Srivastava, H.M., Gupta, K.C. and Goyal, S.P. : The H-Functions of
One and Two Variables with Applications, South Asian Publ., New Delhi
1982.
[72] H. Takayasu, H. : Fractals in the Physical Sciences, Manchester Univ.
Press, Manchester and New York 1990.
[73] Uchaikin, V.V. and Zolotarev, V.M. : Chance and Stability. Stable
Distributions and their Applications, VSP, Utrecht 1999. [Series ”Modern
Probability and Statistics”, No 3]
[74] Zolotarev, V.M. : One-dimensional stable distributions, Amer. Math.
Soc., Providence, R.I. 1986. [English Transl. from the Russian edition,
1982]
	Introduction
	The standard diffusion equation
	The time-fractional diffusion equation
	The Cauchy problem for the time-fractional diffusion equation
	The Signalling problem for the time-fractional diffusion equation
	The Cauchy problem for the symmetric space-fractional diffusion equation
	Conclusions
