Diffusion in heterogeneous media: An iterative scheme for finding approximate solutions to fractional differential equations with time-dependent coefficients

JID:YJCPH AID:5437 /FLA

[m3G; v 1.135; Prn:3/09/2014; 9:18] P.1 (1-15)

Journal of Computational Physics ••• (••••) •••–•••

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

Diffusion in heterogeneous media: An iterative scheme for finding approximate solutions to fractional differential equations with time-dependent coefficients Mauro Bologna a , Adam Svenkeson b,c,∗ , Bruce J. West d , Paolo Grigolini b a

Instituto de Alta Investigación, Universidad de Tarapacá, Casilla 6-D, Arica, Chile Center for Nonlinear Science, University of North Texas, P.O. Box 311427, Denton, TX 76203-1427, USA Army Research Laboratory, 2800 Powder Mill Road, Adelphi, MD 20783, USA d Information Science Directorate, Army Research Office, Research Triangle Park, NC 27709, USA b c

a r t i c l e

i n f o

Article history: Received 29 April 2014 Received in revised form 7 August 2014 Accepted 17 August 2014 Available online xxxx Keywords: Fractional calculus Strange kinetics Anomalous diffusion Fractional index expansion Time-dependent coefficients

a b s t r a c t Diffusion processes in heterogeneous media, and biological systems in particular, are riddled with the difficult theoretical issue of whether the true origin of anomalous behavior is renewal or memory, or a special combination of the two. Accounting for the possible mixture of renewal and memory sources of subdiffusion is challenging from a computational point of view as well. This problem is exacerbated by the limited number of techniques available for solving fractional diffusion equations with time-dependent coefficients. We propose an iterative scheme for solving fractional differential equations with time-dependent coefficients that is based on a parametric expansion in the fractional index. We demonstrate how this method can be used to predict the long-time behavior of nonautonomous fractional differential equations by studying the anomalous diffusion process arising from a mixture of renewal and memory sources. © 2014 Elsevier Inc. All rights reserved.

1. Introduction Anomalous diffusion has been explained using strange kinetics [1,2]; a term adopted from nonlinear dynamics to emphasize the lack of a characteristic scale in space, time or both in the phenomena of interest. This concept has been used to describe complex physical [3], social [4] and biological [5,6] processes. In the presence of strange kinetics the traditional diffusion equation is replaced with a fractional diffusion equation in order to adequately describe subdiffusion and superdiffusion. An assumption that was made until fairly recently, with little or no discussion in the construction of the equation of motion for classical diffusion was the homogeneity of the ambient fluid in which the diffusing particle (tracer) is embedded. In the case considered herein inhomogeneity results in the diffusion coefficient being time dependent, resulting in a scaled Brownian motion process [7–9]. However before we focus on this technical problem of solving a nonautonomous fractional differential equation in time we provide a general context for fractional diffusion equations. Höfling and Franosch [6] point out that in cell biology the diffusion of macromolecules and organelles is anomalous with effects such as time-dependent diffusion coefficients, persistent correlations in time, and non-Gaussian spatial displacements occurring as manifestations of macromolecular crowding in cells and in cellular membranes. Leptos et al. [10] conducted

*

Corresponding author. E-mail addresses: [email protected] (M. Bologna), [email protected] (A. Svenkeson), [email protected] (B.J. West), [email protected] (P. Grigolini). http://dx.doi.org/10.1016/j.jcp.2014.08.027 0021-9991/© 2014 Elsevier Inc. All rights reserved.

JID:YJCPH AID:5437 /FLA

[m3G; v 1.135; Prn:3/09/2014; 9:18] P.2 (1-15)

M. Bologna et al. / Journal of Computational Physics ••• (••••) •••–•••

2

experiments on the motion of tracers suspended in a fluid of swimming Eukaryotic microorganisms of varying concentrations. The interplay between the inanimate tracer particles and the advection by flows from the swimming microorganisms resulted in their displacement having a self-similar probability density function (pdf ) with a Gaussian core and exponential tails. A reanalysis of their data by Eckhardt and Zammert [11] shows an excellent fit to a Mittag-Leffler pdf derived using the CTRW model of Montroll and Weiss [12]. Zaid et al. [13] present a theoretical study of a simplified tracer–swimmer interaction showing that the non-Gaussian effect of the tails of the pdf can be modeled using a fractional diffusion equation. However most complex biological phenomena do not have fractional diffusion equations that give rise to analytic solutions and therefore they must be solved numerically. There is the rub, most fractional diffusion equations cannot be numerically integrated in a straightforward manner. There is a great deal of subtlety in their integration depending on the fractional indices in space and time and therefore herein we adopt an approach that relies on a new method for obtaining an analytic approximation to the solution of a nonautonomous fractional differential equation. The equivalence between fractional diffusion in time and subordination [14] makes it possible to bypass the need of designing an efficient numerical algorithm to solve a problem described by fractional derivatives in time [15], insofar as it leads to a procedure resting on ordinary derivative numerical algorithms, supplemented by the assumption that dynamics are frozen for extended time intervals. However, there are interesting cases where this procedure cannot be applied. Bologna et al. [16] recently studied a diffusion process generated by the joint action of correlated fluctuations and trapping process that by itself would lead to subdiffusion whereas the correlated fluctuations may generate either subdiffusion or superdiffusion. The trapping process, as in Refs. [14] and [15], is simulated by fractional derivatives in time. In this case the subordination procedure is invalidated [16] and recourse to an efficient method of solving the fractional differential equation is required. The present paper is the sequel to that analysis and herein, using iterative techniques, we address the question of how to construct the approximate solution to the fractional diffusion equation describing the joint action of trapping and correlated fluctuations. The results of this study are expected to be of interest to the research work aiming to establish whether the source of anomalous behavior is renewal or memory [5,17]. 1.1. Fractional kinetic equations The general expression for renewal anomalous diffusion is given by the fractional diffusion equation, see for example [3, 4,18,19],

∂α ∂β p ( x , t ) = D p (x, t ), ∂tα ∂ xβ

p (x, t = 0) = δ(x),

(1)

with α ≤ 1 and β ≤ 2. The condition α = 1 and β = 2 corresponds to the ordinary diffusion equation, which only applies to homogeneous media. Real diffusion processes in complex media are called anomalous since they do not scale in the traditional way and they obey the more general fractional diffusion equation, with α = 1 and β = 2. The renewal nature of this picture is easily explained metaphorically for both α < 1 and β < 2 using Aesop’s fable of the race between the tortoise and the hare [20]. The random walker is the hare who takes long naps intermittently disrupted by very long jumps. In fact, the condition α < 1 corresponds to the hare resting at the site with coordinate x for an extended time drawn from a waiting-time pdf ψ(τ ) with the time asymptotic inverse power-law structure

ψ(τ ) ∝

1

τ 1 +α

.

(2)

At the end of the n-th rest time interval of length τn a new time interval τn+1 is selected, with no correlation with the lengths of the earlier time intervals. Thus, at the end of a very long time interval with no event, the system is renewed through the choice of a new time interval that does not have any correlation with the system’s past. The source of anomaly is not given by the memory of the past but by the fact that for μ = 1 + α < 2 the mean waiting time

∞ τ  =

t ψ(t )dt

(3)

0

is infinite. Consequently, longer and longer time intervals are expected to appear between events so as to make τ  = ∞ possible. To avoid ambiguity in our definition of memory we note that while it is true that an individual walker has no memory of the past when making the next step, an increasing fraction of the walkers in an ensemble appears to be immobile as time advances because of the infinite mean waiting time between jumps [21]. In short, there is memory in the evolution of the pdf, as made evident by the fractional derivative in time in Eq. (1), but this definition of memory is different from the slowly-decaying correlation between steps of the random walk that we are referring to. We use the Caputo fractional derivative in time defined in terms of the Laplace transform

L

 α  ∂ = sα p ( x , t ); s p (x, s) − sα −1 p (x, 0) ∂tα

(4)

JID:YJCPH AID:5437 /FLA

[m3G; v 1.135; Prn:3/09/2014; 9:18] P.3 (1-15)

M. Bologna et al. / Journal of Computational Physics ••• (••••) •••–•••

3

where  p (x, s) denotes the Laplace transform of the pdf in time. See Pramukkul et al. [14] for recent work on the issues connecting a fractional derivative in time with the underlying stochastic processes, in particular the connection of the fractional time derivative with the inverse power-law waiting-time pdf. We adopt the Riesz fractional derivative in space defined in terms of the Fourier transform in the symmetric case



 ∂β F p (x, t ); k = −|k|β  p (k, t ) ∂ xβ

(5)

where  p (k, t ) is the Fourier transform of the pdf in space and is known as the characteristic function. As far as the condition β < 2 is concerned, it is important to stress that the meaning of the Riesz fractional derivative is as follows. When the hare makes a jump forward or backward she selects the length of the jump from the pdf





p |ξ | ∝

1

|ξ |μξ

,

(6)

with 1 < μξ < 3. As a consequence of this inverse power-law pdf the second moment of the step length diverges, ξ 2  = ∞, thereby violating the condition for the central limit theorem and resulting in the emergence of non-Gaussian diffusion from the random walk. The process is renewal because the choice of the n-th jump ξn does not have any memory of earlier jumps. We refer the reader to the excellent review article [4] and to the original paper on the fractional calculus and anomalous diffusion [22] for discussions of Eq. (5). 1.2. Time-dependent coefficients In sharp contrast with the renewal character of Eq. (1) anomalous diffusion may emerge from a process with infinite memory as described by the equation [16,23]

∂ ∂2 p (x, t ) = D H t 2H −1 2 p (x, t ). ∂t ∂x

(7)

The phase space equation of evolution for the pdf given by Eq. (7) is equivalent to the Langevin equation

d dt

x(t ) = ξ(t ),

(8)

with ξ(t ) being fractional Gaussian noise [24], namely, a process with Gaussian statistics and power-law correlation function. We notice that the Gaussian hypothesis is fundamental for the dynamical derivation of a diffusion process with time dependent diffusivity, as discussed in [25]. Eq. (7) is the time asymptotic limit of the equation originally found in [25], and it is interesting to notice that these authors addressed also the important problem of studying the influence of a containing potential that is currently being discussed by Jeon, Chechkin and Metzler [26]. Eq. (7) also arises in the context of homogeneous turbulence where the time-dependent coefficient is the dispersion of the velocity fluctuations [27] and Richardson’s empirical mean-square displacement between two tracer particles increasing as t 3 [28,29] entails a scaling parameter in the present context of H = 3/2. This latter case was discussed by Thiel and Sokolov [8], with the simple Langevin equation replaced with scaled Brownian motion such that ξ(t ) → D H t 2H −1 ζ (t ) on the right-hand side of Eq. (8) and ζ (t ) is white Gaussian noise. The distinction between scaled Brownian motion and fractional Brownian motion (FBM) is discussed in detail by Lim and Muniandy [7], as well as by Jeon et al. [9]. In a biological context Jeon et al. [30] combine single particle tracking microscopy data of endogenous lipid granules in living fission yeast cells and establish that at short times the data are best described by CTRW subdiffusion and at long times by FBM subdiffusion. The steps in CTRW have finite variance but inverse power-law waiting times whereas in FBM the steps have power-law correlations. Consequently, there is weak ergodicity breaking [21] at early times so that the time averages and ensemble averages do not coincide in such complex biological systems [31]. Weigel et al. [32] also establish using single-molecule tracking that ergodic and non-ergodic processes coexist in plasma membranes. Tabei et al. [33] find that intracellular transport of insulin granules cannot be described by either the CTRW or FBM alone, and propose a model that mixes the two. These remarks are intended to emphasize the need to advance the theoretical understanding of diffusion processes with mixed origins (see for instance the recent progress made by Schulz et al. [34]). Also, they are intended to caution the reader regarding the interpretation of averages using the solution to fractional diffusion equations when the underlying physical process contains a non-ergodic component. Here we introduce dimensional analysis and adopt the notation that [ y ] denotes the dimension of the variable y. The scaling of the characteristic function is determined by the condition that the dimension of the dynamic variable scales as

[x] = [t ]η where

(9)

η denotes the scaling coefficient, which for Eq. (1) was found to be [16,23]

η=

α β

.

(10)

JID:YJCPH AID:5437 /FLA

[m3G; v 1.135; Prn:3/09/2014; 9:18] P.4 (1-15)

M. Bologna et al. / Journal of Computational Physics ••• (••••) •••–•••

4

In the case of the renewal anomalous diffusion of Eq. (1) this scaling index is found to be equivalent to the Hurst exponent [16,23]

η = H,

(11)

also in the case of the process with infinite memory of Eq. (7). Note that the adoption of the symbol H to denote the scaling coefficient of the anomalous diffusion with infinite memory goes back to the foundation of fractional Brownian motion [35] explaining the persistent and anti-persistent properties of droughts and floods of the Nile river observed by Hurst [36]. The symbol H was adopted by Mandelbrot to honor Hurst’s discoveries. A theoretical problem of conceptual and practical interest is to assess the origin of the deviation from ordinary diffusion, namely, the origin of the deviation from classical scaling

1

η = .

(12)

2

One strategy for determining the source of the deviation from classical scaling is given by [17]. However an even more challenging problem results from the assumption that real diffusive phenomena may not be separately renewal or have infinite memory but comprise a combination of the two [16]. Such a blending of effects would give rise to the fractional diffusion equation with a time-dependent diffusion coefficient

∂α ∂β p (x, t ) = D H t 2H −1 β p (x, t ). α ∂t ∂x

(13)

It is worth noting here that Eq. (13) differs from the fractional Fokker–Planck equation (FFPE)

∂ β ∂ 1 −α ∂ p (x, t ) = K H (t ) β 1−α p (x, t ) ∂t ∂ x ∂t

(14)

that under the condition of assuming the Caputo fractional derivative of this article can be exchanged with the Riemann– Liouville fractional derivative, would be compatible with the form given by Magdziarz et al. [37] (see also Ref. [38]). This is because the fractional time derivative on the right does not commute with the time-dependent diffusion coefficient K H (t ) or more generally it does not commute with a Fokker–Planck operator having time-dependent coefficients such as a time-dependent force. The intuitive physical meaning of Eq. (13) is that the particle diffuses under the influence of correlated fluctuations, a process that in the anti-dispersive condition may lead to diffusion coefficient D (t ) vanishing for t → ∞. In addition to that the particle is trapped for extended times and perceives these correlated fluctuations only when it jumps. It is easy to prove [16] that in this case the scaling η of Eq. (10) becomes

2H − 1 + α

η=

β

.

(15)

As a consequence of the joint action of trapping and anti-dispersive fluctuation ( H < 0.5) the particle diffusion is completely quenched, rather than gradually moving towards D (t ) = 0 for t → ∞. Adopting the Fourier representation, we transform Eq. (13) to obtain [4,22] ∂α

∂tα

 p (k, t ) = −|k|β D H t 2H −1 p (k, t ),

(16)

thereby initiating the study of fractional differential equations with time-dependent coefficients, which we denote for simplicity as dα

dt α

Y (t ) = − f (t )Y (t ),

(17)

where

Y (t ) ≡  p (k, t ).

(18)

Note that in this article we adopt the following notation

f (t ) = 2H χ t 2H −1 . Consequently

χ=

(19)

χ is defined as

|k|β D H 2H

.

(20)

Rivero et al. [39] present a general approach to solving Eq. (17) based on a generalization of the Frobenius regular

α -singular point method. When applicable the method provides an exact series solution around the regular α -singular point. However, in the case of Eq. (17) this method is only applicable when α = 1 − 2H , causing a regular α -singular point at t = 0. Herein we provide an iterative technique for solving Eq. (17) when α = 1 − 2H that is valid for all t, and provide analytic expressions for the asymptotic time limit t → ∞.

JID:YJCPH AID:5437 /FLA

[m3G; v 1.135; Prn:3/09/2014; 9:18] P.5 (1-15)

M. Bologna et al. / Journal of Computational Physics ••• (••••) •••–•••

5

1.3. Preview Since the only systematic way to solve such fractional differential equations as Eq. (17) analytically is in restricted time domains [39] one has recourse to numerical integration. But the numerical solutions of fractional differential equations present a large menu of unsolved problems. Setting these issues aside our purpose herein is to present an iterative scheme for solving Eq. (17). In Section 2 we show that the equivalence between the fractional calculus representation and subordination representation of a dynamic process is violated by the memory-induced time dependence of the diffusion coefficient. This forces us to adopt an iterative scheme to obtain approximate analytic solutions. After a short illustration of the popular predictor–corrector algorithm in Section 3 we move to a discussion of a new iterative scheme in Section 4, based on a parametric expansion of α around α = 1 if α is close to 1, or on the expansion of α around the vanishing value if α  1. Here we also find approximate solutions to Eq. (16) that result from applying the first steps of the iterative process. Section 5 is devoted to concluding remarks and, finally, in Appendix A we shortly review some basic properties of the transformation from operational to chronological or clock time. 2. Breakdown of the equivalence between subordination and fractional calculus When integrating a system of ordinary differential equations there are many well-established numerical schemes to choose from. This is not yet true when moving to the numerical integration of fractional differential equations. Fortunately, a way to bypass the direct integration of fractional differential equations exists through subordination. The natural connection between fractional calculus and subordination [14,40] allows the solution of an ordinary differential equation to be transformed into the solution of its fractional counterpart [15]. However, unsolved problems with this approach arise when the fractional differential equations are not autonomous, as in Eq. (17). To establish the breakdown of the subordination approach to fractional calculus induced by long-term memory, let us focus our attention on the Caputo fractional equation with a time-dependent rate

dα dt α

Y (t ) ≡

t

1

Γ (1 − α )

1

d

(t − τ )α dτ

0

Y (τ )dτ = −2H χ t 2H −1 Y (t ).

The adoption of the subordination approach implies that we first obtain the solution to Eq. (21) with



Y (τ ) = exp −χτ

2H



Y 0.

(21)

α = 1: (22)

Using the theoretical perspective of subordination [41] we interpret τ as the operational time and shift perspectives from the solution in operational time τ to the solution in chronological or observational time t. This shift is accomplished by implementing the transformation

∞ Y (t ) =

∞

dτ P

(S)

0





dτ P ( S ) (t , τ ) exp −χτ 2H Y 0 .

(t , τ )Y (τ ) =

(23)

0

Of course the entire theory is dependent on the choice of the kernel in Eq. (23) that subordinates chronological time to operational time. One strategy for determining the subordination kernel in Eq. (23) is to use the Continuous Time Random Walk (CTRW) of Montroll and Weiss [12], where the random walker makes jumps at discrete times indexed by 1, 2, . . . , n, . . . , corresponding to the discrete experimental times t 1 , t 2 , . . . , tn , . . . . The experimental times can be used to construct the sequence of time intervals tn+1 − tn = ηn . The time intervals {ηn } are assumed to be derived from a waiting-time pdf for the chronological time

ψα (t ) = −



d dt



E α −(λt )α ,

(24)

where





E α −(λt )α =

∞

(−1)k (λt )kα k =0

Γ (kα + 1)

(25)

denotes the Mittag-Leffler function. The discrete time n becomes indistinguishable from the continuous time τ when n 1. Note that the connection between chronological t and operational τ is established by setting η1 + η2 + · · · + ηn ≤ t < η1 + η2 + · · · + ηn + ηn+1 . This can be interpreted as an asymmetric Lévy process with P ( S ) (t , τ ) denoting the distribution density of the variable t at time τ . More details on the function P ( S ) (t , τ ) can be found in Appendix A. In the case α = 0.5 the function P ( S ) (t , τ ) has the analytical expression 2

P

(S)

(t , τ ) =

exp[− 4τλ t ]

√

0

π λ0 t

.

(26)

JID:YJCPH AID:5437 /FLA

[m3G; v 1.135; Prn:3/09/2014; 9:18] P.6 (1-15)

M. Bologna et al. / Journal of Computational Physics ••• (••••) •••–•••

6

By inserting Eq. (26) into Eq. (23) we obtain

Y (t ) = √

∞

Y0

π λ0 t







exp −χτ 2H exp − 0

τ2 4 λ0 t

 dτ

(27)

as the quadrature solution to Eq. (21). 2.1. Example of equivalence Here we double-check the equivalence between fractional calculus and subordination in the case where there is no memory-induced time dependence of the diffusion coefficient. To eliminate the time dependence of the rate we set H = 1/2 in both Eq. (21) and Eq. (27). To establish an analytical comparison between the two theoretical perspectives we refer to Eq. (27), which is based on α = 0.5. Thus we must set α = 0.5 in Eq. (21) as well. Since we have set α = H = 0.5 in our example we can write the Laplace transform of Eq. (21) as

s1/2 Y (s) − s−1/2 Y 0 = −χ  Y (s),

(28)

and rearrange terms to obtain

1  Y (s) = 1/2 s

Y0 s1/2 + χ

(29)

.

It is possible to invert exactly the Laplace transform. Indeed for Eq. (21) we have that

Y (t ) = e χ

2

t

erfc[χ

√

t ]Y 0 ≈

Y0

√

(30)

χ πt

where erfc[ z] = 1 − erf[ z] and erf[ z] is the error function. On the other hand, the exact evaluation of the integral of Eq. (27) with H = 0.5 yields

Y (t ) = e λ0 χ

2

t



erfc[ λ0 t χ ]Y 0 ,

(31)

and, after setting λ0 = 1, Eq. (31) coincides with Eq. (30). In summary when the memory-induced time dependence is suppressed, fractional calculus and subordination yield the same result, as previously determined [14,15]. 2.2. Equivalence breakdown Let us now return to the integral-differential equation with long-term memory

dα dt α

Y (t ) = −2H χ t 2H −1 Y (t ),

(32)

which can be solved using a series expansion. Introducing the series representation for the solution into the fractional differential equation with a time-dependent rate

Y (t ) =

∞

cn t nγ ,

(33)

n =0

we obtain by direct substitution ∞

cn t nγ −(α +2H −1)

n =1

∞

Γ (nγ + 1) cn t nγ . = −2H χ Γ (nγ − α + 1)

(34)

n =0

Note that on the left-hand side of Eq. (34) the fractional derivative of a monomial was introduced just as would be done when using the ordinary calculus. Equating the powers of time on both sides of the equation imposes the equality among scaling indices

γ = α + 2H − 1

(35)

enabling us to write the recurrence relation between expansion coefficients

cn+1 = −2H χ

Γ [(n + 1)γ − α + 1] cn . Γ [(n + 1)γ + 1]

(36)

JID:YJCPH AID:5437 /FLA

[m3G; v 1.135; Prn:3/09/2014; 9:18] P.7 (1-15)

M. Bologna et al. / Journal of Computational Physics ••• (••••) •••–•••

7

Eq. (36) implies that the coefficients can be expressed as the product of the ratio of gamma functions, a procedure familiar from the analysis of ordinary partial differential equations. Thus, we obtain the expression for the expansion coefficients

cn (k) = c 0 (−2H χ )n

−1 n i =0

Γ [(i + 1)γ + 1 − α ] , Γ [(i + 1)γ + 1]

n ≥ 1.

(37)

Note that the series Eq. (33) with the coefficient given by Eq. (37) has an infinite radius of convergence. The problem therefore has been formally solved exactly. Thus, we have two exact solutions: the fractional calculus solution expressed as an infinite series and the subordination solution expressed as a quadrature. Are the two equivalent? That is not an easy question to answer and so we do it in stages, first by comparing the two for specific values of the coefficients and then more generally, but in both cases the solutions are compared asymptotically in time. We now consider the particular case H = 1 and rewrite Eq. (32) as dα

dt α

Y (t ) = −2χ tY (t ).

(38)

Taking the Laplace transform of this equation we obtain the differential equation in the Laplace representation

sα  Y ( s ) − s α −1 Y 0 = 2χ

d ds

 Y (s).

(39)

Solving this latter equation under the condition that  Y (s) has an inverse Laplace transform we obtain by straightforward integration

 Y (s) = Y 0 exp



s α +1

 s

 exp −

2(1 + α )χ

0

u α +1

 α −1 u

2(1 + α )χ

2χ

du .

For s → 0 we obtain, after scaling the integration variable by s, sα Y

 Y (s) ≈

0

2αχ

,

(40)

(41)

that using the Tauberian theorem can be shown to imply

Y (t ) ∼ t −α −1 .

(42)

Note that this asymptotic result is for the case H = 1. Consequently, to compare the exact result of Eq. (27) we set the fractional index α = 1/2 in Eq. (42) to get

Y (t ) ∼ t −3/2

(43)

as the asymptotic scaling of the exact solution to the fractional differential equation with a rate that increases linearly in time. The subordination solution for α = 1/2 is still given by Eq. (27). We can compare the exact subordination solution to the exact series solution asymptotically by setting H = 1. This yields

1 1 Y (t ) = √ ≈ √ . 2 λ0 χ t 1 + 4t λ0 χ

(44)

Upon comparison with Eq. (43), it is clear that the asymptotic solutions from the two approaches do not coincide. More generally the integral Eq. (27) can be evaluated when t → ∞. Indeed, for H = 0 we may write the following approximate expression

Y (t ) = √

∞

1

π λ0 t

≈√



0

∞

1

π λ0 t





exp −χτ 2H exp −



τ2 4 λ0 t

 dτ



exp −χτ 2H dτ

(45)

(46)

0

and the latter integral can be evaluated to yield

Y (t ) ≈

1 Γ ( 2H + 1)

χ

1 2H

√

1

π λ0 t

.

Note that setting H = 1 recovers the result of Eq. (44). Consequently, the two approaches are not equivalent for

(47)

α = 1/2.

JID:YJCPH AID:5437 /FLA

[m3G; v 1.135; Prn:3/09/2014; 9:18] P.8 (1-15)

M. Bologna et al. / Journal of Computational Physics ••• (••••) •••–•••

8

2.3. General condition It is interesting to notice that P ( S ) (t , τ ) can also be interpreted as the pdf of the operational time time t. West and Grigolini [42] show, using the CTRW, that the scaling solution yields a pdf of the form

1 P ( S ) (τ , t ) = α F t



τ at chronological



τ

(48)

.

tα

Inserting Eq. (48) into the quadrature solution Eq. (23) for subordination leads to

1

∞ 

Y (t ) = α t

F

τ







exp −χτ 2H dτ .

tα

(49)

0

Note that the function P ( S ) (τ , t ) for t → ∞ is expected to become flat as any diffusion process does. Thus, in the asymptotic limit we have



lim F

t →∞

τ



tα

= c,

(50)

where c is a constant, so that the asymptotic form of the solution becomes

c

∞

lim Y (t ) = α t

t →∞





exp −χτ 2H dτ ∝ t −α .

(51)

0

In summary, it is seen that subordination always leads to the asymptotic condition t −α , whereas fractional derivatives only generate this asymptotic condition when H = 0.5. In the absence of memory, the result for fractional derivatives is in agreement with the general result obtained in Ref. [14]. The adoption of the arguments of Ref. [14], in fact, shows that in this case fractional calculus and subordination are equivalent. From a physical point of view the breakdown of this equivalence is due to the fact that the time-dependent diffusion coefficient D (t ) can be interpreted as a perturbation exerted at chronological time on a process whose roots are in the operational time. This is a challenging issue that has recently been the subject of debate among some groups [43–45]. 3. Discussion of a current algorithm Developing algorithms capable of accurately simulating the time evolution of fractional trajectories is an area of active research. The nonlocal memory effects induced by the fractional derivative operator require an extension of the traditional integration schemes, and present the main challenge in generating efficient algorithms. The iterative scheme we propose in Section 4 is fundamentally different from existing methods, and as such, is intended to be useful for assessing the validity of existing numerical algorithms. In order to fully appreciate the different angle of approach, in this section we discuss for comparison purposes the predictor–corrector algorithm [46], a commonly used method of generating solution trajectories belonging to a wide class of fractional differential equations. We stress that this section is a brief review of earlier work [46, 47] included in order to make the present article self-contained and these references ought to be consulted for additional details on the predictor–corrector algorithm. The algorithm, developed by Diethelm, Freed, and Ford [46,47], relies on the standard predictor–corrector integration technique for ordinary differential equations, adapted to the purpose of integrating fractional differential equations. As in the ordinary case, the fractional differential equation is turned into an integral equation. The next unknown point along the solution trajectory is predicted through a basic lowest-order approximation of the integrand, and then corrected using a higher-order approximation. To illustrate this algorithm in more detail, consider a fractional differential equation fitting the form of Eq. (17), which we rewrite here in a slightly more general way as

dα dt α





Y (t ) = G t , Y (t ) .

(52)

We would like to find the solution of Eq. (52) on the interval [0, T ]. We divide the interval into discrete time steps of equal length h. Since we are considering the case where 0 < α ≤ 1, only one initial condition is needed to define the trajectory, that is Y (t = 0) = Y 0 . The initial value problem is then turned into the Volterra integral equation

Y (t ) = Y 0 +

1

t

Γ (α )

  (t − τ )α −1 G τ , Y (τ ) dτ .

0

The equivalence between Eq. (52) and Eq. (53) can be checked through their respective Laplace transforms.

(53)

JID:YJCPH AID:5437 /FLA

[m3G; v 1.135; Prn:3/09/2014; 9:18] P.9 (1-15)

M. Bologna et al. / Journal of Computational Physics ••• (••••) •••–•••

9

The task, then, is to find a good approximation to the integral in Eq. (53). For this purpose, following Diethelm and Freed [47], we adopt the product trapezoidal quadrature formula so that t n +1 t n +1   α −1 (tn+1 − τ ) G (τ )dτ ≈ (tn+1 − τ )α −1 G n+1 (τ )dτ . 0

(54)

0

The exact function G (τ ) is replaced by the approximate piecewise function G n+1 (τ ). Note here we are using a concise notation where G (τ ) = G (τ , Y (τ )). The integral on the right-hand side of Eq. (54) can instead be written as a sum, t n +1  n +1

(tn+1 − τ )α −1 G n+1 (τ )dτ = a j ,n+1 G (t j ).

(55)

j =0

0

The coefficients are given by

⎧ α +1 − (n − α )(n + 1)α ] ⎪ j = 0, ⎨ [n α + 1 α + 1 α + 1 a j ,n+1 = [(n − j + 2) − 2(n − j + 1) + (n − j ) ] 1 ≤ j ≤ n, α (α + 1) ⎪ ⎩ 1 j = n + 1. hα

This approximation leads to the corrector formula

Y n +1 = Y 0 +

1

Γ (α )



n

j =0



(56)



(P )  a j ,n+1 G (t j , Y j ) + an+1,n+1 G tn+1 , Y n+1

.

(57)

Eq. (57) by itself does not yet give us an algorithm to generate the fractional trajectory, as the right-hand side depends (P ) on the predicted value Y n+1 . Rather, assuming the predicted value is known, Eq. (57) uses the predicted value to move closer to the real value. To find the predicted value we must make an approximation even simpler than the trapezoidal rule. Again following Diethelm and Freed [47], the predictor formula is derived from a product rectangle rule, t n +1  n

(tn+1 − τ )α −1 G (τ )dτ ≈ b j ,n+1 G (t j ),

(58)

j =0

0

where

b j ,n+1 =

hα 

α

(n + 1 − j )α − (n − j )α .

(59)

This simpler approximation yields the predictor equation (P )

Y n +1 = Y 0 +

1

Γ (α )

n

b j ,n+1 G (t j , Y j ).

(60)

j =0

Eq. (57) and Eq. (60) together form the predictor–corrector algorithm, and can be used to generate the approximate values Y 1 , Y 2 , . . . , Y ( T ) along the fractional trajectory starting from the initial condition Y 0 . 4. Beyond the limits of the current algorithms The equivalence between subordination and fractional calculus makes it possible to adopt subordination for numerical purposes, thereby bypassing the difficulties and numerical inaccuracies that might be generated by the direct integration of fractional differential equations. However, when the fractional differential equations are not autonomous, that equivalence breaks down, and we have to find another method to accurately calculate the fractional trajectory. It is of fundamental importance to rest on efficient algorithms as close as possible to analytical procedures. For that reason, in this section we propose an iterative scheme for finding approximate solutions to fractional differential equations with time-dependent coefficients. Let us study the integral-differential equation dα

dt α

Y (t ) = − f (t )Y (t ),

(61)

where the minus sign has been introduced for convenience. Here f (t ) is a sufficiently smooth arbitrary function. The specific form of f (t ) given by Eq. (19) for the example of anomalous diffusion will be considered later, after the general arguments are illustrated. Eq. (61) may be written explicitly in terms of the Caputo fractional derivative

−

1 f (t )Γ (1 − α )

t 0

1

d

(t − τ )α dτ

Y (τ )dτ = Y (t ).

(62)

JID:YJCPH AID:5437 /FLA

[m3G; v 1.135; Prn:3/09/2014; 9:18] P.10 (1-15)

10

M. Bologna et al. / Journal of Computational Physics ••• (••••) •••–•••

Performing the change of variable

τ = tu, after a little algebra Eq. (62) becomes

1 Y (t ) = − α t f (t )Γ (1 − α )

1

1

d

(1 − u )α du

0

Y (ut )du .

(63)

Using the fact that 0 < α < 1 we may write the following expansion in terms of the for 0 ≤ u ≤ 1, i.e.

α parameter of the factor (1 − u )−α ,

∞ 

(−1)k αk k (1 − u )−α = exp −α ln(1 − u ) = ln (1 − u ), k!

(64)

k =0

and inserting the expansion into Eq. (63) we may write ∞

(−1)k αk 1 Y (t ) = − α t f (t )Γ (1 − α ) k!

1

k =0

lnk (1 − u )

d du

Y (ut )du .

(65)

0

Note that the terms t α and Γ (1 − α ) do not need to be expanded since they are exact and so they include all

α -order.

Assuming a solution of the form

Y (t ) =

∞

αk F k (t ),

(66)

k =0

at zero-order we have from Eq. (65)

1 F 0 (t ) = − α t f (t )Γ (1 − α )

1

 1 F 0 (t ) − Y 0 F 0 (ut )du = − α du t f (t )Γ (1 − α ) d

0

(67)

where F 0 (0) = Y 0 . It follows that at zero-order we have after some algebra

F 0 (t ) =

Y0 1 + Γ (1 − α )t α f (t )

(68)

.

Iterating the procedure and equating the terms at the same example, for F 1 (t ) we obtain

1

1

F 0 + α F 1 (t ) = − α t f (t )Γ (1 − α )

d du

0

α + Γ (1 − α )t α f (t )

F 1 (t ) = − α t f (t )Γ (1 − α )

1

d

1

F 0 (ut )du − α α t f (t )Γ (1 − α )

ln(1 − u )

d du

d du

F 1 (ut )du

0

F 0 (ut )du .

(69)

0

α we find 1

du 0

1

1

Taking into account Eq. (67), at first-order in

1

αk -order we find the formal solution to the problem. For

F 1 (ut )du + α t f (t )Γ (1 − α )

1 ln(1 − u )

d du

F 0 (ut )du .

(70)

0

Solving Eq. (70) with respect to F 1 (t ) and imposing the condition F 1 (0) = 0 we find that

F 1 (t ) =

1 1 + Γ (1 − α )t α f (t )

1 ln(1 − u )

d du

F 0 (ut )du .

(71)

0

Thus, we have the first two lowest-order contributions to the general solution of fractional differential equations fitting the form of Eq. (61). Given that the method of solution is recursive (see Eq. (65) and following), the iterative strategy has the potential to be used as a numerical algorithm. As an example of an application of the formal perturbative solution just described, let us go back to the fractional diffusion equation with a time-dependent coefficient, Eq. (13), and its Fourier representation Eq. (16). Now we have f (t ) =

JID:YJCPH AID:5437 /FLA

[m3G; v 1.135; Prn:3/09/2014; 9:18] P.11 (1-15)

M. Bologna et al. / Journal of Computational Physics ••• (••••) •••–•••

2H χ t 2H −1 as in Eq. (19). We can use the zeroth-order solution in to Eq. (68), the lowest-order contribution F 0 (t ) is

F 0 (t ) =

Y0 1 + 2χ Γ (1 − α ) Ht α +2H −1

11

α to approximately solve Eq. (16). In this case, according (72)

.

Note, in particular, choosing H = 1 and taking the asymptotic limit, that we find

Y (t ) ∼ F 0 (t ) ∼ t −α −1 ,

(73)

which coincides with the previous result of Eq. (42). The second-order contribution to the example can be obtained by inserting Eq. (72) into Eq. (71) and integrating, which can be done analytically for t → ∞. Staying at the lowest order, after substituting Eq. (20) into Eq. (72) we arrive at the approximate solution to the fractional diffusion equation in Fourier space, Eq. (16), with the condition  p (k, 0) = 1,

 p (k, t ) ≈

1 1 + Γ (1 − α ) D H |k|β t α +2H −1

(74)

.

For β given by a rational number, the above Fourier transform can be inverted via standard methods. A simple example of the pdf is obtained by setting β = 2,





|x| 1 , p (x, t ) = exp − 2 Γ (1 − α ) D H t α +2H −1 Γ (1 − α ) D H t α +2H −1

(75)

which is a Laplace distribution with a time-dependent variance. We now consider the case when α approaches one. We note that 1/Γ (1 − α ) → 0 as α → 1 while (1 − u )−α approaches a non-integrable function at u = 1. To overcome this difficulty we perform an integration by parts and rewrite Eq. (63) as



1

Y (t ) = − α t Y˙ (0) + t f (t )Γ (2 − α ) where Y˙ (0) ≡

dY | . dt t =0

1

1 −α

(1 − u )



k =0 k !

t Y˙ (0)δk,0 +

Y (t ) = − α t f (t )Γ (2 − α )

du 2

Y (ut )du

(76)

0

Performing a Taylor expansion for

∞ εk



d2

α in the neighborhood of α = 1, that is ε = 1 − α  1, we have

1

d2

k

ln (1 − u )

du 2

 Y (ut )du .

(77)

0

Assume the solution to this equation has the series expansion form

Y (t ) =

∞

εk F k (t ), ε = 1 − α ,

(78)

k =0

ε0 -order contribution to the solution,  1

which when inserted into Eq. (77) yields the



1

F 0 (t ) = − α t Y˙ (0) + t f (t )Γ (2 − α )

d2

du 2

F 0 (ut )du ,

(79)

0

and consequently

1 F 0 (t ) = − α −1 F˙ 0 (t ). t f (t )Γ (2 − α )

(80)

The solution for F 0 (t ) is promptly found by direct integration to be



t

F 0 (t ) = exp −Γ (2 − α )

 zα −1 f ( z)dz Y 0 ,

(81)

0

reducing the lowest order term to quadratures. The next order contribution in ε to the solution is

1



F 1 (t ) = − α t F˙ 1 (t ) + t f (t )Γ (2 − α )

1 ln(1 − u ) 0

d2 du 2

 F 0 (ut )du .

(82)

JID:YJCPH AID:5437 /FLA

[m3G; v 1.135; Prn:3/09/2014; 9:18] P.12 (1-15)

M. Bologna et al. / Journal of Computational Physics ••• (••••) •••–•••

12

Imposing F 1 (0) = 0 we obtain after some algebra



t F 1 (t ) = −

t

exp −Γ (2 − α )

 zα −1 f ( z)dz dτ

τ

0

1

ln(1 − u ) d2

τ

du 2

F 0 (u τ )du .

(83)

0

Let us apply these results to Eq. (17). It follows that



F 0 (t ) = exp −

 Γ (2 − α ) D H α +2H −1 β |k| Y 0 , t α + 2H − 1

(84)

i.e. a Lévy distribution in x space is obtained by inverse Fourier transforming this expression. 5. Concluding remarks We have argued that the fractional calculus provides a convenient way to represent the physics of anomalous diffusion, and in particular, when describing the diffusion of tracer particles in a heterogeneous ambient fluid such as arises in the study of biological fluids. These fractional diffusion equations have analytic solutions, such as the Lévy distribution and the Mittag-Leffler function, in the absence of memory. However in the more general case when memory is present and the diffusion coefficient is a polynomial in time no analytic solutions have been found. To clarify the connection with other approaches to solving this problem we have demonstrated that those analyses, including our own, that rely on subordination lead to different results than those obtained from a straightforward integration of the fractional diffusion equation with time-dependent coefficients. Two integration strategies have been put forward to solve the anomalous diffusion equations with time-dependent coefficients. Both integration techniques rely on iterative perturbation expansions of the solution in the scaling index of the fractional time derivative, one for α  1 and the other for ε = 1 − α  1, and the utility of the techniques will be determined by the rate of convergence that is obtained for given time dependencies. Herein we have shown that reasonable perturbation solutions are obtained for polynomial time dependencies. The recursive nature of the analytical procedure opens the door for potential numerical implementation. We leave this problem open for future work. Note that the standard numerical methods for solving ordinary differential equations involve calculating approximate values to successive points along the solution trajectory in an iterative fashion starting from an initial value. This approach, as made clear by Section 3, has been carried over to the numerical methods for solving fractional differential equations, but the nonlocality of the fractional operators makes it computationally expensive to proceed this way when seeking asymptotic time behavior. In contrast, the iterative scheme we present readily yields the asymptotic behavior. Further iterations then improve the accuracy of this full solution over a more extended range of time. This alternative approach to the integration of fractional differential equations may be a better fit to the nonlocal structure of the operators. Acknowledgements A.S. and P.G. warmly thank ARO and Welch Foundation for their support through Grants No. W911NF-11-1-0478 and No. B-1577, respectively. M.B. acknowledges financial support from FONDECYT project No. 1120344. Appendix A. The transformation from operational to experimental time This appendix provides details for the derivation of the transformation function P ( S ) (t , τ ). A.1. Stochastic central limit theorem Following Pramukkul et al. [14] we write the Laplace transform of P ( S ) (t , τ ) as a Mittag-Leffler (ML) function in the Laplace variable,

Eα



 −(λt )α u =

∞ dτ P ( S ) (t , τ ) exp(−u τ ),

(A.1)

0

namely

   P ( S ) (t , u ) = E α −(λt )α u .

(A.2)

The derivation of Eq. (A.1) contained in [14] was based on applying the subordination perspective with an inverse powerlaw subordination function ψ(τ ) ∝ τ1μ with μ < 2, not necessarily identical to waiting time pdf corresponding to the ML function. The authors of Ref. [14] made the assumption that the probability of detecting the events generated by this prescription is very small, thereby generating a much slower survival probability. By rescaling time, so as to recover the time

JID:YJCPH AID:5437 /FLA

[m3G; v 1.135; Prn:3/09/2014; 9:18] P.13 (1-15)

M. Bologna et al. / Journal of Computational Physics ••• (••••) •••–•••

13

scale of the original process, the higher-order corrections to the Laplace transform of the ML function vanish. The resulting expression is Eq. (A.1) with the transformation function P ( S ) (t , τ ) generated using as a waiting-time pdf ψ(τ ) exactly the one associated with the ML function. It is well known that

    L E α −(λt )α ; s =

1 s + λα s 1 −α

(A.3)

.

Thus, using Eq. (A.1) we obtain the double Laplace transform

P

( S )∗

(s, u ) =

1 s + λα s 1 − α u

=

s α −1 λα . sα λα + u

(A.4)

Note that we have implemented separate Laplace transforms for the operational and chronological times. A.2. First approach The first approach to determining P ( S ) (t , τ ) relies on inverse Laplace transforming Eq. (A.4) with respect to the operational Laplace variable u followed by the inverse Laplace transform with respect to the chronological Laplace variable s. Eq. (A.4) makes it easy to inverse Laplace transform with respect to u yielding sα s α −1  P ( S ) ( s , τ ) = α e − λα τ . λ

(A.5)

Let us now inverse Laplace transform of this last equation with respect to s,

P ( S ) (t , τ ) =



1

dse st

2π i

sα −1 − ταs e λ , α

(A.6)

λ

Br

which is a contour integral along a Bromwich contour. This is the approach to obtaining P ( S ) (t , τ ) adopted by Fogedby [48] and Stanislavsky [40]. A.3. Second approach The second approach to determining P ( S ) (t , τ ) is realized by noticing that Eq. (A.1) can be written in the form of Eq. (A.2), thereby yielding the inverse Laplace transform with respect to u,

P ( S ) (t , τ ) =



1





due u τ E α −u (λt )α .

2π i

(A.7)

Br

A.4. General case Finding the analytical expression of P ( S ) (t , τ ) is difficult using either of the Laplace inversions. However, we can derive its expression in the asymptotic limit which agrees with the scaling arguments of Section 2.3. In Eq. (A.6) we make the long-time s → 0 approximation

  exp −τ sα sα −1 ≈ sα −1 .

(A.8)

In this limit the contour integral Eq. (A.6) yields

P ( S ) (t , τ ) ≈

sin(πα )Γ (α )

π (λt )α

(A.9)

,

in agreement with Section 2.3. We now turn our attention to the second Laplace inversion Eq. (A.7) and make use of the fact that in the long-time limit the ML function is proportional to 1/t α ,





E α −u (λt )α ≈

1

Γ (1 − α )u (λt )α

=

sin(πα ) Γ (α )

π

u (λt )α

.

(A.10)

Since 1/u is the Laplace transform of 1, we immediately see that the inverse Laplace transform of Eq. (A.10) yields Eq. (A.9) and the two expressions are asymptotically identical.

JID:YJCPH AID:5437 /FLA

14

[m3G; v 1.135; Prn:3/09/2014; 9:18] P.14 (1-15)

M. Bologna et al. / Journal of Computational Physics ••• (••••) •••–•••

References [1] M.F. Shlesinger, G.M. Zaslavsky, J. Klafter, Strange kinetics, Nature 363 (1993) 31–37, http://dx.doi.org/10.1038/363031a0. [2] M. Bologna, P. Grigolini, B.J. West, Strange kinetics: conflict between density and trajectory description, Chem. Phys. 284 (2002) 115–128, http://dx.doi.org/10.1016/S0301-0104(02)00543-8. [3] C. Ingo, R.L. Magin, L. Colon-Perez, W. Triplett, T.H. Mareci, On random walks and entropy in diffusion-weighted magnetic resonance imaging studies of neural tissue, Magn. Reson. Med. 71 (2014) 617–627, http://dx.doi.org/10.1002/mrm.24706. [4] A.A. Dubkov, B. Spagnolo, V.V. Uchaikin, Lévy flight superdiffusion: an introduction, Int. J. Bifurc. Chaos 18 (9) (2008) 2649–2672. [5] E. Barkai, Y. Garini, R. Metzler, Strange kinetics of single molecules in living cells, Phys. Today 65 (2012) 29, http://dx.doi.org/10.1063/PT.3.1677. [6] F. Höfling, T. Franosch, Anomalous transport in the crowded world of biological cells, Rep. Prog. Phys. 76 (4) (2013) 046602, http://dx.doi.org/ 10.1088/0034-4885/76/4/046602. [7] S.C. Lim, S.V. Muniandy, Self-similar Gaussian processes for modeling anomalous diffusion, Phys. Rev. E 66 (2002) 021114, http://dx.doi.org/ 10.1103/PhysRevE.66.021114. [8] F. Thiel, I.M. Sokolov, Scaled Brownian motion as a mean-field model for continuous-time random walks, Phys. Rev. E 89 (2014) 012115, http://dx.doi.org/10.1103/PhysRevE.89.012115. [9] J. Jeon, A.V. Chechkin, R. Metzler, Scaled Brownian motion: a paradoxical process with a time dependent diffusivity for the description of anomalous diffusion, Phys. Chem. Chem. Phys. 16 (30) (2014) 15811–15817, http://dx.doi.org/10.1039/C4CP02019G. [10] K.C. Leptos, J.S. Guasto, J.P. Gollub, A.I. Pesei, R.E. Goldstein, Dynamics of enhanced tracer diffusion in suspensions of swimming eukaryotic microorganisms, Phys. Rev. Lett. 103 (2009) 198103, http://dx.doi.org/10.1103/PhysRevLett.103.198103. [11] B. Eckhardt, S. Zammert, Non-normal tracer diffusion from stirring by swimming microorganisms, Eur. Phys. J. E 35 (2012) 96, http://dx.doi.org/ 10.1140/epje/i2012-12096-7. [12] E.W. Montroll, G.H. Weiss, Random walks on lattices. II, J. Math. Phys. 6 (1965) 167, http://dx.doi.org/10.1063/1.1704269. [13] I.W. Zaid, J. Dunkel, J.M. Yeomans, Lévy fluctuations and mixing in dilute suspensions of algae and bacteria, J. R. Soc. Interface 8 (2011) 1314–1331, http://dx.doi.org/10.1098/rsif.2010.0545. [14] P. Pramukkul, A. Svenkeson, P. Grigolini, M. Bologna, B.J. West, Complexity and the fractional calculus, Adv. Math. Phys. 2013 (2013) 498789, http://dx.doi.org/10.1155/2013/498789. [15] A. Svenkeson, M.T. Beig, M. Turalska, B.J. West, P. Grigolini, Fractional trajectories: decorrelation versus friction, Physica A 392 (2013) 5663–5672, http://dx.doi.org/10.1016/j.physa.2013.07.028. [16] M. Bologna, B.J. West, P. Grigolini, Renewal and memory origin of anomalous diffusion: a discussion of their joint action, Phys. Rev. E 88 (2013) 062106, http://dx.doi.org/10.1103/PhysRevE.88.062106. [17] Y. Meroz, I.M. Sokolov, J. Klafter, Test for determining a subdiffusive model in ergodic systems from single trajectories, Phys. Rev. Lett. 110 (2013) 090601, http://dx.doi.org/10.1103/PhysRevLett.110.090601. [18] F. Mainardi, Y. Luchko, G. Pagnini, The fundamental solution of the space–time fractional diffusion equation, Fract. Calc. Appl. Anal. 4 (2) (2001) 153–192. [19] R. Gorenflo, F. Mainardi, D. Moretti, G. Pagnini, P. Paradisi, Discrete random walk models for space–time fractional diffusion, Chem. Phys. 284 (2002) 521–541, http://dx.doi.org/10.1016/S0301-0104(02)00714-0. [20] M. Bologna, Y. Ahat, B.J. West, P. Grigolini, Can intermittent long-range jumps of a random walker compensate for lethargy?, J. Phys. A, Math. Theor. 44 (15) (2011) 152003, http://dx.doi.org/10.1088/1751-8113/44/15/152003. [21] J.H.P. Schulz, E. Barkai, R. Metzler, Aging effects and population splitting in single-particle trajectory averages, Phys. Rev. Lett. 110 (2013) 020602, http://dx.doi.org/10.1103/PhysRevLett.110.020602. [22] B.J. West, V. Seshadri, Linear systems with Lévy fluctuations, Physica A 113 (1982) 203–216, http://dx.doi.org/10.1016/0378-4371(82)90015-2. [23] R. Cakir, P. Grigolini, A.A. Krokhin, Dynamical origin of memory and renewal, Phys. Rev. E 74 (2006) 021108, http://dx.doi.org/10.1103/ PhysRevE.74.021108. [24] W. Deng, E. Barkai, Ergodic properties of fractional Brownian–Langevin motion, Phys. Rev. E 79 (2009) 011112, http://dx.doi.org/10.1103/ PhysRevE.79.011112. [25] M. Annunziato, P. Grigolini, J. Riccardi, Fluctuation–dissipation process without a time scale, Phys. Rev. E 61 (2000) 4801, http://dx.doi.org/ 10.1103/PhysRevE.61.4801. [26] J.-H. Jeon, A.V. Chechkin, R. Metzler, Scaled Brownian motion: a paradoxical process with a time dependent diffusivity for the description of anomalous diffusion, arXiv:1405.2193 [cond-mat.stat-mech], 2014. [27] G.K. Batchelor, Diffusion in a field of homogeneous turbulence, Math. Proc. Camb. Philos. Soc. 48 (1952) 345–362, http://dx.doi.org/10.1017/ S0305004100027687. [28] L.F. Richardson, Atmospheric diffusion shown on a distance-neighbor graph, Proc. R. Soc. Lond. A 110 (1926) 709–737, http://dx.doi.org/ 10.1098/rspa.1926.0043. [29] M.F. Shlesinger, B.J. West, J. Klafter, Lévy dynamics of enhanced diffusion: application to turbulence, Phys. Rev. Lett. 58 (1987) 1100, http://dx.doi.org/ 10.1103/PhysRevLett.58.1100. [30] J. Jeon, V. Tejedor, S. Burov, E. Barkai, C. Selhuber-Unkel, K. Berg-Sørensen, L. Oddershede, R. Metzler, In vivo anomalous diffusion and weak ergodicity breaking of lipid granules, Phys. Rev. Lett. 106 (2011) 048103, http://dx.doi.org/10.1103/PhysRevLett.106.048103. [31] S. Burov, J. Jeon, R. Metzler, E. Barkai, Single particle tracking in systems showing anomalous diffusion: the role of weak ergodicity breaking, Phys. Chem. Chem. Phys. 13 (2011) 1800–1812, http://dx.doi.org/10.1039/C0CP01879A. [32] A.V. Weigel, B. Simon, M.M. Tamkun, D. Krapf, Ergodic and nonergodic processes coexist in the plasma membrane as observed by single-molecule tracking, Proc. Natl. Acad. Sci. USA 108 (2011) 6438, http://dx.doi.org/10.1073/pnas.1016325108. [33] S.M.A. Tabei, S. Burov, H.Y. Kim, A. Kuznetsov, T. Huynh, J. Jureller, L.H. Philipson, A.R. Dinner, N.F. Scherer, Intracellular transport of insulin granules is a subordinated random walk, Proc. Natl. Acad. Sci. USA 110 (2013) 4911–4916, http://dx.doi.org/10.1073/pnas.1221962110. [34] J.H.P. Schulz, E. Barkai, R. Metzler, Aging renewal theory and application to random walks, Phys. Rev. X 4 (2014) 011028, http://dx.doi.org/ 10.1103/PhysRevX.4.011028. [35] B.B. Mandelbrot, J.W. van Ness, Fractional Brownian motions, fractional noises and applications, SIAM Rev. 10 (4) (1968) 422–437. [36] H.E. Hurst, Long term storage capacity of reservoirs, Trans. Am. Soc. Civ. Eng. 116 (1951) 770–779. [37] M. Magdziarz, A. Weron, J. Klafter, Equivalence of the fractional Fokker–Planck and subordinated Langevin equations: the case of a time-dependent force, Phys. Rev. Lett. 101 (2008) 210601, http://dx.doi.org/10.1103/PhysRevLett.101.210601. [38] M. Magdziarz, J. Gajda, T. Zorawik, Comment on fractional Fokker–Planck equation with space and time dependent drift and diffusion, J. Stat. Phys. 154 (2014) 1241–1250, http://dx.doi.org/10.1007/s10955-014-0919-9. [39] M. Rivero, L. Rodríguez-Germá, J.J. Trujillo, Linear fractional differential equations with variable coefficients, Appl. Math. Lett. 21 (2008) 892–897, http://dx.doi.org/10.1016/j.aml.2007.09.010. [40] A.A. Stanislavsky, Probability interpretation of the integral of fractional order, Theor. Math. Phys. 138 (2004) 418–431, http://dx.doi.org/10.1023/ B:TAMP.0000018457.70786.36.

JID:YJCPH AID:5437 /FLA

[m3G; v 1.135; Prn:3/09/2014; 9:18] P.15 (1-15)

M. Bologna et al. / Journal of Computational Physics ••• (••••) •••–•••

15

[41] I.M. Sokolov, Lévy flights from a continuous-time process, Phys. Rev. E 63 (2000) 011104, http://dx.doi.org/10.1103/PhysRevE.63.011104. [42] B.J. West, P. Grigolini, Complex Webs: Anticipating the Improbable, Cambridge University Press, Cambridge, UK, 2011. [43] B.I. Henry, T.A.M. Langlands, P. Straka, Fractional Fokker–Planck equations for subdiffusion with space- and time-dependent forces, Phys. Rev. Lett. 105 (2010) 170602, http://dx.doi.org/10.1103/PhysRevLett.105.170602. [44] E. Heinsalu, M. Patriarca, I. Goychuk, P. Hanggi, Use and abuse of a fractional Fokker–Planck dynamics for time-dependent driving, Phys. Rev. Lett. 99 (2007) 120602, http://dx.doi.org/10.1103/PhysRevLett.99.120602. [45] A. Stanislavsky, K. Weron, A study of diffusion under a time-dependent force by time subordination, J. Stat. Mech. Theory Exp. 2012 (07) (2012) P07020, http://dx.doi.org/10.1088/1742-5468/2012/07/P07020. [46] K. Diethelm, N.J. Ford, A.D. Freed, A predictor–corrector approach for the numerical solution of fractional differential equations, Nonlinear Dyn. 29 (2002) 3–22, http://dx.doi.org/10.1023/A:1016592219341. [47] K. Diethelm, A.D. Freed, The FracPECE subroutine for the numerical solution of differential equations of fractional order, Forsch. Wissenschaftlich. Rechn. 1999 (1998) 57–71. [48] H.C. Fogedby, Langevin equations for continuous time Lévy flights, Phys. Rev. E 50 (1994) 1657, http://dx.doi.org/10.1103/PhysRevE.50.1657.