alexa
Reach Us +44-7482-875032
High-order Accurate Numerical Methods for Solving the Space Fractional Advection-dispersion Equation | OMICS International
ISSN: 2168-9679
Journal of Applied & Computational Mathematics
Make the best use of Scientific Research and information from our 700+ peer reviewed, Open Access Journals that operates with the help of 50,000+ Editorial Board Members and esteemed reviewers and 1000+ Scientific associations in Medical, Clinical, Pharmaceutical, Engineering, Technology and Management Fields.
Meet Inspiring Speakers and Experts at our 3000+ Global Conferenceseries Events with over 600+ Conferences, 1200+ Symposiums and 1200+ Workshops on Medical, Pharma, Engineering, Science, Technology and Business
All submissions of the EM system will be redirected to Online Manuscript Submission System. Authors are requested to submit articles directly to Online Manuscript Submission System of respective journal.

High-order Accurate Numerical Methods for Solving the Space Fractional Advection-dispersion Equation

Fenga L1, Zhuangb P2, Liu F1* and Turnera I1

1School of Mathematical Sciences, Queensland University of Technology, GPO Box 2434, Brisbane, Qld 4001, Australia

2School of Mathematical Sciences, Xiamen University, Xiamen 361005, China

*Corresponding Author:
Liu F
School of Mathematical Sciences
Queensland University of Technology
GPO Box 2434, Brisbane, Qld 4001, Australia
Tel: +61731381329
E-mail: [email protected]

Received date: November 14, 2015; Accepted date: December 14, 2015; Published date: December 18, 2015

Citation: Fenga L, Zhuangb P, Liu F, Turnera I (2015) High-order Accurate Numerical Methods for Solving the Space Fractional Advection-dispersion Equation. J Appl Computat Math 4:279. doi:10.4172/2168-9679.1000279

Copyright: © 2015 Fenga L, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Visit for more related articles at Journal of Applied & Computational Mathematics

Abstract

In this paper, we consider a type of space fractional advection-dispersion equation, which is obtained from the classical advection-diffusion equation by replacing the spatial derivatives with a generalized derivative of fractional order. Firstly, we utilize the modified weighted and shifted Grunwald difference operators to approximate the Riemann-Liouville fractional derivatives and present the finite volume method. Specifically, we discuss the Crank-Nicolson scheme and solve it in matrix form. Secondly, we prove that the scheme is unconditionally stable and convergent with the accuracy of O(τ2 + h2). Furthermore, we apply an extrapolation method to improve the convergence order, which can be O(τ4 + h4). Finally, two numerical examples are given to show the effectiveness of the numerical method, and the results are in excellent agreement with the theoretical analysis.

Keywords

Finite volume method; Riemann-Liouville fractional derivative; Fractional advection-dispersion equation; Crank-Nicolson scheme; Extrapolation method

Introduction

The fractional advection-dispersion equation (FADE) is a generalization of the classical advection-dispersion equation (ADE). It provides a useful description of transport dynamics in complex systems which are governed by anomalous diffusion and nonexponential relaxation [1]. The FADE was firstly proposed by Chaves [2] to investigate the mechanism of super diffusion and with the goal of having a model able to generate the L´evy distribution and was later generalized by Benson et al. [3,4] and has since been treated by numerous authors. Many numerical methods have been proposed for solving the FADE.

Meerschaert and Tadjeran [5] developed practical numerical methods to solve the one-dimensional space FADE with variable coefficients. Liu et al. [6] transformed the space fractional Fokker-Planck equation into a system of ordinary differential equations (method of lines), which was then solved using backward differentiation formulas. Liu et al. [7] proposed an implicit difference method (IDM) and an explicit difference method (EDM) to solve a space-time FADE. Liu et al. [8] presented a random walk model for approximating a L´evy-Feller advection-dispersion process and proposed an explicit finite difference approximation (EFDA). Momani and Odibat [9] developed two reliable algorithms; the Adomian decomposition method and variational iteration method, to construct numerical solutions of the space-time FADE in the form of a rapidly convergent series with easily computable components. Zhuang et al. [10] discussed a variable-order FADE with a nonlinear source term on a finite domain. Wang et al. [11] developed a fast characteristic finite difference method for the efficient solution of space-fractional transient ADE in one space dimension, which require less storage and computation cost while retaining the same accuracy and approximation property. Shen et al. [12] presented an explicit difference approximation and an implicit difference approximation for a space-time Riesz-Caputo FADE with initial and boundary conditions. Chen et al. [13] considered the variable-order Galilei ADE with a nonlinear source term. Sousa [14] developed a numerical method for fractional advection diffusion problems with source terms by using a Lax-Wendroff-type time discretization procedure, which was explicit and second order accurate. Zhang [15] proposed a novel numerical scheme for solving the space fractional Fokker-Planck equation with the help of the Pad´e approximation. In addition, finite element method, finite volume method, meshless method and spectral method are also employed to approximate the FADE [16-25].

In this paper, we consider the following space FADE with constant coefficients and a source term on the interval [a, b] [19,20]:

equation (1)

subject to the initial condition:

equation (2)

and the zero Dirichlet boundary conditions:

equation (3)

where 1<β<2 and 0 ≤ θ ≤ 1. The variable u(x, t) represents the concentration; V0>0 and K0>0 are the drift and the anomalous diffusion coefficients respectively; θ is the skewness; f(x, t) is a source (or absorbent) term. The left and right Riemann-Liouville fractional derivatives of β order on the finite domain [a, b] are defined as

equation (4)

equation (5)

Applying the finite volume method to solve fractional diffusion equation (FDE) has received less attention. Except Hejazi et al., Zhang et al., Hejazi et al. [19-21] and Liu et al., Yang et al., Feng et al. [26-28], to our knowledge, there are no other reports in the literature about utilizing the finite volume method for the FDE. Finite volume method for Equation.(1) was first proposed by Zhang et al. [20] based on the Riemann-Liouville definition for the fractional derivative, but without theoretical analysis. Hejazi et al. [19] derived the finite volume method for Equation.(1) by utilizing fractionally-shifted Grunwald formulae for the fractional derivative and proved the stability and convergence of the scheme, however, whose accuracy is O(τ + h).

Inspired by the shifted Grunwald difference operators proposed by Tadjeran et al. [29] and multi-step method, Tian et al. [30] presented the weighted and shifted Grunwald difference (WSGD) operators, which was of second-order accuracy; however, they only prove the case of 1<α<2. In this paper, firstly we modify the WSGD operators to satisfy the case of 0<α<1. Then, based on the modified weighted shifted Grunwald difference operators to approximate the Riemann - Liouville fractional derivatives, we derive the Crank - Nicolson scheme and finite volume method for Equation.(1) with accuracy of order O(τ2+h2); we then establish its stability and convergence. Moreover, we apply the extrapolation method to improve the convergence order, which is relatively new and can be O(τ4 + h4).

The outline of the paper is as follows. In Section 2, we modify the WSGD operators and give some lemmas and the second-order approximation for the Riemann-Liouville fractional derivatives. In Section 3, we present the finite volume method for the FADE and derive the Crank-Nicolson scheme. We proceed with the proof of the stability and convergence of the Crank-Nicolson scheme in Section 4 and apply the extrapolation method in Section 5. In order to verify our theoretical analysis, two numerical examples are carried out and the results are compared with the exact solution in Section 6. Finally, some conclusions are drawn in Section 7.

The Second-Order Approximation for the Riemann- Liouville Fractional Derivative

In the following, we suppose the symbol C ia a generic positive constant, which may take different values at different places. Firstly, in the interval [a, b], we take the mesh points xi=a + ih, i=0, 1, · · ·, m, and tn=nτ, n=0, 1, · · ·, N, where h=(b − a)/M, τ=T / N, i.e., h and τ are the uniform spatial step size and temporal step size. Now, we give some definitions about the fractional derivative.

Definition 1: Riemann-Liouville fractional derivatives [31]. The γ (n − 1<γ<n) order left and right Riemann-Liouville fractional derivatives of function v(x) on [a, b] are given by:

(6)

equation (6)

equation (7)

Definition 2: The γ (n − 1<γ<n) order left and right Grunwald- Letnikov fractional derivatives of function v(x) on [a, b] are given by [19],

equation

equation

Definition 3:

equation

equation

Generally the standard Grunwald - Letnikov difference formula is applied to approximate the Riemann-Liouville fractional derivative [19]. Meerschaert and Tadjeran [5] showed that the standard Grunwald-Letnikov difference formula was often unstable for time dependent problems and they proposed the shifted Grunwald difference operators:

equation

equation

whose accuracies are first order, i.e.,

equation (8)

equation (9)

where p is an integer and equation. In fact, the coefficients equation are the coefficients of the power series of function (1 - z)γ.

equation

for all |z| ≤ 1, and they can be evaluated recursively

equation

Tadjeran et al. [29] prove Equation.(8) under the additional hypotheses that 1<γ<2 and p is an integer. In fact, with very minor modifications, their proof holds for the more general case. Based on this point, Hejazi et al. [19] prove Equation.(8) and Equation.(9) hold for 0<γ<1 and p ia a positive real number.

Lemma 1: Let α and p be positive numbers, and suppose thatequation, n ∈ N* and all derivatives of v(x) up to order [α] + n + 2 belong to equation Then if a=-∞ and b=+∞ in Equation. (6) and Equation.(7) respectively, there exist constants cl, p independent of h, f, x such that [19],

equation

and

equation

uniformly for equation, where c0, p=1, c1, p=p – α/2.

Lemma 2: Supposing that 0<α<1, the coefficients equation satisfy [7,19],

equation

Inspired by the shifted Grunwald difference operators and multistep method, Tian et al. [30] derive the WSGD operators:

equation

equation

Lemma 3: Supposing that 1<γ<2, v(x) ∈ equation, equation and equation and their Fourier transforms belong to equation then the WSGD operators satisfy [30],

equation

equation

uniformly for equation, where p, q are integers and p ≠ q.

Tian et al. prove Lemma 3 under the additional hypotheses that 1<γ<2 and p, q are integers. In fact, with very minor modifications, their proof holds for the more general case. Similar to Lemma 1, we can conclude that Lemma 3 holds for 0<γ<1 and p, q are positive real numbers as well.

Lemma 4: Let 0<α<1 and p, q be positive numbers, and suppose that equation, n ∈ N* and all derivatives of v(x) up to order [α] + n + 2 belong to equation Then if a=- ∞ and b=+ ∞ in Equation.(6) and Equation.(7) respectively, there exist constants cl independent of h, f, x such that

equation

and

equation

uniformly for equation, where equation and p ≠ q.

Proof: The proof proceeds the same as Lemma 3, we will not repeat it here. It is also noted that, this lemma is just a special case of Theorem 4, in Kincaid and Cheney [32] with λ=0, l=m=2.

Remark 1: Considering a well-defined function v(x) on the bounded interval [a, b], if v(a)=0 and v(b)=0, the function v(x) can be zero extended for x<a and x>b. And then the α order left and right Riemann-Liouville fractional derivatives of v(x) at each point x can be approximated by the modified WSGD operators with second order accuracy.

equation

equation

To evaluate the fractional derivatives at x ± h/2 , we modify the WSGD operators by adopting (p, q)=( 1/2 , −1/2 ) (the role of p and q is symmetric). When (p, q) =(1/2, −1/2) and 0<α<1, the secondorder discrete approximations with error terms at xi+1/2 and xi−1/2 for the Riemann-Liouville fractional derivatives on the domain [a, b] are:

equation (10)

equation (11)

and

equation (12)

equation (13)

where

equation (14)

Now, we discuss the properties of the coefficients equation, k=0, 1, 2 . . .

Lemma 5: Supposing that 0<α<1, the coefficients equation satisfy

equation

Proof: Combining the definition of equationand the properties of equation, it is easy to derive the value ofequation and equation. When k ≥ 2, by Equation.(14) and Lemma 2.(2)

equation

we have equation.

Moreover,

equation

Recalling Lemma 2.(2), equation when k ≥ 1, we obtain when k ≥ 2

equation

i.e.,

equation

For the sum equation, recalling Lemma 2.(4), we have

equation

Since equation when equation for m ≥ 1.

Recalling Lemma 2.(3), we obtain

equation

Finite Volume Method

In this section, we utilize the finite volume method to approximate Equation.(1) and derive the Crank-Nicolson scheme. We begin by formulating Equation.(1) in conservative form. Using Equations (4) and (5) we have

equation (15)

We set α=β − 1, then 0<α<1. This allow us to write Equation.(15) as

equation

By the definition of the Riemann-Liouville fractional derivatives, we arrive at

equation (16)

Comparing (16) with the general transport equation

equation (17)

we define the total flux as

equation

with advective component V0u(x, t) and diffusive component

equation (18)

A finite volume discretization is applied by integrating Equation.(17) over the i th control volume [xi−1/2, xi+1/2]:

equation

Interchanging the order of differentiation and integration on the left, and performing the integration on the first term on the right, we have

equation

which leads to the standard finite volume discretization

equation (19)

where equation and equation. By the standard volume approximation, we get

equation

equation

The flux P(x, t) has both advective and diffusive components. The main feature of our method is the approximation of the diffusive flux p(xi±1/2, t) by the modified weighted and shifted Grunwald difference operators. For the advective flux, we discretize it by a standard averaging scheme

equation

equation

Thus

equation

By virtue of the mean value theorem, we know the term equation is O(h), which leads to

equation

For the diffusive flux (18), we approximate it by Equations.(10 - 13), yielding

equation

By virtue of the mean value theorem, we obtain the terms equation and equation are O(h), which leads to

equation

Then Equation.(19) may be written as

equation (20)

Now we consider the discretization of time. It is easy to conclude that,

equation

and

equation

Therefore, the Crank-Nicolson scheme of Equation.(20) at (xi, tn−1/2) gives

equation

where

equation

and

equation

Let equation be the approximation solution of u(xi, tn), then we obtain

equation (21)

The boundary and initial conditions are discretized as

equation (22)

Define equation

equation

and

equation

Setting

equation

then we can rewrite Equation.(21) in matrix form

equation (23)

Theoretical Analysis

Stability

Here, we consider the stability of the scheme (23). For simplicity, we just discuss the case of θ=1/2. To begin with, we give the following theorem concerning the coefficients of matrix Q.

Theorem 1: Suppose that 0<α<1, K0>0, V0>0. When θ=1/2 and

equation (24)

the coefficients Qij satisfy

equation (25)

i.e., Q is strictly diagonally dominant.

Proof: When θ=1/2, it is easy to obtain

equation

where r0=− τ/2h<0. For a given i, we consider the sum

equation

Recalling Lemma 5.(2), Lemma 5.(3) and (24), we can drop the absolute value

signs:

equation

The two telescoping sums have the form

equation

which by virtue of Lemma 5.(5) equals equation. Thus we have

equation

due to Lemma 5.(1), i.e.,

equation.

Thus, the proof is completed.

Remark 2: The condition (24) can be satisfied by using a sufficiently small=spatial step h.

Corollary 1: Let equation be any eigenvalue of the matrix Q, then the real part of λ satisfies

c>0. (26)

Proof: By the Gerschgorin’s circle theorem [30], we have

equation

Then

equation

therefore,

equation

Using (25) and Qii>0, it is easy to derive (26).

Corollary 2: The matrix I + Q is strictly diagonally dominant. Therefore, I + Q is invertible and Equation.(23) is solvable.

Theorem 2: The scheme (23) is unconditionally stable.

Proof: Since the eigenvalues of matrix Q satisfy Corollary 1, then the eigenvalues of the matrix (I + Q)−1(I − Q) satisfy

equation

Hence, the spectral radius of the matrix (I + Q)−1(I − Q) is less than one. Therefore, the scheme (23) is unconditionally stable.

Convergence

By Equation.(21), we notice that the local truncation error of the Crank- Nicolson scheme gives,

equation

Theorem 3: Let un be the exact solution vector of the problem (1 - 3), where un=[u(x1, tn), u(x2, tn), . . . , u(xm−1, tn)]T. Then the numerical solution Un unconditionally converges to the exact solution un as h and τ tend to zero, and

equation

Proof: Let equation denote the error at grid points (xi, tn) and equation. Combining Equations.(19) to (21), yields

equation

where

equation

Using the conditions (2), (3) and (22), we obtain the errors equation and equation for i=1, 2, · · · , m−1 and n=0, 1, · · · , N. We can write the system in matrix-vector form

equation

or

equation

Where equation and b=O(τ3 + τh2)(I + Q)−1. By iterating and noting that E0=0, we obtain

equation

Now, from Corollary 1, Corollary 2 and Theorem 2, we have ρ((I + Q) −1)<1 and ρ(M)<1. Therefore, we can choose a vector norm and induced matrix norm || · || such that ||M||<1 and ||(I + Q) −1||<1. Then upon taking norms,

equation

Thus,

equation

which completes the proof.

Improving the Order of Convergence

Here we use the Richardson extrapolation method (REM) to improve the convergence order. Suppose that imageis the approximation solution of function f(x) with an asymptotic expansion

equation

Then we have

equation

Eliminating the middle terms C1h2 on the right, we find

equation (27)

which means the approximation order of f(x) has been improved from O(h2) to O(h3). Moreover, we can improve the approximation order of f(x) from O(h3) to O(h4). Suppose that equationis the approximation solution of function f(x) with an asymptotic expansion

equation

Then we have

equation

Eliminating the middle terms C2h3 on the right, we find

equation (28)

According to Theorem 3, we know that the numerical method converges at the rate of O(τ2 + h2). In order to improve the order of convergence, we apply the numerical method on a (coarse) grid τ=h and then on a finer grid of size τ/2=h/2 and τ/4=h/4. By applying the Richardson extrapolation formulae (27) and (28) consecutively, the convergence order can be improved from O(τ2 + h2) to O(τ4 + h4) [33].

Numerical Examples

In order to demonstrate the effectiveness of the finite volume method, two examples are presented.

Example 1: At first, we consider the following space FADE

equation

where 0<α<1,

equation

equation

and the exact solution is u(x, t)=x6(1 − x)6e−t.

Table 1 describes the L2 error and convergence order of the Crank Nicolson scheme at t=1 with τ=h, where α and θ are corresponding to three distinct values, respectively. It can be seen that the numerical results are in excellent agreement with the exact solution.

θ = 0.2 α = 0.1   α = 0.5   α = 0.9  
τ = h ||E(h,τ)||2 Order ||E(h,τ)||2 Order ||E(h,τ)||2 Order
1/16 1.4643E-06   1.3445E-06   8.7375E-07  
1/32 3.4234E-07 2.10 3.2644E-07 2.04 2.1494E-07 2.02
1/64 8.4533E-08 2.02 8.1135E-08 2.01 5.3505E-08 2.01
1/128 2.1079E-08 2.00 2.0269E-08 2.00 1.3361E-08 2.00
1/256 5.2665E-09 2.00 5.0679E-09 2.00 3.3393E-09 2.00
θ = 0.5 α = 0.1   α = 0.5   α = 0.9  
τ = h ||E(h,τ)||2 Order ||E(h,τ)||2 Order ||E(h,τ)||2 Order
1/16 1.5025E-06   1.3012E-06   8.6108E-07  
1/32 3.4148E-07 2.14 3.2429E-07 2.00 2.1280E-07 2.02
1/64 8.4264E-08 2.02 8.1517E-08 1.99 5.3079E-08 2.00
1/128 2.1007E-08 2.00 2.0474E-08 1.99 1.3268E-08 2.00
1/256 5.2486E-09 2.00 5.1327E-09 2.00 3.3175E-09 2.00
θ = 0.8 α = 0.1   α = 0.5   α = 0.9  
τ = h ||E(h,τ)||2 Order ||E(h,τ)||2 Order ||E(h,τ)||2 Order
1/16 2.0424E-06   1.2834E-06   9.8391E-07  
1/32 4.9460E-07 2.05 3.2118E-07 2.00 2.4430E-07 2.01
1/64 1.2330E-07 2.00 8.0868E-08 1.99 6.1133E-08 2.00
1/128 3.0888E-08 2.00 2.0325E-08 1.99 1.5308E-08 2.00
1/256 7.7366E-09 2.00 5.0972E-09 2.00 3.8311E-09 2.00

Table 1: The error and convergence order for different α and θ.

What follow are the numerical results by applying the extrapolation method (EM). Table 2 shows the error and convergence order of the Crank-Nicolson scheme at t=1 with τ=h, θ=0.8 for the different α by applying Equation.(27) and (28) consecutively, whose accuracy can be O(τ4 + h4). The results are very encouraging and show that the extrapolation method is efficient and reliable.

equation (29)

τ = h α= 0.1 α = 0.5 α = 0.9
||E(h, τ)||2 Order ||E(h, τ)||2 Order ||E(h, τ)||2 Order
1/8 7.3032E-08   1.2766E-08   2.3311E-08  
1/16 3.4313E-09 4.41 8.8236E-10 3.85 3.0979E-10 6.23
1/32 2.0569E-10 4.06 6.1021E-11 3.85 2.2932E-11 3.76
1/64 1.3658E-11 3.91 4.0824E-12 3.90 1.4286E-12 4.00
1/128 8.8210E-13 3.95 2.6629E-13 3.94 8.9184E-14 4.00

Table 2: The error and convergence order for EM with θ=0.8.

where the δ(·) stands for the Dirac delta function. The Fourier transform of the analytical solution to (29) on an infinite domain is given by Benson et al. [4]:

equation

where 0<α<1. Provided our finite computational domain is sufficiently large that no significant concentration reaches the boundaries throughout the simulation, we may take the infinite domain solution as our benchmark.

We first discuss the asymmetric dispersion (θ=0) for different α and the solutions at time t=100 with V=0, K=1 are plotted in Figure 1. The numerical solutions were obtained by using τ=h=1/400, and they agree well with the analytical solutions. The greater asymmetry in the solution for smaller α is readily apparent and the solution for the nonfractional case (α=1) exhibits no asymmetry at all.

computational-mathematics-benchmark-numerical-solutions

Figure 1: Comparison of the benchmark and numerical solutions at time t=100 with V=0, K=1, θ=0 and varying α.

Figure 2 displays the impact of the skewness θ on the dispersion with α=0.4, V=0 and K=1. It can be seen from the figure, when θ>0.5, the dispersion is skewed forward; while when θ<0.5, the dispersion is skewed backward. The numerical solutions were obtained by using τ=h=1/400, and again agree well with the analytical solution.`

computational-mathematics-benchmark-numerical-solutions

Figure 2: Comparison of the benchmark and numerical solutions at time t=100 with V=0, K=1, α=0.4 and varying θ.

In Figure 3, we plots the analytical and numerical solutions at times t=10, 20, 40, 80 with θ=0.4, α=0.4, V=0.5 and K=1. We see that the transport process is biased to the right and at time t=80 a significant concentration has reached the right edge of the graphed region (though not the computational domain, which extends to x=400).

computational-mathematics-benchmark-numerical-solutions

Figure 3: Comparison of the benchmark and numerical solutions at different time with V=0.5, K=1, α=0.4 and θ=0.5.

Conclusions

In this paper, we have developed and demonstrated a second order finite volume method for solving a class of space FADE. Firstly, based on the modified WSGD operators to approximate the Riemann-Liouville fractional derivatives, applying the finite volume method, we derive the Crank-Nicolson scheme of the problem and rewrite it in matrix form. Subsequently, we prove that the scheme is unconditionally stable and convergent with the accuracy of O(τ2 + h2). Moreover, to improve the convergence order, we employ the extrapolation method at both the time step and space step. Finally, some numerical results are given to show the stability, consistency and convergence of our computational approach.

This technique could be extending to two-dimensional or threedimensional problems with complex regions. In the future, we would like to investigate finite volume method for the space FADE in high dimensions.

References

Select your language of interest to view the total content in your interested language
Post your comment

Share This Article

Article Usage

  • Total views: 8971
  • [From(publication date):
    December-2015 - Aug 20, 2019]
  • Breakdown by view type
  • HTML page views : 8764
  • PDF downloads : 207
Top