A High-Accuracy MOC/FD Method for Solving Fractional Advection-Diffusion Equations
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.
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 , physics [2, 3], engineering , mathematics , and science [6, 7], such as diffusion , 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 , 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  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 . Liu et al.  used L2 scheme form  to discretize the fractional Fokker-Planck equation. Deng  proposed finite element method for the space and time fractional Fokker-Planck equation. Meerschaert and Tadjeran  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 . For example, their mathematical models are basically advection dominated. Fractional derivatives play an important role in modelling particle transport in anomalous diffusion . The space-fractional advection-dominated diffusion equation  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 .
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  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
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  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 : where is a positive integer, , and , which can be calculated by
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 .
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
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.
3.2. Convergence Analysis
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 :
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.
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.
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.
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.
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.
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.
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.
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.
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.View at: Google Scholar
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.View at: MathSciNet
I. Podlubny, Fractional Differential Equations, Academic Press, New York, NY, USA, 1999.View at: MathSciNet
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.View at: Publisher Site | Google Scholar | MathSciNet
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.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
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.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
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.View at: Google Scholar
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.View at: Google Scholar
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.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
K. W. Morton and D. F. Mayers, Numerical Solution of Partial Differential Equations, Cambridge University Press, Cambridge, UK, 1994.View at: MathSciNet