Numerical Solution to Anomalous Diffusion Equations for Levy Walks

Abstract

The process of Levy random walks is considered in view of the constant velocity of a particle. A kinetic equation is obtained that describes the process of walks, and fractional differential equations are obtained that describe the asymptotic behavior of the process. It is shown that, in the case of finite and infinite mathematical expectation of paths, these equations have a completely different form. To solve the obtained equations, the method of local estimation of the Monte Carlo method is described. The solution algorithm is described and the advantages and disadvantages of the considered method are indicated.

1. Introduction

At present, the theory of anomalous diffusion is rarely used to describe the combustion processes of a substance, although there are all the prerequisites for this. For example, when a substance is burning, various kinds of turbulence are developing. As a result, the diffusion packet width $\Delta \left(t\right)$ stops obeying the law $\Delta \left(t\right)\propto t^{\gamma }$ with an exponent $\gamma =1/2$ and starts growing with time by the law with an exponent $\gamma \ne 1/2$, which testifies to the appearance of anomalous diffusion. Signs of the appearance of anomalous diffusion at thermal transport in a low-dimensional system are indicated in the papers [1,2]. This assumption is also supported by experimental data. For example, in the work [3], the thermal radiation in the combustion chamber during the combustion of natural gas and acetylene was studied, depending on the level of fuel enrichment with oxygen. Time series analysis showed that the combustion process at any oxygen concentration is subdiffusive in nature.

The assumption about the formation of anomalous diffusion during the combustion of a substance allows us to introduce into consideration the fractional differential equations of anomalous diffusion. An effective coefficient of heat conductivity for the Levy–Fokker–Plank equation was obtained in the papers [4,5]. In the papers [6,7,8], to describe the combustion process, it is proposed to use the fractional differential equation of anomalous diffusion:

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

(1)

where ${}_0{D}_t^{\beta }$ is the fractional Riemann–Liouville derivative [9] of the order $0<\beta <1$ by time and ${\partial }_{xx}$ is the classical particular derivative of the second order by coordinate. In the paper [10], a two-dimensional combustion model with a fractional time derivative was studied. To solve the obtained equation of diffusion, the authors develop an adaptive finite-difference discontinuous Galerkin method. A modification of the Equation (1) in the case of the dependence of the diffusivity on temperature and relaxation time has been considered in the article [11]. In the papers [8,12], a fractional-differential combustion model with the first derivative with respect to time and a fractional derivative with respect to the spatial variable is considered:

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

(2)

Here, ${\partial }_t$ is the partial time derivative, and $\frac{{\partial }^{\alpha }}{{\partial \left|x\right|}^{\alpha }}$ is the fractional-differential Riesz operator [9]. To solve this equation, a finite-difference scheme using an adaptive strategy is described in the papers [8,12]. In the paper [13], a fractional-differential generalization of the kinetic equation was obtained that describes the relationship between the radius of the ball and time in the model of the combustion of a fireball, theoretically predicted by the Soviet physicist Ya.B. Zeldovich [14].

However, when using the Equations (1) and (2), it should be kept in mind that they are based on the process of Continuous Time Random Walk (CTRW) [15,16,17,18,19,20,21,22]. This process assumes that a particle instantly moves from one point in space to another at a random distance $R_i,i=1,2,3,\dots$, and then rests at this point for a random time $T_i,i=1,2,3,\dots$ All these random values $R_i$ and $T_i$ are independent of each other and between one another and they are distributed by the laws $p_R\left(x\right)\propto x^{-\alpha -1},x\to \infty$ and $q_T\left(t\right)\propto t^{-\beta -1},t\to \infty$, respectively.

Depending on the value $\alpha$, the path distribution has different properties. If $\alpha =2$, then the distribution has finite mathematical expectation and variance, at $1<\alpha <2$, mathematical expectation is finite and variance is infinite, at $0<\alpha ⩽1$ both mathematical expectation and variance are infinite. The instantaneousness of jumps means that, for an arbitrarily small time interval from the initial one, the particle can be at an arbitrarily large distance from the source. In some cases, this non-physical behavior of a particle does not lead to any contradiction with the experiment. For example, in the case of the normal diffusion ($\alpha =2,\beta =1$), when the distribution of paths and rest times have a finite variance, the process is described by the classical diffusion equation, the solution to which is expressed in terms of the normal distribution. It is well known that the normal distribution is nonzero on the entire number axis, which indicates the infinite velocity of the particle in the inherent walk model. However, in view of the finite variance of the distribution of paths, the infinite velocity of motion is compensated for by a small value of the paths. The situation is completely different in the case of anomalous diffusion. As the exponent $\alpha$ decreases, the probability of greater paths increases and at values $\alpha <1$ this probability turns out to be significant. Therefore, it is necessary to use the anomalous diffusion model with a certain degree of caution, especially when considering problems with limited spatial geometry or processes limited in time, for example, to describe combustion processes in furnaces.

One of the ways to eliminate the difficulty described above is to introduce a constant final velocity of the particle. One of the first works in which a constant velocity of particle motion was introduced into the model of anomalous diffusion is the work [23]. In this paper, the authors called this model Levy walks. Further study of this model was carried out in the works [24,25,26]. The work [27] is devoted to the study of Levy walks in bounded and semi-bounded spaces. In the work [28], kinetic equations of anomalous diffusion with a finite velocity are obtained, the root-mean-square deviation is investigated, and an exact solution to the kinetic equation in the one-dimensional case is obtained. The work [29] is devoted to the study of statistical moments for the case of Levy random walks with a random finite velocity without traps, and the case of multidimensional walks with traps of arbitrary type with constant velocity is considered in the works and the case of multidimensional walks with traps of arbitrary type with constant velocity is considered in the works [30,31,32,33]. The work [34] examines the influence of the final velocity on the spatial distribution of particles in Levy walks with exponential traps. In this work, it was found that, in the case $1<\alpha ⩽2$, taking account of the constant velocity of particle motion is reduced to a decrease in the diffusion coefficient in the equation of anomalous diffusion (2). In the case $\alpha <1$ distributions take W-like or U-like form and cannot be described by the Equation (2). The papers [35,36,37] study the process of Levy random walks without traps and it was shown that in the case of infinite mathematical expectation of the distribution of paths, the asymptotic distributions have U-like and W-like form. In the paper [38] an expression for current within framework of the Levy walks model was obtained. The obtained expression generalize the Fourier’s law to the case of anomalous thermal transport for the Levy diffusion model.

The work [39] succeeded in addressing the problem of describing Levy random walks in the case $\alpha <1$. In this work, it was shown that, in the case of an infinite mathematical expectation of the distribution of paths, to take account of the finite velocity, it is necessary to replace the fractional Laplacian in the equation of anomalous diffusion by a material derivative of a fractional order. Later, the authors of the works [40,41,42,43,44] come to the same conclusion. The solution to the equation of anomalous diffusion with a material derivative of the fractional order was obtained in the works [44,45,46]. In these works, it was shown that the solution to this equation is expressed in terms of the Lamperti distribution.

In this paper, we consider a method for the numerical solution to the equations of anomalous diffusion taking account of the constant finite velocity of the particle motion between collisions. The work is organized as follows. In Section 2, a kinetic equation describing the considered process of walks is derived. In Section 3, fractional differential equations and the solutions of these equations are obtained, describing the asymptotic (at $t\to \infty$ and $x\to \infty$) distribution of particles. It was shown that, in the case of an infinite mathematical expectation of the distribution of paths ($0<\alpha <1$) and finite mathematical expectation $1<\alpha <2$, the behavior of the process is described by completely different equations. In Section 4, a numerical method for solving the kinetic equation based on the method of local density estimation is considered.

2. Kinetic Equation of the Random Walk Process

To obtain the kinetic equation, we will use the approach proposed in the paper [28], which was developed in the paper [47]. We will consider the density of collisions $f(\mathbf{r},\mathbf{p},t)$, where $\mathbf{r}$ is the radius-vector of the particle, $\mathbf{p}$—is the particle momentum, 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 $d\mathbf{r}$ of the vicinity of the point $\mathbf{r}$ for the interval of time $dt$, at which the momentum of the particle takes on a value from $\mathbf{p}$ to $\mathbf{p}+d\mathbf{p}$. We will consider the nonrelativistic case $\mathbf{p}=m\mathbf{v}$. Without loss of generality, we assume that $m=1$. The paper [47] shows that, with the presence n of discrete states, the value $f(\mathbf{r},\mathbf{v},t)$ can be represented in the form:

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

(3)

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}({\mathbf{\Omega }}^{\prime },\mathbf{\Omega })d{\mathbf{\Omega }}^{\prime }\times \hfill \\ \hfill \int f_i\left(\mathbf{r}-v^{\prime }{\mathbf{\Omega }}^{\prime }\tau ,v^{\prime }{\mathbf{\Omega }}^{\prime },t-\tau \right)h_{ij}(v^{\prime },v)d{v}^{\prime }.\end{array}$$

(4)

Here, $k_i\left(\tau \right)$ is the probability density distribution of the residence time in the state i, $c_{ij}$—the probabilities of transition from the state i into the state j, $W_{ij}({\mathbf{\Omega }}^{\prime },\mathbf{\Omega })$ is the probability density that before collision the velocity had the direction ${\mathbf{\Omega }}^{\prime }$, after collision the direction took the value $\mathbf{\Omega }$, $h_{ij}(v^{\prime },v)$ is the density of the probability of the change in velocity from the value $v^{\prime }$ to v, $s_j(\mathbf{r},\mathbf{v},t)$ is the density of new particle sources in the state j, $\mathbf{v}=v\mathbf{\Omega }$, $v=\left|\mathbf{v}\right|$, $d\mathbf{v}=dvd\mathbf{\Omega }$, the summation is carried out over all possible previous states. The values $c_{ij}$, $W_{ij}({\mathbf{\Omega }}^{\prime },\mathbf{\Omega })$ and $h_{ij}(v^{\prime },v)$ are normalized:

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

(5)

The transition from the density of collisions $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\mathbf{\Omega }\tau ,v\mathbf{\Omega },t-\tau \right)d\tau ,$$

(6)

where $K\left(t\right)={\int }_t^{\infty }k\left(\tau \right)d\tau$. Substituting (3) in (6) we get 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),$$

(7)

where

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

(8)

where $K_j\left(t\right)={\int }_t^{\infty }{k}_j\left(\tau \right)d\tau$. The physical interpretation of the last expression is simple. To detect the 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 the moment of time from t to $t+dt$ the particle must pass to this state at the point $\mathbf{r}-v\mathbf{\Omega }\tau$ at the moment of time $t-\tau$ and stay in this state during the time greater than $\tau$. Transition to the density of particles $\rho (\mathbf{r},t)$ is carried out with the help of the integral,

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

(9)

The system of Equations (4), (7) and (8), together with conditions (5) describes practically any process of random walks with n discrete states under fairly general assumptions about the scattering indicatrix $W_{ij}({\mathbf{\Omega }}^{\prime },\mathbf{\Omega })$ and the law of redistribution of speed $h_{ij}(v^{\prime },v)$. In this work, using these equations, we obtain kinetic equations describing Levy walks with a constant velocity of motion between two successive scatterings of a particle.

We define the process of walks as follows. There is only this state—the state of motion ($n=1$). A particle moves at a constant velocity v between two successive collisions. After the collision, the particle changes its direction, which is determined by the scattering indicatrix $W\left(\mathbf{\Omega }\right)$. After which, the particle continues to move in a new direction with the same constant velocity v. Random times between two successive collisions of a particle $T_i,i=1,2,3,\dots$ are independent random values. Since the motion occurs with a finite speed, then for times $T_i$ the particle covers the path $R_i=v{T}_i,i=1,2,3,\dots$. The values $R_i$ are the paths of particles.

Since there is only one state, then in Equations (4), (7) and (8), we need to put $n=1$, $c_{11}=1$. We will also assume that the source can be represented in the form $s(\mathbf{r},\mathbf{v},t)=s(\mathbf{r},t)h_s\left(v\right)W_s\left(\mathbf{\Omega }\right)$ and, for brevity, we omit the subscript indicating the status number. As a result, we obtain:

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

(10)

$$\begin{array}{cc}\hfill f(\mathbf{r},\mathbf{v},t)& =s(\mathbf{r},t)h_s\left(v\right)W_s\left(\mathbf{\Omega }\right)\hfill \\ & +{\int }_0^{t}k\left(\tau \right)d\tau \int W({\mathbf{\Omega }}^{\prime },\mathbf{\Omega })d{\mathbf{\Omega }}^{\prime }\int f\left(\mathbf{r}-v^{\prime }{\mathbf{\Omega }}^{\prime }\tau ,v^{\prime }{\mathbf{\Omega }}^{\prime },t-\tau \right)h(v^{\prime },v)d{v}^{\prime }.\hfill \end{array}$$

(11)

The motion velocity between two successive collisions will be denoted by $v_0$. Let us also assume that the direction of motion after each collision does not depend on the previous direction of motion. In view of the foregoing, we get:

$$h(v^{\prime },v)=h_s\left(v\right)=\delta (v-v_0),W({\mathbf{\Omega }}^{\prime },\mathbf{\Omega })=W_s\left(\mathbf{\Omega }\right)=W\left(\mathbf{\Omega }\right).$$

Now, substituting these expressions in Equation (11), we obtain:

$$f(\mathbf{r},\mathbf{v},t)=W\left(\mathbf{\Omega }\right)\delta (v-v_0)\left(s(\mathbf{r},t)+{\int }_0^{t}k\left(\tau \right)d\tau \int d{\mathbf{\Omega }}^{\prime }\int f\left(\mathbf{r}-v^{\prime }{\mathbf{\Omega }}^{\prime }\tau ,v^{\prime }{\mathbf{\Omega }}^{\prime },t-\tau \right)d{v}^{\prime }\right).$$

(12)

From this relation, it is clear that, in cases when the densities of transition probabilities $h(v^{\prime },v)$ and $W_{ij}({\mathbf{\Omega }}^{\prime },\mathbf{\Omega })$ do not depend on the previous value $v^{\prime }$ and ${\mathbf{\Omega }}^{\prime }$, then the density of the collision can be represented in the form of the product $f(\mathbf{r},v\mathbf{\Omega },t)=W\left(\mathbf{\Omega }\right)h\left(v\right)F(\mathbf{r},\mathbf{t})$. Thus, we obtain:

$$f(\mathbf{r},v\mathbf{\Omega },t)=W\left(\mathbf{\Omega }\right)\delta (v-v_0)F(\mathbf{r},\mathbf{t}).$$

(13)

Now substituting this relation in (12) and by integrating over $dvd\mathbf{\Omega }$ we get the equation for $F(\mathbf{r},t)$:

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

(14)

Here, we used the normalization condition $\int W\left({\mathbf{\Omega }}^{\prime }\right)d{\mathbf{\Omega }}^{\prime }=1$. The physical meaning of the quantity $F(\mathbf{r},t)$ is quite simple. This is the density of collisions in the volume element $d\mathbf{r}$ of the vicinity of the point $\mathbf{r}$.

Let us now pass from the collision density to the phase density $\psi (\mathbf{r},\mathbf{v},t)$, and then to the density of particles $\rho (\mathbf{r},t)$. To this end, we put the expression (13) in (10). As result, we obtain:

$$\psi (\mathbf{r},\mathbf{v},t)={\int }_0^{t}K\left(\tau \right)W\left(\mathbf{\Omega }\right)\delta (v-v_0)F(\mathbf{r}-v\mathbf{\Omega }\tau ,t-\tau )d\tau .$$

Now integrating this expression over $dvd\mathbf{\Omega }$ and taking account of the ratio (9), we get the equation for the density of particles:

$$\rho (\mathbf{r},t)={\int }_0^{t}K\left(\tau \right)d\tau \int W\left(\mathbf{\Omega }\right)F(\mathbf{r}-v_0\mathbf{\Omega }\tau ,t-\tau )d\mathbf{\Omega }.$$

(15)

For the further solution to the obtained equations, it turns out to be convenient to pass from the time $\tau$ to the particle path $\xi$. By substituting the integration variable $\tau =\xi /v_0$, in Equations (14) and (15), we get:

$$\begin{array}{cc}\hfill \rho (\mathbf{r},t)& =\frac{1}{v_0}{\int }_0^{v_{0}t}P\left(\xi \right)d\xi \int F(\mathbf{r}-\xi \mathbf{\Omega },t-\xi /v_0)W\left(\mathbf{\Omega }\right)d\mathbf{\Omega },\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill F(\mathbf{r},t)& =s(\mathbf{r},t)+{\int }_0^{v_{0}t}p\left(\xi \right)d\xi \int F(\mathbf{r}-\xi {\mathbf{\Omega }}^{\prime },t-\xi /v_0)W\left({\mathbf{\Omega }}^{\prime }\right)d{\mathbf{\Omega }}^{\prime }..\hfill \end{array}$$

(17)

Here, the following notation was introduced: $p\left(\xi \right)=\frac{1}{v_0}k(\xi /v_0)$ is the density of probability of the path distribution and $K\left(\frac{\xi }{v_0}\right)\equiv P\left(\xi \right)={\int }_{\xi }^{\infty }p\left(y\right)dy$.

From this system of equations, it is possible to exclude the equation for $F(\mathbf{r},t)$. To do this, we will substitute Equation (17) in (16) and change the order of integration in the second summand. As a result, we obtain the equation for density,

$$\begin{array}{c}\rho (\mathbf{r},t)=\frac{1}{v_0}{\int }_0^{v_{0}t}P\left(\xi \right)d\xi \int s(\mathbf{r}-\xi \mathbf{\Omega },t-\xi /v_0)W\left(\mathbf{\Omega }\right)d\mathbf{\Omega }\hfill \\ \hfill +{\int }_0^{v_{0}t}p\left(\xi \right)d\xi \int \rho (\mathbf{r}-\xi {\mathbf{\Omega }}^{\prime },t-\xi /v_0)W\left({\mathbf{\Omega }}^{\prime }\right)d{\mathbf{\Omega }}^{\prime }.\end{array}$$

(18)

The first summand in this equation describes unscattered radiation. The second summand describes multiply scattered radiation. This equation describes the random walk of a particle with constant velocity in three-dimensional space with an arbitrary distribution of paths.

Let us simplify the problem and consider one-dimensional particle walks. Let the random walk process occur along the axis x. In this case, the function $W\left(\mathbf{\Omega }\right)$ takes the form:

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

(19)

where ${\omega }_1$ and ${\omega }_2$ are the probabilities of motion in the positive and negative directions of the axis $Ox$ respectively and ${\omega }_1+{\omega }_2=1$. Now substituting (19) in Equation (18) and considering that $\mathbf{\Omega }=(sin\theta cos\phi ,sin\theta sin\phi ,cos\theta )$, $\rho (\mathbf{r},t)=\rho (x,y,z,t)$, $s(\mathbf{r},t)=s(x,y,z,t)$, $d\mathbf{\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)=\frac{1}{v_0}{\int }_0^{v_{0}t}P\left(\xi \right)\left({\omega }_{1}s(x-\xi ,y,z,t-\xi /v_0)+{\omega }_{2}s(x+\xi ,y,z,t-\xi /v_0)\right)d\xi \hfill \\ \hfill +{\int }_0^{v_{0}t}p\left(\xi \right)\left({\omega }_1\rho (x-\xi ,y,z,t-\xi /v_0)+{\omega }_2\rho (x+\xi ,y,z,t-\xi /v_0)\right)d\xi .\end{array}$$

Since random walks along the axis $Ox$ are considered, then,

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

Now substituting these expressions into the previous equation and integrating over the variables y and z, we finally obtain:

$$\begin{array}{c}\rho (x,t)=\frac{1}{v_0}{\int }_0^{v_{0}t}P\left(\xi \right)\left({\omega }_{1}s(x-\xi ,t-\xi /v_0)+{\omega }_{2}s(x+\xi ,t-\xi /v_0)\right)d\xi \hfill \\ \hfill +{\int }_0^{v_{0}t}p\left(\xi \right)\left({\omega }_1\rho (x-\xi ,t-\xi /v_0)+{\omega }_2\rho (x+\xi ,t-\xi /v_0)\right)d\xi .\end{array}$$

(20)

The first component in this equation describes unscattered particles that, after escaping from the source, did not have a single collision and move in positive and negative directions, respectively. The second component describes multiply scattered particles, which at the moment of time $t-\xi /v_0$ had a collision and after that they began their motion in positive and negative directions, respectively.

3. Asymptotic Solution to a Kinetic Equation

An asymptotic solution to this equation can be found. To do this, we perform the Fourier–Laplace transform:

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

(21)

of the Equation (20). As a result, we obtain:

$$\begin{array}{c}\widehat{\rho }(k,\lambda )=\frac{1}{v_0}{\widehat{s}}_2(k,\lambda )\left({\omega }_1\widehat{P}(\lambda /v_0-ik)+{\omega }_2\widehat{P}(\lambda /v_0+ik)\right)\hfill \\ \hfill +\widehat{\rho }(k,\lambda )\left({\omega }_1\widehat{p}(\lambda /v_0-ik)+{\omega }_2\widehat{p}(\lambda /v_0+ik)\right).\end{array}$$

As a result, we obtained a simple algebraic equation, the solution to this equation has the form:

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

Let there be a point instantaneous source $s(x,t)=\delta \left(x\right)\delta \left(t\right)$. This means that $\widehat{s}(k,\lambda )=1$. Considering that

$$\widehat{P}\left(\lambda \right)={\int }_0^{\infty }{e}^{-\lambda \xi }P\left(\xi \right)d\xi =\frac{1-\widehat{p}\left(\lambda \right)}{\lambda },$$

we obtain

$$\widehat{\rho }(k,\lambda )=\frac{1/v_0}{1-W(k,v_0,\lambda )}\left({\omega }_1\frac{1-\widehat{p}\left(\frac{\lambda }{v_0}-ik\right)}{\frac{\lambda }{v_0}-ik}+{\omega }_2\frac{1-\widehat{p}\left(\frac{\lambda }{v_0}+ik\right)}{\frac{\lambda }{v_0}+ik}\right),$$

(22)

where $W(k,v,\lambda )={\omega }_1\widehat{p}(\lambda /v-ik)-{\omega }_2\widehat{p}(\lambda /v+ik)$. This solution describes the spatial distribution of particles with random walks of a particle at a constant velocity, with an arbitrary distribution of paths. This solution is not new and was obtained earlier (see, for example, [40,42,48,49,50,51,52]). In this paper, we will consider the asymptotic solution to this equation in the case when the distribution of paths has asymptotics of the form:

$$p\left(x\right)\propto \alpha x_0^{\alpha }{x}^{-\alpha -1},0<\alpha <2,x\to \infty ,$$

(23)

where $x_0={(\Gamma (1-\alpha )sin\left(\frac{\pi }{2}(1-\alpha )\right))}^{-1/\alpha }$. The cases of other distributions of paths are considered in the work [52].

Let us consider the case $0<\alpha <1$. In this case, the Laplace transform of the density (23) has the form (see [33]):

$$\widehat{p}\left(\lambda \right)\approx 1-{\left(\lambda x_0\right)}^{\alpha }\Gamma (1-\alpha ),\lambda \to 0.$$

Now substituting this expression into the solution (22), we obtain:

$$\widehat{\rho }(k,\lambda )=\frac{{\omega }_1{\left(\lambda -ik{v}_0\right)}^{\alpha -1}+{\omega }_2{\left(\lambda +ik{v}_0\right)}^{\alpha -1}}{{\omega }_1{\left(\lambda -ik{v}_0\right)}^{\alpha }+{\omega }_2{\left(\lambda +ik{v}_0\right)}^{\alpha }}.$$

(24)

This expression completely coincides with the result obtained in [37,49,53]. Taking account of the fact that the multiplier ${(\lambda \pm ikv)}^{\alpha }$ is the Fourier–Laplace transform of the fractional material derivative [40,54],

$${(\lambda \pm ıkv)}^{\alpha }\widehat{f}(k,\lambda )={\int }_0^{\infty }dt{\int }_{-\infty }^{\infty }{e}^{ıkx-\lambda t}f(x,t){\left(\frac{\partial }{\partial t}\pm v\frac{\partial }{\partial x}\right)}^{\alpha }dx,$$

we obtain that the process is described by the fractional differential equation [44,53,54,55]:

$$\begin{array}{c}\left[{\omega }_1{\left(\frac{\partial }{\partial t}-v_0\frac{\partial }{\partial x}\right)}^{\alpha }+{\omega }_2{\left(\frac{\partial }{\partial t}+v_0\frac{\partial }{\partial x}\right)}^{\alpha }\right]\rho (x,t)=\hfill \\ \hfill =\left[{\omega }_1{\left(\frac{\partial }{\partial t}-v_0\frac{\partial }{\partial x}\right)}^{\alpha -1}+{\omega }_2{\left(\frac{\partial }{\partial t}+v_0\frac{\partial }{\partial x}\right)}^{\alpha -1}\right]\delta \left(x\right)\delta \left(t\right).\end{array}$$

(25)

By changing the variable $s=ik$ in (24), one can perform the inverse transformation using the method described in the work [56] (see also [37]). As result, we have:

$$\rho (x,t)=\frac{2sin\left(\pi \alpha \right)}{\pi v_{0}t}\frac{{\omega }_1{\omega }_2{\left(1-{(x/v_{0}t)}^2\right)}^{\alpha -1}}{{\omega }_1^2{(1-x/v_{0}t)}^{2\alpha }+{\omega }_2^2{(1+x/v_{0}t)}^{2\alpha }+2{\omega }_1{\omega }_2\left(1-{(x/v_{0}t)}^2\right)cos\left(\pi \alpha \right)}.$$

(26)

In the case $1<\alpha <2$, the distribution (23) has a mathematical expectation. In view of the fact that the Laplace transform of density (23) takes the form (see [33]):

$$\widehat{p}\left(\lambda \right)\approx 1-m_1\lambda +m_2{\lambda }^{\alpha },\lambda \to 0,$$

where $m_1=\frac{\alpha x_0}{\alpha -1}$, $m_2=\frac{x_0^{\alpha }\Gamma (2-\alpha )}{\alpha -1}$. Now substituting this expression in (22) and considering asymptotics $\lambda \to 0$, $k\to 0$, $|\lambda /kv|\to 0$, we get:

$$\rho (k,\lambda )\approx \frac{1}{\lambda +ik{v}_0({\omega }_2-{\omega }_1)-\frac{m_2}{m_1}\left({\omega }_1{(-ik{v}_0)}^{\alpha }+{\omega }_2{\left(ik{v}_0\right)}^{\alpha }\right)}.$$

(27)

To obtain an equation describing the process of random walks, we rewrite this expression in the form:

$$\lambda \widehat{\rho }(k,\lambda )-1=\widehat{\rho }(k,\lambda )\left(ik{v}_0({\omega }_1-{\omega }_2)+\frac{m_2}{m_1}\left({\omega }_1{(-ik{v}_0)}^{\alpha }+{\omega }_2{\left(ik{v}_0\right)}^{\alpha }\right)\right).$$

Let us simplify the problem and consider symmetric random walks ${\omega }_1={\omega }_2=1/2$, we obtain:

$$\lambda \widehat{\rho }(k,\lambda )-1=\widehat{\rho }(k,\lambda )\frac{m_2}{2m_1}{v}_0^{\alpha }{\left|k\right|}^{\alpha }\left({(-i)}^{\alpha }+i^{\alpha }\right).$$

In view of the fact that $\frac{1}{2}\left({(-i)}^{\alpha }+i^{\alpha }\right)=cos\left(\frac{\pi }{2}\alpha \right)$ we get:

$$\lambda \widehat{\rho }(k,\lambda )-1=-D_{\alpha }\widehat{\rho }(k,\lambda ){\left|k\right|}^{\alpha }.$$

(28)

Here, $D_{\alpha }=\frac{m_2}{m_1}{v}_0^{\alpha }sin\left(\frac{\pi }{2}(\alpha -1)\right)$. To perform the inverse Fourier–Laplace transform of this equation, we need the Riesz fractional differentiation operator:

$$\frac{{\partial }^{\nu }}{{\partial \left|x\right|}^{\nu }}=\frac{-1}{2cos(\pi \nu /2)}\left({}_{-\infty }{D}_x^{\nu }+{}^{\infty }{D}_x^{\nu }\right).$$

The Fourier transform of this operator has the form (see, for example, [9]):

$${\int }_{-\infty }^{\infty }{e}^{ikx}\frac{{\partial }^{\nu }f\left(x\right)}{{\partial \left|x\right|}^{\nu }}dx=-{\left|k\right|}^{\nu }\widehat{f}\left(k\right).$$

Taking this relation into account, it is possible to perform the inverse Fourier–Laplace transform of Equation (28). As a result, we obtain the equation of anomalous diffusion:

$$\frac{\partial \rho (x,t)}{\partial t}=D_{\alpha }\frac{{\partial }^{\alpha }\rho (x,t)}{{\partial \left|x\right|}^{\alpha }},$$

(29)

with the initial condition $\rho (x,0)=\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 (27). It has already been done in the works [52,53]. As a result, the inverse Fourier–Laplace transform has the form:

$$\rho (x,t)=\frac{1}{{\left(D_{\alpha }t\right)}^{1/\alpha }}g\left((x-\omega vt){\left(D_{\alpha }t\right)}^{-1/\alpha };\alpha ,\omega ,0,1\right),$$

where $g(x;\alpha ,\omega ,0,1)$ is the density of the standard stable law and $\omega ={\omega }_1-{\omega }_2$. In the case of symmetric walks ${\omega }_1={\omega }_2=1/2$, this expression takes the form:

$$\rho (x,t)=\frac{1}{{\left(D_{\alpha }t\right)}^{1/\alpha }}g\left(x{\left(D_{\alpha }t\right)}^{-1/\alpha };\alpha ,0,0,1\right).$$

(30)

As we can see, the equation solution (29) is expressed in terms of a symmetrical stable law [20,57,58].

Thus, taking account of the constant velocity of the particle in the Levy walk model leads to the need to consider two cases $0<\alpha <1$ and $1<\alpha <2$. The principal difference between these two cases is that, at $0<\alpha <1$, the mathematical expectation of the distribution (23) is equal to infinity and in the case $1<\alpha <2$, the mathematical expectation is finite. It is this fact that leads to different asymptotic distributions of particles. In the case of an infinite mathematical $0<\alpha <1$ the probability of the appearance of large paths turns out to be significant, and leads to the fact that the overwhelming majority of particles move along the front of the distribution $x=\pm vt$. This is especially clear at values $\alpha <0.6$. In this case, the U—like form of the asymptotic distribution of particles is formed. As the parameter $\alpha$ increases, the share of large paths decreases. This leads to the appearance of multiply scattered particles, which form a hump in the central part of the distributions. As a result, the W—like asymptotic distribution of particles is formed. These solutions are not new and were obtained earlier when considering random walks with constant velocity [49,52,53], and when considering similar models of random walks [37,42,48,59,60] (see also the overview work [50]). In the case of the finite mathematical expectation ($1<\alpha <2$), the process of random walks of a particle with a constant velocity falls under the action of the generalized central limit theorem. As a result, the asymptotic distribution of the particle coordinate is described by a stable law. In this case, the influence of the final velocity is reduced to a decrease in the diffusion coefficient $D\to D_{\alpha }$ [34].

As we can see, in the case of the finite mathematical expectation and infinite mathematical expectation, not only are the asymptotic distributions different, but so are the equations describing the process of random walks. In the case of the finite mathematical expectation, the random walk process is described by the anomalous diffusion equation with the first time derivative (29), and the influence of the finite velocity is reduced to replacing the diffusion coefficient in this equation. In the case of the infinite mathematical expectation for taking account of the finite velocity of motion, as noted in the work [34], it is no longer sufficient to simply replace the diffusion coefficient. For this, it is necessary to replace the operator of the fractional Laplacian with the material derivative of the fractional order, which was first noted in the work [39]. Later, the authors of the works [40,41,42,43,44] come to the conclusion like this. Thus, in the case $0<\alpha <1$ in asymptotics $t\to \infty$ the random walk process is described by Equation (25).

4. Numerical Solution to a Kinetic Equation

It should be noted that the numerical methods for solving Equation (29) are well investigated (see, for example, [8,12]). For a more detailed familiarization with the methods of numerical solutions of equations with fractional derivatives we refer the reader to the reviews [61,62]. However, there are no numerical methods for solving Equation (25) or similar equations with the material derivative of the fractional order. At least, the authors are not familiar with works devoted to methods of the numerical solution of such equations.

We consider a method for a numerical solution to Equations (25) and (29), which is based on the Monte Carlo method. From Section 2, it is clear that the random walk process of a particle with a constant velocity is described by an integral transport equation, which in the one-dimensional case takes the form (20). This allows the use of the Monte Carlo method 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.

From Section 2, the simplest method of numerical solution immediately follows based on modeling trajectories and histogram estimates of the particle distribution density. Each trajectory starts at a moment in time $t=0$ from the origin of coordinates $x=0$. We will consider one-dimensional symmetric walks ${\omega }_1={\omega }_2=1/2$. Therefore, from the origin of coordinates, a particle with equal probability moves to the right or to the left at a random distance $R_1$, spending time $T_1=R_1/v_0$ on it. After that, the direction of motion is again modeled and with equal probability the particle continues to move to the right or to the left, moving a random distance $R_2$, spending time $T_2=R_2/v_0$ on it. After that, 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^{*},$$

(31)

where $T^{*}$ is a given moment in time at which it is necessary to find a solution. As soon as this condition is no longer met, the trajectory is terminated and a new trajectory begins. Random paths $R_k$, $k=1,2,\dots$ are distributed with a density,

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

Thus, random quantities $R_k$, $k=1,2,\dots$ are modeled according to the formula $R=x_0{\zeta }^{-1/\alpha }$, where $\zeta$ is a uniformly distributed random variable on a segment $(0,1]$.

To construct the simplest histogram estimate for the solution to the kinetic equation, the entire area of interest $\Delta =[a,b]$ is broken down into non-overlapping 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 (31). 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 the moment of time $T^{*}$. As a result, the estimate of the density for the interval ${\Delta }_i$ is given by the expression:

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

(32)

where summation is performed over the ensemble N of 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 like $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 estimate the probability density of detecting a particle at the point $x^{*}$ at the 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 reaching a point is calculated $(x^{*},T^{*})$, provided that the next collision is the final one.

To find this probability, let us consider the process of a random walk. Suppose that the density needs to be estimated at the point $x^{*}$ at the moment of time $T^{*}$ (see Figure 1). We will consider symmetrical random walks ${\omega }_1={\omega }_2=1/2$. The particle trajectory begins at the point $x=0$ at the moment of time $t=0$. Let the particle, as a result of a random walk at the moment of time $t^{\prime }$, reach the point $x^{\prime }$. To estimate the density at the point $x^{*}$, it is necessary to calculate the probability of transition from the point ($x^{\prime },t^{\prime }$) to the point $(x^{*},T^{*})$. It is obvious that, to get from the point ($x^{\prime },t^{\prime }$) to the point ($x^{*},T^{*}$), there are two possible alternatives of the trajectory continuation (see, Figure 1). The first option is to cover the distance ${\xi }_1$ in the positive direction of the axis $Ox$. Then, to scatter at the point ($x_1^{\prime },t_1^{\prime }$) and to cover the distance greater than ${\xi }_2$ in the negative direction of the axis $Ox$. The second option is to cover the distance ${\xi }_2$ first in the negative direction of the axis $Ox$. Then, to scatter at the point ($x_2^{\prime },t_2^{\prime }$) and to cover the distance greater than ${\xi }_1$ in the positive direction of the axis $Ox$. Since the particle is moving with a constant velocity $v_0$ all the time, then the paths ${\xi }_1$ and ${\xi }_2$ will take time ${\tau }_1={\xi }_1/v_0$ and ${\tau }_2={\xi }_2/v_0$. The values of the paths ${\xi }_1$ and ${\xi }_2$ are determined by the formulas:

$${\xi }_1=\frac{1}{2}(v_0(T^{*}-t^{\prime })+(x^{*}-x^{\prime })),{\xi }_2=\frac{1}{2}(v_0(T^{*}-t^{\prime })-(x^{*}-x^{\prime })).$$

Thus, the probability of transition from the point ($x^{\prime },t^{\prime }$) to the point ($x^{*},T^{*}$) is determined by the formula:

$$\psi (x^{*},x^{\prime },T^{*},t^{\prime })=\frac{1}{4}(p\left({\xi }_1\right)P\left({\xi }_2\right)+p\left({\xi }_2\right)P\left({\xi }_1\right)),$$

(33)

where $P\left(\xi \right)={\int }_{\xi }^{\infty }p\left({\xi }^{\prime }\right)d{\xi }^{\prime }$, and the multiplier $1/4$ appears as a result of the fact that, on the considered section of the trajectory, the particle changes its direction of motion twice. This transition probability is calculated after each particle scattering.

Theoretically, the contribution to the sought density can be made at the point $x^{*}$ by those particles that can get to the point $x^{*}$ from the current point without having a single scatter in the remaining time. Suppose that, as a result of random walks, the particle is scattered at the point $x_1^{\prime }$ at the moment of time $t_1^{\prime }$ (see Figure 1). Thus, staying at the point $x_1^{\prime }$ at the moment of time $t_1^{\prime }$, the particle can be found at the point $x^{*}$ at the moment of time $T^{*}$ if its path $\xi$ satisfies the inequality $\xi ⩾|x^{*}-x_1^{\prime }|$. In this case, the contribution to the density estimate has the form,

$${\psi }_{unsc}(x^{*},x^{\prime },T^{*},t^{\prime })=\frac{1}{2}P\left(\right|x^{*}-x^{\prime }\left|\right),$$

where we need to substitute $x^{\prime }=x_1^{\prime }$, $t^{\prime }=t_1^{\prime }$. However, it should be noted that, for such a situation to be realized, during the random walk, the particle must undergo scattering at some point exactly lying on the segments $OA$ or $OB$. Considering that the probability that the coordinate of a particle will take exactly a given value is equal to zero, then the contribution to the sought density from such particles will be zero. Nevertheless, unscattered particles make a significant contribution to the points lying on the half lines $x=\pm vt$. These half lines determine the front of the distribution. Beyond this front, the probability of detecting a particle is zero. In view of the fact that the source is at the point $x=0$, then it is obvious that the probability to detect a particle at the point $x=0$ at the moment of time $t=0$ is equal to one. Therefore, if after its appearance, the particle did not have any scattering and its path $\xi ⩾v_0{T}^{*}$, then such particles reach the points $x=\pm v_0{T}^{*}$ and form the front of the distribution. The contribution of such particles at the points $x=\pm v_0{T}^{*}$ is calculated by the formula:

$${\psi }_{unsc}(\pm v_0{T}^{*},0,T^{*},0)=\frac{1}{2}P\left(v_0{T}^{*}\right).$$

The particles will also form the front of the distribution which, after their appearance during scattering, did not change their original direction of motion. Density contribution at points $x=\pm v_0{T}^{*}$ from such particles after each collision is calculated by a similar formula:

$$\psi (\pm v_0{T}^{*},x^{\prime },T^{*},t^{\prime })=\frac{1}{2}P\left(\right|\pm v_0{T}^{*}-x^{\prime }\left|\right).$$

(34)

Here, $x^{\prime }$, $t^{\prime }$ is the coordinate and moment of time at which the particle scattered. It should be noted that if, after its appearance, the particle began its motion in the positive direction, then unscattered or multiply scattered in one direction particles contribute only at the point $v_0{T}^{*}$. Particles that began their motion in the negative direction—to the point $-v_0{T}^{*}$. 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 (x^{*},x_i,T^{*},t_i),$$

where $K\left(T^{*}\right)$ is the number of scatterings that occurred during the time interval $(0,T^{*})$. To estimate the density at the points $x\ne \pm v_0{T}^{*}$, the function $\psi (x^{*},x_i,T^{*},t_i)$ has the form (33), to estimate the density at the points $x=\pm v_0{T}^{*}$ the function $\psi (x^{*},x_i,T^{*},t_i)$ takes the form (34). We finally obtain that the density at the point $x^{*}$ at the moment of time $T^{*}$, is estimated by the formula:

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

(35)

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

In Figure 2, the results of the numerical solution to the Equation (25) are given for the values $\alpha =0.5$, $v=1,T^{*}=1000$. Circles show the results of the histogram’s estimation of the Monte Carlo method (32), asterisks are the results of the local estimate (35) and the solid line is the solution (26). The calculation results are transformed for the variable $-1⩽\zeta ⩽1$ with the help of transformation $\rho (\zeta ,t)=vt\rho (x/vt,t)$, where $-vt⩽x⩽vt$. As we can see from the figure, at the value $\alpha =0.5$, the results of all three solutions coincide. The figure also shows the advantages of the considered solution method. The results of the local estimate (35) do not contain the horizontal component of the error, which is present in the histogram estimate of the Monte Carlo method (32) and is connected with a finite quantity of the interval ${\Delta }_i$. In view of the fact that the contribution from one trajectory to the point $x^{*}$ is calculated $K\left(T^{*}\right)$ times, then this leads to a decrease in the statistical (vertical) error. In the figure, the magnitude of the statistical error does not exceed the size of the symbol. It should be noted that the solution (26) is an asymptotic solution and describes the distribution of particles at $t\to \infty$. As we can see from the figure, in the case $\alpha =0.5$ and time $T^{*}=1000$, the results of all three methods of solution coincide. This means that the process of the random walk of a particle has already entered the asymptotic regime. However, the moment when the random walk process enters the asymptotic regime for different values $\alpha$ is different.

Figure 3 contains the solution results of a kinetic Equation (20) using the method of the histogram estimation (circles), the method of the local estimation (asterisks). Solution results are given for three times $T^{*}={10}^2,{10}^3,{10}^4$. The figure also shows the solution to the Equation (25) (solid curve), which describes the asymptotic distribution of particles. As we can see in the figure, for the value $\alpha =0.8$ at the indicated times, the random walk process does not yet reach the asymptotic regime. However, when time $T^{*}$ increases, the solutions to the kinetic equation gradually approach the asymptotic distribution. It can also be seen from the figure that, at times $T^{*}={10}^3$ and $T^{*}={10}^4$, the solutions to the histogram estimate of the Monte Carlo method and the local estimate coincide, which confirms the validity of the results of the local estimation method. It can be seen from this figure that, at large times $T^{*}={10}^3$ and $T^{*}={10}^4$, there is an increase in the magnitude of the statistical error in the results of the local estimation. This increase in the calculation error is due to the fact that, at greater times, the value of the contribution (33) turns out to be a small magnitude. As a result, the calculation error can already affect it.

The solution to the kinetic equation for the case $\alpha >1$ and specified values $T^{*}$ is given in Figure 4. In this figure, the points are the results of solving the kinetic Equation (20) with the use of the local estimation method, the solid curve is the solution (30). As in the previous cases, the calculation results are given for the variable $\zeta =x/vt$, where $-vt⩽x⩽vt$. The figure shows that, at times $T^{*}=100$ and $T^{*}=1000$, the random walk process has not yet reached the asymptotic regime. The asymptotic regime of random walks is reached at $T^{*}={10}^4$. As one can see, for a given time, the asymptotic solution (30) and the local estimation results are in good agreement with each other.

One obvious shortcoming of the local estimation method should also be pointed out. In the considered method, it is necessary to calculate the contribution (33) after each scattering of a particle. This leads to an increase in the computational operations of the processor, which in turn leads to an increase in the calculation time in comparison with the histogram method for estimating the density. However, this increase in calculation time is a necessary price to pay for the clear advantages of the local estimation method. As noted earlier, the main benefit of the local estimation (35) is the possibility to find a solution at a given point $x^{*}$. This means that the results of the local estimation do not contain a systematic (horizontal) component of the error. If we give a set of points $x_i^{*},i=1,\dots ,n$, then one trajectory will contribute at once to all points of the set $x_i^{*}$. This allows the solution to be constructed as a smooth function of the coordinate x. If we also take into account that the contribution (33) is calculated after each particle scattering, then the contribution from one trajectory is calculated $K\left(T^{*}\right)$ times, where $K\left(T^{*}\right)$ is the number of scatterings of a particle in the interval $[0,T^{*}]$. This leads to a decrease in the statistical component of the error.

5. Conclusions

The use of the theory of anomalous diffusion to describe combustion processes is only at the initial stage of development. At present, there are a few works devoted to the studies undertaken in this area of research [3,6,7,8,10,12]. In all the works, to describe combustion processes, the use of the equations of anomalous diffusion (1) or (2) is proposed. As is known, these equations describe the asymptotic ($t\to \infty ,x\to \infty$) distribution of particles in the CTRW process, which is based on the assumption of the instant travel of particles from one point of space to another. This non-physical behavior of the particle leads to the fact that, at a time instant that is arbitrarily close to the initial one, the particle can be at an arbitrarily large distance from the source. Therefore, it is necessary to use these equations to describe the processes occurring in a limited area of space and develop them in time with a certain degree of caution.

An obvious way to eliminate the instantaneous movement of a particle is to introduce into consideration a constant final velocity of the movement of a particle, which was done in this work. Taking account of the constant velocity of motion shows that, depending on whether the mathematical expectation of the travel is finite or infinite, the asymptotic distribution of particles is described by completely different equations. In the case of the finite mathematical expectation of the path value ($1<\alpha <2$), the asymptotic process is described by the anomalous diffusion Equation (29), and the consideration of the finite velocity is reduced to the substitution of the anomalous diffusion coefficient $D\to D_{\alpha }$. In the case of the infinite mathematical expectation ($0<\alpha <1$), the consideration of the finite velocity leads to a completely different Equation (25) containing the fractional-order material derivative operator.

The main difficulty in using equations in fractional derivatives is to find solutions to these equations. Analytical methods for solving these equations are only at the stage of development. Therefore, the main method of solution is to apply numerical methods. In this paper, a numerical solution method is considered which is based on a local estimate of the Monte Carlo method. This method is based on the idea of modeling the trajectory of a particle’s random walk. The idea of the proposed method consists of the following. Now, after each scattering of a particle, one should calculate the probability of transition from the scattering point to a given point at which it is necessary to estimate the particle density. By giving a set of points, the transition probability should be calculated for each point from the given set. The advantages of this method over the standard histogram estimation of the Monte Carlo method are obvious. Firstly, since the solution is estimated at a specified point, then the estimation results do not contain a systematic (horizontal) component of the error. Secondly, each trajectory will contribute to all points of the given set at once. Taking account of the fact that, before the termination of the simulation of the trajectory, the particle undergoes $K\left(T^{*}\right)$ of scatterings, then one trajectory will contribute $K\left(T^{*}\right)$ times to each point of the given set. This leads to a decrease in the statistical error. However, since now the transition probability is calculated after each particle scattering, this leads to an increase in the arithmetic operations of the processor. As a result, in comparison with the histogram estimate of the Monte Carlo method, more time is spent on modeling the same number of trajectories. However, an increase in the calculation time should be considered a necessary payment for the complete absence of horizontal error.

The calculations made show that the results of the local estimation of the Monte Carlo method are in good agreement with both the results of the histogram estimation of the Monte Carlo method and the results of solving Equations (25) and (29), describing the asymptotic behavior of the process. From the calculations presented it also follows that, at different values of the parameter $\alpha$, the process of random walks reaches the asymptotic regime at different times. This indicates another advantage of using the Monte Carlo method. In fact, Equations (25) and (29) describe the asymptotic behavior of the process at $t\to \infty$, while the Monte Carlo method allows us to find a solution to the kinetic equation at any arbitrary moment in time and, thus, trace the evolution of the distribution of particles in time.