#### 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.

#### 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 where with ; here, the is the Caputo fractional derivative of order with respect to as defined as

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 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.

#### 2. Stochastic Representation

In this section, we will give the stochastic representation of the multiterm time-fractional diffusion equation.

Let be the increasing Lévy motion with Laplace transform , then we get the following theorem.

Theorem 1. *The subordinated process is the stochastic representation of the multiterm time-fractional advection-diffusion equation (1), where the parent process and subordinator of is defined as
respectively. Here, is independent of and , are also independent with each other for different .*

*Proof. *Following the same procedure shown in [33], we first establish the relation between the PDF of and the *PDF* of . From the definition of (see Equation (4)), we have [34], therefore
So, the Laplace transform of can be expressed as
Here, the is obtained from the definition of . Then, using the total probability formula and the independence between and , we can get the PDF of , given by
where is the PDF of the parent process . So the Laplace transform of the above equation yields
Thus, the following relation between and holds
Since the process is given by the Itô stochastic differential equation (3), its PDF obeys the classical advection-dispersion equation [35]
The Laplace transform of the above equation with respect to yields
By changing the variable to and using the relation (9), the above equation yields
Since the Laplace transform of the Caputo fractional derivative is given by , by comparing the above equation with the Laplace transform of Equation (3), we get . So the subordinated process 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., .

Corollary 2. *The subordinated process governed by the multiterm time-fractional diffusion equation (1) is subdiffusive.*

*Proof. *Taking advantage of the relation between and , we can get the following evaluation formula for the mean square displacement
where is the PDF of and its Laplace transformation is given by Equation (4). By inverting the Laplace transform, we can get
Here, , the PDF of inverse time -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 Equation (4), and inverting the Laplace transform, we have
where is the Mittag-Leffler function [36]. Here, we have used the equality
From Equation (16), we can know the model resembles a subdiffusion for , and since , the model resembles a subdiffusion for . This result can be got in another way; to see this, note that in distribution and grows faster than for , so the -stable subordinator dominates as and the -stable subordinator dominates as .

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 of Equation (1) in the presence of a uniform force field ,
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
The time series of length (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 . Typically, is considered in the limit to obtain good statistics. It is easy to show that for Brownian motion, 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. signifies weak ergodicity breaking [38–40].

Corollary 4. *The subordinated process governed by the multiterm time-fractional diffusion equation (1) is weak ergodicity breaking.*

*Proof. *From Equation (16), we have
According to the definition of time-averaged MSD, we calculate the of .
The disparity between the ensemble and the -averaged MSD exists. Even in the limit of long measurement times , £¬, 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 and 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 process is divided into two steps. is the horizon and .

*Step 1. *This step is aimed at simulating the subordinator (see Equation (4)). Since the is the strictly increasing -stable Levy motion with independent increments, then the process on the mesh () can be simulated as follows:
where is the i.i.d. strictly increasing -stable Levy noise [42, 43], generated by
Here, is uniform distribution on , and is exponential distribution with mean 1. Then, we can get the simulation of the process . From the definition of (4), we know the subordinator is the first passage time of . So, for every element , we only need to find the element such that , then . Since is a pure-jump process. For every jump of , there is a corresponding flat period of its inverse These heavy-tailed flat periods of represent long waiting times in which the subdiffusive particle gets immobilized in the trap. The sample path of can be found in Figure 1, and the corresponding waiting time in Figure 2. From the figures, we find the subordinator stop at the same time, which leads the waiting time to be fluctuant.

*Step 2. *This step is aimed at simulating the process . Since the parent process is driven by Brown motion, then we employ the Euler scheme to simulate the process , given by
where is the i.i.d. standard normal noise, . The sample path of 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 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.

#### 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.

#### 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.

#### Acknowledgments

We are grateful to Professor Zhongdi Cen for useful comments and suggestions. This work is supported by the National Natural Science Foundation of China No. 11801288, the Natural Science Foundation of Ningbo No. 2017A610130, and the Natural Science Foundation of Zhejiang Province No. LY17A010020.