# A Fast Implicit Finite Difference Method for Fractional Advection-Dispersion Equations with Fractional Derivative Boundary Conditions

Authors: Taohua Liu; Muzhou Hou
Journal: Advances in Mathematical Physics
 (2017)

---

## Abstract

Fractional advection-dispersion equations, as generalizations of classical integer-order advection-dispersion equations, are used to model the transport of passive tracers carried by fluid flow in a porous medium. In this paper, we develop an implicit finite difference method for fractional advection-dispersion equations with fractional derivative boundary conditions. First-order consistency, solvability, unconditional stability, and first-order convergence of the method are proven. Then, we present a fast iterative method for the implicit finite difference scheme, which only requires storage ofO(K) and computational cost of O(Klog⁡K). Traditionally, the Gaussian elimination method requires storage of O(K2) and computational cost of O(K3). Finally, the accuracy and efficiency of the method are checked with a numerical example.


---

## Body

## 1. Introduction

Fractional derivatives are almost as old as their more familiar integer-order counterparts [1, 2]. Fractional advection-dispersion equations are used to model many problems in physics, biology, and finance [3–5]. Fractional differential equations have attracted fantastic attention of many authors in recent years [6–11].The fractional space derivatives are used to establish anomalous diffusion and dispersion models, where the propagation velocity of the particle beam is different from the classical Brown motion model. When the space second derivative in a diffusion or dispersion model is replaced with a fractional space derivative, it leads to superdiffusion. In this paper, we consider a one-dimensional fractional advection-dispersion equation with fractional derivative boundary conditions:(1)∂pr,t∂t=-vr∂pr,t∂r+dr∂αpr,t∂rα+fr,t,0<r<R,0<t,(2)p0,t=0,βpR,t+dr∂α-1pr,t∂rα-1r=R=wt,0<t,(3)pr,0=ur,0≤r≤R,where 1<α≤2; p(r,t) is the distribution function of the diffusion quantity; v(r) is flow velocity; d(r) is diffusion coefficient; f(r,t) is the source term. We assume v(r)≥0 and d(r)≥0, so that the flow is from left to right. The space fractional partial derivatives in (1) and (2) are left-sided Riemann-Liouville fractional derivatives over the x domain. For a function g(x) over the interval L<x<R, the left-sided Riemann-Liouville fractional derivative of order μ is defined by (4)dμgxdxμ=1Γn-μdndxn∫Lxgξx-ξμ+1-ndξ,where n is an integer, such that n-1<μ≤n. See [12] for details. If μ is an integer, then (4) gives the standard integer derivative. β=0 corresponds to a fractional Neumann boundary condition and β>0 corresponds to a fractional Robin boundary condition. In this paper, we only discuss β>0.For a one-dimensional advection-dispersion equation with constant coefficients, Benson et al. [13] got available analytical solutions using Fourier transform method. However, many problems require a model with variable coefficients. Meerschaert and Tadjeran [14] developed finite difference approximations for one-dimensional fractional advection-dispersion equations with Dirichlet boundary conditions. And, recently, some authors have discussed numerical method for fractional equations with fractional derivative boundary conditions. Jia and Wang [15] developed fast implicit finite difference methods for two-sided space fractional diffusion equations with fractional derivative boundary conditions. Guo et al. [16] developed implicit finite difference method for fractional percolation equation with Dirichlet and fractional boundary conditions. As far as we know, the study on numerical method for one-dimensional advection-dispersion equations with fractional derivative boundary conditions is still limited. This motivates us to examine a fast implicit finite difference method for fractional advection-dispersion equations with fractional derivative boundary conditions in this paper.The rest of the paper is organized as follows. In Section2, we develop an implicit finite difference method and analyze its consistency. In Section 3, the solvability, unconditional stability, and convergence of the method are proven. In Section 4, we present a fast iterative method for the implicit finite difference scheme, which only requires storage of O(K) and computational cost of O(Klog⁡K). In Section 5, a numerical experiment is presented to verify the accuracy and efficiency of the proposed scheme.


## 2. The Implicit Finite Difference Method for the Advection-Dispersion Equations with Fractional Derivative Boundary Conditions and Its Consistency

First, leth=R/K and Δt be the spatial mesh width and the time step, where K is a positive integer. We define a spatial partition by ri=ih for i=0,1,…,K and a temporal partition by tn=nΔt for n=0,1,…,N. Let vi=v(ri), di=d(ri), fin=f(ri,tn), wn=w(tn), ui=u(ri). Denote the exact and numerical solutions at the mesh point (ri,tn) by Pin=p(ri,tn) and pin, respectively. Pi0=u0(ri) and p0n=0.Then, we employ the backward Euler scheme to approximate the first-order time derivative; a right-shifted Grünwald formula to approximate left-sidedα-order Riemann-Liouville fractional derivative and standard Grünwald formula to approximate α-1-order Riemann-Liouville fractional derivative are defined as follows [12]:(5)∂p∂tri,tn+1=Pin+1-PinΔt+OΔt,∂p∂rri,tn+1=Pin+1-Pi-1n+1h+Oh,∂αp∂rαri,tn+1=1hα∑j=0i+1gα,jPi-j+1n+1+Oh,∂α-1p∂rα-1R,tn+1=1hα-1∑j=0Kgα-1,jPK-jn+1+Oh,where gδ,l,(δ=α-1;α) is normalized Grünwald weights:(6)gδ,l=Γl-δΓ-δΓl+1.Moreover, the normalized Grünwald weights satisfy the properties in the following Lemma 1.Lemma 1.
Letμ be a positive real number and the integer n≥1. We have(i)
gμ,0=1, gμ,j=1-μ+1/jgμ,j-1 for j≥1;(ii)
gμ,1<gμ,2<⋯<0, ∑j=0n‍gμ,j>0 for 0<μ<1;(iii)
gμ,2>gμ,3>⋯>0, ∑j=0n‍gμ,j<0 for 1<μ<2.Therefore implicit finite difference method for the advection-dispersion equations with fractional derivative boundary conditions (1)–(3) can be written as follows:(7)pin+1-pinΔt=-vipin+1-pi-1n+1h+dihα∑j=0i+1gα,kpi-j+1n+1+fin,fori=1,2,…,K-1,p0n=0,βpKn+1+dKhα-1∑j=0Kgα-1,jpK-jn+1=wn+1,pi0=ui.Denote the local truncation error byRin+1 for 1≤i≤K. Then, from (7), we have (8)Rin+1=Pin+1-PinΔt--viPin+1-Pi-1n+1h+dihα∑k=0i+1gα,kPi-k+1n+1-fin=OΔt+Oh,1≤i≤K-1,RNn+1=βPKn+1-dKhα-1∑j=0Kgα-1,jPK-jn+1-wn+1=Oh.From (8), this implicit finite difference scheme is consistent.Equations (7) may be rearranged as follows:(9)-gα,0Bipi+1n+1+1+Ei-gα,1Bipin+1-Ei+gα,2Bipi-1n+1-Bi∑k=3i+1gα,kpi-k+1n+1=pin+Δtfin+1,βhα-1+dKgα-1,0pKn+1+dK∑j=1K-1gα-1,jpK-jn+1=hα-1wn+1,where Bi=diΔt/hα, Ei=viΔt/h. Equations (9) can be written in matrix form: (10)Apn+1=qn+fn+1,where(11)pn+1=p1n+1,p2n+1,…,pKn+1T,qn=p1n,p2n,…,pK-1n,0T,fn+1=Δtf1n+1,Δtf2n+2,…,ΔtfK-1n+1,hα-1wn+1T,where the coefficient matrix A=[ai,j]K×K and its entries are (12)ai,j=-gα,i+1-jBi,for1≤j≤i-2,1≤i≤K-1;-Ei-gα,2Bi,for1≤j=i-1≤K-1;1+Ei-gα,1Bi,forj=i,1≤i≤K-1;-gα,0Bi,forj=i+1,1≤i≤K-1;0,forj=i+2,1≤i≤K-1;gα-1,K-jdK,for1≤j≤K-1,i=K;βhα-1+dKgα-1,0,fori=j=K.


## 3. Solvability, Stability, and Convergence Analysis

Having developed a numerical scheme and shown that it is consistent, in the following theorems, we show this method is solvable, unconditionally stable, and convergent.Theorem 2 (solvability).
Ifβ>0, (10) is uniquely solvable.Proof.
Letri be the sum of absolute values of elements along the ith row excluding the diagonal elements Ai,i. According to Lemma 1, we can get(13)ri=∑j=1,j≠iKAi,j=∑j=1i-2-Bigi+1-jα+Ei+Big2α+Big0α=∑j=0igjα-g1αBi+Ei<-g1αBi<Ai,i1≤i≤K-1,rK=∑j=1K-1AK,j=∑j=1K-1dKgK-jα-1=-∑j=1K-1dKgjα-1=-dK∑j=0K-1gjα-1-g0α-1<dKg0α-1<AK,K.Matrix A is strictly diagonal, which guarantees the invertibility of matrix A. Hence, (10) is uniquely solvable.To discuss the stability of the numerical method, we denote byp~in(1≤i≤K) the approximate solution of the difference scheme with the initial condition p~i0(1≤i≤K) and define εin=pin-p~in, εn=(ε1n,ε2n,…,εKn)T, εn∞=max1≤i≤K⁡εin.Theorem 3.
The implicit finite difference method for the advection-dispersion equations with fractional derivative boundary conditions defined by (7) is unconditionally stable.Proof.
From (9) and the definition of εin, we have (14)-gα,0Biεi+1n+1+1+Ei-gα,1Biεin+1-Ei+gα,2Biεi-1n+1-Bi∑k=3i+1gα,kεi-k+1n+1=εin,(15)βhα-1+dKgα-1,0εKn+1+dK∑j=1K-1gα-1,jεK-jn+1=0.By Lemma 1, we have (16)∑j=1K-1gα-1,j=-∑j=1K-1gα-1,j<gα-1,0.It follows from (15) that (17)εKn+1=dK∑j=1K-1gα-1,jεK-jn+1hα-1β+gα-1,0dK≤dK∑j=1K-1gα-1,jεK-jn+1hα-1β+gα-1,0dK<max1≤i≤K-1⁡εin.Therefore, we suppose that εn∞=εi0n, where 1≤i0≤N-1. By Lemmas 1(i) and (iii), then (18)εn+1∞=εi0n+1<1+Ei0-Ei0-∑k=0i0+1gα,kBi0εi0n+1=1+Ei0-gα,1Bi0εi0n+1-gα,0Bi0εi0n+1-Ei0+gα,2Bi0εi0n+1-∑k=3i0+1gα,kBi0εi0n+1≤1+Ei0-gα,1Bi0εi0n+1-gα,0Bi0εi0+1n+1-Ei0+gα,2Bi0εi0-1n+1-∑k=3i0+1gα,kBi0εi0-k+1n+1≤1+Ei0-gα,1Bi0εi0n+1-gα,0Bi0εi0+1n+1-Ei0+gα,2Bi0εi0-1n+1-∑k=3i0+1gα,kBi0εi0-k+1n+1=εi0n≤εn∞,which implies that εn∞<ε0∞.
Thus, the implicit finite difference method defined by (7) is unconditionally stable.For the convergence analysis, we defineein=Pin-pin(1≤i≤K) and en=(e1n,e2n,…,eKn)T, en∞=max1≤i≤K⁡ein.Theorem 4.
There exists a positive constantC, independent of h and Δt, such that (19)en+1∞≤CΔt+h,1≤n≤N.Proof.
From (9), we have (20)-gα,0Biei+1n+1+1+Ei-gα,1Biein+1-Ei+gα,2Biei-1n+1-Bi∑k=3i+1gα,kei-k+1n+1=ein+ΔtRin+1,βhα-1+dKgα-1,0eKn+1+dK∑j=0K-1gα-1,jeK-jn+1=hα-1RKn+1,with e0n+1=0, and ei0=0, i=0,1,…,K. If en+1∞=eKn+1≥max1≤i≤K-1⁡ein+1, we have (21)βhα-1+dKgα-1,0eKn+1+dK∑j=0K-1gα-1,jeK-jn+1=hα-1RKn+1.By Lemma 1(ii), then (22)βhα-1+dK∑j=0K-1gα-1,jeKn+1=βhα-1+dKgα-1,0eKn+1-dK∑j=1K-1gα-1,jeKn+1≤βhα-1+dKgα-1,0eKn+1-dK∑j=1K-1gα-1,jeK-jn+1≤βhα-1+dKgα-1,0eKn+1-dK∑j=1K-1gα-1,jeK-jn+1≤βhα-1+dKgα-1,0eKn+1+dK∑j=1K-1gα-1,jeK-jn+1=hα-1RKn+1.With Lemma 1 and the Stirling formula for the Gamma function [17], we have (23)∑j=0K-1gα-1,j-1=-1K-1α-2K-1-1=ΓKΓ2-αΓK-α+1=OKα-1,as K→∞. Combining (22) and (23), we get (24)eKn+1≤dK∑j=0K-1gα-1,j-1hα-1RKn+1≤C1h.If en+1∞=ei0n+1>eKn+1, 1≤i0≤K-1, then (25)en+1∞=ei0n+1<1+Ei0-Ei0-∑k=0i0+1gα,kBi0ei0n+1=1+Ei0-gα,1Bi0ei0n+1-gα,0Bi0ei0n+1-Ei0+gα,2Bi0ei0n+1-∑k=3i0+1gα,kBi0ei0n+1≤1+Ei0-gα,1Bi0ei0n+1-gα,0Bi0ei0+1n+1-Ei0+gα,2Bi0ei0-1n+1-∑k=3i0+1gα,kBi0ei0-k+1n+1≤1+Ei0-gα,1Bi0ei0n+1-gα,0Bi0ei0+1n+1-Ei0+gα,2Bi0ei0-1n+1-∑k=3i0+1gα,kBi0ei0-k+1n+1=ei0n+ΔtRi0n+1≤max1≤i≤K-1⁡ein+ΔtC2Δt+h.By induction, we can finally obtain (26)en∞≤C1h+nΔtC2Δt+h≤CΔt+h.This completes the proof.


## 4. Storage and Fast Krylov Subspace Method of the Coefficient Matrix

The following theorem shows that the storage of matrixA can be stored in O(K) memories, instead of O(K2) memories.Theorem 5.
Coefficient matrixA can be stored in 4K+1 memories.Proof.
We express coefficient matrixA defined by (12) in block form as follows: (27)A=AK-1,K-1AK-1,KAK,K-1TaK,K,where AK-1,K-1 is the (K-1)×(K-1) matrix that consists of the first K-1 rows and the first K-1 columns. AK,K-1T is an (K-1)-dimensional row vector that consists of the first K-1 entries in the Kth row of matrix A, and AK-1,K is an (K-1)-dimensional column vector that consists of the first K-1 entries in the Kth column of matrix A. Matrix AK-1,K-1 has the similar form as that in the finite difference method for (1)–(3) with the fractional derivative boundary conditions [15]. AK-1,K-1 can be shown to have the following decomposition: (28)AK-1,K-1=I-diag⁡Bii=1K-1·Aα+diag⁡Eii=1K-1·Aβ,where I is the identity matrix of the order K-1. And A(α) and A(β) are (N-1)th-order Toeplitz matrix given as follows:(29)Aα=gα,1gα,00⋯00gα,2gα,1gα,0⋱00gα,3gα,2gα,1⋱00⋮⋱⋱⋱⋱⋮gα,K-2gα,K-3gα,K-4⋱gα,1gα,0gα,K-1gα,K-2gα,K-3⋯gα,2gα,1,Aβ=100⋯00-110⋱000-11⋱00⋮⋱⋱⋱⋱⋮000⋱10000⋯-11.According to the definition of (12), we can write AK-1,K and AK,K-1T as follows: (30)AK-1,K=-BK-1·gα,0·IK-1,AK,K-1T=dK·gα,where g(α)=(gα,K-1,gα,K-2,…,gα,1), and IK-1 is the (K-1)th column of I.Thus, in order to storeA, we need only to store (B1,B2,…,BK-1)T,(E1,E2,…,EK-1)T, (gα,0,gα,1,…,gα,K-1)T, g(β)=(1,-1,…,0)T, AK,K-1=(aK,1,aK,2,…,aK,K-1)T, I, and aK,K, which totally have 4K+1 parameters, instead of storing the super triangular matrix A in a straightforward but naive way that required K2 memories.Theorem 6.
For any vectorμK∈RK and TK to be a Toeplitz matrix of order K, the matrix-vector multiplication TKμK can be carried out in O(Klog⁡K) operations without any lossy compression.Proof.
A Toeplitz matrix is a matrix in which each descending diagonal from left to right is constant. LetTK be a Toeplitz matrix of order K of the following form: (31)TK=t0t-1t-2⋯t2-Kt1-Kt1t0t-1⋱t3-Kt2-Kt2t1t0⋱⋱t3-K⋮⋱⋱⋱⋱⋮tK-2tK-3⋱⋱t0t-1tK-1tK-2tK-3⋯t1t0.The matrix-vector multiplication TKμK can be carried out via the fast Fourier transform (FFT). First, we embed Toeplitz matrix TK into a 2K×2K circulant matrix C2N as follows [18, 19]: (32)C2K=TKBKBKTK,BK=0t1-K⋯t-2t-1tK-10t1-K⋱t-2⋮tK-10⋱⋯t2⋱⋱⋱t1-Kt1t2⋯tK-10.It is known that a circulant matrix C2K can be decomposed as follows [19, 20]: (33)C2K=F2K-1diag⁡F2Kc2KF2K,where c2K is the first column vector of C2K. F2K is the 2K×2K discrete Fourier transform matrix in which (j,k)-entry (F2K)j,k of the matrix F2K is given by (34)F2Kj,k=12Kexp⁡2πijk2K,i=-1,0≤j,k≤2K-2.For any vector μK∈RK, we extend the vector by 0 to a 2K-dimensional vector μ2K=μK0. It is well known that the matrix-vector multiplication F2Kμ2K can be carried out in O(2Klog⁡2K)=O(Klog⁡K) operations via the fast Fourier transform. Equation (33) shows that F2Kμ2K can be evaluated in O(Klog⁡K) operations. We observe (35)C2KμK0=F2K-1diag⁡F2Kc2KF2KμK0.At the same time, (36)C2KμK0=TKBKBKTKμK0=TKμKBKμK.Thus, the upper half of the vector multiplication F2Kμ2K gives the matrix-vector multiplication TKμK. So we can evaluate TKμK in O(Klog⁡K) operations. Then, by (28) we evaluate AK-1,K-1μK-1 in O(Klog⁡K) operations. Then, from (27), for any vector μK∈RK, the matrix-vector multiplication TKμK is evaluated in O(Klog⁡K) operations without any lossy compression.Remark 7.
If the system matrixA from (10) is a dense and full coefficient matrix, the Gaussian elimination method requires storage of O(K2) and computational cost of O(K3). In this paper, the system matrix A is supertriangular; the traditional Gaussian elimination method requires storage of O(K2) and computational cost of O(K2). We present a fast iterative method for the implicit finite difference scheme, which only requires storage of O(K) and computational cost of O(Klog⁡K).


## 5. A Numerical Example

The following fractional differential equation was considered:(37)∂pr,t∂t=-vr∂pr,t∂r+dr∂αpr,t∂rα+fr,tOn a finite domain[0,1]. The flow velocity is v(r)=α/4-2α. Diffusion coefficient is d(r)=Γ(2-α)·rα-1. The source term is f(r,t)=5e-t(r2+3α-4/4-2α). The initial function is u(r)=5r(1-r). And the boundary conditions are w(t)=5e-t(1/2-α-2/(3-α)(2-α)).Note that the exact solution to this problem isp(r,t)=5e-tr(1-r). In the numerical experiments, we consider three different α in the case, respectively. Table 1 shows ehN∞ the maximum error of the numerical solution of at the time T=1 and β=1.0. The stability and its order O(Δt+h) of convergence are proved.Table 1
Error behavior for the implicit finite difference scheme (7), for example, at time T=1 and β=1.0.
Δ
t
=
h
α
=
1.9
α
=
1.8
α
=
1.6
e
h
N
∞
Error
-
rate
e
h
N
∞
Error
-
rate
e
h
N
∞
Error
-
rate
2
-
4
0.0299
—
0.0351
—
0.0565
—
2
-
5
0.0155
1.9290
0.0185
1.8973
0.0308
1.8344
2
-
6
0.0079
1.9620
0.0096
1.9271
0.0164
1.8781
2
-
7
0.0040
1.9750
0.0049
1.9592
0.0087
1.8851
2
-
8
0.0020
2.0000
0.0025
1.9600
0.0047
1.8510Figures1–3 show the numerical solutions obtained by applying the Euler method (10) discussed above, with α=1.9, α=1.8, and α=1.6 respectively, with T=1, β=1.0 and Δt=h=2-6. The numerical solution compares well with the exact analytic solution to the fractional partial differential equation in these test cases.Figure 1
The behavior solution atT=1 with Δt=h=2-6, α=1.9.Figure 2
The behavior solution atT=1 with Δt=h=2-6, α=1.8.Figure 3
The behavior solution atT=1 with Δt=h=2-6, α=1.6.



---

*Source: 8716752-2017-09-24.xml*