Discrete Dynamics in Nature and Society
Volume 2015 (2015), Article ID 303857, 9 pages
http://dx.doi.org/10.1155/2015/303857

Research Article

A Discrete Method Based on the CE-SE Formulation for the Fractional Advection-Dispersion Equation

Silvia Jerez and Ivan Dzib

Department of Applied Mathematics, CIMAT, Jalisco s/n, 36240 Guanajuato, GTO, Mexico

Received 16 November 2014; Revised 15 January 2015; Accepted 17 January 2015

Academic Editor: Manuel De la Sen

Copyright © 2015 Silvia Jerez and Ivan Dzib. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

We obtain a numerical algorithm by using the space-time conservation element and solution element (CE-SE) method for the fractional advection-dispersion equation. The fractional derivative is defined by the Riemann-Liouville formula. We prove that the CE-SE approximation is conditionally stable under mild requirements. A numerical simulation is performed for the one-dimensional case by considering a benchmark with a discontinuous initial condition in order to compare the results with the analytical solution.

1. Introduction

Recently, the analytical and numerical framework of fractional differential equations has attracted much attention because of its extended use in the modeling of nonlocal transport phenomena which appear in many fields like magnetic plasma, turbulence, and flows in porous media [1–5]. A nonlocal transport can lead an anomalous diffusion because of large tracer displacements (super-diffusion) or trapping structure vortices (sub-diffusion). The standard local diffusion description is linked to Gaussian processes and diffusion differential models are obtained. However, a nonlocal transport flow presents a non-Gaussian particle motion and in this case fractional diffusion models are derived [6–8]. One of the most considered proposals for the modeling of the mass flow in media with anomalous diffusion is the fractional advection-dispersion equation (FADE); see, for example, [7–10].

Different definitions for the fractional derivative can be used to describe the FADE, and the most common ones are the Caputo, Riemann-Liouville, and Grünwald-Letnikov derivatives [11–13]. Here, we consider the one-dimensional FADE given by

and use the Riemann-Liouville formula to describe the fractional derivatives:

where and represent the space and the time, respectively, is the solute concentration, is the fractional dispersion, is the skewness of the diffusive transport taking values in the interval , and is the fractional order. In our case we consider and by definition of the ceil function we get . We can rewrite the fractional equation (1) in a more compact form as follows:

by defining the fractional operator as

The analytical solution of (3) is only known under particular initial and boundary conditions, so numerical solutions are calculated and studied for most of the problems. In the last years, a lot of works have been focused on the development of numerical algorithms for fractional differential equations. The first efforts are based on extensions of standard methods used for the nonfractional case like the finite difference method (FDM) [14–17], the finite volume method (FVM) [18, 19], and the finite element method (FEM) [20, 21]. More recently other methods have been developed, for example, the spectral method [22], the collocation method [23], or a combination of both [24]. Generally, the FDM is constructed using the shifted Grünwald formula for the discretization of the fractional derivative, the FVM under a conservative formulation is based on the discretization of the integral using the Riemann-Liouville definition or the Caputo definition of the fractional derivative, and the FEM and collocation methods are variational techniques based on the discretization of the integral on some interpolating polynomial that approximates the solution. Therefore, all these techniques generate algorithms with dense coefficient matrices for which it is complicated to analyze their stability and a high computational cost to solve them is required.

As an alternative to avoid the discretization of the integral associated with the fractional derivative, in this paper we propose a numerical algorithm based on the space-time conservation element and solution element (CE-SE) method developed by Chang and To [25]. For this method, the solution is approximated by a first-order Taylor expansion which must fulfill the integral form of the equation in the conservation elements and the differential form in the solution elements. Therefore, the space-time domain is discretized two times by: a rhomboid mesh and a rectangular mesh. The theoretical framework of the CE-SE method for the nonfractional case is well established and it has been applied to a large number of real problems with successful results [26–33]. An important property of the fractional differential operators is the nonlocality which is achieved with the discretization of the integral associated with the fractional derivative. In our case we seek to keep the nonlocality by providing information of the numerical solution at nearby points in both meshes and using nonlocal approximations of the standard derivatives. It is important to remark that we consider this work as a first step to open a venue for a different way to approximate the fractional differential equations like the FADE equation.

The paper is organized as follows. In Section 2, we give the discretization of the space-time domain for the CE-SE method, its extension to the fractional advection-dispersion equation (4), and the von Neumann stability analysis of the proposed algorithm. In Section 3, we present numerical experiments of the fractional CE-SE method for the FADE with an initial condition type shock and it is compared with the analytical solution of the problem. Finally, this work is completed with a section of conclusions and an appendix with the calculation of the fractional Riemann-Liouville of a linear polynomial of two variables.

2. CE-SE Method for the 1D FADE

2.1. Discretization of the Space-Time Domain

An important feature of the CE-SE method is how the space-time domain is discretized. Let us start with the usual discretization by taking a set of points such that for and for each we take such that with being the left end of the spatial interval and the spatial and time step sizes are defined by and with being a constant of proportionality. Next, the space-time domain is partitioned by two types of mesh: nonoverlapping opened rhombi centered in the solution node which are named solution elements (SEs) and nonoverlapping closed rectangles which are called conservation element (CEs); see Figure 1.

Figure 1: Mesh discretization for the CE-SE method.

By definition, is the interior of the space-time region bounded by the opened rhombus in Figure 2. In each , the flow variables are assumed to be continuous with possible discontinuities at the interface between the two rhombuses. The second mesh is given by nonoverlapping rectangular regions, , which fill the computational domain, see Figure 2. In each , a local space-time flux balance is enforced so a global flux conservation relation is satisfied over the entire space-time domain. Note that each consists of three different s and if we compute the solution in a integer time step then the spatial index is running in the half-integer points and vice versa, see Figure 1.

(a)

(b)

Figure 2: The graphic to the left shows a conservation element and the graphic to the right shows the associated solution element .

2.2. The Fractional CE-SE Method

Now, we construct an explicit algorithm based on the CE-SE method [25, 34] for the one-dimension FADE. First, the solution of (3) is approximated in each by a first-order Taylor expansion of two variables which is denoted by and is given by

Denoting by , , and for all and integers, we need to know these coefficients to obtain the approximation (5). Assuming that they are known in the time step , we are going to calculate the coefficients and imposing that the numerical approximation (5) verifies the integral form of (3) in each . We split the rectangle in two subrectangles and solve the integral form in each subrectangle; see Figure 3. Consider

for . Applying the Green theorem to the integral of the left hand side of (6) we have that

where is the parameterization of given by:

(a)

(b)

Figure 3: In both graphics (a) and (b) different subdivisions of the conservation element, , are shown. In (a) the two regions CEs that the is divided into are shown whereas in (b) the four SEs regions are shown.

Substituting the values of and for each and replacing and its partial derivatives by the coefficients ’s, ’s, and ’s in (7), we obtain

Now, we rewrite the integral of the right hand side of (6) as a function of the coefficients ’s and ’s. For that, it is necessary to calculate the Riemann-Liouville fractional derivative of the function defined in (5). This calculation is carry out in the Appendix. Let be a spatial interval, and for any we write and ( the right-end index of ) and substituting them in (A.5), we get

where

Since the approximation defined in (5) must fulfill the FADE (1) in each , then the coefficient can be obtained by

for . Denoting by the Courant number and by the anomalous diffusion number, we can rewrite the previous formula for the coefficient as

Using (15) for and substituting (9) and (11) into (6) for the case of and for substituting (10) and (12) into (6), we have the following discrete system:

The system (16) can be written in matrix form as follows:

where , , and are the coefficient matrices that are given by

The CE-SE method is a family of schemes [35], and the developed algorithm (16) is called the a-scheme which is the core algorithm of the CE-SE family. Next, we are going to prove the convergence of this fractional a-scheme verifying its consistency with the linear FADE (1) and obtaining its stability conditions. For the stability analysis, we use the von Neumann technique [36] which is a local analysis. Thus, the boundary effects are not considered which is a reasonable assumption since in most cases the numerical instability appears in the interior nodes of the discretization of the domain.

2.3. Stability and Convergence Analysis

We start the von Neumann stability assuming that the solution and its partial derivative with respect to can be expressed as a sum of eigenmodes at each node of the mesh:

where is the spatial wave number and is a complex number. Now, we substitute (19) in the discrete system (17) and obtain

and multiplying both sides by yields

with

Taking the square of Euclidean norm in both sides of (21), we have

where is the amplification factor of each mode in the expansions (19). By calculating we have

where , , , and . The stability of the method is achieved if the magnitude of the amplification factor is less than 1. By definition, the values of the positive coefficients and are close; thus the following condition is fulfilled

By definition of the coefficient , we get ; thus if we take small enough.

Now, we analyze the local truncation error to obtain the consistency of the fractional CE-SE method (17) with the FADE (1). Given that we assume the flow variables , , and to be continuous in each and considering that the numerical solution is a first-order Taylor approximation (5), we conclude that the fractional a-scheme (17) is second-order accurate in space and time. Given the above results, we have the following theorem using the Lax Equivalence theorem [36] which states that For a uniformly solvable linear finite difference scheme which is consistent with a well-posed linear evolution problem, the stability is a necessary and sufficient condition for its convergence.

Theorem 1. The fractional CE-SE scheme (17) is convergent for a small enough.

It is well known that the a-scheme of the CE-SE family is a nondissipative algorithm. We are interested in flow problems with discontinuities, so we need to introduce some numerical dissipation in the approximation. In the next section, we modify the a-scheme using another way to approximate the derivative , and the modified algorithm is known as a--scheme.

2.4. An Approximation of the Derivative for Discontinuities

In order to introduce information about the behavior of the solution in the calculation of the partial derivative of with respect to , we propose another approximation for the coefficient . By adding both equations of the system (16), we can rewrite the fractional CE-SE algorithm as follows:

As a first approximation, we consider the coefficient given in [25, 28]:

where the ratio , with a positive constant, can be considered as a sloper limiter coefficient which measures the smoothness of the function by using the forward difference and the backward difference:

with . As the points are not solution nodes, we can obtain as follows:

The estimation (27) for the partial spatial derivative has been used in nonfractional differential problems with nonsmooth solutions obtaining successful results [25, 26, 28]. In the boundary nodes, the coefficient is calculated from the forward difference and the coefficient is obtained with the backward difference .

Remark 2. Notice that we have used a weighted finite difference approximation to approximate the coefficient in the step , but other approximations can be used. We think that this can be a way to introduce nonlocal information in the calculation of the numerical solution.

3. Numerical Results

For the fractional case it is important to remark that there are few benchmark problems with known analytical solution. So in this section, we consider (3) with a jump condition of the shock type:

with a positive constant function and analyze the performance of the fractional CE-SE method (26)-(27) by comparing with the exact solution, given by Sousa [37]. So, taking and the boundary conditions and , the exact solution is

where and the cumulative probability function is defined by

with

For all and , the cumulative probability function is obtained using the identity . Next, we show the numerical results of the fractional a--scheme for the FADE (1) with initial condition (30) using the same test problem given in [37].

Example 3. We consider the solution of (1) with (30) in the interval for , , and and the boudary conditions are and . In Figures 4 and 5, we show the profile of the analytical solution (31) and the CE-SE approximation (26)-(27) at for several values of the parameter . Moreover, we compare the global error of the fractional a- scheme versus some finite difference methods that appear in [37] for a fixed time using the norm given by

where is the difference between the exact and the numerical solutions; see Table 1. Notice that for the three schemes as is closer to one the numerical error is larger than in the cases as is closer to two, since for smaller the solution loses differentiability. This fact is visible at the extremes of the shock in Figure 4, where we display the CE-SE behavior for and 1.4. Although the three methods present similar numerical errors, the CE-SE method has the advantage of holding a low computational cost compared to the other two methods. Finally in Figure 5, we show the accuracy of the fractional CE-SE method as the step size is diminished.

Table 1: Global -Error in (34) with at time .

Schemes

CESE
6.
8.
1.
1.

Upwind
7.
1.
1.
2.

Lax-Wendroff
6.
7.
1.
1.

(a)

(b)

Figure 4: Numerical solution (red line) of methods (26)-(27) taking and the exact solution (blue line) of the FADE. The graphic to the left shows the results for and the graphic to the right shows the results for , for both cases .

(a)

(b)

Figure 5: Numerical solution (red line) of the method (26)-(27) and the exact solution (blue line) for the FADE with . The graphic to the left shows the results for and the graphic to the right shows the results for .

4. Conclusions

In this work we have constructed a discrete method under the CE-SE formulation for the fractional advection-diffusion equation described by the Riemann-Liouville derivative. The new fractional method avoids the fact that the calculation at the new point depends on the function values at points of further regions at that point. We have also proved the convergence of the fractional CE-SE method and show its performance for a benchmark problem with a shock initial condition. Future work will focus on the analysis of different proposals to approximate the partial spatial derivative of the solution to include the nonlocal information of the flow. We also consider extending the fractional CE-SE method to two dimensions in order to cover a wider range of flow problems important.

Appendix

Fractional Derivative of a Two-Variable Linear Polynomial

Let be a first-degree polynomial function of two variables as follows:

where and are the independent variables with as the indices of the -discretization and -discretization, respectively, and are constant coefficients. Using the linearity property, we rewrite the Riemann-Liouville fractional derivative (4) of the polynomial function (A.1) as follows:

for . Next, we calculate the fractional derivatives that appear in (A.2) using the formulas (2) and obtain

Substituting (A.3) in (A.2), we have

Later on algebraic manipulations, the Riemann-Liouville fractional derivativative of the function defined in (A.1) is given by

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

P. J. Torvik and R. L. Bagley, “On the appearance of the fractional derivative in the behavior of real materials,” Journal of Applied Mechanics, vol. 51, no. 2, pp. 294–298, 1984.

T. R. Ellsworth, P. J. Shouse, T. H. Skaggs, J. A. Jobes, and J. Fargerlund, “Solute transport in unsaturated soil: experimental design, parameter estimation, and model discrimination,” Soil Science Society of America Journal, vol. 60, no. 2, pp. 397–407, 1996.

Y. Pachepsky, D. Benson, and W. Rawls, “Simulating scale-dependent solute transport in soils with the fractional advective-dispersive equation,” Soil Science Society of America Journal, vol. 64, no. 4, pp. 1234–1243, 2000.

D. del-Castillo-Negrete, B. A. Carreras, and V. E. Lynch, “Fractional diffusion in plasma turbulence,” Physics of Plasmas, vol. 11, no. 8, pp. 3854–3864, 2004.

B. M. Schulz and M. Schulz, “Numerical investigations of anomalous diffusion effects in glasses,” Journal of Non-Crystalline Solids, vol. 352, no. 42–49, pp. 4884–4887, 2006.

A. A. Greenenko, A. V. Chechkin, and N. F. Shul’ga, “Anomalous diffusion and Lévy flights in channeling,” Physics Letters, A, vol. 324, no. 1, pp. 82–85, 2004.

D. del-Castillo-Negrete, “Fractional diffusion models of nonlocal transport,” Physics of Plasmas, vol. 13, no. 8, Article ID 082308, 2006.

W. Chen, H. Sun, X. Zhang, and D. Korošak, “Anomalous diffusion modeling by fractal and fractional derivatives,” Computers & Mathematics with Applications, vol. 59, no. 5, pp. 1754–1758, 2010.

D. A. Benson, S. W. Wheatcraft, and M. M. Meerschaert, “The fractional-order governing equation of Lévy motion,” Water Resources Research, vol. 36, no. 6, pp. 1413–1423, 2000.

Y. Zhang, D. A. Benson, M. M. Meerschaert, and E. M. LaBolle, “Space-fractional advection-dispersion equations with variable parameters: diverse formulas, numerical solutions, and application to the Macrodispersion Experiment site data,” Water Resources Research, vol. 43, no. 5, Article ID W05439, 2007.

M. Weilbeer, Efficient Numerical Methods for Fractional Differential Equations and Their Analytical Background, Papierflieger, 2005.

D. Baleanu, K. Diethelm, E. Scalas, and J. J. Trujillo, Fractional Calculus Models Andnumerical Methods, vol. 3, World Scientific, Singapore, 2009.

C. Li, D. Qian, and Y. Chen, “On Riemann-Liouville and Caputo derivatives,” Discrete Dynamics in Nature and Society, vol. 2011, Article ID 562494, 15 pages, 2011.

M. M. Meerschaert and C. Tadjeran, “Finite difference approximations for fractional advection-dispersion flow equations,” Journal of Computational and Applied Mathematics, vol. 172, no. 1, pp. 65–77, 2004.

C. Li and F. Zeng, “Finite difference methods for fractional differential equations,” International Journal of Bifurcation and Chaos in Applied Sciences and Engineering, vol. 22, no. 4, Article ID 1230014, 2012.

F. Liu, P. Zhuang, V. Anh, I. Turner, and K. Burrage, “Stability and convergence of the difference methods for the space-time fractional advection-diffusion equation,” Applied Mathematics and Computation, vol. 191, no. 1, pp. 12–20, 2007.

A. Ashyralyev and F. Dal, “Finite difference and iteration methods for fractional hyperbolic partial differential equations with the Neumann condition,” Discrete Dynamics in Nature and Society, vol. 2012, Article ID 434976, 15 pages, 2012.

X. Zhang, J. W. Crawford, L. K. Deeks, M. I. Stutter, A. G. Bengough, and I. M. Young, “A mass balance based numerical method for the fractional advection-dispersion equation: theory and application,” Water Resources Research, vol. 41, no. 7, Article ID W07029, 2005.

H. Hejazi, T. Moroney, and F. Liu, “Stability and convergence of a finite volume method for the space fractional advection-dispersion equation,” Journal of Computational and Applied Mathematics, vol. 255, pp. 684–697, 2014.

Q. Huang, G. Huang, and H. Zhan, “A finite element solution for the fractional advection-dispersion equation,” Advances in Water Resources, vol. 31, no. 12, pp. 1578–1589, 2008.

Y. Zheng, C. Li, and Z. Zhao, “A note on the finite element method for the space-fractional advection diffusion equation,” Computers & Mathematics with Applications, vol. 59, no. 5, pp. 1718–1726, 2010.

X. Li and C. Xu, “A space-time spectral method for the time fractional diffusion equation,” SIAM Journal on Numerical Analysis, vol. 47, no. 3, pp. 2108–2131, 2009.

A. Saadatmandi, M. Dehghan, and M.-R. Azizi, “The Sinc-Legendre collocation method for a class of fractional convection-diffusion equations with variable coefficients,” Communications in Nonlinear Science and Numerical Simulation, vol. 17, no. 11, pp. 4125–4136, 2012.

M. Zayernouri and G. E. Karniadakis, “Fractional spectral collocation method,” SIAM Journal on Scientific Computing, vol. 36, no. 1, pp. A40–A62, 2014.

S. C. Chang and W. M. To, “A new numerical framework for solving conservation laws—the method of space-time conservation element and solution element,” Tech. Rep. 104495, NASA, Cleveland, Ohio, USA, 1991.

X. Y. Wang, S. C. Chang, and P. C. E. Jorgenson, “Prediction of sound waves propagatingthrough a nozzle withoutlwilh a shock wave using the space-time CE/SE method,” NASA 209937, 2000.

Z.-C. Zhang, S. T. J. Yu, and S.-C. Chang, “A space-time conservation element and solution element method for solving the two- and three-dimensional unsteady Euler equations using quadrilateral and hexahedral meshes,” Journal of Computational Physics, vol. 175, no. 1, pp. 168–199, 2002.

S. Jerez, J. V. Romero, M. D. Rosell{\'o}, and F. J. Arnau, “A semi-implicit space-time CESE method to improve mass conservation through tapered ducts in internal combustion engines,” Mathematical and Computer Modelling, vol. 40, no. 9-10, pp. 941–951, 2004.

J. Wang, K. Liu, and D. Zhang, “An improved CE/SE scheme for multi-material elastic-plastic flows and its applications,” Computers & Fluids, vol. 38, no. 3, pp. 544–551, 2009.

S. Qamar and S. Mudasser, “On the application of a variant CE/SE method for solving two-dimensional ideal MHD equations,” Applied Numerical Mathematics, vol. 60, no. 6, pp. 587–606, 2010.

S. C. Chang, A New Approach for Constructing Highly Stable High Order CESE Schemes, National Aeronautics and Space Administration, Glenn Research Center, 2010.

S. Sadighi, A. Ahmad, and M. Shirvani, “Dynamic simulation of a pilot scale vacuum gas oilhydrocracking unit by the spacetime CE/SE method,” Chemical Engineering and Technology, vol. 35, no. 5, pp. 919–928, 2012.

S. Qamar and M. Yousaf, “The space-time CESE method for solving special relativistic hydrodynamic equations,” Journal of Computational Physics, vol. 231, no. 10, pp. 3928–3945, 2012.

D. Yang, S. Yu, and J. Zhao, “Convergence and error bound analysis for the space-time CESE method,” Numerical Methods for Partial Differential Equations. An International Journal, vol. 17, no. 1, pp. 64–78, 2001.

S.-C. Chang, “Courant number insensitive CE/SE schemes,” in Proceedings of the 38th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, Indianapolis, Ind, USA, July 2002.

L. N. Trefethen, Finite Difference and Spectral Methods for Ordinary and Partial Differential Equations, Cornell University, Ithaca, NY, USA, 1996.

E. Sousa, “Finite difference approximations for a fractional advection diffusion problem,” Journal of Computational Physics, vol. 228, no. 11, pp. 4038–4054, 2009.