Numerical Method for Solving of the Anomalous Diffusion Equation Based on a Local Estimate of the Monte Carlo Method

Abstract

This paper considers a method of stochastic solution to the anomalous diffusion equation with a fractional derivative with respect to both time and coordinates. To this end, the process of a random walk of a particle is considered, and a master equation describing the distribution of particles is obtained. It has been shown that in the asymptotics of large times, this process is described by the equation of anomalous diffusion, with a fractional derivative in both time and coordinates. The method has been proposed for local estimation of the solution to the anomalous diffusion equation based on the simulation of random walk trajectories of a particle. The advantage of the proposed method is the opportunity to estimate the solution directly at a given point. This excludes the systematic component of the error from the calculation results and allows constructing the solution as a smooth function of the coordinate.

1. Introduction

Anomalous diffusion processes are characterized by a power-law dependence of the width of the diffusion packet on time $\Delta \left(t\right)\propto D_{\alpha }{t}^{\gamma }$, where $D_{\alpha }$ is the diffusion coefficient [1,2,3,4]. Depending on the value of the exponent $\gamma$, different diffusion regimes are distinguished: $\gamma <1/2$ (a sub-diffusion), $\gamma =1/2$ (a normal diffusion), $\gamma >1/2$ (a super-diffusion). If $\gamma =1$, then the quasi-ballistic regime is established, if $\gamma >1$, then this regime bears the name super-ballistic. More detailed information about various regimes can be found in the works [5,6,7,8].

The Continuous Time Random Walk (CTRW) model underlies the model of anomalous diffusion [9]. This model was first introduced in the work [10] and was further developed in the works [4,11,12,13] (see also survey articles [9,14,15]). The CTRW model describes the random walk of a particle using a hopping-trap mechanism. The random walk of a particle represents a sequence of instantaneous random jumps of the value $R_i$, $i=1,2,3,\dots$ and rest states with random rest times $T_i$, $i=1,2,3,\dots$, which successively change one another. With a power-law distribution of the free path value $p\left(x\right)\propto x^{-\alpha -1},x\to \infty ,0<\alpha ⩽2$ and rest time $q\left(t\right)\propto t^{-\beta -1},t\to \infty ,0<\beta ⩽1$, the width of the diffusion packet will grow according to the power law $\Delta \left(t\right)\propto D_{\alpha }{t}^{\gamma }$. Within this model, normal diffusion is obtained if the probability density of the path distribution $p\left(x\right)$ has a finite second moment $\langle R^2\rangle <\infty$ [6], and the probability density of the rest time distribution $q\left(t\right)$ has the finite mathematical expectation $\langle T\rangle <\infty$ [7].

The asymptotic distribution of particles in the CTRW process was first obtained by M. Kotulski in the article [16]. Later, independently from M. Kotulski, these distributions were obtained in the article [17] in which these distribution were called fractional stable distributions. It is well known that the anomalous diffusion equation is an asymptotic description of the CTRW process. The articles [18,19,20] are devoted to the solution of the anomalous diffusion equation. These papers show that the solution of the anomalous diffusion equation is expressed through the classes of the fractional stable and stable distributions.

Anomalous diffusion generalizes normal diffusion to the case of considering transport processes in inhomogeneous and turbulent media. The model of anomalous diffusion was first used to describe charge transfer in amorphous semiconductors [10,11,12,13]. Later, the model of anomalous diffusion became widespread in the description of transport processes in turbulent plasma [21,22,23,24,25,26,27], propagation of cosmic rays in the galaxy [28,29,30,31,32], studying the diffusion of microRNA in the cell [33], the fluctuation of prices on exchanges and currency exchange rates [34]. The theory of combustion is one of the few areas where anomalous diffusion has not become widespread yet. However, recently, this direction has been intensively developing. The authors of [35,36] point out that anomalous diffusion occurs during heat transfer in low-dimensional systems. In the paper [37] devoted to the study of thermal radiation during the combustion of natural gas and acetylene, it was found that the combustion process was of a subdiffusion nature.

The assumption about the formation of anomalous diffusion in the combustion process allows introducing fractional differential equations of diffusion into consideration. In the papers [38,39,40], the effective thermal conductivity coefficient was obtained for the Levy–Fokker–Planck and fractional Boltzmann equations. The authors of [41,42,43] propose to use the fractional differential equation of diffusion to describe the combustion process

$${}_0{D}_t^{\beta }u={\partial }_{xx}u+s\left(u\right),t>0,0<x<L,$$

where ${}_0{D}_t^{\beta }$ is the fractional Riemann-Liouville derivative (Appendix A, A2) of the order $0<\beta <1$ in terms of time and ${\partial }_{xx}$ is the classical second-order partial derivative with respect to the coordinate. The paper considers the problem when the source $f\left(u\right)$ is singular, and the initial and boundary conditions are chosen in the form $u(x,0)=u_0\left(x\right)$, $u(0,t)=u(L,t)=0$. The paper [41] examines the damping effect in the framework of the investigated model, and the paper [42] explores the phenomenon of explosion and the possibility of describing this phenomenon with the help of the passage to the limit $\beta \to 0$. The authors of [44] examine a two-dimensional combustion model with a fractional time derivative. To solve the fractional-differential diffusion equation, the authors develop an adaptive finite-difference discontinuous Galerkin method. In the paper [45], the authors consider a fractional-differential combustion model with the first derivative with respect to time and a fractional derivative with respect to the spatial variable.

$${\partial }_{t}u=\frac{{\partial }^{\alpha }}{{\partial \left|x\right|}^{\alpha }}u+s\left(u\right),u(x,0)=u_0\left(x\right),u(a,t)=u(b,t)=0.$$

Here, ${\partial }_t$ is the partial time derivative, and $\frac{{\partial }^{\alpha }}{{\partial \left|x\right|}^{\alpha }}$ is the fractional-differential Riesz operator (A4). Using this fractional-differential model of combustion, the paper investigates the effect on the damping phenomenon of a quantity of the order of the fractional derivative $\alpha$, the spatial size of the area under study, and the initial conditions. To solve the fractional-differential combustion equation, finite-difference methods of solving the equations are used. In the article [46], a fractional-differential generalization of the kinetic equation was obtained that describes the dependence of the radius of the ball on time in the model of combustion of a fireball, which was theoretically predicted by the Soviet physicist Ya.B. Zeldovich [47].

As we can see, the use of the anomalous diffusion model and the fractional-differential generalization of the diffusion equation for modeling combustion processes is only in the initial state. One of the reasons for this is the complexity of solving such kind of equations. Finite-difference methods were used to solve fractional-differential diffusion equations in the works considered above [41,42,43,45]. However, in these works, fractional differential equations are studied in which only one of the derivatives has a fractional order: the time derivative or the coordinate derivative. In this paper, we will propose a method for the numerical solution to the anomalous diffusion equation with a fractional derivative in both time and coordinate and with a source of a special type.

$${}_0{D}_t^{\beta }\rho (x,t)=D\frac{{\partial }^{\alpha }}{{\partial \left|x\right|}^{\alpha }}\rho (x,t)+\frac{t^{-\beta }}{\mathrm{\Gamma }(1-\beta )}\delta \left(x\right),$$

with boundary conditions $\rho (x,t)\to 0$ if $x\to \pm \infty$ and $\rho (x,t)=0$, if $t<0$. Here, $0<\beta ⩽1$, $0<\alpha ⩽2$, D is the generalized diffusion coefficient.

2. Master Equation of the CTRW Process

To obtain the master equation of the CTRW process, we use the approach proposed in the paper [48] and further developed in the paper [49]. Consider the collision density $f(\mathbf{r},\mathbf{p},t)$, where $\mathbf{r}$ is the radius of the particle vector, $\mathbf{p}$ is the particle impulse, and t is time. The value $f(\mathbf{r},\mathbf{p},t)d\mathbf{r}d\mathbf{p}dt$ is the number of collisions in the volume element, and $d\mathbf{r}$ is the vicinity of the point $\mathbf{r}$ for the time interval $dt$, at which the particle impulse takes the value from $\mathbf{p}$ to $\mathbf{p}+d\mathbf{p}$. We will consider the nonrelativistic case $\mathbf{p}=m\mathbf{v}$. Without losing generality, we will assume that $m=1$. It was shown in the paper [49] that if there are n discrete states, the value $f(\mathbf{r},\mathbf{v},t)$ can be presented in the form

$$f(\mathbf{r},\mathbf{v},t)=\sum _{j=1}^n{f}_j(\mathbf{r},\mathbf{v},t),$$

(1)

where

$$\begin{array}{c}{f}_j(\mathbf{r},\mathbf{v},t)=s_j(\mathbf{r},\mathbf{v},t)+\sum _{i=1}^n{c}_{ij}{\int }_0^t{k}_i\left(\tau \right)d\tau \int W_{ij}({\Omega }^{\prime },\Omega )d{\Omega }^{\prime }\times \hfill \\ \hfill \int f_i\left(\mathbf{r}-v^{\prime }{\Omega }^{\prime }\tau ,v^{\prime }{\Omega }^{\prime },t-\tau \right)h_{ij}(v^{\prime },v)d{v}^{\prime }.\end{array}$$

(2)

Here, $k_i\left(\tau \right)$ is the probability density distribution of residence time in the state i and $c_{ij}$ are the probabilities of passing from the state i into the state j, $W_{ij}({\Omega }^{\prime },\Omega )$ is the density of probability of the fact that before the collision, the velocity had its direction ${\Omega }^{\prime }$; after the collision, the direction obtained the value $\Omega$, $h_{ij}(v^{\prime },v)$ is the probability density of changing the velocity from the value $v^{\prime }$ to v, $s_j(\mathbf{r},\mathbf{v},t)$ is the density of the sources of new particles in the state j, $\mathbf{v}=v\Omega$, $v=\left|\mathbf{v}\right|$, $d\mathbf{v}=dvd\Omega$, the summation is carried out over all possible previous states. The quantities $c_{ij}$, $W_{ij}({\Omega }^{\prime },\Omega )$, and $h_{ij}(v^{\prime },v)$ have normalization

$$\sum _{j=1}^n{c}_{ij}=1,\int W_{ij}({\Omega }^{\prime },\Omega )d\Omega =1,\int h_{ij}(v^{\prime },v)dv=1.$$

(3)

The passage from the collision density $f(\mathbf{r},\mathbf{v},t)$ to the phase density $\psi (\mathbf{r},\mathbf{v},t)$ is carried out with the help of the integral

$$\psi (\mathbf{r},\mathbf{v},t)={\int }_0^{t}K\left(\tau \right)f\left(\mathbf{r}-v\Omega \tau ,v\Omega ,t-\tau \right)d\tau ,$$

(4)

where $K\left(t\right)={\int }_t^{\infty }k\left(\tau \right)d\tau$. Substituting (1) in (4), we obtain that the phase density has the form of the sum

$$\psi (\mathbf{r},\mathbf{v},t)=\sum _{j=1}^n{\psi }_j(\mathbf{r},\mathbf{v},t),$$

(5)

where

$${\psi }_j(\mathbf{r},\mathbf{v},t)={\int }_0^t{K}_j\left(\tau \right)f_j\left(\mathbf{r}-v\Omega \tau ,v\Omega ,t-\tau \right)d\tau ,$$

(6)

where $K_j\left(t\right)={\int }_t^{\infty }{k}_j\left(\tau \right)d\tau$. The physical interpretation of the latter expression is simple. To detect a particle in the state j of the vicinity $d\mathbf{r}$ of the point $\mathbf{r}$ with the velocity in the interval from $\mathbf{v}$ to $\mathbf{v}+d\mathbf{v}$ at a moment of time from t to $t+dt$, the particle is supposed to pass into this state in the point $\mathbf{r}-v\Omega \tau$ at the moment of time $t-\tau$ and stay in this state some time more than $\tau$. The passage to the particle density $\rho (\mathbf{r},t)$ is carried out with the use of the integral

$$\rho (\mathbf{r},t)=\int \psi (\mathbf{r},\mathbf{v},t)d\mathbf{v}.$$

(7)

The system of Equations (2), (5) and (6) together with conditions (3) describe practically any process of random walk with n discrete states under fairly general assumptions about the scattering indicatrix $W_{ij}({\Omega }^{\prime },\Omega )$ and the velocity redistribution law $h_{ij}(v^{\prime },v)$. In this study, using these equations, we will obtain master equations that describe the CTRW process.

The CTRW process is determined in the following way. At random times $T_1$, $T_1+T_2$, $T_1+T_2+\dots +T_j$, the particle makes instantaneous jumps with the value $R_1,R_2,\dots ,R_j$. Random jumps $R_i$ and random times $T_i$ are independent between one another as well as between each other. Thus, in the CTRW process, there are two states: $j=1$ is the state of rest and $j=2$ is the state of motion. By definition, in the CTRW process, after each jump, the particle enters a state of rest, and after each state of rest, the particle makes a jump. This means that the probabilities of passing $c_{ij}$ from the state i into the state j have the following values

$$c_{11}=0,c_{12}=1,c_{21}=1,c_{22}=0.$$

(8)

The state of rest is consistent with the velocity $v=0$, and infinite velocity $v_2=\infty$ corresponds to the state of motion (instant jump). To write the equation for $f_2(\mathbf{r},\mathbf{v},t)$, we first assume that in state 2, the particle moves with some constant velocity $v_2$, and then, we will carry out the passage $v_2\to \infty$. Taking account of the foregoing, we have

$$h_{21}(v^{\prime },v)\equiv h_1\left(v\right)=\delta \left(v\right),h_{12}(v^{\prime },v)\equiv h_2\left(v\right)=\delta (v-v_2),$$

(9)

$\delta \left(x\right)$ is the Dirac delta function. Furthermore, we assume that the direction of motion of the velocity after each collision does not depend on the previous direction of motion

$$W_{12}({\Omega }^{\prime },\Omega )=W_{21}({\Omega }^{\prime },\Omega )=W(\Omega ).$$

(10)

We represent the density of sources for each of the states in the form

$$s_j(\mathbf{r},\mathbf{v},t)={\sigma }_j{s}_j(\mathbf{r},t)h_j^s\left(v\right)W_1(\Omega ),j=1,2.$$

(11)

Here, ${\sigma }_j$ is the probability of the birth of a particle in the state j$({\sigma }_1+{\sigma }_2=1)$, $h_j^s\left(v\right)$ is the initial velocity modulus distribution, $W_j(\Omega )$ is the velocity direction distribution, and $s_j(\mathbf{r},t)$ is the spatio-temporal density of a source distribution. Without loss of generality, let us assume that a particle begins its history from a state of rest. Thus, without loss of generality, we assume that the particle starts its history from the state of rest. Thus,

$${\sigma }_1=1,{\sigma }_2=0,h_1^s\left(v\right)=\delta \left(v\right),h_2^s\left(v\right)=\delta (v-v_2),W_1(\Omega )=W_2(\Omega )=W(\Omega ).$$

(12)

Taking account of (8)–(12) for the collision density $f_j(\mathbf{r},\mathbf{p},t)$, we obtain the system of two equations

$$\begin{array}{cc}\hfill f_1(\mathbf{r},\mathbf{v},t)& =W(\Omega )\delta \left(v\right)\left(s_1(\mathbf{r},t)+{\int }_0^t{k}_2\left(\tau \right)d\tau \int d{\Omega }^{\prime }\int f_2(\mathbf{r}-v^{\prime }{\Omega }^{\prime }\tau ,v^{\prime }{\Omega }^{\prime },t-\tau )d{v}^{\prime }\right),\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill f_2(\mathbf{r},\mathbf{v},t)& =W(\Omega )\delta (v-v_2){\int }_0^t{k}_1\left(\tau \right)d\tau \int d{\Omega }^{\prime }\int f_1(\mathbf{r}-v^{\prime }{\Omega }^{\prime }\tau ,v^{\prime }{\Omega }^{\prime },t-\tau )d{v}^{\prime }.\hfill \end{array}$$

(14)

From this relation, we can see that in the case when the transition probability densities $h_{ij}(v^{\prime },v)$ and $W_{ij}({\Omega }^{\prime },\Omega )$ do not depend on the previous value $v^{\prime }$ and ${\Omega }^{\prime }$, then the collision density can be represented in the form of a product $f_j(\mathbf{r},v\Omega ,t)=W(\Omega )h_j\left(v\right)F_j(\mathbf{r},\mathbf{t})$, $j=1,2$, where $h_1\left(v\right)$ and $h_2\left(v\right)$ have the form (9). Thus, we obtain

$$f_1(\mathbf{r},v\Omega ,t)=W(\Omega )\delta \left(v\right)F_1(\mathbf{r},\mathbf{t}),f_2(\mathbf{r},v\Omega ,t)=W(\Omega )\delta (v-v_2)F_2(\mathbf{r},\mathbf{t}).$$

(15)

Now, substituting these relations in (13) and () performing the integration over $dvd\Omega$, for $F_j(\mathbf{r},t),j=1,2$ we get a system of equations

$$\begin{array}{c}{F}_1(\mathbf{r},t)=s_1(\mathbf{r},t)+{\int }_0^t{k}_2\left(\tau \right)d\tau \int W\left({\Omega }^{\prime }\right)d{\Omega }^{\prime }\int F_2(\mathbf{r}-v^{\prime }{\Omega }^{\prime }\tau ,t-\tau )\delta (v^{\prime }-v_2)d{v}^{\prime }=\hfill \\ \hfill s_1(\mathbf{r},t)+{\int }_0^t{k}_2\left(\tau \right)d\tau \int F_2(\mathbf{r}-v_2{\Omega }^{\prime }\tau ,t-\tau )W\left({\Omega }^{\prime }\right)d{\Omega }^{\prime },\end{array}$$

(16)

$$F_2(\mathbf{r},t)={\int }_0^t{k}_1\left(\tau \right)d\tau \int W\left({\Omega }^{\prime }\right)d{\Omega }^{\prime }\int F_1(\mathbf{r}-v^{\prime }{\Omega }^{\prime }\tau ,t-\tau )\delta \left(v^{\prime }\right)d{v}^{\prime }={\int }_0^t{k}_1\left(\tau \right)F_1(\mathbf{r},t-\tau )d\tau .$$

(17)

In the equation for $F_2(\mathbf{r},t)$, the condition of normalization $\int W\left({\Omega }^{\prime }\right)d{\Omega }^{\prime }=1$ was used. The physical meaning of quantities $F_j(\mathbf{r},t)$ is simple enough. This is the density of collisions in a volume element $d\mathbf{r}$ of the point vicinity $\mathbf{r}$, as a result of which the particle passes into the state j for a period of time $t,t+dt$ for any value of the particle velocity.

Let us now pass from the collision density to the phase density ${\psi }_j(\mathbf{r},\mathbf{v},t)$ and then to the density of particles ${\rho }_j(\mathbf{r},t)$. For this purpose, we substitute the expressions (15) in (6)

$$\begin{array}{cc}\hfill {\psi }_1(\mathbf{r},\mathbf{v},t)& ={\int }_0^t{K}_1\left(\tau \right)W(\Omega )\delta \left(v\right)F_1(\mathbf{r}-v\Omega \tau ,t-\tau )d\tau ,\hfill \\ \hfill {\psi }_2(\mathbf{r},\mathbf{v},t)& ={\int }_0^t{K}_2\left(\tau \right)W(\Omega )\delta (v-v_2)F_2(\mathbf{r}-v\Omega \tau ,t-\tau )d\tau .\hfill \end{array}$$

Now, integrating each of these expressions over $dvd\Omega$ in view of (5) and (7) for the particle density, we obtain the system of equations

$$\begin{array}{cc}\hfill \rho (\mathbf{r},t)& ={\rho }_1(\mathbf{r},t)+{\rho }_2(\mathbf{r},t),\hfill \\ \hfill {\rho }_1(\mathbf{r},t)& ={\int }_0^t{K}_1\left(\tau \right)F_1(\mathbf{r},t-\tau )d\tau ,\hfill \\ \hfill {\rho }_2(\mathbf{r},t)& ={\int }_0^t{K}_2\left(\tau \right)d\tau \int W(\Omega )F_2(\mathbf{r}-v_2\Omega \tau ,t-\tau )d\Omega ,\hfill \end{array}$$

(18)

where $F_1(\mathbf{r},t)$ and $F_2(\mathbf{r},t)$ are determined by Equations (16) and (17). As one can see, the total particle density is the sum of the particle density in each of the states. It should be noted that this system of equations is not yet a system of master equations to describe the CTRW process. This system describes the random walk of a particle with two states: state 1 is rest, and state 2 is motion with a constant final velocity $v_2$.

To obtain the system of equations for the CTRW process, it is necessary to go to the limit $v_2\to \infty$. It should be pointed out that in Equations (16) and (18), $\tau$ is an independent variable, which has the meaning of the time in the state of motion. Therefore, if we pass now to the limit $v_2\to \infty$, then it means that $\tau \to 0$ and the probability is $k_2\left(\tau \right)d\tau =0$. To avoid this, it is necessary to pass to a new variable $\xi =v\tau$, which has the meaning of the free path of a particle. As a result, we obtain the system of equations

$$\begin{array}{cc}\hfill {\rho }_1(\mathbf{r},t)& ={\int }_0^{t}Q\left(\tau \right)F_1(\mathbf{r},t-\tau )d\tau ,\hfill \\ \hfill {\rho }_2(\mathbf{r},t)& =\frac{1}{v_2}{\int }_0^{v_{2}t}P\left(\xi \right)d\xi \int W(\Omega )F_2(\mathbf{r}-\xi \Omega ,t-\xi /v_2)d\Omega ,\hfill \\ \hfill F_1(\mathbf{r},t)& =s_1(\mathbf{r},t)+{\int }_0^{v_{2}t}p\left(\xi \right)d\xi \int F_2(\mathbf{r}-\xi {\Omega }^{\prime },t-\xi /v_2)W\left({\Omega }^{\prime }\right)d{\Omega }^{\prime },\hfill \\ \hfill F_2(\mathbf{r},t)& ={\int }_0^{t}q\left(\tau \right)F_1(\mathbf{r},t-\tau )d\tau .\hfill \end{array}$$

Here, for convenience, the following notation was introduced $p\left(\xi \right)=\frac{1}{v_2}{k}_2\left(\frac{\xi }{v_2}\right)$, $P\left(\xi \right)\equiv K_2\left(\frac{\xi }{v_2}\right)={\int }_{\xi }^{\infty }p\left(y\right)dy$, $q\left(\tau \right)\equiv k_1\left(\tau \right)$, $Q\left(\tau \right)\equiv K_1\left(\tau \right)$.

Let us now pass in this system of equations to the limit $v_2\to \infty$. From the equation for ${\rho }_2(\mathbf{r},t)$, it is seen that due to the presence of the multiplier $1/v_2$ before the sign of the integral ${\rho }_2(\mathbf{r},t)\to 0$ at $v_2\to \infty$. There is a simple explanation for this fact. Since now the particle is moving with infinite velocity, it instantly moves from one point in space to another. As a result, the probability of detecting a particle in a state of motion in the time interval $t+dt$ is equal to zero and as a sequence, ${\rho }_2(\mathbf{r},t)=0$. Thus, for the CTRW process, we arrive at the system of equations

$$\begin{array}{cc}\hfill \rho (\mathbf{r},t)& ={\rho }_1(\mathbf{r},t)={\int }_0^{t}Q\left(\tau \right)F_1(\mathbf{r},t-\tau )d\tau ,\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill F_1(\mathbf{r},t)& =s_1(\mathbf{r},t)+{\int }_0^{\infty }p\left(\xi \right)d\xi \int F_2(\mathbf{r}-\xi {\Omega }^{\prime },t)W\left({\Omega }^{\prime }\right)d{\Omega }^{\prime },\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill F_2(\mathbf{r},t)& ={\int }_0^{t}q\left(\tau \right)F_1(\mathbf{r},t-\tau )d\tau .\hfill \end{array}$$

(21)

As one can see, for the CTRW process, the particle distribution density $\rho (\mathbf{r},t)$ only the density of the spatial distribution of particles at rest is determined.

The obtained system of equations can be simplified, and it is possible to pass to the equation of the particle density only $\rho (\mathbf{r},t)$. For this purpose, we need to put Equation () into Equation (). Then, we need to put the obtained equation for the collision density $F_1(\mathbf{r},t)$ into Equation (19) and change the integration order in terms of time. As a result, we get the equation for density

$$\rho (\mathbf{r},t)={\int }_0^{t}Q\left(\tau \right)s_1(\mathbf{r},t)d\tau +{\int }_0^{t}q\left(\tau \right)d\tau {\int }_0^{\infty }p\left(\xi \right)d\xi \int W(\Omega )\rho (\mathbf{r}-\xi \Omega ,t-\tau )d\Omega .$$

(22)

This equation is the master equation of the CTRW process in three-dimensional space. The result obtained coincides with similar results obtained in the works [4,14,48].

Let us simplify the problem and consider the one-dimensional case. Let the random walk process occur along the axis x. In this case, the function $W(\Omega )$ takes the form

$$W(\Omega )=W(\theta ,\phi )=\frac{1}{sin\theta }({\omega }_1\delta \left(\phi \right)+{\omega }_2\delta (\phi -\pi ))\delta (\theta -\pi /2),$$

(23)

where ${\omega }_1$ and ${\omega }_2$ are probabilities of motion in the positive and negative direction of the axis $Ox$ correspondingly and ${\omega }_1+{\omega }_2=1$. Substituting now (23) in (22) and taking account that $\Omega =(sin\theta cos\phi ,sin\theta sin\phi ,cos\theta )$, $\rho (\mathbf{r},t)=\rho (x,y,z,t)$, $s_1(\mathbf{r},t)=s_1(x,y,z,t)$, $d\Omega =sin\theta d\theta d\phi$ and integrating the resulting equation over the angular variables, we obtain

$$\begin{array}{c}\rho (x,y,z,t)={\int }_0^{t}Q\left(\tau \right)s_1(x,y,z,t)d\tau \hfill \\ \hfill +{\int }_0^{t}q\left(\tau \right)d\tau {\int }_0^{\infty }p\left(\xi \right)({\omega }_1\rho (x-\xi ,y,z,t-\tau )+{\omega }_2\rho (x+\xi ,y,z,t-\tau ))d\xi .\end{array}$$

Since we are considering random walks along the axis $Ox$, then

$$\rho (x,y,z,t)=\rho (x,t)\delta \left(y\right)\delta \left(z\right),s_1(x,y,z,t)=s_1(x,t)\delta \left(y\right)\delta \left(z\right).$$

Substituting now these expressions into the previous equation and integrating over $dydz$, we finally get

$$\rho (x,t)={\int }_0^{t}Q\left(\tau \right)s_1(x,t)d\tau +{\int }_0^{t}q\left(\tau \right)d\tau {\int }_0^{\infty }p\left(\xi \right)({\omega }_1\rho (x-\xi ,t-\tau )+{\omega }_2\rho (x+\xi ,t-\tau ))d\xi .$$

(24)

This equation describes one-dimensional walks of a particle in the CTRW process. The solution to this equation will be sought by the standard method of the Fourier–Laplace transform. Performing the Fourier–Laplace transform

$$\widehat{\rho }(k,\lambda )={\int }_0^{\infty }dt{\int }_{-\infty }^{\infty }{e}^{ikx-\lambda t}\rho (x,t)dx$$

we get that the solution to Equation (24) has the form

$$\widehat{\tilde{\rho }}(k,\lambda )=\frac{1-\tilde{q}\left(\lambda \right)}{\lambda }\frac{{\widehat{\tilde{s}}}_1(k,\lambda )}{1-\tilde{q}\left(\lambda \right)({\omega }_1\widehat{p}\left(k\right)+{\omega }_2\widehat{p}(-k))}.$$

(25)

Here, $\widehat{q}\left(\lambda \right)$ is the Laplace transform of density $q\left(t\right)$, $\widehat{p}\left(k\right)$ is the Fourier transform of density $p\left(\xi \right)$, ${\widehat{s}}_1(k,\lambda )$ is the Fourier–Laplace transform of source function $s_1(x,t)$ and

$${\int }_0^{\infty }{e}^{-\lambda t}Q\left(t\right)dt=\frac{1-\widehat{q}\left(\lambda \right)}{\lambda }.$$

3. The Equations of Anomalous Diffusion

The obtained Equation (24) is an exact description of the walk process, and its solution in Fourier–Laplace images (25) is its exact solution. However, it is possible to perform the inverse Fourier–Laplace transform of the solution (25) and to write down an equation describing the process of a walk, only if to consider the asymptotic solution at $t\to \infty$ and $x\to \infty$. The form of this equation is determined by two characteristics of the distribution of free path and rest time: mathematical expectation and variance [6,7,14]. If the variance of free paths $\langle R^2\rangle$ and rest times $\langle T^2\rangle$ are finite, then the random walk process is described by a standard diffusion equation. If the variance of rest time is finite $\langle T^2\rangle <\infty$, the mathematical expectation of free paths is finite ($\langle R\rangle <\infty$ and $\langle R^2\rangle =\infty$) or infinite ($\langle R\rangle =\infty$), then the random walk process is described with the anomalous diffusion equation with the first derivative with respect to time and fractional derivative with respect to the coordinate. In the case $\langle T\rangle =\infty$ and $\langle R^2\rangle <\infty$, we obtain the equation of anomalous diffusion with a fractional time derivative and a second coordinate derivative. In the case $\langle T\rangle =\infty$ and $\langle R^2\rangle =\infty$, then the random walk process is described with a fractional derivative equation in both time and coordinate. Based on the foregoing, we will consider all these cases separately.

At first, we consider the simplest case $\langle T\rangle <\infty$ and $\langle R^2\rangle <\infty$. Let rest times and free paths have the exponential distribution

$$\begin{array}{cc}\hfill q\left(t\right)& =\nu e^{-\nu t},\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill p\left(x\right)& =\mu e^{-\mu x}.\hfill \end{array}$$

(27)

Performing the Laplace transform of the density $q\left(t\right)$ and the Fourier transform of the density $p\left(x\right)$, we obtain

$$\tilde{q}\left(\lambda \right)=\frac{\nu }{\nu +\lambda },\widehat{p}\left(k\right)=\frac{\mu }{\mu +ik}.$$

The asymptotic solution at $t\to \infty$ and $x\to \infty$ is of interest to us. According to Tauberian theorems, the behavior of the function at $t\to \infty$ or $x\to \infty$ corresponds to the behavior of the transformant at $\lambda \to 0$ or $k\to 0$. Expanding images $\tilde{q}\left(\lambda \right)$ and $\widehat{p}\left(k\right)$ in a series and leaving one summand in the expansion $\tilde{q}\left(\lambda \right)$ and two summands in the expansion $\widehat{p}\left(k\right)$, we obtain

$$\begin{array}{cc}\hfill \tilde{q}\left(\lambda \right)& \approx 1-\langle T\rangle \lambda ,\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill \widehat{p}\left(k\right)& \approx 1+i\langle R\rangle k-\frac{\langle R^2\rangle }{2}{k}^2.\hfill \end{array}$$

(29)

Here, we use the fact that $1/\nu =\langle T\rangle$, $1/\mu =\langle R\rangle$ and $1/{\mu }^2=\langle R^2\rangle$. We will consider symmetrical random walks (${\omega }_1={\omega }_2=1/2$) with the point instantaneous source $s_1(x,t)=\delta \left(x\right)\delta \left(t\right)$. We substitute the expansions (28) and (29) in the solution (25). As a result, we obtain

$$\widehat{\tilde{\rho }}(k,\lambda )=\frac{1}{\lambda +D_A{k}^2},D_A=\frac{\langle R^2\rangle }{2\langle T\rangle },$$

where $D_A$ is the diffusion coefficient. If we write this expression in the form

$$\lambda \widehat{\tilde{\rho }}(k,\lambda )-1=-D_A{k}^2\widehat{\tilde{\rho }}(k,\lambda )$$

(30)

and take into account that

$$\begin{array}{cc}\hfill {\int }_0^{\infty }{e}^{-\lambda t}\frac{\partial f\left(t\right)}{\partial t}dt& =\lambda \tilde{f}\left(\lambda \right)-{\tilde{f}}_0,\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill {\int }_{-\infty }^{\infty }{e}^{ikx}\frac{{\partial }^{2}f\left(x\right)}{\partial x^2}dx& =-k^2\widehat{f}\left(k\right),\hfill \end{array}$$

(32)

where ${\tilde{f}}_0$ is the Laplace transform of the initial condition, then it is clear that (30) is nothing but the Fourier–Laplace transform of the diffusion equation

$$\frac{\partial \rho (x,t)}{\partial t}=D_A\frac{{\partial }^2\rho (x,t)}{\partial x^2}$$

(33)

with zero boundary conditions at infinity and the initial condition $\rho (x,0)=\delta \left(x\right)$. As it is well known, the solution to this equation is expressed in terms of the normal distribution and has the form

$$\rho (x,t)={\left(4\pi D_{A}t\right)}^{-1/2}exp\{-x^2/\left(4D_{A}t\right)\}.$$

(34)

An important property of this solution is the self-semilarity of the density profile [9]. Self-similarity is understood as a special form of solution symmetry, which means that a change in the scale of some variables can be compensated by a change in the scale of other variables. In the case under consideration, the solution (34) can be represented in the form $\rho (x,t)={\left(2D_{A}t\right)}^{-1/2}g\left(x{\left(2D_{A}t\right)}^{-1/2}\right)$, where $g\left(x\right)={\left(2\pi \right)}^{-1/2}exp\{-x^2/2\}$ is the density of the normal law (or Gauss distribution). As we can see, the solution $\rho (x,t)$ at an arbitrary moment of time can be obtained from the density of normal law $g\left(x\right)$ by scaling transformation of coordinate and density. Thus, in the case under consideration, the width of diffusion packet grows according to the law $\Delta \left(t\right)\propto {\left(2D_{A}t\right)}^{1/2}$.

We consider the case $\langle T\rangle =\infty$ and $\langle R^2\rangle <\infty$. Let us assume that free paths have the distribution (27), and the rest times are distributed according to the law

$$q\left(t\right)=\left\{\begin{array}{cc}0,\hfill & t<t_0\hfill \\ \beta t_0^{\beta }{t}^{-\beta -1},\hfill & t⩾t_0,\hfill \end{array}\right.,0<\beta <1.$$

(35)

A characteristic property of this distribution is that this distribution has moments of order $k<\beta$, i.e., $\langle T^k\rangle <\infty$, if $k<\beta$, and $\langle T^k\rangle =\infty$, if $k⩾\beta$. Considering this circumstance, in the work [50], it was shown that the Laplace transform of density (35) had the form

$$\tilde{q}\left(\lambda \right)\approx 1-{\left(\lambda t_0\right)}^{\beta }\mathrm{\Gamma }(1-\beta ),\lambda \to 0,$$

(36)

where $\mathrm{\Gamma }\left(x\right)$ is Euler’s gamma function. Now, let us substitute this expansion as well as the expansion (29) in Equation (25). As in the previous case, we will consider symmetrical random walks $({\omega }_1={\omega }_2=1/2)$ with the point instantaneous source $s_1(x,t)=\delta \left(x\right)\delta \left(t\right)$. As a result, from (25), we obtain

$$\widehat{\tilde{\rho }}(k,\lambda )=\frac{{\lambda }^{\beta -1}}{{\lambda }^{\beta }+D_B{k}^2},D_B=\frac{\frac{1}{2}\langle R^2\rangle }{t_0^{\beta }\mathrm{\Gamma }(1-\beta )}.$$

We represent this equation in the form

$${\lambda }^{\beta }\widehat{\tilde{\rho }}(k,\lambda )=-D_B{k}^2\widehat{\tilde{\rho }}(k,\lambda )+{\lambda }^{\beta -1}.$$

To perform the inverse Fourier–Laplace transform of this expression, we will use the formulas (A3) and (32). As a result, the random walk process in the case under consideration is described with the fractional differential equation

$${}_0{D}_t^{\beta }\rho (x,t)=D_B\frac{{\partial }^2\rho (x,t)}{\partial x^2}+\frac{t^{-\beta }}{\mathrm{\Gamma }(1-\beta )}\delta \left(x\right)$$

with a diffusion coefficient $D_B$. In the work [18], it was shown that the solution to this equation has the form

$$\rho (x,t)={\left(D_B{t}^{\beta }\right)}^{-1/2}q\left(x{\left(D_B{t}^{\beta }\right)}^{-1/2};2,\beta \right),$$

where $q(x;\alpha ,\beta )$ is the density of a fractional-stable law [17,51,52]

$$q(x,\alpha ,\beta )={\int }_0^{\infty }g\left(x{y}^{\beta /\alpha };\alpha ,0\right)g(y;\beta ,1)y^{\beta /\alpha }dy.$$

(37)

Here, $g(x;\alpha ,0)$ and $g(y,\beta ,1)$ are the densities of symmetric and one-sided strictly stable laws [53,54] and $0<\alpha ⩽2$, $0<\beta ⩽1$.

We consider the case $\langle T\rangle <\infty$ and $\langle R^2\rangle =\infty$. Let the rest times have the exponential distribution (26) and free paths have the power distribution

$$p\left(x\right)=\left\{\begin{array}{cc}0,\hfill & x<x_0\hfill \\ \alpha x_0^{\alpha }{x}^{-\alpha -1},\hfill & x⩾x_0,\hfill \end{array}\right.,0<\alpha <2.$$

(38)

As mentioned earlier, distributions of this kind have moments not exceeding $\alpha$, i.e., $\langle X^k\rangle <\infty$, if $k<\alpha$ and $\langle X^k\rangle =\infty$, if $k⩾\alpha$. Since the parameter $\alpha$ takes values from the interval $0<\alpha ⩽2$, then at values $0<\alpha ⩽1$, the mathematical expectation of free paths $\langle R\rangle$ is infinite, and at values $1<\alpha ⩽2$, the mathematical expectation is finite. In this regard, it is necessary to consider these two cases separately.

At first, we consider the case $0<\alpha <1$. Let us perform the Fourier transform of the density (38)

$$\widehat{p}\left(k\right)=\alpha x_0^{\alpha }{\int }_{x_0}^{\infty }{x}^{-\alpha -1}{e}^{ikx}dx.$$

(39)

Integrating this expression once by parts, we obtain

$$\widehat{p}\left(k\right)=e^{ik{x}_0}+ik{x}_0^{\alpha }{\int }_{x_0}^{\infty }{x}^{-\alpha }{e}^{ikx}dx.$$

In this integral, we change the integration variable $ikx=-t$. Then, we will pass to the redistribution of $k\to 0$ and keep the summands that do not exceed $x^{\alpha }$. As a result, we obtain

$$\widehat{p}\left(k\right)\approx 1-{(-ik{x}_0)}^{\alpha }{\int }_0^{-i\infty }{t}^{-\alpha }{e}^{-t}dt.$$

(40)

Next, we will use the well-known result (see §1.5, the formula (31) in [55])

$${\int }_0^{\infty }{t}^{\gamma -1}{e}^{-ctcos\psi -ictsin\psi }dt=\mathrm{\Gamma }\left(\gamma \right)c^{-\gamma }{e}^{-i\gamma \psi },$$

where $-\frac{\pi }{2}<\psi <\frac{\pi }{2},\Re \gamma >0$, or $\psi =\pm \frac{\pi }{2},0<\Re \gamma <1$. If we consider that $z=t{e}^{i\psi }$, then this formula can be represented in the form

$${\int }_0^{\infty e^{i\psi }}{z}^{\gamma -1}{e}^{-cz}dz=\mathrm{\Gamma }\left(\gamma \right)c^{-\gamma },$$

(41)

where $-\frac{\pi }{2}<\psi <\frac{\pi }{2},\Re \gamma >0$, or $\psi =\pm \frac{\pi }{2},0<\Re \gamma <1$. Comparing this formula with the integral on the right-hand side (40), we obtain

$$\widehat{p}\left(k\right)\approx 1-{(-ik{x}_0)}^{\alpha }\mathrm{\Gamma }(1-\alpha ),0<\alpha <1.$$

(42)

Now, we consider the case $1<\alpha <2$. Integrating twice in parts the Fourier transform of the density (39), we get

$$\widehat{p}\left(k\right)=e^{ik{x}_0}+\frac{ik{x}_0}{\alpha -1}{e}^{ik{x}_0}+\frac{k^2{x}_0^{\alpha }}{1-\alpha }{\int }_{x_0}^{\infty }{x}^{-\alpha +1}{e}^{ikx}dx.$$

Now we will change the integration variable $ikx=-t$ in this expression, and then, we will pass to the limit $k\to 0$. When passing to the limit, we keep only summands with a degree not exceeding $\alpha$. As a result, we have

$$\widehat{p}\left(k\right)\approx 1+\frac{ik{x}_0}{\alpha -1}+\frac{{(-ik{x}_0)}^{\alpha }}{\alpha -1}{\int }_0^{-i\infty }{t}^{1-\alpha }{e}^{-t}dt.$$

If we use the formula (41) now, then we finally get

$$\widehat{p}\left(k\right)\approx 1+\frac{ik{x}_0}{\alpha -1}+\frac{{(-ik{x}_0)}^{\alpha }}{\alpha -1}\mathrm{\Gamma }(2-\alpha ),1<\alpha <2.$$

(43)

Thus, we have obtained asymptotic expressions for the Fourier transform of the free path distribution (38) for two cases: $0<\alpha <1$ and $1<\alpha <2$.

Now, we get back to the expression (25) and consider the multiplier ${\omega }_1\widehat{p}\left(k\right)+{\omega }_2\widehat{p}(-k)$. As in the previous cases, we will consider the symmetrical random walks (${\omega }_1={\omega }_2=1/2$). Using the expression (42) and taking account of the relation ${(-i)}^{\alpha }+i^{\alpha }=cos(\pi \alpha /2)$ for the case $0<\alpha <1$, we obtain

$${\omega }_1\widehat{p}\left(k\right)+{\omega }_2\widehat{p}(-k)=1-C{k}^{\alpha },C=x_0^{\alpha }\mathrm{\Gamma }(1-\alpha )cos(\pi \alpha /2).$$

Using the doubling formula for the Gamma function $\mathrm{\Gamma }\left(2z\right)={\left(2\pi \right)}^{-1/2}{2}^{2z-1/2}\mathrm{\Gamma }\left(z\right)$$\mathrm{\Gamma }(z+1/2)$ and the symmetry formula $\mathrm{\Gamma }(\frac{1}{2}+z)\mathrm{\Gamma }(\frac{1}{2}-z)=\frac{\pi }{cos\left(\pi z\right)}$, the coefficient C can be represented in the form

$$C=2^{-\alpha }\frac{\sqrt{\pi }{x}_0^{\alpha }\mathrm{\Gamma }\left(1-\frac{\alpha }{2}\right)}{\mathrm{\Gamma }\left(\frac{1}{2}+\frac{\alpha }{2}\right)}.$$

In the case $1<\alpha <2$, it is necessary to use the expression (43). As a result, we obtain

$${\omega }_1\widehat{p}\left(k\right)+{\omega }_2\widehat{p}(-k)=1-C^{\prime }{k}^{\alpha },C^{\prime }=x_0^{\alpha }\frac{\mathrm{\Gamma }(2-\alpha )}{\alpha -1}sin\left(\frac{\pi }{2}(\alpha -1)\right),$$

where we also assume symmetrical random walks $({\omega }_1={\omega }_2=1/2)$. Using the doubling formula for the Gamma function and the symmetry formula $\mathrm{\Gamma }\left(z\right)\mathrm{\Gamma }(1-z)=\frac{\pi }{sin\left(\pi z\right)}$, it is possible to show that $C^{\prime }=C$. Thus,

$$\frac{1}{2}(\widehat{p}\left(k\right)+\widehat{p}(-k))=1-C{k}^{\alpha },\alpha \in (0,1)\cup (1,2).$$

(44)

Now, we substitute this expression and the expression (28) in the solution (25). As a result, we obtain

$$\widehat{\tilde{\rho }}(k,\lambda )=\frac{1}{\lambda +D_C{k}^{\alpha }},D_C=\frac{C}{\langle T\rangle }.$$

(45)

We rewrite this expression in the form

$$\widehat{\tilde{\rho }}(k,\lambda )\lambda -1=-D_C{k}^{\alpha }\widehat{\tilde{\rho }}(k,\lambda ).$$

If we use now the relations (31) and (A5), then it is possible to perform the inverse Fourier–Laplace transform of this equation easily. As a result, we obtain

$$\frac{\partial \rho (x,t)}{\partial t}=D_C\frac{{\partial }^{\alpha }\rho (x,t)}{{\partial \left|x\right|}^{\alpha }}+\delta \left(x\right)\delta \left(t\right).$$

The solution to this equation can be obtained by performing the inverse Fourier–Laplace transform of the solution (45). It was done in the work [18], which showed that the solution to this equation has the form

$$\rho (x,t)={\left(D_{C}t\right)}^{-1/\alpha }q\left(x{\left(D_{C}t\right)}^{-1/\alpha };\alpha ,1\right),$$

where the density $q(x;\alpha ,\beta )$ is determined by the expression (37).

It remains to consider the case $\langle T\rangle =\infty$ and $\langle R^2\rangle =\infty$. Suppose the rest times have the distribution (35), and free paths are distributed according to the law (38). Since the parameter $\alpha$ takes the values in the range $0<\alpha <2$, then, as it was mentioned earlier, there is a necessity to consider two cases: $0<\alpha <1$ and $1<\alpha <2$. This is determined by the fact that when passing to the limit $k\to 0$, it necessary to take account of different summands in the expansion of the function image $p\left(x\right)$. However, such a passage to the limit has already been performed when considering the previous case. It was shown that in the case of symmetric random walks, the multiplier ${\omega }_1\widehat{p}\left(k\right)+{\omega }_2\widehat{p}(-k)$,which is a component of the expression (25) has the form (44). The Laplace transform of the density (35) was also obtained by us earlier. It has the form (36). Now, we substitute the expansions (36) and (44) in the solution (25) and will keep summands in the obtained expression that do not exceed $k^{\alpha }$ and ${\lambda }^{\beta }$. As a result, we get

$$\widehat{\tilde{\rho }}(k,\lambda )=\frac{{\lambda }^{\beta -1}}{{\lambda }^{\beta }+D_D{k}^{\alpha }},D_D=2^{-\alpha }\frac{x_0^{\alpha }\sqrt{\pi }\mathrm{\Gamma }(1-\frac{\alpha }{2})}{t_0^{\beta }\mathrm{\Gamma }(1-\beta )\mathrm{\Gamma }\left(\frac{1}{2}(1+\alpha )\right)}.$$

Here, $D_D$ is the diffusion coefficient. If we use the relations (A3) and (A5), it is possible to show easily that the obtained expression is the solution to the equation

$${}_0{D}_t^{\beta }\rho (x,t)=D_D\frac{{\partial }^{\alpha }\rho (x,t)}{{\partial \left|x\right|}^{\alpha }}+\frac{t^{-\beta }}{\mathrm{\Gamma }(1-\beta )}\delta \left(x\right).$$

The paper [18] shows that the solution of this equation has the form

$$\rho (x,t)={\left(D_D{t}^{\beta }\right)}^{-1/\alpha }q\left(x{\left(D_D{t}^{\beta }\right)}^{-1/\alpha };\alpha ,\beta \right),$$

where $q(x,\alpha ,\beta )$ is the density of a fractional-stable law (37).

As we can see, in all considered cases, the asymptotic distribution of particles at $t\to \infty$ and $x\to \infty$ is described with equations of the same type. The difference between these equations lies only in the order of the derivative with respect to time or coordinate and in the diffusion coefficient. If we assume that the case $\langle T\rangle <\infty$ corresponds to the parameter value $\beta =1$, and the case $\langle R^2\rangle <\infty$ corresponds to the parameter value $\alpha =2$, then the random walk process is described with the anomalous diffusion equation

$${}_0{D}_t^{\beta }\rho (x,t)=D\frac{{\partial }^{\alpha }\rho (x,t)}{{\partial \left|x\right|}^{\alpha }}+\frac{t^{-\beta }}{\mathrm{\Gamma }(1-\beta )}\delta \left(x\right),$$

(46)

where $0<\beta ⩽1$, $0<\alpha ⩽2$, and the generalized diffusion coefficient D for different parameter values takes different values

$$D=\left\{\begin{array}{ccc}{D}_A,& \beta =1,& \alpha =2,\\ D_B,& 0<\beta <1,& \alpha =2,\\ D_C,& \beta =1,& 0<\alpha <2,\\ D_D,& 0<\beta <1,& 0<\alpha <2.\end{array}\right.$$

(47)

As the paper [18] shows, the solution of this equation has the form

$$\rho (x,t)={\left(D{t}^{\beta }\right)}^{-1/\alpha }q\left(x{\left(D{t}^{\beta }\right)}^{-1/\alpha };\alpha ,\beta \right),$$

(48)

where $q(x;\alpha ,\beta )$ is the density of a fractional-stable law (37). According to the properties of fractional-stable laws [51] in case of parameter values $\beta =1,\alpha =2$, the density $q(x;2,1)$ becomes the density of the normal law, and in case $\beta =1$, the density $q(x;\alpha ,1)$ is the density of a stable law. Thus, Equation (46) describes the random walk process in all considered cases. The transition from one case to another is carried out only by replacing the value of the generalized diffusion coefficient D. Different types of the generalized diffusion coefficient are determined by different distributions of the rest times and the value of a free path. It is seen from (48) that this solution possesses the self-semilarity property. Therefore, the diffusion packet expands with time according to the law with exponent $\gamma =\beta /\alpha$, i.e., $\Delta \left(t\right)\propto D^{1/\alpha }{t}^{\beta /\alpha }$. As we can see in the case of normal diffusion, we obtain well-known result $\Delta \left(t\right)\propto \sqrt{D_A}{t}^{1/2}$.

4. The Solution to the Equation of Anomalous Diffusion

As we can see from Section 2, the CTRW process is described by the integral transport equation, which in the one-dimensional case takes the form (24). Therefore, the Monte Carlo method can be used to find a solution to this equation. The advantage of Monte Carlo methods is that they allow one to find a solution in multidimensional problems, as well as for various boundary and initial conditions. In this paper, we consider the solution to the anomalous diffusion Equation (46) under the condition $\rho (x,t)=0$ if $t<0$ and $\rho (x,t)\to 0$ if $x\to \pm \infty$.

From Section 2, the simplest method of stochastic solution to the anomalous diffusion equation immediately follows (46), based on trajectory modeling and histogram density estimation. Each trajectory begins at the moment of time $t=0$ from the origin of coordinates $x=0$ from the state of rest. In the state of rest, the particle stays for a random time $T_1$. Then, with equal probability, the particle jumps to the right or to the left at a distance $R_1$. After that, the particles will enter a state of rest. Then, the process continues in the same way. The construction of the trajectory continues as long as the condition is met

$$\sum _{k=1}^{N\left(T^{*}\right)}{T}_k⩽T^{*},$$

(49)

where $T^{*}$ is a given moment in time at which a solution is to be found. As soon as this condition is no longer met, the trajectory is terminated. Then, a new trajectory begins.

Depending on the parameter values $\alpha$ and $\beta$, the rest times $T_k$ and free path values $R_k$ have different distributions. As noted earlier, the determining parameters that influence the form of the differential equation are the mathematical expectation of the rest time $\langle R^2\rangle$. If the mathematical expectation of the rest time has a finite value, then this corresponds to the case $\beta =1$. If the variance of a free path has a finite value, then this corresponds to the case $\alpha =2$. In this paper, exponential distributions are chosen as such distributions. Therefore, in the case $\beta =1$, the rest times have the distribution (26). In the case $\alpha =2$, free paths have the distribution (27). The values $0<\beta <1$ correspond to the case of the infinite mathematical expectation of the rest time. In this case, the rest times are distributed with the density (35). The values $0<\alpha <2$ correspond to the case of the infinite variance of free paths. In this case, free paths are distributed with the density (38). Thus, with the value $\beta =1$, random times $T_k$, $k=1,2,\dots$ are modeled according to the formula $T=-{\nu }^{-1}log\left(\zeta \right)$, and in the case $0<\beta <1$ according to the formula $T=t_0{\zeta }^{-1/\beta }$. If $\alpha =2$ random free paths $R_k$, $k=1,2,\dots$ are modeled according to the formula $R=-{\mu }^{-1}log\left(\zeta \right)$. If $0<\alpha <2$, then $R=x_0{\zeta }^{-1/\alpha }$. Here, $\zeta$ represents equally distributed random values on the segment $(0,1]$.

To construct the simplest histogram estimate of the solution (46), all the region of interest $\Delta =[a,b]$ is broken down into disjoint intervals ${\Delta }_i=(x_i,x_{i+1}]$, $i=1,2,\dots M-1$, $x_1=a,x_M=b$. To construct a histogram, the trajectory is modeled until the condition is met (49). As soon as this condition ceases to be met, the trajectory is terminated, and the contribution from this trajectory is calculated

$$h_j\left({\Delta }_i\right)=\frac{I\left({\Delta }_i\right)}{{\Delta }_i}.$$

where $I\left({\Delta }_i\right)$ is the interval indicator ${\Delta }_i$

$$I\left({\Delta }_i\right)=\left\{\begin{array}{cc}1,& x^{*}\in {\Delta }_i,\\ 0,& x^{*}\notin {\Delta }_i,\end{array}\right.$$

where $x^{*}$ is the coordinate of a particle at a moment of time $T^{*}$. As a result, the density estimate for the interval ${\Delta }_i$ is given with the expression

$$\tilde{\rho }({\Delta }_i,t)=\frac{1}{N}\sum _{j=1}^N{h}_j\left({\Delta }_i\right),$$

(50)

where the summation is performed over an ensemble of N independent trajectories.

Despite the simplicity of this estimate, it has several disadvantages. Firstly, the estimate of the solution $\tilde{\rho }({\Delta }_i,t)$ is sought for the interval ${\Delta }_i$. This is the source of the systematic (horizontal) component of the error ${\delta }_x$. Secondly, this estimate also contains the statistical component of the error $\widehat{\delta }$, which decreases as $N^{-1/2}$ at $N\to \infty$. It is impossible to eliminate these errors completely; one can only reduce their value. However, a decrease in one of these values leads to an increase in the other value or to an increase in the calculation time.

It is possible to get rid of the systematic component of the error completely ${\delta }_x$ if to consider one of the varieties of a local estimate. As in the case of the histogram estimate, the problem is to determine the probability density of detecting a particle at the point $x^{*}$ at a moment of time $T^{*}$. The main element of solving the problem of transport theory by the Monte Carlo method, trajectory modeling, remains unchanged. The difference lies in the estimation method. In the case of a local estimate, the probability of a particle hitting a point $(x^{*},T^{*})$ is calculated assuming that the next collision is the final one. This probability is calculated after each collision (state change) of the particle. As a result, for the CTRW process, this probability is given by the value

$$\psi (\xi ,\tau )=\frac{1}{2}p\left(\xi \right)Q\left(\tau \right),$$

(51)

where $Q\left(\tau \right)={\int }_{\tau }^{\infty }q\left(t\right)dt$. It should be pointed out that for the CTRW process, this probability is calculated only for the transition <<rest>>→<<jump>>. For transitions <<jump>>→<<rest>>, this probability will be equal to zero. As a result, the contribution to the density estimate from each individual trajectory has the form

$$h_j(x^{*},T^{*})=\sum _{i=1}^{K\left(T^{*}\right)}\psi \left(\right|x^{*}-x_i|,T^{*}-t_i),$$

where $K\left(T^{*}\right)$ is the number of state change acts <<rest>>→<<jump>> that occurred for the interval of time $(0,T^{*})$. The density estimate takes the form

$$\tilde{\rho }(x^{*},T^{*})=\frac{1}{N}\sum _{j=1}^N{h}_j(x^{*},T^{*}),$$

(52)

where the summation is performed over an ensemble of N independent trajectories.

As we can see, the local estimate evaluates the density at a given point $x^{*}$. This means that this estimate does not contain a systematic component of the error ${\delta }_x$, which is connected with the finite value of the interval ${\Delta }_i$, as it was in the histogram estimate. Moreover, since each individual trajectory contributes more than once, as was the case with the histogram estimate and $K\left(T^{*}\right)$ times, then this leads to a decrease in the statistical error.

The results of solving Equation (46) are given in Figure 1, Figure 2 and Figure 3. In these figures, the points correspond to the results of the local estimation (52), the circles correspond to the results of the histogram estimation (50), and the solid curves correspond to the solution (48) with the corresponding diffusion coefficient determined from the relation (47). The solution results are given for different points in time. From Figure 1, it is clear that for the parameter value $\beta =0.3$, the results of the local and histogram estimation coincide with the solution (48) at time $T^{*}=10$. This means that by this time, the walk process has already entered the asymptotic regime. Thus, for the given parameter values, the estimate (52) can be used to solve Equation (46) at times $T^{*}⩾10$. In the case $\beta =0.6$ with the time values $T^{*}=10$ and 100, it is clear (see Figure 1 on the right) that the results of the local estimate and solution (48) differ. This means that at such times, the random walk process has not yet reached the asymptotic regime. However, at times $T^{*}={10}^3$ and ${10}^4$, the random walk process already reaches the asymptotic regime. This means that at times $T^{*}⩾{10}^3$, the estimate (52) can be used to solve numerically Equation (46) for the given parameter values. Similar conclusions can be drawn for other presented solution results. For an exponential distribution of rest times and $\alpha =0.7$ and $\alpha =1.4$ (Figure 2), it is clear that the random walk process reaches the asymptotic regime at time $T^{*}=10$. Thus, for the indicated values of the parameters, the estimate (52) can be used to solve Equation (46) at times $T^{*}⩾10$. In the case $\beta =0.7$ and $\alpha =0.9$ (see Figure 3 on the left), the random walk process becomes asymptotic at time $T^{*}=1000$, and for the values $\beta =0.7$, $\alpha =1.4$ (see Figure 3 on the right) at time $T^{*}=100$.

5. Conclusions

The use of the theory of anomalous diffusion and equations in fractional derivatives to describe combustion processes is only at the initial stage of research. At the moment, there are not many papers in which this approach is used to describe combustion processes. However, existing experimental studies (see, for example [37]) indicate the legitimacy of this approach. One of the main difficulties in using equations in fractional derivatives is finding their solutions. Analytical methods for solving equations of this kind are only at the stage of development. Therefore, one of the main methods of solving equations in fractional derivatives is finite difference methods. In the papers [41,42,43,44,45], the numerical solution of the equation of anomalous diffusion is investigated, which is expressed in terms of derivatives of fractional order only in the case when one of the derivatives (or derivative with respect to time or derivative with respect to coordinate) is not of integer order. In this paper, we propose a numerical method for solving the anomalous diffusion equation in which both the time derivative and the coordinate derivative can be of non-integer order.

This paper considers the numerical method for solving the anomalous diffusion equation based on the use of a local estimate. This method is based on the idea of modeling random realizations of particle trajectories. However, unlike the histogram method for estimating the distribution density, each individual trajectory contributes to the estimate not once, but several times. In the proposed approach, after each act of a change in the state of a particle, the probability is calculated that to get to a given point with the coordinate $x^{*}$ in a moment of time $t^{*}$, the following collision will turn out to be the final. For the considered model of walks, this probability has the form (51). This local estimate has several advantages over the histogram estimate. Since the density is estimated at a given point $x^{*}$, then this means that the result of the estimate does not contain a systematic component of the error ${\delta }_x$. Specifying a set of points $x_i^{*}$, $i=1,2,\dots ,M$, it becomes possible to estimate the solution at several points at once. Moreover, one trajectory will contribute to all points at once $x_i^{*}$. As a result, the probability $\psi (\xi ,\tau )$, determined by (51), is calculated after each act of changing the state of a particle; then, this leads to a considerable decrease in statistical error. Taking account of the fact that from one trajectory, the contribution can be calculated at once to all points $x_i^{*}$, $i=1,2,\dots ,M$ of a given set, this means that the desired solution can be constructed as a smooth function of the coordinate x.

In conclusion, one more important point should be noted. As shown in this study, the equation of anomalous diffusion (46) describes the asymptotic distribution (at $t\to \infty$) of particles in the CTRW process. A characteristic feature of this process is that the random rest times of a particle are characterized by a distribution with an infinite mathematical expectation $\langle T\rangle =\infty$, and the random free paths in the case of the exponent $\alpha \in (1,2)$ are characterized by an infinite second moment $\langle R^2\rangle =\infty$, and in the case of the exponent value $\alpha \in (0,1]$, then by the infinite mathematical expectation $\langle R\rangle =\infty$. Taking into account that in the process of CTRW, a particle instantly moves from one point in space to another (instantaneous jumps), this leads to a non-physical result: in an arbitrarily small time interval, a particle can be at an arbitrarily large distance from the point of its previous position. It should be noted that the random walk process is also characterized by instantaneous jumps that leads to the diffusion Equation (33). However, as it was shown at the beginning of Section 3, to obtain normal diffusion, it is necessary to assume that the distributions of the rest time and the distribution of the jump value of a particle have a finite mathematical expectation and a finite variance. As a result, the instantaneous motion of a particle from one point in space to another is compensated by the small value of these jumps. In the case of a power-law distribution of the jump value (38), the probability of large jumps remains significant for any jump value. As the value decreases of $\alpha$, this probability increases. This property is characteristic of power distributions. Therefore, when using the equation of anomalous diffusion (46) to describe combustion processes, especially in furnaces where the geometry is given, a certain amount of care must be taken. It should be understood that the solution to the anomalous diffusion equation decreases according to a power law at $x\to \infty$ and is different from zero in the entire space. The latter means that the probability of detecting a particle at an arbitrarily large distance from the source at an arbitrarily close moment of time to the initial time is different from zero. Taking account of the fact that this probability decreases according to a power law, then this probability will be significant, and it can no longer be neglected.