Zeitschrift für Analysis und ihre Anwendungen

Journal for Analysis and its Applications Volume 18 (1999), No. 2, 231–246

Approximation of Lévy-Feller Diffusion by Random Walk

R. Gorenflo and F. Mainardi

Dedicated to Prof. L. von Wolfersdorf on occasion of his 65th birthday

Abstract. After an outline of W. Feller’s inversion of the (later so called) Feller potential operators and the presentation of the semigroups thus generated, we interpret the two-level difference scheme resulting from the Grünwald-Letnikov discretization of fractional derivatives as a random walk model discrete in space and time. We show that by properly scaled transition to vanishing space and time steps this model converges to the continuous Markov process that we view as a generalized diffusion process. By re-interpretation of the proof we get a discrete probability distribution that lies in the domain of attraction of the corresponding stable Lévy distribution. By letting only the time-step tend to zero we get a random walk model discrete in space but continuous in time. Finally, we present a random walk model for the timeparametrized Cauchy probability density. Keywords: Stable probability distributions, Riesz-Feller potentials, pseudo-differential equa-

tions, Markov processes, random walks

AMS subject classification: 26 A 33, 44 A 20, 45 K 05, 60 E 07, 60 J 15, 60 J 60

1. Introduction

Let

0 < α ≤ 2  and  |θ| ≤ \{ \begin{matrix}α & if 0 < α ≤ 1 \\ 2 − α & if 1 < α ≤ 2 \end{matrix}  (1.1)

(θ real) and denote by p_{α}(x; θ) for x ∈ ℝ the stable probability density whose characteristic function (Fourier transform) is

p̂_{α}(κ; θ) = exp ( − |κ|^{α}e^{i(sign κ) \frac{θπ}{2}} )  (κ ∈ ℝ)  (1.2)

(see, e.g., [4], [17], [19] for the general theory of stable probability distributions). In particular we recommend [4], Feller’s parametrization being close to ours. For a generic function f on ℝ we denote by f̂ its Fourier transform

f̂(κ) = ∫_{−∞}^{+∞}e^{iκx}f(x) dx  (κ ∈ ℝ)  (1.3)

R. Gorenflo: Free Univ. of Berlin, Dept. Math. & Comp. Sci., Arnimallee 2-6, D-14195 Berlin e-mail: gorenflo@math.fu-berlin.de F. Mainardi: Univ. of Bologna, Dept. Phys., Via Irnerio 46, I-40126 Bologna, Italy e-mail: mainardi@bo.infn.it

ISSN 0232-2064 / $ 2.50 \begin{matrix}c \\ ◯ \end{matrix} Heldermann Verlag Berlin

and we then have, in the case of ∫_{−∞}^{+∞}|f̂(κ)| dκ < ∞,

f(x) = \frac{1}{2π} ∫_{−∞}^{+∞}e^{−iκx}f̂(κ) dκ (x ∈ ℝ).

For t > 0 we rescale p_{α} by the similarity variable x t^{− \frac{1}{α}} to obtain

g_{α}(x, t; θ) := t^{− \frac{1}{α}} p_{α}(xt^{− \frac{1}{α}}; θ)  (x ∈ ℝ, t > 0).  (1.4)

This function g_{α}(·, t; θ) again is a stable probability density, and by interpreting x as space and t as time variable we have in g_{α} a description of a Markov process that can be considered as a generalized diffusion process. In fact, we have in

g_{2}(x, t; 0) = \frac{1}{2\sqrt{π}} t^{− \frac{1}{2}} exp ( − \frac{x^{2}}{4t} )

the classical Gauss process and in

g_{1}(x, t; 0) = \frac{1}{π} \frac{t}{x^{2} + t^{2}} ,

the Cauchy process. For a few other pairs (α, θ) leading to elementary or well-investi-gated special functions, see [19]. A general representation of all stable probability densities in terms of Fox H functions has been only recently achieved (see [18]). The Fourier transform of g_{α} being

ĝ_{α}(κ, t; θ) = exp ( − t|κ|^{α}e^{i(sign κ) \frac{θπ}{2}} )  (κ ∈ ℝ)  (1.5)

we recognize g_{α}(x, t; θ) as the fundamental solution (Green function for the Cauchy problem) of the pseudo-differential equation

\frac{∂u(x, t)}{∂t} = D_{θ}^{α}u(x, t)  (x ∈ ℝ, t > 0)  (1.6)

where the pseudo-differential operator D_{θ}^{α} has the symbol −|κ|^{α}e^{i(sign κ) \frac{θπ}{2}} . For initial values

u(x, 0) = f(x)  (x ∈ ℝ, f ∈ L_{1}(ℝ))  (1.7)

we then have as solution to (1.6)

u(x, t) = ∫_{−∞}^{+∞}g_{α}(x − ξ, t; θ)f(ξ) dξ  (1.8)

and for all t > 0 then

u(·, t) ∈ C^{∞} ∩ L_{1}(ℝ)  and ∫_{−∞}^{+∞}u(x, t) dx = ∫_{−∞}^{+∞}f(x) dx.

William Feller in his pioneering paper [3] has shown that the pseudo-differential operator D_{θ}^{α} can be viewed as the operator inverse to the Feller potential operator (the name ”Feller potential” is used in [16]) which is a linear combination of two Weyl integrals. Honouring both Lévy and Feller for their essential contributions [11], [12] and [3] we call the process described by (1.6) Lévy-Feller diffusion.

We now give, in our notation, a formal account of the essentials of Feller’s theory (for more details see [8]). With the Weyl integrals

\begin{matrix}(I_{+}^{α}φ)(x) & = \frac{1}{Γ(α)} ∫_{−∞}^{x} (x − ξ)^{α−1}φ(ξ) dξ \\ (I_{−}^{α}φ)(x) & = \frac{1}{Γ(α)} ∫_{x}^{+∞}(ξ − x)^{α−1}φ(ξ) dξ \end{matrix} \}  (x ∈ ℝ)  (1.9)

and (for 0 < α < 2 but α ≠ 1) the coefficients

c_{+} = c_{+}(α; θ) = \frac{sin(\frac{π}{2} (α − θ))}{sin(απ)},  c_{−} = c_{−}(α; θ) = \frac{sin(\frac{π}{2} (α + θ))}{sin(απ)}  (1.10)

and (by passing to the limit α = 2)

c_{+}(2, 0) = c_{−}(2, 0) = −\frac{1}{2},  (1.11)

the Feller potentials are given as

(I_{θ}^{α}φ)(x) = c_{−}(α, θ)(I_{+}^{α}φ)(x) + c_{+}(α, θ)(I_{−}^{α}φ)(x) .  (1.12)

Note that in accordance with [16] we omit the singular case α = 1.

Feller [3] has shown the operator I_{θ}^{α} to possess the semigroup property

I_{θ}^{α}I_{θ}^{β} = I_{θ}^{α+β} for 0 < α, β < 1 with α + β < 1,

and so analytic continuation to negative exponents can be justified to obtain the operator

D_{θ}^{α} := −I_{θ}^{−α} = −{c_{+}(α, θ)I_{+}^{−α} + c_{−}(α, θ)I_{−}^{−α}}  (1.13)

for 0 < α ≤ 2 but α ≠ 1, the parameter θ restricted as in (1.1), with (see [16])

I_{±}^{−α} = \{ \begin{matrix}±\frac{d}{dx}I_{±}^{1−α} & if 0 < α < 1 \\ \frac{d^{2}}{dx^{2}} I_{±}^{2−α} & if 1 < α ≤ 2. \end{matrix}  (1.14)

From [3], equating −\frac{θπ}{2α} to Feller’s parameter δ, we take the symbol of the pseudodifferential operator D_{θ}^{α} as D̂_{θ}^{α} = −|κ|^{α}e^{i(sign κ)θπ/2}. In particular, we have D_{0}^{2} = \frac{d^{2}}{dx^{2}} but D_{0}^{1} ≠ \frac{d}{dx}.

For the rest of this paper we always keep in mind the distinction of the following two cases:

(a) 0 < α < 1 and |θ| ≤ α.
(b) 1 < α ≤ 2 and |θ| ≤ 2 − α.

Henceforth, for ease of notation, we shall omit the arguments of the coefficients c_{+} = c_{+}(α, θ) and c_{−} = c_{−}(α, θ). We have

c_{±} \{ \begin{matrix}≥ 0 & in the case (a) \\ ≤ 0 & in the case (b) \end{matrix}  (1.15)

and

c_{+} + c_{−} = \begin{matrix}cos & \frac{θπ}{2} \\ cos & \frac{απ}{2} \end{matrix} \{ \begin{matrix}> 0 & in the case (a) \\ < 0 & in the case (b). \end{matrix}  (1.16)

The reader is asked not to worry about the foregoing purely formal description of Feller’s considerations. It will merely serve us as a motivation for constructing a difference scheme via the Grünwald-Letnikov discretization of fractional derivatives, a difference scheme which by interpretation as a random walk model will be shown to converge (in a sense to be specified in Section 3).

2. Random walks, discrete in space and time

In this section we define a random variable Y assuming only integers as values, its probability distribution depending on three parameters α, θ and μ. By aid of this random variable we define a random walk on an equidistant grid {jh| j ∈ ℤ} with a space-step h > 0. We show that after introduction of a time-step τ > 0 this random walk admits an interpretation as an explicit difference scheme for the Cauchy problem (1.6) - (1.7), namely for

\begin{matrix}\frac{∂u(x, t)}{∂t} & = D_{θ}^{α}u(x, t) (x ∈ ℝ, t > 0) \\ u(x, 0) & = f(x). \end{matrix} \}  (2.1)

In the next section we shall show that the probability distribution of the discrete random variable Y belongs to the domain of attraction of the Lévy distribution with the parameters α and θ, proceeding in a way which simultaneously proves ”convergence” of the random walk (if τ = μh^{α} → 0) to the corresponding Lévy-Feller diffusion characterized by (1.4).

Let Y be a random variable assuming its values in ℤ, P (Y = k) = p_{k} for k ∈ ℤ, with probabilities p_{k} defined as follows. With a parameter μ, restricted by

0 < μ ≤ \{ \begin{matrix}\begin{matrix}cos & \frac{απ}{2} \\ cos & \frac{θπ}{2} \end{matrix} & in the case (a) \\ \frac{1}{α} \frac{| cos \frac{απ}{2} |}{cos \frac{θπ}{2}} & in the case (b) \end{matrix}  (2.2)

put in the case (a)

\begin{matrix}p_{0} & = 1 − μ(c_{+} + c_{−}) \\ p_{k} & = (−1)^{k+1}μc_{+}({α \atop k}) \\ p_{−k} & = (−1)^{k+1}μc_{−}({α \atop k}) (k ∈ ℕ) \end{matrix} \}  (2.3)

and in the case (b)

\begin{matrix}p_{0} & = 1 + μα(c_{+} + c_{−}) \\ p_{1} & = −μ[c_{+}({α \atop 2}) + c_{−}], p_{−1} = −μ[c_{−}({α \atop 2}) + c_{+}] \\ p_{k} & = (−1)^{k}μc_{+}({α \atop k + 1}), p_{−k} = (−1)^{k}μc_{−}({α \atop k + 1}) (k ≥ 2). \end{matrix} \}  (2.4)

One sees that all p_{k} ≥ 0, and by rearrangement it turns out that

∑_{k∈ℤ} p_{k} = 1 − μ(c_{+} + c_{−}) ∑_{j=0}^{+∞}(−1)^{j}({α \atop j}) = 1 − 0.

Remark 2.1. It is worthwhile here to observe the fact which will also be useful in Section 3 that for all α > 0 the series ∑_{k=0}^{∞}(−1)^{k}\binom{α}{k}z^{k} for (1 − z)^{α} converges absolutely and uniformly on the closed unit disk |z| ≤ 1, due to the asymptotics |\binom{α}{k}| ∼ Γ(α + 1)\frac{| sin(απ)|}{π}k^{−(α+1)} for k → ∞, valid for non-integer α > 0. This asymptotics can be deduced by use of the reflection formula for the gamma function and Stirling’s asymptotics.

We obtain a random walk on the grid {jh| j ∈ ℤ} starting at the point 0, by defining random variables

S_{n} = hY_{1} + hY_{2} + … + hY_{n}  (n ∈ ℕ)  (2.5)

with the Y_{j} as independent identically distributed random variables, all having the same probability distribution as the random variable Y.

Let us write our random walk in an alternative way. Discretizing the space variable x and the time variable t by grid points x_{j} = jh and instants t_{n} = nτ, with h > 0, τ > 0, j ∈ ℤ, n ∈ ℕ_{0} and denoting by y_{j}(t_{n}) the probability of sojourn of the random walker in point x_{j} at instant t_{n}, the recursion S_{n+1} = S_{n} + hY_{n+1} (following from (2.5)) means

y_{j}(t_{n+1}) = ∑_{k∈ℤ} p_{k}y_{j−k}(t_{n})  (j ∈ ℤ, n ∈ ℕ_{0}),  (2.6)

and the random walker starting at point x_{0} = 0 means y_{0}(0) = 1 and y_{j}(0) = 0 for j ≠ 0. However, in the recursion scheme (2.6) it is legitimate to use a more general initial sojourn probability distribution {y_{j}(0)| j ∈ ℤ}. There is yet another possible interpretation of (2.6), namely as a redistribution scheme of an extensive quantity (e.g. mass, charge, or may be probability), y_{j}(t_{n}) being imagined as a clump of this extensive quantity, sitting in point x_{j} at instant t_{n}. Then (2.6) is a conservative and non-negativity preserving redistribution scheme. In fact, from all p_{k} ≥ 0 and ∑_{k∈Z} p_{k} = 1 it follows immediately for all n ∈ ℕ that

∑_{j∈ℤ} y_{j}(t_{n}) = ∑_{j∈ℤ} y_{j}(0)  if ∑_{j∈ℤ} |y_{j}(0)| < ∞
all y_{j}(t_{n}) ≥ 0  if all y_{j}(0) ≥ 0.

Such redistribution schemes have been shown to be useful for discretization of diffusion processes modelled by second order linear parabolic differential equations (see, e.g., [6], [7], [9]) as they discretely imitate essential properties of the continuous process.

To come nearer to the Cauchy problem (2.1) we relate the time step τ to the space step h by the scaling relation

τ = μh^{α}  (2.7)

and remark that the y_{j}(t_{n}) are then intended as approximations to

∫_{x_{j}− \frac{h}{2}}^{x_{j}+ \frac{h}{2}}u(x, t_{n}) dx

which, if u(·, t_{n}) is continuous, is also ≈ hu(x_{j}, t_{n}). It is again a matter of rearrangement to show that (2.6) is equivalent to the explicit difference scheme

\frac{y_{j}(t_{n+1}) − y_{j}(t_{n})}{τ} = _{h}D_{θ}^{α}y_{j}(t_{n})  (2.8)

where (in analogy to (1.13)) _{h}D_{θ}^{α} = −{c_{+} _{h}I_{+}^{−α} + c_{−} _{h}I_{−}^{−α}} with the Grünwald-Letnikov discretization (see [16]) of the fractional derivatives in the form

hI_{±}^{−α}y_{j} = \{ \begin{matrix}h^{−α} ∑_{k=0}^{∞}(−1)^{k}({α \atop k})y_{j∓k} & in the case (a) \\ h^{−α} ∑_{k=0}^{∞}(−1)^{k}({α \atop k})y_{j±1∓k} & in the case (b). \end{matrix}  (2.9)

Notice the shift of index in the case (b) which among other things has the effect that in the special case α = 2 (the classical diffusion equation) we obtain the standard symmetric three-point difference scheme. For more details and discussions see [8].

Instead of trying to work out a convergence proof for the difference scheme (2.8), thereby using the Lax-Richtmyer theory of consistency, stability and convergence (in effect the Lax equivalence theorem, see [5] or [14]) we prefer to present in the next section a proof in the true spirit of random walks. We leave the numerical analysis aspect to a forthcoming paper.

3. Convergence and domain of attraction

We will show that for fixed t = nτ > 0 the discrete distribution of the sojourn probabilities y_{j}(t_{n}) (j ∈ ℤ) with initial condition y_{j}(0) = δ_{j0} (Kronecker symbol) converges completely to the probability distribution with density

g_{α}(x, t; θ) = t^{− \frac{1}{α}} p_{α}(xt^{− \frac{1}{α}}; θ)  (x ∈ ℝ)  (3.1)

as n → +∞. Let us remind that this probability distribution has the characteristic function

ĝ_{α}(κ, t; θ) = ∫_{−∞}^{+∞}g_{α}(x, t; θ) e^{iκx}dx = exp ( − t|κ|^{α}e^{i(sign κ) \frac{θπ}{2}} ).  (3.2)

To avoid confusion of language one meets in probability theory let us agree to use the terminology adopted in [10]. From this source we take Definitions 3.1 - 3.4, Remark 3.1 and Theorem 3.1.

Definition 3.1. Let (F_{n}) be a sequence of uniformly bounded, non-decreasing right-continuous functions defined on ℝ. We say that F_{n} converges weakly to a bounded non-decreasing right-continuous function F on ℝ if F_{n}(x) → F (x) at all continuity points of F. In this case we write F_{n}→^{w}F.

Definition 3.2. Let (F_{n}) be as in Definition 3.1. Then (F_{n}) is said to converge completely to F if

(i) F_{n}→^{w}F and

(ii) F_{n}(∓∞) → F (∓∞) as n → ∞.

In this case we write F_{n}→^{c}F.

Theorem 3.1 (Continuity theorem). Let (F_{n}) be a sequence of probability distribution functions, and let (φ_{n}) be the sequence of the corresponding characteristic functions,

φ_{n}(κ) = ∫_{−∞}^{+∞}e^{iκx} dF_{n}(x) (κ ∈ ℝ).

Then (F_{n})converges completely to a probability distribution function F if and only if φ_{n}(κ) → φ(κ) for all κ ∈ ℝ as n → ∞, where φ(κ) is continuous at κ = 0. In this case the limit function φ is the characteristic function of the limit distribution function F,

φ(κ) = ∫_{−∞}^{+∞}e^{iκx} dF (x) (κ ∈ ℝ).

Definition 3.3. In the cases where the functions F_{n} and F are probability distribution functions such that F_{n}→^{c}F, let X_{n} and X be random variables corresponding to F_{n} and F, respectively. Then we say that X_{n} converges in law to X.

Definition 3.4. Let (X_{n}) be a sequence of independent identically distributed random variables with common probability distribution function F. Suppose there exist sequences (a_{n}) and (b_{n}) of constants, with b_{n} > 0, such that the sequence of sums b_{n}^{−1} ∑_{k=1}^{n} X_{k}−a_{n} converges in law to some random variable with probability distribution function G. Then we say that F is attracted to G. The set of all probability distribution functions attracted to G is called the domain of attraction of the distribution function G.

Remark 3.1. A stable probability distribution is characterized by having a domain of attraction.

Let us now state, with the notations of Sections 1 and 2, our

Theorem 3.2. Let the independent identically distributed random variables Y_{1}, Y_{2}, Y_{3}, … have the common probability distribution function F of the random variable Y with P (Y = k) = p_{k}, k ∈ ℤ, and p_{k} given by (2.3), (2.4), respectively. Let X(t) with t > 0 be the random variable with probability density g_{α}(x, t; θ) and let G_{α}(·, t; θ) be the corresponding distribution function,

G_{α}(x, t; θ) = ∫_{−∞}^{x} g_{α}(ξ, t; θ) dξ.

Then F is attracted by G_{α}(·, t; θ), indeed: for n → ∞ the distribution function of the random variable

X_{n} = (\frac{t}{μn})^{\frac{1}{α}} \{Y_{1} + Y_{2} + … + Y_{n}\}  (3.3)

converges completely to G_{α}(·, t; θ), the distribution function of the random variable X(t).

Proof. Using the scaling relation (2.7), namely τ = μh^{α}, and the substitution t = nτ of time, we get

h = (t/(μn))^{\frac{1}{α}} ,  (3.4)

and comparing (3.3) and (2.5) we see that X_{n} = S_{n}, the random variable taking values in the grid {jh| j ∈ ℤ} at the fixed instant t_{n} = nτ = t. In view of Theorem 3.1 and the fact that ĝ_{α}(κ, t; θ) is continuous at κ = 0, it only remains to prove that the characteristic function of the sojourn probabilities y_{j}(t_{n}), namely the function

ŷ(κ, t; h) = ∑_{j∈ℤ} y_{j}(t_{n}) e^{ijκh},  with t_{n} = t,  (3.5)

tends for all κ ∈ ℝ (as h → 0) to ĝ_{α}(κ, t; θ) = exp ( − t|κ|^{α}e^{i(sign κ) \frac{θπ}{2}} ). Let us calculate ŷ(κ, t; h) for ease of notation via the generating functions

p̃(z) = ∑_{j∈ℤ} p_{j}z^{j}  and  ỹ(z, t_{n}) = ∑_{j∈ℤ} y_{j}(t_{n}) z^{j}  (3.6)

of the transition probabilities and the sojourn probabilities. The series in (3.6) converge absolutely and uniformly on the periphery |z| = 1 of the unit circle, representing there a continuous function, and due to the fact that the random walk occurs on the grid {jh| j ∈ ℤ} change to characteristic functions p̂(κ) and ŷ(κ, t_{n}) is accomplished via z = e^{iκh} (κ ∈ ℝ). Using the binomial series for (1 − z)^{α}, absolutely convergent on |z| = 1 if α > 0 (see Remark 2.1), we readily verify the identities

p̃(z) = \{ \begin{matrix}1 − μ\{c_{+}(1 − z)^{α} + c_{−}(1 − z^{−1})^{α}\} & in the case (a) \\ 1 − μ\{c_{+}z^{−1}(1 − z)^{α} + c_{−}z(1 − z^{−1})^{α}\} & in the case (b). \end{matrix}  (3.7)

From the discrete convolution (2.6) we deduce ỹ(z, t_{n}) = ỹ(z, 0) (p̃(z))^{n}, and the special initial condition y_{j}(0) = δ_{j0} for j ∈ ℤ gives ỹ(z, 0) ≡ 1, hence

ỹ(z, t_{n}) = (p̃(z))^{n}.  (3.8)

In view of (3.4), (3.8) and the fixation t = t_{n} = nτ we have to show that, with z = e^{iκh},

(p̃(z))^{n} → exp ( − t|κ|^{α}e^{i(sign κ) \frac{θπ}{2}} ) = ĝ_{α}(κ, t; θ)

as n → ∞. More clearly, using (2.7) and t = t_{n} = nτ, we have to show that the function

ŷ(κ, t; h) = (p̃(e^{iκh}))^{\frac{t}{μh^{α}}}  (κ ∈ ℝ)  (3.9)

has the property

lim_{h→0} ŷ(κ, t; h) = ĝ_{α}(κ, t; θ).  (3.10)

Let us first treat the case (a): 0 < α < 1 and |θ| ≤ α. Then κ = |κ| sign κ and

p̃(e^{iκh}) = 1 − μ{c_{+}(1 − e^{i|κ|h sign κ})^{α} + c_{−}(1 − e^{−i|κ|h sign κ})^{α}}.

We see that p̃(e^{i0h}) = 1, whereas we can get the result for κ < 0 by complex conjugation of that for κ > 0. So, for notational ease, we treat in detail the case κ > 0. In this case

p̃(e^{iκh}) = 1 − μ{c_{+}(1 − e^{iκh})^{α} + c_{−}(1 − e^{−iκh})^{α}}  (3.11)

and for small h by Taylor

(1 − e^{iκh})^{α} = (−iκh + O(h^{2}))^{α}
= (−i)^{α}(κh)^{α}(1 + O(h))^{α}
= e^{−i \frac{πα}{2}} (κh)^{α} + O(h^{α+1})

and

(1 − e^{−iκh})^{α} = e^{i \frac{πα}{2}} (κh)^{α} + O(h^{α+1}).

Inserting this into (3.7) we find

p̃(e^{iκh}) = 1 − μκ^{α}h^{α}\{c_{+}e^{−i \frac{πα}{2}} + c_{−}e^{i \frac{πα}{2}} \} + O(h^{α+1}).

By use of (1.10) for c_{+} and c_{−} and the complex representation of sin \frac{π(α∓θ)}{2} a straightforward calculation yields for (fixed) κ > 0

p̃(e^{iκh}) = 1 − μκ^{α}h^{α}e^{i \frac{πθ}{2}} + O(h^{α+1}) = 1 − μ|κ|^{α}h^{α}e^{i \frac{πθ}{2}} + O(h^{α+1}),

for κ < 0 by complex conjugation

p̃(e^{iκh}) = 1 − μ|κ|^{α}h^{α}e^{−i \frac{πθ}{2}} + O(h^{α+1}).

So, finally, for all κ ∈ ℝ,

p̃(e^{iκh}) = 1 − μ|κ|^{α}h^{α}e^{i(sign κ) \frac{πθ}{2}} + O(h^{α+1})  (3.12)

and by (3.9)

log ŷ(κ, t; h) = \frac{t}{μh^{α}} \{ − μ|κ|^{α}h^{α}e^{i(sign κ) \frac{θπ}{2}} + O(h^{α+1})\}
= −t|κ|^{α}e^{i(sign κ) \frac{θπ}{2}} + O(h),

hence, as desired, (3.10).

In the case (b): 1 < α ≤ 2 and |θ| ≤ 2 − α, we have by (3.7)

p̃(e^{iκh}) = 1 − μ{c_{+}e^{−iκh}(1 − e^{iκh})^{α} + c_{−}e^{iκh}(1 − e^{−iκh})^{α}}

and in comparison to the case (a) we have because of e^{∓iκh} = 1 + O(h) within {…} the additional asymptotic term

O(h)(1 − e^{iκh})^{α} + O(h)(1 − e^{−iκh})^{α} = O(h^{α+1}),

hence again (3.12) for all κ ∈ ℝ, and again we arrive at (3.10)

It is instructive to take a look at the very special case α = 2 and θ = 0 (the classical diffusion equation u_{t} = u_{xx}). In this case

\begin{aligned}p̃(z) & = 1 + μ\{\frac{1}{z} − 2 + z\} \\ p̃(e^{iκh}) & = 1 + μ\{e^{i \frac{κh}{2}} − e^{−i \frac{κh}{2}} \}^{2} = 1 − 4μ sin^{2} \frac{κh}{2} \end{aligned}

and one finds (p̃(e^{iκh}))^{\frac{t}{μh^{2}}} → exp(−tκ^{2}) as h → 0.

4. A random walk model, discrete in space, continuous in time

Consider the difference scheme (2.8) which is equivalent to the redistribution scheme (or random walk model) (2.6) with the coefficients given by (2.3) or (2.4), respectively. By sending the parameter μ → 0 (letting the time step τ tend to 0) we obtain an infinite system of ordinary differential equations

\begin{matrix}y′_{j}(t) & = _{h}D_{θ}^{α}y_{j}(t) \\ y_{j}(0) & given \end{matrix} \}  (j ∈ ℤ)  (4.1)

describing a time-continuous redistribution scheme over the grid {jh| j ∈ ℤ} in time t ≥ 0 of the form

y′_{j}(t) = ∑_{k∈ℤ} q_{k}y_{j−k}(t)  (j ∈ ℤ).  (4.2)

Interpreting y_{j}(t) as a clump of an extensive quantity sitting in point x_{j} = jh at instant t we have, for (4.2) to describe such a redistribution scheme, the balancing conditions

\begin{matrix}q_{0} & < 0 \\ q_{k} & ≥ 0 (0 ≠ k ∈ ℤ) \end{matrix} \}  (4.3)

∑_{k∈ℤ} q_{k} = 0.  (4.4)

In analogy to our redistribution scheme (2.6) of Section 2 system (4.2) also is conservative and non-negativity preserving. In fact, it can be shown (we leave this as an exercise to the reader) that system (4.2) under conditions (4.3) - (4.4) is uniquely solvable if

∑_{j∈ℤ} |y_{j}(0)| < ∞  (4.5)

and that then ∑_{j∈ℤ} |y_{j}(t)| < ∞ for all t > 0. It can further be shown that then

∑_{j∈ℤ} y′_{j}(t) = 0, hence ∑_{j∈ℤ} y_{j}(t) ≡ ∑_{j∈ℤ} y_{j}(0)

for all t > 0. If furthermore y_{j}(0) ≥ 0 (j ∈ ℤ), then y_{j}(t) ≥ 0 (j ∈ ℤ) for all t > 0.

The interpretation of (4.2) with (4.3) and (4.4) as a redistribution scheme means: |q_{0}|y_{j}(t) is the rate of outflow from the point x_{j} = jh being transferred to other points, and this must equal the sum of the rates q_{k}y_{j}(t), received by the points x_{j+k} (k ≠ 0).

Using in (4.1) again the Grünwald-Letnikov discretization (2.9) we find the following. In the case (a) 0 < α < 1 and |θ| ≤ α:

\begin{matrix}q_{0} & = −h^{−α}(c_{+} + c_{−}) = −h^{−α} \begin{matrix}cos & \frac{θπ}{2} \\ cos & \frac{απ}{2} \end{matrix} \\ q_{k} & = h^{−α}(−1)^{k+1}c_{+}({α \atop k}) (k ∈ ℕ) \\ q_{−k} & = h^{−α}(−1)^{k+1}c_{−}({α \atop k}) (k ∈ ℕ). \end{matrix} \}  (4.6)

In the case (b) 1 < α ≤ 2 and |θ| ≤ 2 − α:

\begin{matrix}q_{0} & = h^{−α}(c_{+} + c_{−})α = −αh^{−α} \frac{cos \frac{θπ}{2}}{| cos \frac{απ}{2} |} \\ q_{1} & = −h^{−α}[c_{+}({α \atop 2}) + c_{−}], q_{−1} = −h^{−α}[c_{−}({α \atop 2}) + c_{+}] \\ q_{k} & = (−1)^{k}h^{−α} c_{+}({α \atop k + 1}) (k ≥ 2) \\ q_{−k} & = (−1)^{k}h^{−α} c_{−}({α \atop k + 1}) (k ≥ 2). \end{matrix} \}  (4.7)

By playing again with infinite sums of binomial coefficients it is readily verified that (4.3) and (4.4) are fulfilled.

For solving system (4.2) with initial values y_{j}(0) (j ∈ ℤ) with ∑_{j} |y_{j}(0)| < ∞ and q_{k} given by (4.6) - (4.7), we apply the method of generating functions. The series

q̃(z) = ∑_{k∈ℤ} q_{k}z^{k}  and  ỹ(z, t) = ∑_{k∈ℤ} y_{k}(t)z^{k}  (4.8)

converge absolutely and uniformly on the periphery |z| = 1 of the unit circle, and system (4.2) is then (with |z| = 1) equivalent to

\frac{∂ỹ(z, t)}{∂t} = ∑_{j∈ℤ} y′_{j}(t) z^{j} = ∑_{j∈ℤ} (∑_{k∈ℤ} q_{k} y_{j−k}(t)) z^{j},

hence we have the z-parameterized ordinary differential equation problem

\begin{matrix}\frac{∂ỹ(z, t)}{∂t} & = q̃(z)ỹ(z, t) (t ≥ 0, |z| = 1) \\ ỹ(z, 0) & = ∑_{k∈ℤ} y_{k}(0) z^{k}. \end{matrix} \}

The solution is ỹ(z, t) = ỹ(z, 0) e^{tq̃(z)}, or simply

ỹ(z, t) = e^{tq̃(z)}  (4.9)

in the special case y_{j}(0) = δ_{j0} for j ∈ ℤ which means ỹ(z, 0) ≡ 1.

By inspection (using the binomial series) we see that

q̃(z) = \{ \begin{matrix}−h^{−α}\{c_{+}(1 − z)^{α} + c_{−}(1 − z^{−1})^{α}\} & in the case (a) \\ −h^{−α}\{c_{+}z^{−1}(1 − z)^{α} + c_{−}z(1 − z^{−1})^{α}\} & in the case (b). \end{matrix}

Changing to characteristic functions via z = e^{iκh} (κ ∈ ℝ), we take from our calculations of Section 3 for small h

q̃(e^{iκh}) = −|κ|^{α}e^{i(sign κ) \frac{πθ}{2}} + O(h).

Then with (4.9) we get for ŷ(κ, t; h) := ỹ(e^{iκh}, t) in analogy to (3.10) the limit relation

lim_{h→0} ŷ(κ, t; h) = exp ( − t|κ|^{α}e^{i(sign κ) \frac{θπ}{2}} ) = ĝ_{α}(κ, t; θ) (κ ∈ ℝ).

We have interpreted (4.2) as a time-continuous redistribution scheme. We can interpret it probabilistically as a random walk model discrete in space (over the grid {jh| j ∈ ℤ}), but continuous in time. At any instant t of time the random walker can jump to another grid point. After arriving at a point x_{m} he will remain sitting there for a random time interval whose length is exponentially distributed. More precisely: when we know that at instant t^{∗} he is sitting at point x_{m}, then the conditional sojourn probabilities for sitting at points x_{j} are η_{j}(t^{∗}) = δ_{mj} (j ∈ ℤ) and (4.2) gives by re-conditioning the equation

η′_{m}(t) = −|q_{0}| η_{m}(t),

for the time interval [t^{∗}, t^{∗}+t̃] of sojourn at x_{m}. Its solution is η_{m}(t) = 1 − e^{−|q_{0}|(t−t^{∗})} (t ≥ t^{∗}), from which we deduce that the time t̃the wanderer remains sitting at any point x_{m} is exponentially distributed with parameter |q_{0}|. Hence, the random walker, after arriving at point x_{m} sits there for a random time interval of length t̃ and then jumps to another point x_{j} in instant t = t^{∗} + t̃. The conditional probability of jumping to the point x_{j} (with j ≠ m) is then given as \frac{q_{j−m}}{|q_{0}|}. For general information on time-continuous discrete Markov processes we refer the reader to [15]. It should finally be remarked that the conditional density η_{m}(t) (t ≥ t^{∗}) can also be obtained in the limit of τ → 0 from the conditional geometric probability distribution relevant in the random walk model (2.6) with the transition probabilities of (2.3) - (2.4) and the scaling condition (2.7).

We can now state a theorem analogous to Theorem 3.2, namely

Theorem 4.1. Let a random walker start in point 0 at instant t = 0 and jump over the grid points jh (j ∈ ℤ) with h > 0, the probabilities y_{j}(t) of sojourn in point jh at instant t ≥ 0 evolving according to (4.2) with y_{j}(0) = δ_{j0}. Then for fixed t > 0 the distribution function G_{α}(·, t; θ; h) given as G_{α}(x, t; θ; h) = ∑_{kh≤x} y_{k}(t) converges

completely to G_{α}(·, t; θ) as h → 0.

5. A random walk model for the Cauchy process

For completeness we present a random walk model, discrete in space and time, for the omitted case α = 1, namely for the Cauchy process (α = 1 and θ = 0). This model cannot be obtained via the Grünwald-Letnikov approach, neither directly nor by a passage to the limit α → 1. We have, with the notations of Sections 1 - 3 the process

g_{1}(x, t; 0) = \frac{1}{π} \frac{t}{x^{2} + t^{2}} = t^{−1}p_{1}(\frac{x}{t}; 0) (x ∈ ℝ, t > 0)

with the Cauchy probability density

p_{1}(x, 0) = \frac{1}{π(x^{2} + 1)} (x ∈ ℝ).

The corresponding characteristic functions are

p̂_{1}(κ; 0) = e^{−|κ|}  and  ĝ_{1}(κ, t; 0) = e^{−t|κ|}.

Let Y be a random variable assuming its values in ℤ with P (Y = k) = p_{k} for k ∈ ℤ defined as follows. With a parameter μ restricted to 0 < μ ≤ \frac{π}{2} put

\begin{matrix}p_{0} & = 1 − \frac{2μ}{π} \\ p_{k} & = \frac{μ}{π|k|(|k| + 1)} (0 ≠ k ∈ ℤ). \end{matrix} \}  (5.1)

Then all p_{k} ≥ 0 and ∑_{k∈ℤ} p_{k} = 1. Indeed,

∑_{k=1}^{∞} \frac{1}{k(k + 1)} = ∑_{k=1}^{∞} (\frac{1}{k} − \frac{1}{k + 1}) = 1.

Proceeding as in Sections 2 and 3 we produce a random walk on the grid {jh| j ∈ ℤ}, with h > 0, letting the walker start in point 0 at instant 0. We define

S_{n} = hY_{1} + hY_{2} + … + hY_{n} (n ∈ N)

where the Y_{k} are independent identically distributed random variables all having their distribution common with the random variable Y. By allowing the walker to jump in instants t_{n} = nτ (τ > 0, n ∈ N) and using S_{n} as the random variable assuming values in the grid {jh| j ∈ ℤ} we have a random walk. We relate the steps h and τ of space and time by τ = μh and again define generating functions

p̃(z) = ∑_{j∈ℤ} p_{j}z^{j}  and  ỹ(z, t_{n}) = ∑_{j∈ℤ} y_{j}(t_{n}) z^{j}  (5.2)

with y_{j}(t_{n}) as probability of sojourn in point x_{j} = jh at instant t_{n} = nτ. Then ỹ(z, t_{n}) = (p̃(z))^{n}, and in view of \frac{t}{τ} = \frac{t}{μh} our aim is to show that for all κ ∈ ℝ and t > 0 the limit relation

p̃(e^{iκh})^{\frac{t}{μh}} → e^{−t|κ|} as h → 0 (5.3)

holds. This limit relation then implies that the distribution of the random variable Y lies in the domain of attraction of the Cauchy distribution, more precisely that the distribution of the random variable

X_{n} = \frac{t}{μn}\{Y_{1} + Y_{2} + … + Y_{n}\}

converges completely to the Cauchy distribution G_{1}(·, t; 0) with

G_{1}(x, t; 0) = ∫_{−∞}^{x} \frac{t}{π(t^{2} + ξ^{2})} dξ = \frac{1}{2} + \frac{1}{π} arctan \frac{x}{t}.

Let us prove (5.3). We observe that the series ∑_{j∈ℤ} p_{j} z^{j} converges absolutely and uniformly on the periphery |z| = 1 of the unit circle, hence represents there a continuous function so that the following calculations are legitimate. In fact,

∑_{k=1}^{∞} \frac{z^{k}}{k(k + 1)} = 1 − (1 − z^{−1}) log(1 − z) for |z| ≤ 1.

We tacitly take the limit 1 for z = 1. Then

p̃(z) = 1 − \frac{μ}{π}\{(1 − z^{−1}) log(1 − z) + (1 − z) log(1 − z^{−1})\}.  (5.4)

Passing to the characteristic function via z = e^{iκh} (κ ∈ ℝ) we get

p̃(e^{iκh}) = 1 − \frac{2μ}{π} \{(1 − cos(κh)) log |1 − e^{iκh}| − sin(κh) arctan \frac{sin(κh)}{1 − cos(κh)}\} .

Using lim_{u→±∞} arctan u = ±\frac{π}{2} we obtain for h → 0 the asymptotics

p̃(e^{iκh}) = 1 − \frac{2μ}{π} \{|κ|h\frac{π}{2} + o(|κ|h)\} = 1 − μ|κ|h + o(h)

which implies (5.3).

Approximation of Lévy-Feller Diffusion 245

6. Concluding remarks

It is instructive to observe that in view of Theorem 3.1 our proof of Theorem 3.2 can be re-interpreted as a proof of existence of Lévy’s stable distributions (for α ≠ 1). In fact, assuming to be ignorant of these we can find them as limiting distributions by sending n → ∞ in (3.3). And gratis (the discrete probabilities being all non-negative cannot become negative in the limit) we get that the limiting densities are everywhere non-negative, for all values of the parameter α between 0 and 2, with omission of the value 1. For this we actually need neither the theory of the inversion of the Feller potentials nor the method of positive-definite functions. Thus we have an alternative way of solving a problem that surmounted Cauchy’s capabilities [2] who had considered the functions exp(−|κ|^{α}) as candidates of cosine transforms of probability densities but could only prove them to have this property in the special cases α = 1 and α = 2. Lévy in [11], [12] introduced the whole scale of stable densities, Bochner in [1] has given an elegant proof for the full range 0 < α ≤ 2 that the inverse Fourier transforms of the functions exp(−|κ|^{α}) are non-negative, hence probability densities. He used the theory of positive-definite functions that we can avoid. A well readable account of Bochner’s method can be looked up in [13].

Acknowledgements. We are grateful to the Italian Consiglio Nazionale delle Ricerche and to the Research Commission of Free University of Berlin for supporting joint efforts of our research groups in Berlin and Bologna. This paper is one of the fruits of this collaboration. We appreciate the careful work of the referees, in particular the constructive-critical comments of one of them.

References

[1] Bochner, S.: Stable laws of probability and completely monotone functions. Duke Math. J. 3 (1937), 726 – 728.

[2] Cauchy, A.: Calcul des probabilités. Sur les résultats moyens d’observations de même nature, et sur les résultats les plus probables. Sur la probabilité des erreurs qui effectent des résultats moyens d’observations de même nature. Comptes rendus 37 (1853), 198 – 206 and 264 – 272. Reprinted in: Oeuvres Complètes d’Augustin Cauchy, Ser. I, Tome 12, pp. 94 – 114. Paris: Gauthier-Villars 1900.

[3] Feller, W.: On a generalization of Marcel Riesz’ potentials and the semigroups generated by them. Meddelanden Lunds Universitets Matematiska Seminarium (Comm. Sém. Mathém. Université de Lund), Tome suppl. dédié à M. Riesz. Lund 1952, pp. 73 – 81.

[4] Feller, W.: An Introduction to Probability Theory and its Applications, Vol. 2, 2nd ed. New York: Wiley 1971 (1st ed. 1966).

[5] Goldstein, J. A.: Semigroups of Linear Operators and Applications. Oxford - New York: Oxford Univ. Press 1985.

[6] Gorenflo, R.: Nichtnegativitäts- und substanzerhaltende Differenzenschemata für lineare Diffusionsgleichungen. Numer. Math. 14 (1970), 448 – 467.

[7] Gorenflo, R.: Conservative difference schemes for diffusion problems. In: Intern. Ser. Numer. Math.: Vol. 39). Basel: Birkhäuser-Verlag 1978, pp. 101 – 124.

[8] Gorenflo, R. and F. Mainardi: Random walk models for space-fractional diffusion processes. Fract. Cal. Appl. Anal. 1 (1998), 167 – 191.

[9] Gorenflo, R. and M. Niedack: Conservative difference schemes for diffusion problems with boundary and interface conditions. Computing 25 (1980), 299 – 316.

[10] Laha, R. G. and V. K. Rohatgi: Probability Theory. New York: Wiley 1979.

[11] Lévy, P.: Calcul des probabilités. Paris: Gauthier-Villars 1925.

[12] Lévy, P.: Théorie de l’addition des variables aléatoires, 2nd ed. Paris: Gauthier-Villars 1954 (1st ed. 1937).

[13] Montroll, E. W. and B. J. West: On an enriched collection of stochastic processes. In: Studies in Statistical Mechanics, Vol. VII: Fluctuation Phenomena (eds.: E. W. Montroll and J. L. Lebowitz). Amsterdam: North-Holland 1979, pp. 61 – 175.

[14] Richtmyer, R. D. and K. W. Morton: Difference Methods for Initial-Value Problems, 2nd Edn. New York: Intersci. Publ. (Wiley) 1967 (1st ed. 1957).

[15] Ross, S. M.: Introduction to Probability Models. New York: Acad. Press 1972.

[16] Samko, S. G., Kilbas, A. A. and O. I. Marichev: Fractional Integrals and Derivatives: Theory and Applications (Translated from the Russian). Amsterdam: Gordon and Breach 1993.

[17] Samorodnitsky, G. and M. S. Taqqu: Stable non-Gaussian Random Processes. New York: Chapman and Hall 1994.

[18] Schneider, W. R.: Stable distributions: Fox function representation and generalization. Lect. Notes Phys. 262 (1986), 497 – 511.

[19] Zolotarev, V. M.: One-dimensional Stable Distributions (Translated from the Russian). Providence, R.I.: Amer. Math. Soc. 1986.

Received 18.06.1998