Monte Carlo simulation of uncoupled continuous-time random walks yielding a stochastic solution of the space-time fractional diffusion equation

We present a numerical method for the Monte Carlo simulation of uncoupled continuous-time random walks with a Lévy $\alpha$-stable distribution of jumps in space and a Mittag-Leffler distribution of waiting times, and apply it to the stochastic solution of the Cauchy problem for a partial differential equation with fractional derivatives both in space and in time. The one-parameter Mittag-Leffler function is the natural survival probability leading to time-fractional diffusion equations. Transformation methods for Mittag-Leffler random variables were found later than the well-known transformation method by Chambers, Mallows, and Stuck for Lévy $\alpha$-stable random variables and so far have not received as much attention; nor have they been used together with the latter in spite of their mathematical relationship due to the geometric stability of the Mittag-Leffler distribution. Combining the two methods, we obtain an accurate approximation of space- and time-fractional diffusion processes almost as easy and fast to compute as for standard diffusion processes.

Figure 1

(Color online) Sample paths of CTRWs with scale parameters ${\gamma }_t=0.001$, ${\gamma }_x={\gamma }_t^{\beta ∕\alpha }$, and different choices of $\alpha$ and $\beta$. With smaller $\alpha$ the jumps become larger; with smaller $\beta$ the waiting times become longer.

Figure 2

(Color online) The Mittag-Leffler complementary cumulative distribution function sampled from Eq. (20) (circles) and computed analytically (solid line) [39], as well as its approximations for $t\to 0$ (Weibull function, long dashes) and $t\to \infty$ (power law, short dashes).

Figure 3

(Color online) Decay of the probability density $p_{{\gamma }_x,{\gamma }_t}(x,t;\alpha ,\beta )$ with $\alpha =1.7$, $\beta =0.8$, ${\gamma }_t=0.1$, and ${\gamma }_x={\gamma }_t^{\beta ∕\alpha }$. The crest at $x=0$ is the survival function $\Psi \left(t\right)=E_{\beta }(-{(t∕{\gamma }_t)}^{\beta })=P(0^{+},t)-P(0^{-},t)$, where $P(x,t)={\int }_{-\infty }^{x}p(u,t)du$.

Figure 4

(Color online) Convergence of $t^{\beta ∕\alpha }{p}_{{\gamma }_x,{\gamma }_t}(x,t;\alpha ,\beta )$ to the scaling function $W(x∕t^{\beta ∕\alpha };\alpha ,\beta )$, Eq. (8), at $t=2$ for selected values of $\alpha$ and $\beta$. The curves are shown in a time-independent way as scaling plots and appear in the same order from bottom to top as reported in the legend—i.e., with decreasing ${\gamma }_t$. The curve with the smallest ${\gamma }_t$ is almost indistinguishable from its theoretical limit $W$ (solid black line). However, in spite of the impression that may arise from the few terms and the ranges chosen here, in general the function sequences are not monotonic. The scale parameters ${\gamma }_x$ and ${\gamma }_t$ tend to 0 as ${\gamma }_x^{\alpha }={\gamma }_t^{\beta }$. The central peak decreases when the ratio $t∕{\gamma }_t$ becomes larger, as is evident in Fig. 3.

Figure 5

(Color online) Convergence of ${\max }_{x\ne 0}\mid p_{{\gamma }_x,{\gamma }_t}(x,t;\alpha ,\beta )-u(x,t;\alpha ,\beta )\mid$ for selected values of $\alpha$ and $\beta$ when ${\gamma }_x$, ${\gamma }_t\to 0$ with ${\gamma }_x^{\alpha }={\gamma }_t^{\beta }$.