Journal of Applied Mathematics

Volume 2013 (2013), Article ID 648595, 8 pages

http://dx.doi.org/10.1155/2013/648595

Research Article

A High-Accuracy MOC/FD Method for Solving Fractional Advection-Diffusion Equations

Lijuan Su and Pei Cheng

School of Mathematical Sciences, Anhui University, Hefei, Anhui 230601, China

Received 15 December 2012; Revised 2 February 2013; Accepted 21 February 2013

Academic Editor: Morteza Rafei

Copyright © 2013 Lijuan Su and Pei Cheng. 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.

Fractional-order diffusion equations are viewed as generalizations of classical diffusion equations, treating super-diffusive flow processes. In this paper, in order to solve the fractional advection-diffusion equation, the fractional characteristic finite difference method is presented, which is based on the method of characteristics (MOC) and fractional finite difference (FD) procedures. The stability, consistency, convergence, and error estimate of the method are obtained. An example is also given to illustrate the applicability of theoretical results.

1. Introduction

The history of fractional calculus is almost as long as integer-order calculus. However, because of lack of application background, fractional calculus was developed very slowly. Up to now, fractional calculus has been applied to almost every field of hydrology [1], physics [2, 3], engineering [4], mathematics [5], and science [6, 7], such as diffusion [8], oscillation, and viscoellastic dynamics. The fractional derivatives in space are often used to model anomalous diffusion because particles spread faster than the classical models predict. The main physical purpose for investigating diffusion equation of fractional order is to describe phenomena of anomalous diffusion in transport progresses through complex or disordered systems including fractal media, and fractional kinetic equations have been proved especially useful in the context of anomalous slow diffusion.

There are several different methods to solve differential equations of fractional order [9], for instance, the method of images, the Fourier transform methods, the method of separation of variables, the Mellin transform methods, and the Laplace transform methods. Wyss [10] considered the time fractional diffusion equation and the solution is given in closed form in terms of the Fox functions. In the last decade, different numerical methods have been developed to solve the fractional diffusion equations [11, 12]. Several explicit finite difference schemes were proposed and analyzed in [13]. Liu et al. [14] used L2 scheme form [5] to discretize the fractional Fokker-Planck equation. Deng [15] proposed finite element method for the space and time fractional Fokker-Planck equation. Meerschaert and Tadjeran [16] considered a Grünwald-Letnikov approximation for the two-sided space-fractional partial differential equations.

Many difficult problems arise in the numerical simulation of fluid flow within porous media in petroleum reservoir simulation and in subsurface contaminant transport and remediation [17]. For example, their mathematical models are basically advection dominated. Fractional derivatives play an important role in modelling particle transport in anomalous diffusion [18]. The space-fractional advection-dominated diffusion equation [19] is

where the function is the drift of the process, is the fractional order of the spatial derivative, physical considerations restrict , the functions may be interpreted as the coefficient of dispersion, and is a source/sink term. Moreover, when (1) is convection dominated diffusion equation, we have .

Equation (1) contains a Riemann-Liouville fractional derivative of order of a function , for , , defined by [9]

where is an integer such that , is the gamma function.

The standard finite difference methods (SFDM) for solving (1) were discussed in [16, 19, 20], and their stability was given in those papers, too. However, these methods often generate numerical solutions with severe nonphysical oscillations. Because of the nature of hyperbolic types [17, 21], the convection-diffusion problems involve moving sharp fronts or boundary layers. While upwind finite difference methods eliminate nonphysical oscillations, they tend to generate numerical solutions with excessive numerical diffusion and smear the moving steep fronts. The efficient grid refinements are required if the standard methods are used in computations. But for high dimensional large-scale problems, it may lead to very large linear systems and need long computation time.

In this paper, we will consider a new numerical difference scheme proposed by the characteristic method, which is named the fractional characteristic finite difference method (CFDM). It is proved that the new method is unconditionally stable, consistent, and convergent; therefore, there is no stability limitation on the size of . Moreover, the CFDM allows large time steps to be used without the loss of accuracy, which is much better than the SFDM and other numerical methods. The reason is that, for problems with significant convection, the solution changes much less rapidly in the characteristic direction than in the direction.

The outline of the paper is as follows. In Section 2, we describe the fractional characteristic finite difference schemes for solving (1). In Section 3, the stability, consistency, and convergence are proved, and the maximum error estimate is derived, too. The numerical experiments are given in Section 4. The conclusions of this paper are given in Section 5.

2. Fractional Characteristic Difference Method

In this section, we consider the combination of characteristic methods with finite difference techniques [22] and seek a new numerical method of (1), which is the convection-dominated problem and reflects the almost hyperbolic nature.

2.1. The Idea of Characteristic Methods

Let and denote the characteristic direction associated with the operator by (see [22]), where

Then, the left part of (1) can be rewritten as

Next, in order to formulate the fractional CFDM of (1), we give some notations. Let be the temporal mesh or time step, , , , is the spatial mesh, the coordinates of the mesh points are , , , with and are positive integers. The values of the solution on these grid points are , and the numerical estimate of the exact value of at the point is denoted by . Finally, define ,  , and .

The characteristic derivative is approximated in the following manner:

where .

Then, by (5) the differential operators at the point are discretized to be

where , is calculated by the piecewise-linear interpolation.

Let , and let symbol be the integral part of . Let and , using linear interpolation to , we get [23]

Generally, the quantity is called the Courant (or CFL) number.

2.2. The Discrete Versions of the Fractional Derivative

First, we develop the right-shifted Grünwald discrete approximation to the fractional derivatives by means of improving the standard form of the Grünwald-Letnikov definition [16]:

where is a positive integer, , and , which can be calculated by

Following (5), we can discretize the Riemann-Liouville operator [12]

where is a positive integer, , and is the integral part of .

Here, according to the definition of , we can draw the following conclusion [9, 12]: when , ; when , .

According to the shifted Grünwald-Letnikov formulae (8) and (10), the definition (2) can be discretized by

with .

Inserting (6), (7), and (11) into (1), neglecting the truncation error, after rearranging the terms, we finally get the fractional characteristic difference scheme

where is associated with the diffusion coefficient, and the initial values are .

3. Stability and Convergence Analysis

In this section, before studying the stability and convergence of the CFDM, we give the following lemma, which was proved previously [16].

Lemma 1. According to the definition of in (9), when , the coefficients have the following properties:

, and , , and ,, .

In the following section, we study the stability and convergence of the fractional CFDM for (1) with the given initial and boundary conditions.

3.1. Stability Analysis

Theorem 2. The fractional characteristic finite difference scheme (12) of (1) with , based on the shifted Grünwald approximation (8) to the space fractional derivative, is unconditionally stable.

Proof. Incorporating the boundary conditions, (12) is rewritten as , where , , and are two matrices.

When and , we have

where , so , too. Thus, holds.

For matrix , we can easily get

Therefore,

For matrix , we can easily see that for all , and

At the same time, according to Lemma 1, we have

for , we can obtain

Therefore, is diagonally dominant by rows. According to the matrix theory [24, 25],

Then, from (15) and (19), we can get

Thus, the characteristic finite difference scheme (12) is unconditionally stable.

In addition, taking into account (6) and (11), we have

so this method is consistent with a local truncation error which is .

3.2. Convergence Analysis

Let ; be the exact solution of (1) at mesh point . Define , ; and . Using , substitution into (12) leads to

where , and , ; .

From (6), (11), and (21), we have the following.

Lemma 3. There is a positive constant , such that

where ; .

Let , because is finite, we can obtain the following result.

Theorem 4. Let be the numerical solution computed by the CFDM (12). Then there is a positive constant , such that

Proof. From (22) and , we can get

According to (19) and Lemma 3, using (25), we finally get

Because , from (26), we obtain the error estimate . Therefore, we get the result , where is a positive constant.

Finally, according to Theorem 4, we draw the following conclusion: when and , holds; that is, the characteristic finite difference scheme (12) is convergent. In summary, the fractional CFDM for (1) is unconditionally stable, consistent, and convergent.

4. Numerical Simulations

In this section, we give some numerical results for a special case of (1) to confirm our theoretical analysis discussed above. In addition, we also compare our new CFDM with the SFDM and point out the advantages of the new method.

For (1), if we assume that and are two constant functions, the function , and the initial condition is the dirac delta function. Then, under the zero Dirichlet boundary conditions, the exact solution of (1) could be found by the Fourier transform methods [9]:

In the numerical experiments, the data are chosen as follows: , , , and . For simplicity, we assume that the initial time , and the initial values are calculated from (27). The numerical solutions are obtained from the CFDM scheme (12) discussed above. The errors of numerical solutions are measured in and norms which are defined as follows:

where and are the numerical solution and analytical solution at , respectively. is the space step, and is the time step. The ratios of convergence in time are calculated by

where and are two time steps, and very small spatial step size is taken in computation. Similarly, the ratios of convergence in space can be calculated, too.

4.1. The Properties of the CFDM

In this subsection, we use the CFDM (12) to compute the numerical solutions to (1), and list the computational results in the figures and tables. From Figure 1, we can see that the numerical solutions of the CFDM coincide with the analytic solutions, with the diffusion coefficient , the space step , and the time step at , for and 1.6. It is clearly shown that the CFDM is stable for (1), which allows the large time step to be used and gains a very good numerical approximation to the exact solution.

Figure 1: Exact solutions and numerical solutions obtained by the CFDM with , , and at . The solid lines correspond to the exact solutions. The starred and dotted lines correspond to numerical solutions for and 1.6, respectively.

Next, we examine the ratio of convergence in space for the CFDM. Table 1 shows and errors and ratios in space with the diffusion coefficient and the small time step at time , for two different fractional orders and 1.6, respectively. The second and fourth columns show the and errors, respectively. Column three shows the ratio of the CFDM in norm, and column five shows the ratio of the CFDM in norm. It is clearly shown that the CFDM method is of first-order accuracy in space in both and norms.

Table 1: The errors and ratios in space calculated by CFDM (12) with and at time , for two different fractional orders and 1.6, respectively.

-CFDM
Ratio
-CFDM
Ratio

1/25

—

—

1/30

0.9325

0.9581

1/35

0.9788

1.0001

1/40

1.0207

1.0401

1/45

1.0581

1.0771

1/60

—

—

1/65

1.0612

1.0961

1/70

1.1002

1.1342

1/75

1.1408

1.1720

1/80

1.1813

1.2108

Meanwhile, Table 2 shows and errors and ratios in time with and at time , for two different fractional orders and 1.6, respectively. The table shows that the CFDM method is of first-order accuracy in time in both and norms, too. Furthermore, from the table, we can see that even with the very large time steps size of to 1/30, which corresponds to the Courant numbers of 14.14 to 84.85, the CFDM can still generate the accurate numerical solutions without noticeable numerical errors. In conclusion, according to Tables 1 and 2, the accuracy of the CFDM method is of order , which confirms the result of Theorem 4.

Table 2: The errors and ratios in time calculated by CFDM (12) with and at time , for two different fractional orders and 1.6, respectively.

-CFDM
Ratio
-CFDM
Ratio

1/5

—

—

1/10

1.0558

1.0052

1/15

1.0536

1.0175

1/20

1.0607

1.0268

1/30

1.0749

1.0383

1/10

—

—

1/15

0.8238

0.8661

1/18

0.7322

0.7997

1/20

0.8987

0.8604

1/25

0.7771

0.7327

4.2. Comparison with the SFDM

To gain a better understanding of the CFDM, we use the well-developed and well-received SFDM to simulate the same model problem for comparison, which usually contains the implicit and explicit finite difference schemes. We only consider the case of in this subsection.

For the implicit SFDM, it is also unconditionally stable. Table 3 and Figure 2 both show that the CFDM is much better than the implicit SFDM, if they have the same assumptions. Firstly, Table 3 shows the comparison of errors calculated by the CFDM (12) and the implicit SFDM for two different diffusion coefficients and , with some different space steps and time steps for at time . From the table, we can easily see that the errors generated by the CFDM is obviously much smaller than the errors generated by the implicit SFDM, which has one- even two-order higher accuracy than the latter.

Table 3: Comparison of errors calculated by the CFDM (12) and the implicit SFDM with , different space steps and time steps at time , for two different diffusion coefficients and .

-CFDM
-im-SFDM

0.1

1/60
1/30

1/50
1/40

1/80
1/50

0.02

1/80
1/20

1/90
1/30

1/100
1/45

Figure 2: The errors between exact solutions and numerical solutions calculated by the CFDM (12) and the implicit SFDM at different time levels with , , and for .

Moreover, Figure 2 shows that the errors between exact solutions and numerical solutions calculated by the CFDM and the SFDM at different time levels with , the diffusion coefficient , the space step , and the time step . It is easily to see that the CFDM can generate accurate numerical solution even if relatively large time step is used, and this significantly improves efficiency [17, 26]. However, the implicit SFDM does not have this kind of property.

Next, we compare the CFDM and the implicit SFDM numerical solutions with and at time in Figure 3. To keep similar error, for the CFDM, the grid sizes are chosen as and , and it generates the error of . However, for the implicit SFDM, we have to choose very small grid sizes, which are at least and . In this case, the implicit SFDM generates the error of , which is still larger than that of the CFDM. As a result, it is evidently shown that the CFDM can allow large time steps and space steps without the loss of accuracy.

Figure 3: Numerical solutions and exact solution at with . The solid line corresponds to the exact solution and the dotted line corresponds to numerical solution of the CFDM with , , and . The stared line corresponds to numerical solution of the implicit SFDM with and .

In Section 4 discussed above, we only compare the numerical results of the CFDM with the implicit SFDM. In fact, the SFDM for (1) also contains several other explicit schemes which are conditionally stable. There are the upwind explicit scheme, central explicit scheme, and the Lax-Wendroff explicit scheme, and so forth. Of course, our new fractional CFDM is obviously better than all of them for the space-fractional convection-dominated diffusion problems because the explicit SFDM schemes are conditionally stable and they does not allow large time steps.

Finally, we compare the CFDM with the explicit SFDM schemes simply. Table 4 shows that the errors of three different explicit schemes with different and for and at time . As a result, the much refined time steps and space steps have to be used to improve the accuracy of the numerical solution in Table 4, in order to obtain a numerical solution with comparable accuracy to the numerical solution by the CFDM with the relatively coarse spatial grids and time steps in Tables 2 or 3. Moreover, under the same assumptions of Table 3, the three explicit SFDM schemes are all unstable.

Table 4: The errors calculated by the upwind, central and Lax-Wendroff explicit difference schemes with different space steps and time steps for and at time .

-Upwind
-Central
-Lax-Wendroff

1/50
1/1000

1/80
1/1500

1/100
1/2000

1/120
1/3000

5. Conclusions

In actual applications, we often use the space-fractional derivatives to model convection-diffusion problems, and in many diffusion processes arising in physical problems, convection essentially dominates diffusion. In this paper, we propose a new scheme called the fractional CFDM and prove that it is unconditionally stable, consistent, and convergent. For this new method, even though the time steps are very coarse, the resulting CFDM solutions are more accurate than the SFDM solutions with much finer temporal grid sizes. Thus, the CFDM is more efficient and superior to the SFDM, especially for the convection-dominated problems.

Acknowledgments

This research was supported by the National Natural Science Foundations of China Grant (no. 11126179 and no. 11226247), the 211 project of Anhui University (no. 02303319 and no. KJTD002B), the Scientific Research Award for Excellent Middle-Aged and Young Scientists of Shandong Province no. BS2010HZ012, and the Nature Science Foundation of Anhui Provincial (no. 1308085QA15). The authors acknowledge the anonymous reviewers for their helpful comments.

References

D. Benson, S. W. Wheatcraft, and M. M. Meerschaert, “Application of a fractional advection-dispersion equation,” Water Resources Research, vol. 36, pp. 1403–1413, 2000.

A. Ashyralyev and Z. Cakir, “On the numerical solution of fractional parabolic partial differential equations with the Dirichlet condition,” Discrete Dynamics in Nature and Society, vol. 2012, Article ID 696179, 15 pages, 2012.

J. P. Roop, “Computational aspects of FEM approximation of fractional advection dispersion equations on bounded domains in ,” Journal of Computational and Applied Mathematics, vol. 193, no. 1, pp. 243–268, 2006.

R. L. Magin, Fractional Calculus in Bioengineering, Begell House Publishers, 2006.

K. B. Oldham and J. Spanier, The Fractional Calculus, Academic Press, New York, NY, USA, 1974.

H. Yang, “Existence of mild solutions for a class of fractional evolution equations with compact analytic semigroup,” Abstract and Applied Analysis, vol. 2012, Article ID 903518, 15 pages, 2012.

M. De la Sen, “Positivity and stability of the solutions of Caputo fractional linear time-invariant systems of any order with internal point delays,” Abstract and Applied Analysis, vol. 2011, Article ID 161246, 25 pages, 2011.

V. J. Ervin, N. Heuer, and J. P. Roop, “Numerical approximation of a time dependent, nonlinear, space-fractional diffusion equation,” SIAM Journal on Numerical Analysis, vol. 45, no. 2, pp. 572–591, 2007.

I. Podlubny, Fractional Differential Equations, Academic Press, New York, NY, USA, 1999.

W. Wyss, “The fractional diffusion equation,” Journal of Mathematical Physics, vol. 27, no. 11, pp. 2782–2785, 1986.

R. Gorenflo and E. A. Abdel-Rehim, “Convergence of the Grünwald-Letnikov scheme for time-fractional diffusion,” Journal of Computational and Applied Mathematics, vol. 205, no. 2, pp. 871–881, 2007.

C. Lubich, “Discretized fractional calculus,” SIAM Journal on Mathematical Analysis, vol. 17, no. 3, pp. 704–719, 1986.

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

F. Liu, V. Anh, and I. Turner, “Numerical solution of the space fractional Fokker-Planck equation,” in Proceedings of the International Conference on Boundary and Interior Layers—Computational and Asymptotic Methods (BAIL '02), vol. 166, no. 1, pp. 209–219, 2004.

W. Deng, “Finite element method for the space and time fractional Fokker-Planck equation,” SIAM Journal on Numerical Analysis, vol. 47, no. 1, pp. 204–226, 2008.

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.

H. Wang, R. E. Ewing, G. Qin, S. L. Lyons, M. Al-Lawatia, and S. Man, “A family of Eulerian-Lagrangian localized adjoint methods for multi-dimensional advection-reaction equations,” Journal of Computational Physics, vol. 152, no. 1, pp. 120–163, 1999.

P. Zhuang, F. Liu, V. Anh, and I. Turner, “New solution and analytical techniques of the implicit numerical method for the anomalous subdiffusion equation,” SIAM Journal on Numerical Analysis, vol. 46, no. 2, pp. 1079–1095, 2008.

R. Schumer, D. A. Benson, M. M. Meerschaert, and B. Baeumer, “Multiscaling fractional advection-dispersion equations and their solutions,” Water Resources Research, vol. 39, pp. 1022–1032, 2003.

L. Su, W. Wang, and Z. Yang, “Finite difference approximations for the fractional advection-diffusion equation,” Physics Letters A, vol. 373, no. 48, pp. 4405–4408, 2009.

H. Wang, H. K. Dahle, R. E. Ewing, M. S. Espedal, R. C. Sharpley, and S. Man, “An ELLAM scheme for advection-diffusion equations in two dimensions,” SIAM Journal on Scientific Computing, vol. 20, no. 6, pp. 2160–2194, 1999.

J. Douglas, and T. F. Russell, “Numerical methods for convection-dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures,” SIAM Journal on Numerical Analysis, vol. 19, no. 5, pp. 871–885, 1982.

K. W. Morton and D. F. Mayers, Numerical Solution of Partial Differential Equations, Cambridge University Press, Cambridge, UK, 1994.

J. M. Varah, “A lower bound for the smallest singular value of a matrix,” Linear Algebra and Its Applications, vol. 11, pp. 3–5, 1975.

R. S. Varga, Matrix Iterative Analysis, Springer, Berlin, Germany, 2000.

H. Wang, “An optimal-order error estimate for a family of ELLAM-MFEM approximations to porous medium flow,” SIAM Journal on Numerical Analysis, vol. 46, no. 4, pp. 2133–2152, 2008.