Statistical Modeling of Diffusion Processes with a Fractal Structure Ivan Sokolovskyy1[0000-0001-8477-1678] and Natalia Shakhovska1[0000-0002-6875-8534] 1 Department of Artificial Intelligence Systems, Lviv Polytechnic National University, NULP, Lviv, UKRAINE, Sokolovskyy.vani@gmail.com natalia.b.shakhovska@lpnu.ua Abstract. A discrete continuous-time random walk model for non-Markov dif- fuse processes with fractal structure is presented. On the basis of the apparatus of integro-differentiation of fractional order, finite-difference approximations of diffuse models of fractional order are obtained taking into account the effects of memory and self-organization A modification of the statistical modeling meth- od (Monte Carlo method) was carried out and an algorithm for its implementa- tion was constructed to study the diffusion process of fractional order in time. Keywords: non-integer integro-differentiation apparatus, fractal structures, dif- fusion processes, statistical modeling method, discrete model. 1 Introduction Today, the fractional intero-differential apparatus is well developed and is used to explain and simulate complex systems in nature. The development of the idea of using the fractional integro-differential apparatus to model complex systems is handled by many scientific schools in the world that are associated with the names: F. Mainardi [8], I. Podlubny [5], S. Samko, A. Kilbas [9], V. Uchajkin [1] and others. Such special attention and interest in using non-integer integro-differentiation is ex- plained by the fact that the mathematical apparatus of differentiation and fractional- order integration allows modeling of various processes and systems, which are char- acterized by the effects of memory, spatial nonlocality and self-organization. A particular advantage of fractional-differential models, as opposed to integer ones, is the ability to describe and explore more accurately real-world models with the above characteristics and effects. Fractal integro-differential parameters have been success- fully applied in the fields such as physics, biology, chemistry and biochemistry, hydrology, medicine, technology, finance. Fractional-order differential equations describe the evolution of physical systems with residual memory, which occupy an intermediate position between Markov systems and systems that are characterized by total memory. In particular, the fractionality index indicates the proportion of states of the system that persist throughout the entire process of its functioning. Copyright © 2019 for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0) 2019 IDDM Workshops. 2 It is believed that the presence of a fractional derivative with time in equations is interpreted as a reflection of a special property of the process - memory (eridarity), and in the case of a stochastic process - non-Markovian behavior. Fractional spatial derivatives reflect the self-similar heterogeneity of the structure or the medium in which the process develops. Such structures are called fractal [5]. The use of fraction- al order differential equation apparatus is important for studying the processes of anomalous diffusion in the study of anomalous properties of complex-structured in- homogeneous structures. Such structures have significant effects on memory and spatial nonlocality. 2 Analysis of Research Abnormal diffusion processes  X 2 ( )  ~ 2 K    (1   ) are characterized by a depar- ture from the linear law  X 2 ( )  ~ K   of mean-square displacement and the presence of the fractional index  depending on the time  , where K1 , K - are the usual diffusion coefficients of the dimension cm 2 * sec  1 and the generalized diffusion coefficient cm 2 * sec   , () is Gamma-function. The fractional index   1 characterizes various modes of diffusion processes: 0    1 a slow diffuse process, 1    2 an accelerated diffuse process. The case   2 is described by the wave equation. This approach describes the so-called non- Gaussian processes in dynamic systems. They are characterized by the presence of correlation dependencies for arbitrarily large space-time scales. According to [10, 11], a differential apparatus based on frac- tional-order derivatives can be used to model anomalous diffusion processes, or direct modeling of the dynamics of particles and their collisions in the system. Fractional-order differential equations are characterized by strong nonlocality and spatial correlation and are based on both spatial and evolutionary fractional differen- tial operators. A characteristic feature of fractional operators of differentiation and integration is the absence of an explicit physical and geometric interpretation of such operations [1, 9, 16]. There are several approaches to solving this problem, which can conditionally be divided into three directions: probabilistic, geometric and physical [18, 19]. The presence of different approaches to the determination of fractional derivatives give rise to ambiguity regarding the correctness and physical meaningfulness of the formulation of initial and boundary conditions depending on the type of fractional derivative [5, 7, 16]. It is also important that the fractional derivatives and integrals included in the in- tegro-differential equations and describe a certain process can be used in the sense of the Riemann-Liouville, Caputo, Wright, Weil, Grunwald-Letnikov, Marcho deriva- tives. At present, to solve fractional-order differential equations, both analytical [20, 21, 22] and numerical methods are used [2, 3, 6, 23, 24, 25]. 3 In the works [4, 20, 21, 22], found are analytical solutions to heat conduction prob- lems with boundary conditions of the first kind containing derivatives of fractional order with time and spatial variable. In particular, one-dimensional cases of problems for an infinite straight line, a semi-bounded straight line, and problems without initial conditions are considered. Relaxation processes at the phase boundary are of complex nature, which leads to nonlinear and nonlocal heat-transfer processes. However, one of the analytical methods used to solve fractional-differential equations is the Laplace transformation method [19, 22]. Analytical solutions of boundary- value problems with fractional derivatives often have considerable difficulties; therefore, numerical methods are more efficient and easier to apply. The theory of numerical methods for solving differential equations in fractional partial derivatives is fragmentary and far from being complete [26, 27]. That is why a considerable number of works is devoted to finding optimal numerical methods. 3 Formulation of Problems The fractional-order differential equation takes the form:   u ( x,  )   u ( x,  )  K (1)   x  where u ( x,  ) - function of diffusion in the region   x,   : 0  x  L,0    T . Note that the diffusion equation (1) with a fractional-order differential operator is associated with a random walk process with continuous time if the asymptotic behav- ior of the  ( ) - density of the waiting time is determined by the relation [1, 10]:   ( )   ,0    1  1 According to [12, 16], the Laplace transformation for  ( s )  exp(  s   )  1  ( s ) is characterized by asymptotics, while for 0    1 , the function  ( ) in such Laplace transformation corresponds to the conditions of distribution density. The types of some functions  ( ) are given in [16]. In par- ticular,  ( ) can be selected in the form of Mitag-Leffler functions for which the Laplace transformation has the abovegiven form. By applying the Laplace transformation with respect to the time variable and the Fourier transformation for the spatial variable, one can obtain the fractional-order diffusion equation with the fractional Kaputo differentiation operator [5]. The differential fractional-order operators in equation (1) are defined by Riemann- Liouville formulas:   u ( x, ) 1 d    1 t      u ( x, t ) dt (2)   (t   )     (   1   ) d    1   4    u ( x, ) 1 d    1 x      u ( , t ) dt (3)   (x  t) x    (   1   ) dx    1      u ( x, )   1   1 d    1   (t  x)    u ( , t )dt , (4)    (   1   )    1 x x dx where () - is Gamma function,   - is an integer part. It is known that analytical methods of implementing differential equations with fractional-order derivatives encounter great difficulties. Therefore, only for some cases, exact solutions of equation (1) were obtained mainly with boundary conditions of the first kind. Finite-difference methods are used to obtain the numerical solution. They are based on the approximation of fractional derivatives using the Grunwald-Letnikov formulas. Such fractional derivatives are a direct generalization in terms of finite differences and are determined by the dependencies for function u ( x,  ) on the interval [a, b]  u ( x)    h u ( x) 1 m    lim  lim  ( 1) k  u ( x  kh) , (5) x  h0 h  h0  k 0 h k Where m   x  a  , the value of h  0 correspond to the left-side derivatives, and  h  h  0 to the right-side. Similarly, you can write for a function on R : 1 m   lim  ( 1) k  u ( x  kh),   0 (6) h0  k 0 k h   k    where ( 1) k    k    k  1 . If the function u (x) is continuous, and its derivative u ' ( x ) is integrated оn the interval [a, x] , then the Riemann-Liouville and Grunwald-Letnikov derivatives coin- cide, in particular, with the Caputo derivatives as well [5, 16]. For further studies, we introduce a uniform grid with respect to spatial and temporal variables in the region   x,   : 0  x  L,0    T .     , h    k , x n :  k  k , k  0, K ,    t , x  ( n  1) h, n  1, N , h  K n l   N  1 Then, according to (2)-(4), the fractional derivatives of equation (1) on the grid can be approximated by the following dependencies [9, 17]: 5   u( xn , j) 1   k      u( xn , j  k )     k  1 (7)    h k 0    u( x n , j) 1 n  1 k      u( x n  k  1, j)    k  1 (8) x  h  k  0    u( x n , j) 1 N  n  1 k      u( x n  k  1, j)  k  0   k  1 (9) x  h The above approximations using shifted Grunwald-Letnikov formulas allow ob- taining conditionally stable explicit and stable implicit first-order accuracy schemes for fractional differential equations. In fractal-structured media, it was shown in [12, 28] that the fundamental solutions of fractional differential equations with respect to temporal and spatial variables are characterized by properties which are intrinsic to the distribution densities of random variables. 4 Simulation of Random Walks Discrete models of Markov random walk for the classical Brownian motion are the basis for the application of statistical methods for studying conventional diffusion processes [10, 16]. For a one-dimensional case where displacements are possible only at two nearest points, it is possible to write u ( x,    )  p1 u ( x  h,  )  p 2 u ( x  h,  ) Where p1 , p 2 are the probabilities of particle displacement by one step k th step the p1  p 2  1, x  nh,  k , p n  u ( nh, k ) - is the probability that at the n process is at the k th point. u ( x,  )  a 2  u ( x,  ) , a  const 2  2 x The implementation of the statistical test method involves the use of appropriate difference schemes [4, 14, 16]. The diffusion equation in the region   x,   : 0  x  L,0    T  takes the form: u ( x,  )  u ( x,  ) 2 (10)  a2 , a  const  2 x To use the statistical test method in order to study the model (10), we use the Crank-Nicholson difference scheme [4]. In this case, the algorithm for calculating the probability of a random variable location at a point of time  j takes the form: 6 j 1  1  un 2a 1 2 (11) h     2 n 1 n 1 u u    a j  1 j  1 (1   ) a j u n 1  u  j n  1 1     2(1   ) a  j    u n    h h 2  h 2   where x n  (n  1)h, j  j , , h - uniform splitting steps, n, j - numbers of splitting nodes,   0;1 . According to [13], the value u nj  1 can be interpreted as the probability of a ran- dom variable location at the point x n at the time  j . That is, over a period of time [ i ,  i  1] we can consider the Markov process with corresponding probabilities:  j  p uj (12) un  k  k n  k The coefficients p k are determined from the difference relation (12) and are the probability of uniform random walks of some particle M around the nodes of the difference grid (11) of approximation of the diffusion equation (1). In particular, for a six-point pattern for two time j and j  1 , such probabilities take the form: h  2(1   ) a (1   ) a a 2 p0  , p1  , p 2  , p k  0, k  3,4,... (13) 2  2a 2  2a 2  2a h h h In particular, p  2 in formula (13) corresponds to the value j  1 . In addition, the  transfer coefficients satisfy the condition  p  1 as well as the stability condi- k  k tion [17] of the difference scheme (11) for   0.5 . The relations (12), (13) characterize the standard random walk model for the Gaussian process. It is believed that a random particle location in the internal node of the difference scheme can move to neighboring nodes, that is, to carry out the move- ment of a unit length or remain in the location node, namely, to perform a zero-length step. 5 Simulation of Random Walks for Fractional-Diffusion Process To construct a discrete model of random walk for the differential equation of diffu- sion of fractional order (1), we use [12, 13, 15]. Equally important in this respect is to establish the relationship of the fundamental solution of the fractional-order differen- tial equation (1) with the time variable of the fractionally stable distributions [16]. In particular, the fundamental solution of the equation (*) with the fractional differentia- tion operator Caputo: 7 G ( x,  )    1  1 n x n   n  11    n  n n 2a n  0 a  is a probability density function of the fractionally stable distribution [9, 13]. For   0.5 we get a Gaussian solution [16]. Then, in the equation (11) obtained in terms of the relations in discrete form, you can move from the function u ( x,  ) to the probability density (or relative diffusion concentration). This, in turn, makes it possible to move on to the probability of a par- ticle location at a nodal point of a discrete grid. In this regard, we introduce the desig- nation. Given the relation (7)-(9) according to [13, 17], we can write:    m     n  m   K     m  0   m 1 Vn h   (14)    i    k 1   i    k 1  C  V n  i  1  C  V n  i  1  i  0     i  1 i  0     i  1  k k 1     k  m   K    m   Vn  Vn m  2   m 1   Vn h  (15)   C  i  k 1   i k 1    V n  i  1  C  V n  i  1    i  0    i 1  i  0    i 1     The relation for determining V kn can be considered as a modeling scheme for a random process with discrete time. The value V kn characterizes the probability of the particle location at the point x n at the time  k during a random walk in the differ- ence grid. The coefficients V kn correspond to the probability of transitions in space and time. We denote them by p1i , p 2 m . Then (15) can be written as:   k 1  p  p1i V kn  km (16) Vn  i 2m V n i   m2 Using the relation (16) for each case of random walk, we can get a finite number of values p and p . Іn addition, these values must satisfy the conditions of non- 1i 2m negativity, normalizing, which are typical of probabilistic characteristics, that is     p1i  0, p 2m  0,  p1i   p 2m  1 . The conditions k V n  0,  V ik   V i0 i   m2 i   i   must also be met. The second condition ensures the preservation of the total number of particles. 8 6 Numerical Experiment Numerical experiments were performed for materials with density   460kg / m 3 , K 2  0.5 . For the equation () was specified boundary conditions u( , )  u  b,   . Initial conditions: u( x,0)  x  a , if 0  x  (b  a) / 2 , u( x,0)  b  x , (b  a) / 2  x  b . Fig. 1. Changing the function U over time for   0.7 with different values x . Fig. 2. Changing the function U over time for   0.9 with different values x . Quantities h  0.1 ,   0.1 . Figures 1 and 2 show the dependences of the function change u ( x,  ) , respectively, for the fractal coefficient   0.7 and   0.9 for dif- ferent values x , which are respectively x  0;0.25;0.5;0.75;1. 9 The analysis of graphic dependences indicates the effect of the parameter  on the function change u ( x,  ) . Increasing the parameter  increases the number of maxi- mum values of the curve for different values of change of coordinates. Conclusions On the basis of discrete models of diffusion processes with fractional derivatives shows a modification of the method of random walk on the implementation of the mathematical model of a diffusion process subject to temporal nonlocality. For nu- merical implementation of such a mathematical model, finite-difference schemes are proposed, an algorithm is developed and software is created. The results of numerical experiments for the implementation of a mathematical model of diffusion processes for different values of the fractional index according to time are given by the statisti- cal method. References 1. Uchajkin, V.: Method of fractional derivatives. Ulyanovsk: Publishing house “Artishok”, 512 p. (2008). 2. Sokolowskyi, Y., Shymanskyi, V., Levkovych, M., Yarkun, V.: Mathematical and soft- ware providing of research of deformation and relaxation processes in environments with fractal structure. Proceedings of the 12th International Scientific and Technical Conference on Computer Sciences and Information Technologies, CSIT 2017, Lviv, Vol. 1, pp. 24-27 3. Sokolovskyy, Y., Levkovych, M., Mokrytska O., Atamanvuk, V.: Mathematical Modeling of Two-Dimensional Deformation-Relaxation Processes in Environments with Fractal Structure. 2nd IEEE International Conference on Data Stream Mining and Processing, DSMP 2018, Lviv, pp. 375-380. 4. Beibalaev V.D.: The numerical method of solving regional problems for a two- dimensional weekend. The solution of thermal conductivity with higher order. Vestn. Sam. state tech. un-that. Ser. Phys.-mat. Sciences. №5 (21), pp. 244-251. (2010) 5. Podlubny, I.: Fractional Differential Eguations. Vol. 198 of Mathematics in Science and Engineering, Academic Press, San Diego, Calif, USA. (1999). 6. Sokolovskyy, Y., Shymanskyi, V., Levkovych, M.: Mathematical Modeling of Non- Isothermal Moisture Transfer and Visco-Elastic Deformation in the Materials with Fractal Structure. 11th International Scientific and Technical Conference on Computer Sciences and Information Technologies, CSIT 2016, Lviv, pp. 91-95. 7. Sokolovskyy, Y., Levkovych, M., Mokrytska O., Kaplunskyy, Y.: Mathematical models of biophysical processes taking into account memory effects and self-similarity. 1st Interna- tional Workshop on Informatics and Data-Driven Medicine, IDDM 2018, Lviv, pp. 215- 228. 8. Mainardi F.: The Fundamental Solutions for the Fractional Diffusion-Wave Equation Appl. Math. Lett. Vol. 9, No. 6, pp. 23-28. (1996) 9. Samko S.G., Kilbas A.A., Marichev O.I.: Integrals and Virology and Other Orders and Ac- tions. Minsk. Science and technology. 688 p. (1987). 10. Gorenflo R., Mainardi F., Moretti D., Pagnini G., Paradisi P.: Discrete random walk mod- els for space-time fractional diffusion. Chemical Physics. Vol. 284, pp. 521-544 (2002) 10 11. Gorenflo R., Mainardi F., Moretti D., Paradisi P.: Time fractional diffusion: a discrete random walk approach. Nonlinear Dynamics. Vol. 29, pp. 129-143. (2002) 12. Gorenflo R., Vivoli A., Mainardi F.: Discrete and continuous random walk models for space-time fractional diffusion. Nonlinear Dynamics. Vol 38, pp. 101-116. (2004) 13. Gorenflo R., Mainardi F.: Power laws, random walks, and fractional diffusion processes as well scaled refinement limits. Fractional Differentiation and its Applications. Pp. 497–516. (2005) 14. Bisquert J.: Fractional diffusion in the multiple-trapping regime and revision of the equiv- alence with the continuous-time random walk. Physical Review Letters. Vol. 91, №1. (2003) 15. Bondarenko A. N., Ivaschenko D. S., Seleznev V. A.: Inverse Sommerfeld Problem for Fractal Media. Proceedings of 6th Korea—Russia International Symposium on Science and Technology KORUS'02, June 24 — 30, 2002 at Novosibirsk State Technical Universi- ty. Novosibirsk. Vol. 1, pp. 246—252. (2002) 16. Uchaikin V.V.: Anomalous self–similar diffusion and Levy–stable laws. Physics–Uspekhi. Vol. 173, №8, pp. 847-876 (2003) 17. Reviznikov D. L., Slastushenskiy Y. V.: Numerical simulation of anomalous diffusion in polygonal billiard gas channel. Matem. Mod. Vol. 25:5, pp. 3–14. (2013) 18. Machado J.A.T., Galhano A.M.S.F.: Fractional order inductive phenomena based on the skin effect. Non-linear Dyn. Vol. 68(1–2), pp. 107–115. (2012) 19. Podlubny I.: Fractional Differential Equations. Mathematics in Science and Engineering, Academic Press, San Diego, Calif, USA, Vol. 198. 340 p. (1999) 20. Datsko B. Y., Gafiychuk V.V.: Different types of instabilities and complex dynamics in reaction-diffusion systems with fractional derivatives. Computational and Nonlinear Dy- namics. DOI No: CND-09-1119. (2012) 21. Gafiychuk V., Datsko B.: Mathematical modeling of different types of instabilities in time fractional reaction-diffusion systems. Computers and Mathematics with Applications. Vol. 59, pp. 1101-1107. (2010) 22. Kexue L., Jigen P.: Laplace transform and fractional differential equations. Appl. Math. Lett. Vol. 24, pp. 2019-2023. (2011) 23. Al-Khaled K.: Numerical solution of time-fractional partial differential equations using Sumudu decomposition method. Rom. Journ. Phys. Bueharest. Vol. 60 (1-2), pp. 99-110. (2015) 24. Diethelm K.: An algorithm for the numerical solution of differential equations of fractional order. Electron. Trans. Numer. Anal. Vol. 5, pp. 1-6. (1997) 25. Ford N.J., Simpson A.: Numerical and analytical treatment of differential equations of fractional order. In Proc. Internat. Conf. on Sci. Comp. & Math. Modeling, Inst. Of Ap- plied Science and Computations, Charleston. Pp. 60-65. (2000) 26. Chen Y.Q., Moore K.L.: Analytical stability bounded for a class of delayed fractional- order dynamic systems. Nonlinear Dynam. Vol. 29. Pp. 191-200. (2002) 27. Ford N.J., Simpson A.: The numerical solution of fractional differential equations: speed versus accuracy. Numer. Algorithms, Vol. 26(4). Pp. 333-346. (2001) 28. Machado J.A.T.: Entropy analysis of integer and fractional dynamical systems. Nonlinear Dyn. Vol. 62(1–2). Pp. 371–378. (2010)