^{1}School of Mathematical Sciences, Queensland University of Technology, GPO Box 2434, Brisbane, Qld 4001, Australia
^{2}School of Mathematical Sciences, Xiamen University, Xiamen 361005, China
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
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.
Finite volume method; Riemann-Liouville fractional derivative; Fractional advection-dispersion equation; Crank-Nicolson scheme; Extrapolation method
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]:
(1)
subject to the initial condition:
(2)
and the zero Dirichlet boundary conditions:
(3)
where 1<β<2 and 0 ≤ θ ≤ 1. The variable u(x, t) represents the concentration; V_{0}>0 and K_{0}>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
(4)
(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}+h^{2}); 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} + h^{4}).
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.
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 x_{i}=a + ih, i=0, 1, · · ·, m, and t_{n}=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)
(6)
(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],
Definition 3:
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:
whose accuracies are first order, i.e.,
(8)
(9)
where p is an integer and . In fact, the coefficients are the coefficients of the power series of function (1 - z)^{γ}.
for all |z| ≤ 1, and they can be evaluated recursively
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 that, n ∈ N* and all derivatives of v(x) up to order [α] + n + 2 belong to Then if a=-∞ and b=+∞ in Equation. (6) and Equation.(7) respectively, there exist constants c_{l}, p independent of h, f, x such that [19],
and
uniformly for , where c0, p=1, c1, p=p – α/2.
Lemma 2: Supposing that 0<α<1, the coefficients satisfy [7,19],
Inspired by the shifted Grunwald difference operators and multistep method, Tian et al. [30] derive the WSGD operators:
Lemma 3: Supposing that 1<γ<2, v(x) ∈ , and and their Fourier transforms belong to then the WSGD operators satisfy [30],
uniformly for , 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 , n ∈ N* and all derivatives of v(x) up to order [α] + n + 2 belong to Then if a=- ∞ and b=+ ∞ in Equation.(6) and Equation.(7) respectively, there exist constants cl independent of h, f, x such that
and
uniformly for , where 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.
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 x_{i+1/2} and x_{i−1/2} for the Riemann-Liouville fractional derivatives on the domain [a, b] are:
(10)
(11)
and
(12)
(13)
where
(14)
Now, we discuss the properties of the coefficients , k=0, 1, 2 . . .
Lemma 5: Supposing that 0<α<1, the coefficients satisfy
Proof: Combining the definition of and the properties of , it is easy to derive the value of and . When k ≥ 2, by Equation.(14) and Lemma 2.(2)
we have .
Moreover,
Recalling Lemma 2.(2), when k ≥ 1, we obtain when k ≥ 2
i.e.,
For the sum , recalling Lemma 2.(4), we have
Since when for m ≥ 1.
Recalling Lemma 2.(3), we obtain
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
(15)
We set α=β − 1, then 0<α<1. This allow us to write Equation.(15) as
By the definition of the Riemann-Liouville fractional derivatives, we arrive at
(16)
Comparing (16) with the general transport equation
(17)
we define the total flux as
with advective component V_{0}u(x, t) and diffusive component
(18)
A finite volume discretization is applied by integrating Equation.(17) over the i th control volume [x_{i−1/2}, x_{i+1/2}]:
Interchanging the order of differentiation and integration on the left, and performing the integration on the first term on the right, we have
which leads to the standard finite volume discretization
(19)
where and . By the standard volume approximation, we get
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(x_{i±1/2}, t) by the modified weighted and shifted Grunwald difference operators. For the advective flux, we discretize it by a standard averaging scheme
Thus
By virtue of the mean value theorem, we know the term is O(h), which leads to
For the diffusive flux (18), we approximate it by Equations.(10 - 13), yielding
By virtue of the mean value theorem, we obtain the terms and are O(h), which leads to
Then Equation.(19) may be written as
(20)
Now we consider the discretization of time. It is easy to conclude that,
and
Therefore, the Crank-Nicolson scheme of Equation.(20) at (x_{i}, t_{n−1/2}) gives
where
and
Let be the approximation solution of u(x_{i}, t_{n}), then we obtain
(21)
The boundary and initial conditions are discretized as
(22)
Define
and
Setting
then we can rewrite Equation.(21) in matrix form
(23)
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, K_{0}>0, V_{0}>0. When θ=1/2 and
(24)
the coefficients Q_{ij} satisfy
(25)
i.e., Q is strictly diagonally dominant.
Proof: When θ=1/2, it is easy to obtain
where r_{0}=− τ/2h<0. For a given i, we consider the sum
Recalling Lemma 5.(2), Lemma 5.(3) and (24), we can drop the absolute value
signs:
The two telescoping sums have the form
which by virtue of Lemma 5.(5) equals . Thus we have
due to Lemma 5.(1), i.e.,
.
Thus, the proof is completed.
Remark 2: The condition (24) can be satisfied by using a sufficiently small=spatial step h.
Corollary 1: Let 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
Then
therefore,
Using (25) and Q_{ii}>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
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,
Theorem 3: Let u^{n} be the exact solution vector of the problem (1 - 3), where u^{n}=[u(x_{1}, t_{n}), u(x_{2}, t_{n}), . . . , u(x_{m−1}, t_{n})]^{T}. Then the numerical solution U^{n} unconditionally converges to the exact solution u^{n} as h and τ tend to zero, and
Proof: Let denote the error at grid points (x_{i}, t_{n}) and . Combining Equations.(19) to (21), yields
where
Using the conditions (2), (3) and (22), we obtain the errors and for i=1, 2, · · · , m−1 and n=0, 1, · · · , N. We can write the system in matrix-vector form
or
Where and b=O(τ^{3} + τh^{2})(I + Q)^{−1}. By iterating and noting that E0=0, we obtain
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,
Thus,
which completes the proof.
Here we use the Richardson extrapolation method (REM) to improve the convergence order. Suppose that is the approximation solution of function f(x) with an asymptotic expansion
Then we have
Eliminating the middle terms C_{1}h^{2} on the right, we find
(27)
which means the approximation order of f(x) has been improved from O(h^{2}) to O(h^{3}). Moreover, we can improve the approximation order of f(x) from O(h^{3}) to O(h^{4}). Suppose that is the approximation solution of function f(x) with an asymptotic expansion
Then we have
Eliminating the middle terms C2h3 on the right, we find
(28)
According to Theorem 3, we know that the numerical method converges at the rate of O(τ^{2} + h^{2}). 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} + h^{2}) to O(τ^{4} + h^{4}) [33].
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
where 0<α<1,
and the exact solution is u(x, t)=x^{6}(1 − x)^{6}e^{−t}.
Table 1 describes the L_{2} 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} + h^{4}). The results are very encouraging and show that the extrapolation method is efficient and reliable.
(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]:
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.
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.`
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).
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} + h^{2}). 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.
Make the best use of Scientific Research and information from our 700 + peer reviewed, Open Access Journals