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

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.


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 twoterm 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 where 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][17][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 twosided 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][31][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 Section 2, 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.

Stochastic Representation
In this section, we will give the stochastic representation of the multiterm time-fractional diffusion equation.
Let U α ðτÞ be the increasing Lévy motion with Laplace transform Ee −kU α ðτÞ = e −τk α , then we get the following theorem. Theorem 1. The subordinated process YðtÞ = XðE t Þ is the stochastic representation of the multiterm time-fractional advection-diffusion equation (1), where the parent process and subordinator of Y t is defined as respectively. Here, E t is independent of BðτÞ and A τ = ∑ r i=0 U α i ðl i τÞ, 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 E t and the PDF uðt, τÞ of A τ . From the definition of E t (see Equation (4)), we have PðE t < τÞ = PðA τ > tÞ [34], therefore So, the Laplace transform of gðτ, tÞ can be expressed aŝ Here, the Ee −kA τ = e −τ∑ r i=0 l i k α i is obtained from the definition of AðτÞ. Then, using the total probability formula and the independence between XðτÞ and E t , we can get the PDF pðx, tÞ of YðtÞ, given by where f ðx, τÞ is the PDF of the parent process XðτÞ. So the Laplace transform of the above equation yieldŝ Thus, the following relation between f ðx, τÞ and pðx, τÞ holdsp

Advances in Mathematical Physics
Since the process XðτÞ is given by the Itô stochastic differential equation (3), its PDF f ðx, τÞ obeys the classical advection-dispersion equation [35] ∂f The Laplace transform of the above equation with respect to τ yields By changing the variable k to ∑ r i=1 l i k α i and using the relation (9), the above equation yields Since the Laplace transform of the Caputo fractional derivative is given by LfD α t f ðtÞ, sg = s αf ðsÞ − s α−1 f ð0Þ, by comparing the above equation with the Laplace transform of Equation (3), we get Cðx, tÞ = pðx, tÞ. So the subordinated process YðtÞ = XðE t Þ is called the stochastic representation of the multiterm time-fractional diffusion equation.

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. Proof. Taking advantage of the relation between pðx, tÞ and f ðx, τÞ, we can get the following evaluation formula for the mean square displacement where gðτ, tÞ is the PDF of E t and its Laplace transformation is given by Equation (4). By inverting the Laplace transform, we can get 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., In order to calculate the mean square displacement, we can first compute its Laplace transform. By using the Laplace transform of gðτ, tÞ Equation (4), and inverting the Laplace transform, we have where E α,β ðxÞ = ∑ ∞ k=0 ðx k /ðΓðαk + βÞÞÞ is the Mittag-Leffler function [36]. Here, we have used the equality 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 of XðtÞ of Equation (1) in the presence of a uniform force field VðxÞ = F, 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]) To connect to single-particle tracking experiments, we now turn to the time-averaged MSD of the stochastic process, defined by 3 Advances in Mathematical Physics The time series xðtÞ 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 x ðtÞ. Typically, δ 2 ðΔÞ is considered in the limit Δ ≪ T to obtain good statistics. It is easy to show that for Brownian motion, hδ 2 ðΔÞi = hx 2 ðΔÞi = 2K 1 Δ 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. hδ 2 ðΔÞi ≠ hx 2 ðΔÞi signifies weak ergodicity breaking [38][39][40].  (1) is weak ergodicity breaking.
Proof. From Equation (16), we have According to the definition of time-averaged MSD, we calculate the δ 2 ðΔÞ of YðtÞ.

Stochastic Simulation
The stochastic representation provides two ways to get the solution of the multiterm time-fractional advectiondiffusion equation (1). One way is to get the analytical solution by substituting f ðx, τÞ 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 timefractional diffusion equation.
The algorithm of simulation of the subordinated process XðE t Þ is divided into two steps. T is the horizon and Δt = T/N.
Step 1. This step is aimed at simulating the subordinator E t (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: where Z j is the i.i.d. strictly increasing α i -stable Levy noise [42,43], generated by 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 E t is the first passage time of A τ . So, for every element t i , we only need to find the element τ j such that A τ j−1 < t i ≤ A τ j , then E t i = τ 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 E t : These heavy-tailed flat periods of E t represent long waiting times in which the subdiffusive particle gets immobilized in the trap. The sample path of E t can be found in Figure 1, and the corresponding waiting time in Figure 2. From the

Advances in Mathematical Physics
figures, we find the subordinator E t stop at the same time, which leads the waiting time to be fluctuant.
Step 2. This step is aimed at simulating the process YðtÞ = XðE t Þ. Since the parent process is driven by Brown motion, then we employ the Euler scheme to simulate the process X ðE t Þ, given by where ξ i is the i.i.d. standard normal noise, ξ i~N ð0, 1Þ.
The sample path of XðE t Þ 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 log ðpðx, 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.

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.

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that no conflicts of interest exist regarding this manuscript.   (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.