# Stochastic Representation and Monte Carlo Simulation for Multiterm Time-Fractional Diffusion Equation

Authors: Longjin Lv; Luna Wang
Journal: Advances in Mathematical Physics
 (2020)

---

## Abstract

In this paper, we mainly study the solution and properties of the multiterm time-fractional diffusion equation. First, we obtained the stochastic representation for this equation, which turns to be a subordinated process. Based on the stochastic representation, we calculated the mean square displacement (MSD) and time average mean square displacement, then proved some properties of this model, including subdiffusion, generalized Einstein relationship, and nonergodicity. Finally, a stochastic simulation algorithm was developed for the visualization of sample path of the abnormal diffusion process. The Monte Carlo method was also employed to show the behavior of the solution of this fractional equation.


---

## Body

## 1. Introduction

Recently, the diffusion equations that generalize the usual one have received considerable attention due to the broadness of their physical applications, in particular, to the anomalous diffusion. In fact, fractional diffusion equations and the nonlinear fractional diffusion equations have been successfully applied to several physical situations such as percolation of gases through porous media [1], thin saturated regions in porous media [2], standard solid-on-solid model for surface growth [3], thin liquid films spreading under gravity [4], in the transport of fluid in porous media and in viscous fingering [5], modeling of non-Markovian dynamical processes in protein folding [6], relaxation to equilibrium in a system (such as polymer chains and membranes) with long temporal memory [7], and anomalous transport in disordered systems [8], diffusion on fractals [9], and the multiphysical transport in porous media, such as electroosmosis [10, 11]. Moreover, some underlying processes can be more accurately and flexibly modeled by multiterm FPDEs. For example, the multiterm time-fractional diffusion-wave equation is a satisfying mathematical model for viscoelastic damping [12]. In [13], a two-term fractional diffusion equation has been successfully used for distinguishing different states in solute transport.The multiterm time-fractional advection-diffusion equations are linear integrodifferential equations. They are obtained from their corresponding classical multiterm advection-diffusion equations by replacing the first-order time derivative by fractional derivative, which reads(1)PDtCx,t=−v∂∂xCx,t+D∂2∂x2Cx,t,where PDt=∑i=1rliDtαi with 0<αr<⋯<α1<1; here, the Dtαi is the Caputo fractional derivative of order αi with respect to t as defined as
(2)Dtαift=1Γ1−αi∫0tf′τt−ταidτ,0<αi<1.There are growing interests in studying these equations because of their importance in modeling many physical, biological, medical, chemical, and many other fields. For example, the over diagnostic ultrasound frequencies, acoustic absorption in biological tissue exhibits a power law with a noninteger frequency [14, 15]. Also, in a complex inhomogeneous conducting medium, experimental evidence shows that the sound waves propagate with the power law of noninteger order. For further applications on physics and on real phenomena [16–18], the Caputo time-fractional operator has been widely used instead of the second time derivative to model mathematically such problems in order to discuss the effect of the memory on the studied system [19].Many analytic and numeric methods are employed to solve this equation. Daftardar-Gejji and Bhalekar considered the multiterm time-fractional diffusion-wave equation using the method of separation of variables [20]. Luchko [21] studied the well-posedness of the multiterm time-fractional diffusion equation based on an approximate maximum principle. Jiang et al. studied the multiterm time-space fractional advection-diffusion equation based on the spectral representation of the fractional Laplacian operator [22]. Meshless analysis based on improved moving least-squares approximation was introduced to solve the two-dimensional two-sided space-fractional wave equation in [23]. Using the method of series expansion, Ye et al. [24] studied the multiterm time-space fractional partial differential equations in 2D and 3D domains. An efficient operational formulation of the spectral tau method for a multiterm time-space fractional differential equation with Dirichlet boundary conditions was proposed in [25]. Liu et al. presented numerical approximations for multiterm time-fractional diffusion equations by using the spectral method in [26] and for multiterm time-fractional wave equations by means of FDMs in [27], respectively. In [28], the authors used finite difference rules to get the approximate solutions of the time-fractional multiterm wave equations. Recently, the stochastic representation method was introduced to solve the fractional diffusion equation. In [29], Kolokoltsov built the relation between the stochastic process and time-fractional diffusion equations with Caputo or Riemann-Liouville derivatives. These generalized Caputo derivatives were further extended to nonmonotone processes, yielding two-sided and even multidimensional extensions. Based on the stochastic representation, the mathematical properties of the related fractional equation were discussed in [30–32]. These papers inspired the research of this paper.In this paper, we introduce the stochastic representation method to solve this multiterm time-fractional diffusion equation. This paper is organized as follows. In Section2, we derive a subordinated process, whose PDF is rightly the solution of this equation, where the parent process is a classical diffusion process and the subordinator is the inverse time of the sum of Lévy motions with different parameter. Taking advantage of this result, we study the properties of this multiterm time-fractional diffusion equation in Section 3. We also employ the Monte Carlo method to simulate the solution for this equation in the next section. Section 5 presents our conclusions.


## 2. Stochastic Representation

In this section, we will give the stochastic representation of the multiterm time-fractional diffusion equation.LetUατ be the increasing Lévy motion with Laplace transform Ee−kUατ=e−τkα, then we get the following theorem.Theorem 1.
The subordinated processYt=XEt is the stochastic representation of the multiterm time-fractional advection-diffusion equation (1), where the parent process and subordinator of Yt is defined as
(3)dXτ=vdτ+2DdBτ,(4)Et=infτ>0:Aτ>t,respectively. Here, Et is independent of Bτ and Aτ=∑i=0rUαiliτ, Uαiτ are also independent with each other for different i.Proof.
Following the same procedure shown in [33], we first establish the relation between the PDF gτ,t of Et and the PDFut,τ of Aτ. From the definition of Et (see Equation (4)), we have PEt<τ=PAτ>t [34], therefore
(5)gτ,t=∂∂τ∫t∞uy,τdy=−∂∂τ∫0tuy,τdy.
So, the Laplace transform ofgτ,t can be expressed as
(6)g^τ,k=∫0∞e−ktgτ,tdt=−∂∂τ∫0∞uy,τ∫y∞e−ktdtdy=−1k∂∂τ∫0∞e−kyuy,τdy=−1k∂∂τEe−kAτ=∑i=1rlikαi−1e−τ∑i=1rlikαi.
Here, theEe−kAτ=e−τ∑i=0rlikαi is obtained from the definition of Aτ. Then, using the total probability formula and the independence between Xτ and Et, we can get the PDF px,t of Yt, given by
(7)px,t=∫0∞fx,τgτ,tdτ,where fx,τ is the PDF of the parent process Xτ. So the Laplace transform of the above equation yields
(8)p^x,k=∫0∞fx,τg^τ,kdτ=∑i=1rlikαi−1∫0∞fx,τe−τ∑i=1rlikαidτ.
Thus, the following relation betweenfx,τ and px,τ holds
(9)p^x,k=∑i=1rlikαi−1f^x,∑i=1rlikαi.
Since the processXτ is given by the Itô stochastic differential equation (3), its PDF fx,τ obeys the classical advection-dispersion equation [35]
(10)∂fx,τ∂τ=−∂∂xv+D∂2∂x2fx,τ.
The Laplace transform of the above equation with respect toτ yields
(11)kf^x,k−fx,0=−∂∂xv+D∂2∂x2f^x,k.
By changing the variablek to ∑i=1rlikαi and using the relation (9), the above equation yields
(12)∑i=1rlikαi−1kp^x,k−px,0=−∂∂xv+D∂2∂x2p^x,k.
Since the Laplace transform of the Caputo fractional derivative is given byLDtαft,s=sαf^s−sα−1f0, by comparing the above equation with the Laplace transform of Equation (3), we get Cx,t=px,t. So the subordinated process Yt=XEt is called the stochastic representation of the multiterm time-fractional diffusion equation.


## 3. Some Properties

The superiority of the stochastic representation approach to the fractional differential equation is that it not only helps us to understand the physical process by providing a description of the dynamical system governed by the fractional differential equation but also provides a way to get the properties of the corresponding equations. For simplicity, we retain two terms for the time-fractional operator, i.e.,0<α2<α1≤1.Corollary 2.
The subordinated process governed by the multiterm time-fractional diffusion equation (1) is subdiffusive.Proof.
Taking advantage of the relation betweenpx,t and fx,τ, we can get the following evaluation formula for the mean square displacement
(13)X2Et=∫0∞gτ,t∫0∞x2fx,τdxdτ=2D∫0∞τgτ,tdτ,where gτ,t is the PDF of Et and its Laplace transformation is given by Equation (4). By inverting the Laplace transform, we can get
(14)gτ,t=tα1Γα1+1+tα2Γα2+1∗gα1τ,t∗gα2τ,t.
Here,gαiτ,t, the PDF of inverse time αi-stable Lévy motion, can be expressed in the form of a Fox function, i.e.,
(15)1τH1110τtαi1,11,αi.
In order to calculate the mean square displacement, we can first compute its Laplace transform. By using the Laplace transform ofgτ,t Equation (4), and inverting the Laplace transform, we have
(16)X2Et0=2Dtα1Eα1−α2,α1+1−tα1−α2∝t2Dα1Γα1+1,t≪1,t2Dα2Γα2+1,t≫1,where Eα,βx=∑k=0∞xk/Γαk+β is the Mittag-Leffler function [36]. Here, we have used the equality
(17)∫0∞e−pttαk+β−1Eα,βk±atαdt=k!pα−βpα∓ak+1.
From Equation (16), we can know the model resembles a α1 subdiffusion for t⟶0+, and since 0<α<1, the model resembles a α2 subdiffusion for t⟶∞. This result can be got in another way; to see this, note that Uατ=τ1/αUα1 in distribution and τ1/α2 grows faster than τ1/α1 for 0<α2<α1<1, so the α2-stable subordinator dominates as τ⟶∞ and the α1-stable subordinator dominates as τ⟶0+.Corollary 3.
The generalized Einstein relation holds for the multiterm time-fractional diffusion equation (1).Proof.
With the help of the stochastic representation, we can calculate the first moment ofXt of Equation (1) in the presence of a uniform force field Vx=F,
(18)XEtF=FEt=Ftα1Eα1−α2,α1+1−tα1−α2.
Comparing the above result with Equation (16) shows that the generalized Einstein relation holds (the definition of generalized Einstein relation can be found in [37])
(19)XEtF=const×X2Et0.
To connect to single-particle tracking experiments, we now turn to the time-averaged MSD of the stochastic process, defined by(20)δ2Δ¯=1T−Δ∫0T−Δxt+Δ−xt2dt.
The time seriesxt of length T (the measurement time) is thus evaluated in terms of squared differences of the particle position separated by the so-called lag time Δ, which defines the width of the window slid along the time series xt. Typically, δ2Δ¯ is considered in the limit Δ≪T to obtain good statistics. It is easy to show that for Brownian motion, δ2Δ¯=x2Δ=2K1Δ as long as the measurement is sufficiently long. Therefore, we call the process ergodic: ensemble averages and long-time averages are equivalent in the limit of long measurement times. δ2Δ¯≠x2Δ signifies weak ergodicity breaking [38–40].Corollary 4.
The subordinated processYt=XEt governed by the multiterm time-fractional diffusion equation (1) is weak ergodicity breaking.Proof.
From Equation (16), we have
(21)Y2Δ∝2DΓα1+1Δα1.
According to the definition of time-averaged MSD, we calculate theδ2Δ¯ of Yt.
(22)δ2Δ¯=1T−Δ∫0T−ΔYt+Δ−Yt2dt=2DTα1−1Eα1−α2,α1−Tα1−α2Δ≠Y2Δ.
The disparity between the ensemble and the -averaged MSD exists. Even in the limit of long measurement timesT⟶∞, limT→∞δ2Δ¯∝2DTα2−1/Γα2+1Δ.£¬, and therefore, the disparity still exists, which ends our proof.


## 4. Stochastic Simulation

The stochastic representation provides two ways to get the solution of the multiterm time-fractional advection-diffusion equation (1). One way is to get the analytical solution by substituting fx,τ and gτ,t into Equation (7). The other way is to simulate the stochastic representation, then use Monte Carlo to simulate the solution. The Monte Carlo method is firstly proposed to simulate the solution of fractional order equation [41]. Here, we mainly introduce how to simulate the sample path of the stochastic representation and get the simulated solution of the multiterm time-fractional diffusion equation.The algorithm of simulation of the subordinated processXEt is divided into two steps. T is the horizon and Δt=T/N.Step 1.
This step is aimed at simulating the subordinatorEt (see Equation (4)). Since the Uαiτ is the strictly increasing αi-stable Levy motion with independent increments, then the process Uαiτ on the mesh τj=jΔτ (j=0,1,⋯,n) can be simulated as follows:
(23)Uαiτ0=0,Uαiτj=Uαiτj−1+Δτ1/αiZj,where Zj is the i.i.d. strictly increasing αi-stable Levy noise [42, 43], generated by
(24)Zj=sinαU+π/2cosU1/a×cosU−αU+π/2ω1−a/a.
Here,U is uniform distribution on −π/2,π/2, and ω is exponential distribution with mean 1. Then, we can get the simulation of the process Aτj=Uα1τj+Uα2τj. From the definition of (4), we know the subordinator Et is the first passage time of Aτ. So, for every element ti, we only need to find the element τj such that Aτj−1<ti≤Aτj, then Eti=τj. Since Uαiτ is a pure-jump process. For every jump of Uα1τ+Uα2τ, there is a corresponding flat period of its inverse Et. These heavy-tailed flat periods of Et represent long waiting times in which the subdiffusive particle gets immobilized in the trap. The sample path of Et can be found in Figure 1, and the corresponding waiting time in Figure 2. From the figures, we find the subordinator Et stop at the same time, which leads the waiting time to be fluctuant.Step 2.
This step is aimed at simulating the processYt=XEt. Since the parent process is driven by Brown motion, then we employ the Euler scheme to simulate the process XEt, given by
(25)X0=0,XEti=XEti−1+vEti−Eti−1+2DEti−Eti−11/2ξi,where ξi is the i.i.d. standard normal noise, ξi~N0,1. The sample path of XEt can be found in Figure 3. Then, the Monte Carlo method can be employed to estimate the solution of Equation (1) (see Figure 4). From the figure, we find that the solution of Equation (1) has sharp peak and heavy tails, in contrast with normal distribution, which is called the stretched Gaussian distribution. The figure for logpx,t is plotted to show these results more clearly (see Figure 4). Here, we remark that all the numerical results are obtained by the software Matlab.Figure 1
The sample paths of the stochastic processEt, where α1=0.9 and α2=0.5.Figure 2
The corresponding waiting times withα1=0.9 and α2=0.5.Figure 3
The sample paths of the stochastic processXEt, where α1=0.9, α2=0.5, v=0.05, and D=0.5.Figure 4
The approximated solutions (up) of Equation (1) with different α1 and α2 are obtained for t=1. Noting that we estimate 5000 times to get these solutions. Below is the corresponding Logarithmic solution.


## 5. Conclusions

In this paper, an advection-diffusion equation with multiterm time-fractional derivatives is employed. We obtained its stochastic representation, which is driven by the Brown motion and the inverse time of the sum of Lévy motions with different parameters. Then, the mean square displacement indicates the model is subdiffusive and the generalized Einstein relation is also retained, but weak ergodicity is breaking. At last, an algorithm is constructed to simulate the sample paths of the stochastic process. With the help of stochastic representation, the Monte Carlo method is employed to approximate the solution of the corresponding equation. We find that the solution is heavy tailed and sharp peaked, which is common in statistical physics and finance. So, we expect that the results obtained here may be useful for the discussion of the anomalous diffusion systems.



---

*Source: 1315426-2020-07-14.xml*