A Finite Difference Method on Quasi-uniform Mesh for Time-Fractional Advection-Diffusion Equations with Source Term

The present paper deals with the numerical solution of time-fractional advection-diffusion equations involving the Caputo derivative with source term by means of an unconditionally stable implicit finite difference method on quasi-uniform grids. We use a special quasi-uniform mesh in order to improve the numerical accuracy of the classical discrete fractional formula for the Caputo derivative. The stability and the convergence of the method are discussed. The error estimates established for a quasi-uniform grid and a uniform one are reported to support the theoretical results. Numerical experiments are carried out to demonstrate the effectiveness of the method.


Introduction
The fractional partial differential equations (FPDEs) have become increasingly popular in recent years. The interest in these equations comes from their mathematical structure and from their applications. FPDEs have applied in various areas of engineering, science, finance, applied mathematics, bio-engineering and so on. Though several ways to solve them theoretically have been proposed [34,21,35], including Green function method and Laplace and Fourier transform method [37,31], the Adomian decomposition [8,7], the Homotopy Perturbation Methods [17,32], generally, numerical solution techniques are preferred when dealing with fractional models since the analytical solutions are available for a few simple cases. In recent years, efficient numerical methods have been developed to solve fractional differential equation, including finite difference methods [29,39,28,36,6], the finite volume method [16], the finite element method [11,43,42] and the spectral method [25,9].
As also for the non-fractional differential equations, finite difference methods are one of the most important classes of numerical methods for solving FPDEs. Zhuang and Liu [46] obtained an implicit difference approximation to solve the time-fractional diffusion equations. Lin and Xu [25,26] proposed the numerical solution by finite/difference approximations for a time-fractional diffusion equation. Liu et al. [27] developed an explicit difference method and an implicit difference method for solving a space-time fractional advection dispersion equation on a finite domain. In [24,3,4], high order numerical difference schemes were constructed in order to solve the Caputo-type advection-diffusion equations. Zhang et al. [44] obtained a finite difference method for FPDEs involving the Caputo derivative on a non-uniform mesh. Recently Jannelli et al. [19,20] determined exact and numerical solutions for the time-fractional advection-diffusion differential equations involving Riemann-Liouville derivative with a non linear source term by means the Lie symmetries. They transformed the fractional partial differential equation into a fractional ordinary differential equations, which is then solved using the implicit backward differentiation formulas.
The main goal of this paper is to construct an unconditionally stable implicit finite difference method for solving the time-fractional advection-diffusion equations (FADEs) with a nonhomogeneous source term involving the Caputo fractional derivative on quasi-uniform grids. We choose to use a special quasi-uniform mesh in order to improve the numerical accuracy of the classical discrete fractional formula for the Caputo derivative, since the fractional derivatives are integrals with weakly singular kernel and the discretization on the uniform mesh may lead to poor accuracy. The consistency, stability and convergence of the proposed difference method are investigated. Three numerical examples are given to show the reliability and efficiency of the derived difference method.

The mathematical model
We consider the following linear time-fractional advection-diffusion equation with the initial and boundary conditions given by where u is the field variable that can represent, for example, the solute concentration, K 1 and K 2 are the constant fluid velocity and the dispersion coefficient, respectively. The time fractional derivative ∂ α ∂t α u(x,t) is the α order Caputo fractional derivative defined by The function f (x,t) can be used to represent sources and sinks. φ (x), ϕ(t) and ψ(t) are known smooth functions. We take K 1 = 0 and K 2 > 0 and we assume the problem (1) has a unique and sufficiently smooth solution under the above initial and boundary conditions. The fractional equation (1) has been treated by a number of authors. It is presented as a useful approach for the description of transport dynamics in complex systems which are governed by anomalous diffusion and non-exponential relaxation patterns [30]. The FADE is also used in groundwater hydrology research to model the transport of passive tracers carried by fluid flow in a porous medium [1].
It is to note that, when α = 1, the model (1) reduces to the classical advectiondiffusion equation with source term used in order to describe several phenomena of relevant interest in many fields of applied sciences. In [18] and [12], a fractional step approach with variable time step is used in order to solve numerically mathematical models that describe evolution problems on a three-dimensional domains

Discretization in time on a quasi-uniform mesh
The main goal of this work is to construct an unconditionally stable implicit finite difference method defined on a quasi-uniform mesh. In general, the existence of a weakly singular kernel (t − s) −α , 0 < α < 1, in fractional derivatives makes it more difficult to get a higher-order scheme. Particularly when the solutions are not suitably smooth, numerical methods on uniform meshes seem to have a poor convergent rate. For these reasons, numerical schemes on non-uniform meshes have been developed in the last years. In this Section, first, we construct the quasi-uniform mesh and then, in order to approximate the Caputo derivative on quasi-uniform grid, we define a suitable discrete fractional derivative formula.
We divide the interval [0, T ] into N subintervals [t n−1 ,t n ], for n = 1, · · · , N and with 0 = t 0 < t 1 < · · ·t N = T . We denote the time step ∆t n as and let A sequence of mesh is called quasi-uniform if there exists a finite constant β such that ∆t max /∆t min ≤ β .
In this case, it holds that ∆t max ≤ β T N −1 . When β = 1, it hods that ∆t n = T N −1 for all n = 1, · · · , N and the mesh is reduced to a uniform mesh.
A sequence of meshes is not quasi-uniform if The non-uniform mesh and quasi-uniform mesh methods have been used for solving differential equations by several authors. In [13,14], quasi-uniform meshes were used for solving a class of boundary values problems on infinite domains. In this work, we are interested on the quasi-uniform mesh [44] defined as follows . It is to note that the time steps {∆t n } N n=1 are a monotonically decreasing sequence with ∆t 1 = O(N −1 ) and ∆t N = O(N −2 ). Figure (1) shows a sample of the quasi-uniform grid (3) obtained for N = 10.
To motivate the choice of the quasi-uniform mesh (3), we discretize the Caputo derivative (2) of a function v(t) by means of the following classical approximate formula PSfrag replacements t Figure 1: Quasi-uniform grid for N = 10 that is the well-known so-called L1 formula defined in [34] where R n is the local truncation error. The L1 formula has been used for solving the fractional differential equations with Caputo derivatives (see [46,15]). Moreover, using the relationship between Caputo derivative and Riemann-Liouville fractional derivative, the L1 formula was also applied to time fractional diffusion equation with Riemann-Liouville fractional derivative (see [22,45]). High-order approximations such as compact difference scheme [10,15,45,5] and spectral method [25,26,23] were applied to improve the spatial accuracy of fractional diffusion equations. It is important to note that it is rather difficult to get a high-order time approximation due to the singularity of fractional derivatives. For any temporal mesh, for 0 < α < 1 and v(t) ∈ C 2 [0, T ], it can be verified that [44] For the uniform mesh, i.e., ∆t n = ∆t for all n = 1, 2, · · · , N, then R n = O(∆t 2−α ) (see [38,25]). From the truncation error estimate of the L1 formula, it is clear that the accuracy is dependent on the fractional order α. This is justified since a weakly singular kernel (t − s) −α is contained in the integral. In order to improve the accuracy of the L1 numerical approximation of derivative of fractional order, it is possible to consider non-uniform meshes. Numerical methods developed with the non-uniform meshes have been developed for solving integro-differential equations with weakly singular. Mustapha [33], Yuste and Quintana-Murillo [40,41] proposed an implicit finite-difference time-stepping method for discretizing of the time diffusion equation.
Zhang et al. [44] obtained a numerical integration formula for any α ∈ (0, 1) by employing the special quasi-uniform grid (3) and obtained the following results: for the quasi-uniform mesh (3) Then, we obtain See [44] for more details. We use these estimates to study the consistency, stability and convergence of the method presented in the following sections.

An implicit finite difference method
For the derivation of the implicit difference method, first, we construct a computational uniform grid in the x direction, that is the spatial size of the mesh ∆x = x j − x j−1 is constant, for 1 ≤ j ≤ J, and quasi-uniform in the time direction. We define the mesh points (x j ,t n ) with x j = a + j∆x, j = 0, · · · , J and t n = t n−1 + ∆t n , for n = 1, · · · , N, with ∆t n defined by (3). J and N are positive integers. We denote by U n j the numerical approximation provided by the difference method of the exact solution u(x j ,t n ) at the mesh points (x j ,t n ), for j = 0, · · · , J and n = 0, · · · , N.
As usual, we discretize the first ∂ /∂ x and second order ∂ 2 /∂ x 2 spatial derivatives by means of the second order three-point central difference formula so that According to the discretization for the Caputo derivative (4), we approximate the time fractional derivative in (1) as follows Replacing u(x j ,t n ) with its numerical approximation U n j and neglecting the local truncation errors, the time-fractional advection-diffusion equation (1) is discretized as follows is the source term and The initial and boundary conditions can be rewritten as For any temporal meshes on [0, T ] and for any n = 1, 2, · · · , N, it holds that [46] T n,k > 0 T n,k > T n,k−1 , Taking into account that T n,n = (∆t n ) −α , we can write then, we obtain the following implicit finite difference scheme where we set The equation (11) can be written in vectorial form as where K is a tridiagonal matrix and where we denote with L the following difference operator Here and in the following we assume the convention that the summation is equal to zero if the lower bound is larger that the upper bound. The obtained method (12) is implicit, in order to compute the numerical solution U n a system with the tridiagonal coefficients matrix K has to be solved. It is interesting to note that the operator L is a kind operator with memory, due to the non-local character of the fractional derivative, this means that the effect on U at time t n , U n , depends on all the previous values, U 0 ,U 1 , · · · ,U n−1 , evaluated at all the previous time t 0 ,t 1 , · · · ,t n−1 .
The main difference compared to the non-fractional case is that, in order to evaluate L , the numerical solutions for all the n previous time values t 0 ,t 1 , · · · ,t n−1 are required, while for non-fractional equations, only the solution to the previous value t n−1 is used. The computational cost to compute the solution at the time t n from the solution at the time t n−1 grows as n, that is, grows as the number of terms in the summation that compares in the second term of the (13). This implies that the computational cost to go from t 0 to t n grows as n 2 .

Consistency, Stability and Convergence
In this Section, we discuss the consistency, the stability and the convergence of the implicit finite difference scheme.
Stability. For the stability analysis of the implicit finite difference scheme we rewrite the equations (11) as LetŪ n j be another approximate solution of the difference scheme (11), and let ρ n j = U n j −Ū n j , be the corresponding round-off error. We let ρ n = (ρ n 0 , ρ n 1 , · · · , ρ n J ) T , and we consider the infinity norm The round-off error satisfies the following round-off equations In order to check whether the finite difference scheme is stable, we study how the size of the round-off error ρ n evolves in time.
Convergence. Let u(x j ,t n ) be the exact solution of Eqs.
(1) at mesh point (x j ,t n ) for j = 0, 1, · · · , J and n = 0, 1, · · · , N. Denoting ε n j = u(x j ,t n ) − U n j , we get the error equations with ε 0 j = 0, for j = 0, 1, · · · , J and ε n J = 0, for n = 1, 2, · · · , N. R n j is the local truncation error. We introduce the following norm then, we obtain where R max = max n,i |R n i |. For n = 1, we have Then To obtain the error estimates of the numerical solutions, we need the uniform error bounds on all time levels. Then, applying the first of the estimates (5), we can consider where C R is a positive constant that is dependent on T , α and the exact solution u(x,t), but independent of N and ∆x. Then, we prove that the solution of the difference method (11), with initial and boundary conditions given by (10), is convergent.

Numerical experiments
In this Section, we report some numerical examples of the FPDEs to demonstrate the accuracy and efficiency of the numerical method. The first examples are chosen in such a way the exact solution of FPDE can be evaluated analytically to show that the approach proposed in this paper properly works. The exact solutions allow to verify the accuracy and the order of convergence of the numerical solution. We use the proposed method on the quasi-uniform grid and on a uniform grid and compare the numerical results, observing that the finite difference scheme generates more accurate numerical solutions on the quasi-uniform grid than on the uniform one. In the third test, we solve a FADE of physical interest with the source term chosen as a linear function of the field variable. in order to show that the method is applicable to a wide class of FADEs. By this last test, we illustrate how the changes in the solution behavior arise when the fractional order is varied.
Example 1. We consider the following FADE where the source term is given by The exact solution is found to be u(x,t) = e x t β .
In this example, we take K 1 = K 2 = 1 with α = 0.5 and β = 5. Figure 2 shows the comparison between the numerical solution U n j and the exact solution u(x j ,t n ) at different times computed with N = J = 20. From the Figure 2, it can be seen that the numerical solution U n j is in good agreement with the exact solution u(x j ,t n ). The exact solution is reported with the solid line.
In order to investigate the temporal error and the convergence order of the numerical difference method, we define the maximum error between the exact solution u(x j ,t N ) and the numerical solution U N j at the final time t N e ∞ (N, In order to show the efficiency of the method, we defineŪ n j the numerical solution obtained by the implicit difference method on a uniform grid defined by ∆t n = T /N, for n = 1, · · · , N. We use notations similar to the (22) and (23) to define the maximum errorē ∞ (N, J) between the exact solution u(x j ,t N ) and the numerical solutionŪ N j at the final time t N and the convergence OrderŪ . In this test, we fix J = 100, a value large enough such that the spatial error is negligible as compared with the temporal error. Table 1 shows the values of e ∞ andē ∞ and the corresponding numerical convergence orders for α = 0.1, 0.5 and 0.9. It can be seen that the method is stable and convergent for solving the problem (21) on both the computational grids. By using the L1 formula on quasi-uniform mesh, we improve the order of accuracy in time of the proposed method. In fact, we observe that the solutions are more accurate on the quasi-uniform mesh and the convergence order on the quasi-uniform mesh is greater than the convergence order obtained with the uniform mesh. The numerical results agree well with the theoretical results.  Table 1: e ∞ ,ē ∞ and convergence orders, related to the numerical solutions U n j and U n j respectively, for different values of N and α, with J = 100.
Example 2. We consider the following FADE where the source term is given by and the exact solution is We take K 1 = K 2 = 1 and set α = 0.1. Figure 3 shows the comparison between the numerical solution U n j and the exact solution u(x j ,t n ) at different times computed with N = J = 20. The exact solution is reported with the solid line. The numerical solution U n j is in good agreement with the exact solution u(x j ,t n ). The values of e ∞ andē ∞ and the corresponding numerical convergence orders obtained for α = 0.1, 0.5 and 0.9 are reported in the table 2. The results confirm that the method is stable and convergent for solving the problem (24) on both the computational grids. The numerical solutions computed on quasi-uniform mesh are more accurate and the convergence order on the quasi-uniform mesh is greater than the convergence order obtained with the uniform mesh. The errors e ∞ andē ∞ satisfy the relationships (22). The convergence orders Order and OrderŪ satisfy the relationships (23).
Example 3. In this test, in order to show the efficiency of the method, we solve the following FADE where the source term is chosen as a linear function of the field variable  Table 2: e ∞ ,ē ∞ and convergence orders, related to the numerical solutions U n j and U n j respectively, for different values of N and α, with J = 100.
The model describes the one-dimensional transport problem of a concentration u(x,t) of a chemical or biological species in a flowing medium such as air or water. The species concentration is assumed horizontally and vertically well mixed such that it varies only in the longitudinal or downstream direction. Moreover, a steady and uniform flow field is imposed and the effects of the dispersion are constant in time and space. A reaction where the transformation rate β is proportional to the species concentration is considered; according to the sign of rate, we may have either decay effects or not.
In this example, we take a = 0, b = 5, K 1 = 1, K 2 = 1 and β = 0.2. Figure  4 shows the solution behavior at different times obtained for α = 0.5 and with N = J = 100. The height of the solution profile decreases as the time increases.

Concluding remarks
In this paper, we develop an unconditionally stable finite difference method on quasi-uniform grid for solving time-fractional advection-diffusion equations involving the Caputo fractional derivative. The Caputo time derivative is discretized by means of a direct generalization of the well-known fractional L1 formula (4) to the case of quasi-uniform meshes. The L1 formula takes into account the non-local character of the time-fractional operator and allows to improve the order of accuracy in time of the proposed method by using of the quasi-uniform time discretization obtained by the mesh (3). We prove the stability and convergence of the proposed method. Numerical experiments are carried out to support the theoretical results. The reported numerical experiments point out that the difference method is more accurate on the quasi-uniform grid than on the uniform mesh and the convergence order is greater on the quasi-uniform mesh than on the uniform mesh. Moreover, it is important to note that, in view of its simplicity, the method is applicable to a wide class of FADEs occurring in applied sciences.