Combined compact difference scheme for the time fractional convectionв
Applied Mathematics and Computation 246 (2014) 464–473
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Combined compact difference scheme for the time fractional convection–diffusion equation with variable coefficients q Mingrong Cui School of Mathematics, Shandong University, Jinan, Shandong 250100, China
a r t i c l e
i n f o
Keywords: Compact Combined Fractional differential equation Convection–diffusion Mixed Difference High-order
a b s t r a c t Fourth-order combined compact finite difference scheme is given for solving the time fractional convection–diffusion–reaction equation with variable coefficients. We introduce the flux as a new variable and transform the original equation into a system of two equations. Compact difference is used as a high-order approximation for spatial derivatives of integer order in the coupled partial differential equations. The Caputo fractional derivative is discretized by a high-order approximation. Both Dirichlet and Robin boundary conditions are discussed. Convergence analysis is given for the problem of integer order with constant coefficients under some assumption. Numerical results are provided to verify the accuracy and efficiency of the proposed algorithm. Ó 2014 Elsevier Inc. All rights reserved.
1. Introduction Fractional differential equations(FDEs) can be used to in mathematical models for physical, biological, geological and financial systems, and the review article [1] and monograph [2] have given detailed discussions on fractional differential equations. In the recent years, many numerical methods have been proposed and studied for solving FDEs up to now, e.g., [3–5]. Compact finite difference schemes have high-order of accuracy and the desirable tridiagonal nature of the finite-difference equations (see [6,7]), and one-dimensional fractional sub-diffusion equation was recently solved by the compact finite dif4 ference scheme with convergence order Oðs þ h Þ in [8,9], a higher order one for the temporal variable in [10], and for two dimensional problems, ADI schemes with error analysis were given in papers [11,12], with the compact schemes for the time fractional convection–diffusion equations with constant coefficients considered in [13]. While high order methods for partial different equations with constant coefficients have been discussed by many authors, it is not so easy to give corresponding results for variable coefficients. One way to give high order difference schemes is to start from lower order scheme, then follow a substitution of the lower order truncation error by replacing the derivatives using the successive derivation of the original continuous problem. Thus high order schemes for variable coefficients must be examined carefully as many derivatives appear, and our aim in this paper is to give some easily implemented schemes for problems with variable coefficients. We must admit for problems with variable coefficients, following this way would be a tedious task. As pointed in [14], the idea of posing a boundary-value problem in the form of a first-order system as a governing conservation equation and an associated constitutive equation, and then using a mixed approximation scheme, has been
q This work was partially supported by Natural Science Foundation of Shandong Province, Grant No. BS2010HZ012, and National Natural Science Foundation of China, Grant No. 11171189. E-mail address: [email protected]
http://dx.doi.org/10.1016/j.amc.2014.08.025 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.
M.R. Cui / Applied Mathematics and Computation 246 (2014) 464–473
465
receiving increasing attention. For the steady and unsteady convection–diffusion equations of integer order, many research papers, e.g., [15,16], and recently a generalization work in [17], discussed solving them by using the compact difference scheme combined with mixed method. For two point boundary problem, this method was also discussed in paper [18]. After introducing the first derivative of the unknown solution as a new variable, then we can establish an approximation system between the derivative and the original solution. These schemes are also called combined compact difference(CCD) schemes, see papers [19,20]. The CCD schemes have significantly better spectral representation [21], and they have been proved very useful solving the shallow water equations in spherical geometry [22], the one dimensional convection diffusion equation [23], high-order compact upwind difference schemes has been developed in [24], and the two-dimensional lid-driven cavity problem [25,26], and Navier–Stokes equations [27]. In this paper we solve the time fractional convection–diffusion–reaction problems using the CCD scheme. The problems considered here are convection–diffusion–reaction problems with variable coefficients, thus they can be used for solving the fractional Fokker–Planck equation [28,29]. We give the CCD schemes and the variable coefficients can be handled easily. Besides, both Dirichlet and Robin conditions are discussed, and as we use the flux as a new variable, the solution and the flux can be obtained simultaneously without losing accuracy. We consider the following model equation C c 0 Dt p
þ bðx; tÞ
@p @ @p þ cðx; tÞp ¼ f ðx; tÞ; aðx; tÞ @x @x @x
L1 < x < L2 ; 0 < t 6 T
ð1Þ
where the Caputo fractional derivative C0 Dct v ð0 < c < 1Þ of the function v ðx; tÞ is defined by [2] C c 0 Dt
1 v ðx; tÞ ¼ Cð1 cÞ
Rt 0
ðt sÞc @@vs ðx; sÞds:
The initial condition for (1) is
pðx; 0Þ ¼ wðxÞ;
L1 < x < L2
ð2Þ
with the Dirichlet boundary condition given by
pðL1 ; tÞ ¼ u1 ðtÞ;
pðL2 ; tÞ ¼ u2 ðtÞ; 0 < t 6 T
ð3Þ
or with the Robin boundary condition
@p þ a1 ðx; tÞp jx¼L1 ¼ b1 ðtÞ; @x @p þ a2 ðx; tÞp jx¼L2 ¼ b2 ðtÞ; 0 < t 6 T: @x
ð4Þ
The paper is organized as follows. In Section 2, after introducing the flux as a new variable and then changing the original problem into two coupled equations, we use the high-order compact difference to approximate the spatial derivative and then adopt the high order discretization for the time fractional derivative. Then we get an implicit CCD scheme for the time fractional convection–diffusion equation with the local truncation error discussed. Treatment of the Robin boundary is discussed in Section 3. We give the convergence analysis for the integer order problem with constant coefficients under some assumption in Section 4. Finally, some numerical examples are given in Section 5 to verify the order of convergence of the schemes. This paper closes with a summary in Section 6. Throughout this paper, the symbol C is a generic positive constant, it may take different values at different places. We use P the ‘‘empty sum’’ convention ql¼p Sl ¼ 0 for q < p. 2. CCD scheme for Dirichlet boundary condition 2.1. Partition of the domain For the numerical solution of Eq. (1) we introduce a uniform grid of mesh points ðxj ; t n Þ, with xj ¼ L1 þ jh; j ¼ 0; 1; . . . ; N x þ 1; and t n ¼ ns; n ¼ 0; 1; . . . ; N. Here both N x and N are positive integers, h ¼ ðL2 L1 Þ=ðN x þ 1Þ is the mesh-width in x, and s ¼ T=N is the time step. For any function v ðx; tÞ, we let v nj ¼ v ðxj ; tn Þ, e.g., the theoretical solution u (i.e., the flux, see below) at the mesh point ðxj ; t n Þ is denoted by unj , and U nj stands for the numerical solution of the flux at the same mesh point. 2.2. Derivation of the CCD scheme In order to give the CCD scheme for the model equation (1), we introduce the flux uðx; tÞ ¼ aðx; tÞ @p as a new variable. @x The unknowns p and u may correspond to pressure and velocity in the flow equations in the mathematical modeling of compressible/incompressible miscible displacement problems, respectively. Then Eq. (1) becomes
(
C c 0 Dt p
þ @u bðx;tÞ u þ cðx; tÞp ¼ f ðx; tÞ; @x aðx;tÞ
¼ 0; u þ aðx; tÞ @p @x
L1 < x < L2 ; 0 < t 6 T:
We use high-order compact approximations for the first order derivatives in (5), i. e.,
ð5Þ
466
M.R. Cui / Applied Mathematics and Computation 246 (2014) 464–473
8 c > < C0 Dt pnj þ > : a1n unj þ j
bn
dx
2
1þh6 d2x
n
j
4
dx
2
1þh6 d2x
4
unj anj unj þ cnj pnj ¼ f j þ Oðh Þ;
pnj ¼ Oðh Þ;
ð6Þ
1 6 j 6 Nx ; 1 < n 6 N
where dx ; d2x denote the approximations of the first and second order derivatives, respectively. That is, for a mesh function v nj ,
dx v nj ¼
v njþ1 v nj1 2h
;
d2x v nj ¼
v nj1 2v nj þv njþ1 h2
:
We need to notice that the subscript j runs from 1 to N x , that is, we try to find numerical solutions for the unknowns pnj ð1 6 j 6 N x Þ and unj ð0 6 j 6 N x þ 1Þ. 2 Multiply the operator 1 þ h6 d2x on both sides of Eq. (6), we obtain
8 C n 2 2 2 > > 1 þ h6 d2x Dct pnj þ dx unj 1 þ h6 d2x ba u j þ 1 þ h6 d2x ½cpnj > > > 0 < 2 n 4 ¼ 1 þ h6 d2x f j þ Oðh Þ; 1 6 j 6 Nx ; 1 < n 6 N; > > > > > : dx pn þ 1 þ h2 d2x 1 u n ¼ Oðh4 Þ; 1 6 j 6 Nx ; 1 < n 6 N: j 6 a j
ð7Þ
The Caputo fractional derivative C0 Dct p appeared in (7) needs to be discretized now, and a high-order approximation can be obtained by the following lemma. Lemma 2.1 [30]. Suppose f ðtÞ 2 C 2 ½0; t k . Let
1 Rðf ðt k ÞÞ ¼ Cð1 cÞ
Z
tk
0
" # 0 k1 X f ðsÞ sc ds a f ðt Þ ðdkj1 dkj Þf ðt j Þ dk1 f ðt 0 Þ ; Cð2 cÞ 0 k ðt k sÞc j¼1
then
jRðf ðt k ÞÞj 6
" # 1 1 c 22c 00 þ ð1 þ 2c Þ max jf ðtÞjs2c ; 06t6t k Cð2 cÞ 12 2c 1c
1c
where c 2 ð0; 1Þ and aj ¼ ðj þ 1Þ j . c Using Lemma 2.1 replace the fractional derivative C0 Dt p by the following equality, C c 0 Dt pðxj ; t n Þ
¼
s c Cð2 cÞ
"
a0 pnj
n1 X
#
ðdnl1
dnl Þplj
dn1 p0j
þ Oðs2c Þ:
ð8Þ
l¼1
Using (8) and multiply both sides of (7) by
sc Cð2 cÞ, then an approximation of (1)–(3) can be given by
8 n1 X > 2 2 > > 1 þ h2 d2x Pnj ðdnl1 dnl Þ 1 þ h6 d2x Plj dn1 1 þ h6 d2x P0j > 6 > > > l¼1 > > > > c Cð2 cÞfd U n 1 þ h2 d2 b U n þ 1 þ h2 d2 ½cPn g > > þ s x < j j 6 x 6 x a j n c h2 2 > ¼ s Cð2 cÞ 1 þ 6 dx f j ; 1 6 j 6 Nx ; 1 < n 6 N; > > > > > 2 n > > dx Pnj 1 þ h6 d2x 1a U j ¼ 0; 1 6 j 6 Nx ; 1 < n 6 N; > > > > > : 0 Pj ¼ wðxj Þ; 0 6 j 6 Nx þ 1; Pn0 ¼ u1 ðtn Þ; PnNx þ1 ¼ u2 ðt n Þ; 1 6 n 6 N: With the definitions of the discrete operators, we can write (9) on the nodal points explicitly.
8 h n b n b n i n n 3 c b c > > > h s Cð2 cÞðU jþ1 U j1 Þ s Cð2 cÞ a U j1 þ 4 a U j þ a U jþ1 > > h i > > > þsc Cð2 cÞ ðcPÞn þ 4ðcPÞn þ ðcPÞn þ ðPn þ 4Pn þ Pn Þ > > j1 j jþ1 j1 j jþ1 > > > > n1 > X > < ðdnl1 dnl ÞðPlj1 þ 4Plj þ Pljþ1 Þ dn1 ðP0j1 þ 4P 0j þ P0jþ1 Þ l¼1 > > > n n n > > ¼ sc Cð2 cÞðf j1 þ 4f j þ f jþ1 Þ; 1 6 j 6 Nx ; 1 < n 6 N; > > > h i > 1 n n n >3 n n 1 1 > ðP P Þ þ U þ 4 U þ ð UÞ ¼ 0; 1 6 j 6 Nx ; 1 < n 6 N; > jþ1 j1 > h a a a j1 j jþ1 > > > : 0 n n Pj ¼ wðxj Þ; 0 6 j 6 Nx þ 1; P0 ¼ u1 ðtn Þ; PNx þ1 ¼ u2 ðt n Þ; 1 6 n 6 N:
ð9Þ
M.R. Cui / Applied Mathematics and Computation 246 (2014) 464–473
467
For every time level t n , we have obtained 2N x equations, but the total number of unknowns is 2N x þ 2, that is, we need to seek the values of the variables Pnj ð1 6 j 6 N x Þ and U nj ð0 6 j 6 N x þ 1Þ. Hence we establish two additional boundary relations for the flux. These equations can be given as in paper [7], but we shift the subscript to obtain a tridiagonal coefficient matrix.
n n 1a 0 U n0 3 1a 1 U n1 ¼ 1h 17 Pn þ 32 Pn1 þ 32 Pn2 16 P n3 ; 6 0 n n n P 32 PnNx 32 PnNx 1 þ 16 P nNx 2 : 1a Nx þ1 U nNx þ1 3 1a Nx U nNx ¼ 1h 17 6 N x þ1
ð10Þ
Therefore, we can give the CCD scheme for (1)–(3) as follows. Step 1: As we have fP0 g, we first solve for fU0 g.
8 0 0 0 1 > U þ 3 1a 1 U 01 ¼ 1h ð17 P0 32 P01 32 P02 þ 16 P03 Þ; > a 0 0 6 0 > > > > < 1a 0j1 U 0j1 þ 4 1a 0j U 0j þ 1a 0jþ1 U 0jþ1 ¼ 3h ðP0j1 P 0jþ1 Þ; 1 6 j 6 Nx ; > 310 U 0 þ 10 U 0 > ¼ 1h 17 P0 þ 32 P0Nx þ 32 P0Nx 1 16 P0Nx 2 ; > a Nx Nx a Nx þ1 N x þ1 6 N x þ1 > > > : 0 Pj ¼ wðxj Þ; 0 6 j 6 Nx þ 1:
ð11Þ
Step 2: After we have obtained fPn1 ; Un1 gðn P 1Þ, we get the numerical solutions fPn ; Un gð1 6 n 6 NÞ by solving the following equations.
n 8 1n n U þ 3 1a 1 U n1 þ 1h ð32 Pn1 þ 32 Pn2 16 P n3 Þ ¼ 17 Pn ; > a 0 0 6h 0 > > > n n n > 1 > U n þ 4 1a j U nj þ 1a jþ1 U njþ1 þ 3h ðPnjþ1 Pnj1 Þ ¼ 0; 1 6 j 6 Nx ; > a j1 j1 > > > 1n n n n > n n n > Pn ; > 3 a Nx U Nx þ 1a Nx þ1 U Nx þ1 þ 1h 16 PNx 2 32 PNx 1 32 PNx ¼ 17 > 6h N x þ1 > h i > n n > n > c > s Cð2 cÞ ba j1 þ 3h U nj1 þ 4ðbaÞj U nj þ ba jþ1 3h U njþ1 > > > > < þð1 þ sc Cð2 cÞcnj1 ÞPnj1 þ 4ð1 þ sc Cð2 cÞcnj ÞPnj > > þð1 þ sc Cð2 cÞcn ÞPnjþ1 > jþ1 > > > > X n1 > > > > ðdnl1 dnl ÞðPlj1 þ 4Plj þ Pljþ1 Þ þ dn1 ðP0j1 þ 4P 0j þ P0jþ1 Þ >¼ > > > l¼1 > > > n n n > > þsc Cð2 cÞðf j1 þ 4f j þ f jþ1 Þ; 1 6 j 6 Nx ; > > > : 0 Pj ¼ wðxj Þ; 0 6 j 6 Nx þ 1; P n0 ¼ u1 ðt n Þ; PnNx þ1 ¼ u2 ðtn Þ; 1 6 n 6 N:
ð12Þ
On each time level t n we denote the exact solution vector of dimension N x by pn ¼ pðtn Þ ¼ ðpn1 ; pn2 ; . . . ; pnNx ÞT and the approxT imate solution vector Pn ¼ Pðt n Þ ¼ ðPn1 ; Pn2 ; . . . ; PnNx Þ . Similarly we have un ¼ uðt n Þ ¼ ðun0 ; un1 ; ; unNx þ1 ÞT and T Un ¼ Uðt n Þ ¼ ðU n0 ; U n1 ; . . . ; U nNx þ1 Þ , respectively. Then we can write (11) and (12) in the matrix form,
A1 U0 ¼ F0
ð13Þ
and
A1 B2
B1 A2
Un
Pn
! n1 X O O D11 Ul þ ¼ l D21 O C P l l¼0
D12 D22
Fn n
G
:
ð14Þ
Let sC ¼ sc Cð2 cÞ, the matrices A1 2 RðNx þ2ÞðNx þ2Þ ; A2 2 RNx Nx ; B1 2 RðNx þ2ÞNx ; B2 2 RNx ðNx þ2Þ , and C l 2 RNx Nx ð1 6 l 6 n 1Þ in (14) are given by
0
n
ð1aÞ0 B B 1n Bð Þ B a0 B B B 0 A1 ¼ B B . B . B . B B .. B . @ 0
n
0
n
ð1aÞ2
3ð1aÞ1 4ð1aÞ1 n
n
n
.. . n
ð1aÞ3 .. .
.. . .. . .. .
ð1aÞ1 .. . .. .
4ð1aÞ2 .. . .. .
ð1aÞNx 1
n
4ð1aÞNx
0
3ð1aÞNx
n n
0 .. . .. .
1
C C C C C C C C; C C 0 C C C n ð1aÞNx þ1 C A n ð1aÞNx þ1
468
M.R. Cui / Applied Mathematics and Computation 246 (2014) 464–473
0
3 2
3 2
16
0
0
0
0
0 C C C 0 C C C 0 C C C 0 C C; C 0 C C C C 3 C C C 0 A
B 0 3 0 0 B B B 3 0 3 0 0 B B .. B 0 3 0 3 0 . B .. 1B B 0 0 3 0 3 . B1 ¼ B hB . .. .. .. .. .. B .. . . . . . B B . B . 0 3 0 B . B . B . @ . 0 0 3
0
b1;1
B B 0 B B2 ¼ sC B B .. @ . 0 with b1;j ¼ 3h
0
bn a j1
1 6
b2;1
b3;1
0
b1;2 .. .
b2;2 .. . 0
b3;2 .. .
32
b1;Nx 1
ð1 6 j 6 N x 1Þ; b2;j
4ð1 þ sC cn1 Þ
B B B 1 þ sC cn1 B B A2 ¼ B 0 B B .. B . @ 0
1 þ sC cn2
1
32
.. . .. . b2;Nx
0 0 0
C C C C; C A
b3;Nx þ1
n n ¼ 4 ba j ð1 6 j 6 N x Þ; b3;j ¼ 3h ba jþ1 ðfor 1 6 j 6 N x þ 1Þ.
0
4ð1 þ sC cn2 Þ 1 þ sC cn3 .. .. . . .. .. . . 0
1
.. . .. . .. . 1þs
n C c N x 1
0 .. . 0 1 þ sC cnNx 4ð1 þ sC cnNx Þ
1 C C C C C C; C C C A
and
1 0 B .. .C B1 4 1 . .. C C B C B .. .. .. C l ¼ ðdnl1 dnl ÞB . . . 0C C: B0 C B. . C B. . . 1 4 1A @. 0 0 1 4 0
4
1
0
Remark 2.1. We can use other choices of fourth order approximations for the boundary points to replace (10), such as those two equations given in paper [18].
n n 1 1 1 11 4 1 U n1 U n2 ¼ P n0 4Pn1 þ 6Pn2 Pn3 þ Pn4 ; a 1 a 2 h 12 3 4 n n 1 1 1 1 4 11 n U nNx 1 4 U nNx ¼ PnNx 3 þ PnNx 2 6PnNx 1 þ 4PnNx þ PNx þ1 : a Nx 1 a Nx h 4 3 12
4
Remark 2.2. Using (7), let Rn1 and Rn2 be defined by
" # n1 n X 2 2 2 n sc h2 2 n l 0 þ dx unj 1 þ h6 d2x ba u j þ 1 þ h6 d2x ½cpnj 1 þ h6 d2x f j ; p Rn1 ¼ Cð2 1 þ d ðd d Þp d p n1 j nl1 nl j j cÞ 6 x l¼1
2 n Rn2 ¼ dx pnj þ 1 þ h6 d2x 1a u j ; 4
1 6 j 6 Nx ; 1 < n 6 N: 4
4
We can see that jRn1 j 6 Cðs2c þ h Þ, jRn2 j 6 Ch , hence the local truncation error of the scheme (11) and (12) is Oðs2c Þ þ Oðh Þ.
M.R. Cui / Applied Mathematics and Computation 246 (2014) 464–473
469
3. CCD scheme for Robin boundary condition We need to add P n0 and PnNx þ1 as unknowns now. That is, we seek solution vectors Pn ; Un 2 RNx þ2 now. Using (4), we have
n 1 U n þ an1;0 Pn0 ¼ bn1;0 ; a 0 0
n 1 Un þ an2;Nx þ1 PnNx þ1 ¼ bn2;Nx þ1 : a Nx þ1 Nx þ1
The CCD scheme now becomes the following two steps. Step 1 : We first solve for U 0j ð0 6 j 6 N x þ 1Þ and P00 ; P0Nx þ1 .
8 1n 0 U þ a01;0 P 00 ¼ b01;0 ; > > > a 0 0 > 10 0 17 0 1 3 0 3 0 1 0 > 1 0 0 > > > a 0 U 0 þ 3 a 1 U 1 6h P0 ¼ h 2 P 1 2 P2 þ 6 P3 ; > > > > > 10 U 0 þ 410 U 0 þ 10 U 0 3 P0 ¼ 3 P 0 ; > > a 1 1 a 2 2 h 0 h 2 > < a 0 0 10 0 10 0 0 1 0 3 U þ 4 a j U j þ a jþ1 U jþ1 ¼ h ðP0j1 P 0jþ1 Þ; 2 6 j 6 Nx 1; a j1 j1 > > 10 0 10 > 0 0 > 1 0 3 0 3 0 > > > a Nx 1 U Nx 1 þ 4 a Nx U Nx þ a Nx þ1 U Nx þ1þ h PNx þ1 ¼ h PNx 1 ; > > > 0 0 0 > > P0 ¼ 1h 32 P0Nx þ 32 P0Nx 1 16 P0Nx 2 ; 1a Nx þ1 U 0Nx þ1 þ a02;Nx þ1 PnNx þ1 ¼ b02;Nx þ1 ; 3 1a Nx U 0Nx þ 1a Nx þ1 U 0Nx þ1 þ 17 > 6h Nx þ1 > > > : 0 Pj ¼ wðxj Þ; 1 6 j 6 Nx :
ð15Þ
Step 2 : We get the numerical solutions fPn ; Un gð1 6 n 6 NÞ by solving the following equations.
8 1n n U þ an1;0 P n0 ¼ bn1;0 ; > > > a0n 0 1n n 1 17 n 3 n 3 n 1 n > > n 1 > > > a 0 U 0 þ 3 a 1 U 1 þ h 6 P0 þ 2 P1 þ 2 P2 6 P3 ¼ 0; > n n n n n n > n n > > > 1a j1 U j1 þ 4 1a j U j þ 1a jþ1 U jþ1 þ 3h ðPjþ1 Pj1 Þ ¼ 0; > > h i > > > sc Cð2 cÞ bn þ 3 U n þ 4bn U n þ bn 3 U n > > j1 j jþ1 a h a a h j1 j jþ1 > > > > > > þð1 þ sc Cð2 cÞcnj1 ÞPnj1 þ 4ð1 þ sc Cð2 cÞcnj ÞPnj > > > > n c n > > < þð1 þ s Cð2 cÞcjþ1 ÞPjþ1 n1 X > > ¼ ðdnl1 dnl ÞðPlj1 þ 4Plj þ Pljþ1 Þ þ dn1 ðP0j1 þ 4P 0j þ P0jþ1 Þ > > > > l¼1 > > > c n n n > þs Cð2 cÞðf j1 þ 4f j þ f jþ1 Þ; > > > > n n n > > > 3 1 U þ 1a Nx þ1 U nNx þ1 > > a Nx Nx > > 1 1 n 3 n 3 n 17 n > > > þ h 6 PNx 2 2 PNx 1 2 PNx þ 6 PNx þ1 ¼ 0; > > n > n n n > > 1a Nx þ1 U Nx þ1 þ an2;Nx þ1 PNx þ1 ¼ b2;Nx þ1 ; > > > : 0 Pj ¼ wðxj Þ; for 1 6 j 6 Nx and 1 6 n 6 N:
ð16Þ
Remark 3.1. As we substitute the derivatives appeared on the boundary conditions for the corresponding fluxes directly, and the equations on the interior points are unchanged, thus the local truncation error of the scheme (15) and (16) is also 4 Oðs2c Þ þ Oðh Þ. 4. Convergence analysis for the integer order problem with constant coefficients under some assumption In this section, we give the convergence analysis for the problem of integer order with constant coefficients under the assumption that both the errors of the primitive unknown and its derivative are zero on the boundary. To give the convergence analysis, we need to introduce the inner product and the induced norm. Let
ðf ; gÞ ¼ h
Nx X
f j gj ;
1=2
kf k ¼ ðf ; gÞ
;
j¼1
and we have the following theorem. Theorem 4.1. We assume that c ¼ 1 and the solution of equation (1) is sufficiently smooth, the coefficients appeared in equation (1) are constants, and the errors of both the primitive variable and its derivative are zero on the boundary. Then the convergence order (1) of the compact difference scheme defined by (11) and (12) is given by 4
kpn Pn k þ kun Un k 6 Cðs þ h Þ;
n ¼ 1; 2; . . . ; N:
ð17Þ
470
M.R. Cui / Applied Mathematics and Computation 246 (2014) 464–473
Proof. Denote the errors by nn ¼ pn Pn and fn ¼ un Un , then from (1) and (9) (with c ¼ 1 now), we get
! n ! n ! n1 2 2 2 h 2 nj nj h 2 b h 2 n 1 þ dx þ dx fj 1 þ dx f þ 1 þ dx ½cnnj ¼ Rn1;j ; a j 6 s 6 6 2
dx nnj
h þ 1 þ d2x 6
!
n 1 f ¼ Rn2;j ; a j
4
ð18Þ
1 6 j 6 Nx ; 1 < n 6 N;
ð19Þ
4
where jRn1;j j 6 Cðs þ h Þ; jRn2;j j 6 Ch , for 1 6 j 6 N x ; 1 < n 6 N.
For any real vector functions f ; g defined on the mesh points xj ð1 6 j 6 N x Þ, let A be the matrix so that 2 T 1 þ h6 d2x f ; g ¼ f Ag holds, it follows that
0
4
B B1 B B B 1B B0 A¼ B. 6B . B. B B. B .. @
1
0
4
1
1 .. . .. .
4 .. . .. .
1 0 .. .. .C . . .. C C C .. .. C . .C 1 C C; .. .. C . . 0C C C 1 4 1C A
0
0
1
4
which is symmetric, positive definite, and in fact, we have
! ! 2 h 2 1 þ dx f ; f 6 kf k2 : 6
1 kf k2 6 3
ð20Þ
Hence we can introduce a weighted norm kf kA ¼ T
1=2
ðf Af Þ
ðg T AgÞ
1þ
1=2
6 12 ðkf k2A þ kgk2A Þ. Now we can multiply (18) by hnnj , (19) by hfnj , and sum for 1 6 j 6 N x to get
! ! 2 h 2 nn nn1 n dx ; n þ ðdx fn ; nn Þ 6 s 2
n
1=2 2 2 T 1 þ h6 d2x f ; f . Besides, we have 1 þ h6 d2x f ; g ¼ jf Agj 6
n
ðdx n ; f Þ þ
h 1 þ d2x 6
!
2
1þ
n ! 1 f ; fn ¼ ðRn2 ; fn Þ; a
h 2 d 6 x
! n ! b f ; nn þ a
1þ
! ! 2 h 2 dx ½cnn ; nn ¼ ðRn1 ; nn Þ; 6
1 < n 6 N:
ð21Þ
ð22Þ
As we assume that nn0 ¼ nnNx þ1 ¼ fn0 ¼ fnNx þ1 ¼ 0 holds true, therefore,
ðdx nn ; fn Þ ¼ h
Nx X nnjþ1 nnj1
2h
j¼1
fnj ¼
Nx x þ1 x 1 1 NX 1 NX 1X 1 nnj fnj1 nnj fnjþ1 ¼ nn ðfn fnjþ1 Þ þ ðnnNx þ1 fnNx þ nnNx fnNx þ1 nn0 fn1 nn1 fn0 Þ 2 j¼2 2 j¼0 2 j¼1 j j1 2
n
¼ ðdx f ; nn Þ: Hence, we can add (21) and (22) together to obtain
1þ
! ! 2 h 2 nn nn1 n dx ;n þ 6 s
¼ ðRn1 ; nn Þ þ ðRn2 ; fn Þ; thus
2
1þ
h 2 d 6 x
! ! n 1 f ; fn þ a
1þ
! ! 2 h 2 dx ½cnn ; nn 6
2
1þ
h 2 d 6 x
!
n ! b f ; nn a
1 < n 6 N;
! n ! 2 1 1 h 2 b 1 1 1 2 2 n 2 n1 2 ðkn kA kn kA Þ þ kfkA þ cknkA 1 þ dx f ; nn 6 ðknn k2 þ kfn k2 Þ þ 3 kRn1 k2 þ akRn2 k2 2s a a 3 4a 4 6 1 1 6 knn k2A þ kfn k2A þ 3 kRn1 k2 þ akRn2 k2 ; 1 < n 6 N: 4a 4
Note that
2
h 1 þ d2x 6
! ! n b b n f ; n ¼ jj a a
! ! 2 2 h 2 n n 1 n 2 b 1 þ dx f ; n 6 kf kA þ knn k2A ; 4a 6 a
471
M.R. Cui / Applied Mathematics and Computation 246 (2014) 464–473
consequently, we get 2
1þ2 c1
b a
! !
1 a
s knn k2A þ skfk2A 6 knn1 k2A þ 6s
1 n 2 kR1 k þ akRn2 k2 : 4
ð23Þ
2
2
Based on the sign of c 1 ba , we divide the discussion of (23) in the following two cases. If c P 1 þ ba , then we have
knn k2A þ
n X 1 1 1 l 2 4 2 skfk2A 6 knn1 k2A þ 6s kRn1 k2 þ akRn2 k2 6 6s kR1 k þ akRl2 k2 6 Cðs þ h Þ : a 4 4 l¼0 2
Otherwise, let c0 ¼ j2ðc 1 ba Þj, for
16
1 2
1 þ 2ðc 1 ba Þs
¼
2
s small enough such that both 1 þ 2ðc 1 ba Þs > 0 and
1 6 1 þ 2c0 s 6 2 1 c0 s
hold, note that n0j ¼ 0ð0 6 j 6 N x þ 1Þ, we get
knn k2A þ
n X 2 1 1 1 skfk2A 6 ð1 þ 2c0 sÞknn1 k2A þ 12s kRn1 k2 þ akRn2 k2 6 12s ð1 þ 2c0 sÞnl kRl1 k2 þ akRl2 k2 6 Cðs þ h4 Þ : a 4 4 l¼1
Thus we can obtain
knn k2A þ
2 1 skfk2A 6 Cðs þ h4 Þ a
ð24Þ
again, hence we can prove (17) by using (20) and (24). The proof of the theorem is completed.
h
Remark 4.1. In the proof we have used the assumption nn0 ¼ nnNx þ1 ¼ fn0 ¼ fnNx þ1 ¼ 0. For problems with Dirichlet boundary conditions, we have nn0 ¼ nnNx þ1 ¼ 0, likewise we have fn0 ¼ fnNx þ1 ¼ 0 for Neumann boundary problems. However, we cannot give a condition to make all these conditions hold true at present. Remark 4.2. For the partial differential equations of fractional order c 2 ð0; 1Þ, in the estimate we will have a factor ð1 C sc Þn multiplied on the righthand side, but this term tends to infinity when c < 1, so we only discuss c ¼ 1 here. 5. Numerical experiments We give some computational results in this section, to test the convergence order of our scheme. In the numerical simulations, we use the discrete l2 -error norm which include errors on the boundary points for both test Problems 1 and 2. Table 1 Results for test Problem 1 with c ¼ 0:2; 0:4; 0:6; 0:8 and 1:0.
c
h
s
N
keN pk
order
keN uk
order
0.2 0.2 0.2 0.2
0.1 0.05 0.025 0.0125
0.1 0.0214 0.0046 9.8431e-4
10 47 218 1016
3.6451e-5 2.3881e-6 1.5520e-7 9.9066e-9
– 3.9320 3.9437 3.9696
2.2814e-4 1.3622e-5 8.4283e-7 5.2457e-8
– 4.0659 4.0146 4.0060
0.4 0.4 0.4 0.4
0.1 0.05 0.025 0.0125
0.1 0.0177 0.0031 5.5243e-4
10 57 320 1810
2.2911e-4 1.4357e-5 9.0508e-7 5.6675e-8
3.9962 3.9876 3.9973
0.0014 8.1324e-5 4.8879e-6 2.9865e-7
– 4.1056 4.0564 4.0327
0.6 0.6 0.6 0.6
0.1 0.05 0.025 0.0125
0.1 0.0138 0.0019 2.6287e-4
10 72 525 3804
9.0349e-4 5.6485e-5 3.5221e-6 2.2017e-7
– 3.9996 4.0034 3.9997
0.0056 3.1955e-4 1.9003e-5 1.1593e-6
– 4.1313 4.0717 4.0349
0.8 0.8 0.8 0.8
0.1 0.05 0.025 0.0125
0.1 0.0099 9.8431e-4 9.7656e-5
10 101 1016 10240
0.0029 1.8331e-4 1.1459e-5 7.1620e-7
– 3.9837 3.9997 4.0000
0.0182 0.0010 6.1813e-5 3.7703e-6
– 4.1859 4.0159 4.0352
1.0 1.0 1.0 1.0
0.1 0.05 0.025 0.0125
0.1 0.0063 3.9063e-4 2.4414e-5
10 160 2560 40960
0.0085 5.3033e-4 3.3146e-5 2.0716e-6
– 4.0025 4.0000 4.0000
0.0524 0.0030 1.7877e-4 1.0904e-5
– 4.1265 4.0688 4.0352
472
M.R. Cui / Applied Mathematics and Computation 246 (2014) 464–473
Table 2 Results for test Problem 2 with c ¼ 0:2; 0:4; 0:6; 0:8 and 1:0.
c
h
s
N
keN pk
order
keN uk
order
0.2 0.2 0.2 0.2
0.1 0.05 0.025 0.0125
0.1 0.0214 0.0046 9.8431e-4
10 47 218 1016
0.0082 5.7642e-4 3.9408e-5 2.6200e-6
– 3.8304 3.8706 3.9108
0.0053 3.7669e-4 2.5776e-5 1.7144e-6
– 3.8145 3.8693 3.9103
0.4 0.4 0.4 0.4
0.1 0.05 0.025 0.0125
0.1 0.0177 0.0031 5.5243e-4
10 57 320 1810
0.0554 0.0037 2.4204e-4 1.5378e-5
– 3.9043 3.9342 3.9763
0.0361 0.0024 1.5831e-4 1.0063e-5
– 3.9109 3.9222 3.9756
0.6 0.6 0.6 0.6
0.1 0.05 0.025 0.0125
0.1 0.0138 0.0019 2.6287e-4
10 72 525 3804
0.2006 0.0131 8.2802e-4 5.1725e-5
– 3.9367 3.9838 4.0007
0.1309 0.0086 5.4158e-4 3.3847e-5
– 3.9280 3.9891 4.0001
0.8 0.8 0.8 0.8
0.1 0.05 0.025 0.0125
0.1 0.0099 9.8431e-4 9.7656e-5
10 101 1016 10240
0.5495 0.0351 0.0022 1.3538e-4
– 3.9686 3.9959 4.0224
0.3584 0.0229 0.0014 8.8586e-5
– 3.9682 4.0318 3.9822
1.0 1.0 1.0 1.0
0.1 0.05 0.025 0.0125
0.1 0.0063 3.9063e-4 2.4414e-5
10 160 2560 40960
1.2862 0.0784 0.0048 3.0034e-4
– 4.0361 4.0297 3.9984
0.8390 0.0512 0.0032 1.9653e-4
– 4.0345 4.0000 4.0253
" keNp k ¼ kpN PN k ¼
N x þ1 X
" keNu k
N
N
¼ ku U k ¼
2
#1=2
ðpNj PNj Þ h
j¼0
N x þ1 X
; #1=2
ðuNj
2 U Nj Þ h
:
j¼0 4
As in paper [12], since the local truncation error of the schemes are proved to be Oðs2c Þ þ Oðh Þ, then if we decrease the mesh size of h to h=2 and s to s=24=ð2cÞ , we expect to get a factor 1/16 for the error kei s; h; tN k with index i indicates p or u. If we decrease h to h=4, then we multiply s by s=28=ð2cÞ . The numbers N in Table 1 and Table 2 are calculated by 1s , that is, we use the round number of 1s, hence N s – 1 but N s 1. The experimental order of convergence r i ðs; hÞ is computed by the formula
1 kei 24=ð2cÞ s; 2h; tN k A; ri ðs; hÞ ¼ log2 @ kei s; h; t N k 0
i ¼ p or u:
5.1. Test problem 1 Consider the Dirichlet problem (1)–(3) with the coefficients aðxÞ ¼ 1 þ x; bðxÞ ¼ 1 þ 2x; cðxÞ ¼ 1 þ 5x, the theoretical solution is pðx; tÞ ¼ ex t1þc , and the right hand side function f ðxÞ ¼ tex ðCð2 þ cÞ þ 6xt c Þ. The errors for c ¼ 0:2; 0:4; 0:5; 0:6; 0:8 and 1.0, respectively, together with the experimental convergence order are shown in Table 1. We see that the order is approximately 4, that is, the performance meets our expectations. 5.2. Test problem 2 Inspired by the example in paper [28], we use the Caputo derivative now and consider the Robin problem (1), (2) and (4) for the following equation C c 0 Dt p
¼
@ 1 @2p p þ 2 þ f ðx; tÞ; @x 1 þ x @x
0 < x < 1; 0 < t 6 1;
1þc 1 1 with bðx; tÞ ¼ 1þx ; cðx; tÞ ¼ ð1þxÞ . 2 and f ðx; tÞ ¼ 8ðx þ 1ÞCð2 þ cÞt now, and the exact solution is pðx; tÞ ¼ 8ðx þ 1Þt
For this test problem, the errors for p and u with c ¼ 0:2; 0:4; . . . ; 1:0, together with the experimental convergence order are shown in Table 2. The results show that the order is also approximately 4, that is, the l2 norm of both keNp k and keNu k is 4
Oðs2c Þ þ Oðh Þ.
M.R. Cui / Applied Mathematics and Computation 246 (2014) 464–473
473
6. Conclusion In this paper we propose a CCD scheme for solving the one dimensional time fractional convection–diffusion–reaction equation with variable coefficients. Both Dirichlet and Robin boundary conditions are considered. We solve both the unknown and its gradient simultaneously. The proposed algorithm has high accuracy and good efficiency. The local trunca4 tion error of the CCD scheme is Oðs2c Þ þ Oðh Þ. We give the convergence analysis for the problem of integer order with constant coefficients under some assumption.The conclusion is verified by some numerical experiments. Acknowledgement The author thanks the referees for their valuable comments and suggestions, which help to improve this manuscript significantly. References [1] R. Metzler, J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Phys. Rep. 339 (2000) 1–77. [2] I. Podlubny, Fractional Differential Equations, An Introduction to Fractional Derivatives, Fractional Differential Equations, Some Methods of Their Solution and Some of their Applications, Academic Press, San Diego, 1999. [3] C.-M. Chen, F. Liu, I. Turner, V. Anh, A Fourier method for the fractional diffusion equation describing sub-diffusion, J. Comput. phys. 227 (2007) 886– 897. [4] P. Zhuang, F. Liu, V. Anh, I. Turner, New solution and analytical techniques of the implicit numerical method for the anomalous subdiffusion equation, SIAM J. Numer. Anal. 46 (2008) 1079–1095. [5] P. Zhuang, F. Liu, V. Anh, I. Turner, Numerical methods for the variable-order fractional advection-diffusion equation with a nonlinear source term, SIAM J. Numer. Anal. 47 (2009) 1760–1781. [6] R.S. Hirsch, Higher order accurate difference solutions of fluid mechanics problems by a compact difference technique, J. Comput. Phys. 24 (1975) 90– 109. [7] S.K. Lele, Compact finite difference schemes with spectral-like resolution, J. Comput. Phys. 103 (1992) 16–42. [8] M.R. Cui, Compact finite difference method for the fractional diffusion equation, J. Comput. Phys. 228 (2009) 7792–7804. [9] C.-M. Chen, F. Liu, V. Anh, I. Turner, Numerical schemes with high spatial accuracy for a variable-order anomalous subdiffusion equation, SIAM J. Sci. Comput. 32 (2010) 1740–1760. [10] G.-H. Gao, Z.-Z. Sun, A compact finite difference scheme for the fractional sub-diffusion equations, J. Comput. Phys. 230 (2011) 586–595. [11] M.R. Cui, Compact alternating direction implicit method for two-dimensional time fractional diffusion equation, J. Comput. Phys. 231 (2012) 2621– 2633. [12] M.R. Cui, Convergence analysis of high-order compact alternating direction implicit schemes for the two-dimensional time fractional diffusion equation, Numer. Algor. 62 (2013) 383–409. [13] M.R. Cui, A high-order compact exponential scheme for the fractional convection-diffusion equation, J. Comput. Appl. Math. 255 (2014) 404–416. [14] G.F. Carey, W.F. Spotz, Higher-order compact mixed methods, Commun. Numer. Methods Eng. 13 (1997) 553–564. [15] W.F. Spotz, G.F. Carey, High-order compact scheme for the steady stream-function vorticity equations, Int. J. Numer. Methods Eng. 38 (1995) 3497– 3512. [16] W.F. Spotz, G.F. Carey, Extension of high-order compact schemes to time-dependent problems, Numer. Methods Partial Differ. Eqs. 17 (2001) 657–672. [17] S. Abide, X. Chesneau, B. Zeghmati, Compact mixed methods for convection/diffusion problems, Appl. Math. Comput. 218 (2012) 5867–5876. [18] J. Zhao, Highly accurate compact mixed methods for two point boundary value problems, Appl. Math. Comput. 188 (2007) 1402–1408. [19] F.C. Chu, C. Fan, A three-point combined compact difference scheme, J. Comput. Phys. 140 (1998) 370–399. [20] F.C. Chu, C. Fan, A three-point sixth-order nonuniform combined compact difference scheme, J. Comput. Phys. 148 (1999) 663–674. [21] K. Mahesh, A family of high order finite difference schemes with good spectral resolution, J. Comput. Phys. 145 (1998) 332–358. [22] T. Nihei, K. Ishii, A fast solver of the shallow water equations on a sphere using a combined compact difference scheme, J. Comput. Phys. 187 (2003) 639–659. [23] J. Zhang, J.J. Zhao, Truncation error and oscillation property of the combined compact difference scheme, Appl. Math. Comput. 161 (2005) 241–251. [24] Q. Zhou, Z.H. Yao, F. He, M.Y. Shen, A new family of high-order compact upwind difference schemes with good spectral resolution, J. Comput. Phys. 227 (2007) 1306–1339. [25] T.K. Sengupta, V. Lakshmanan, V.V.S.N. Vijay, A new combined stable and dispersion relation preserving compact scheme for non-periodic problems, J. Comput. Phys. 228 (2009) 3048–3071. [26] T.K. Sengupta, V.V.S.N. Vijay, S. Bhaumik, Further improvement and analysis of CCD scheme: dissipation discretization and de-aliasing properties, J. Comput. Phys. 228 (2009) 6150–6168. [27] W.J. Chen, J.C. Chen, Combined compact difference method for solving the incompressible Navier–Stokes equations, Int. J. Numer. Methods Fluids 68 (2012) 1234–1256. [28] S. Chen, F. Liu, P. Zhuang, V. Anh, Finite difference approximations for the fractional Fokker–Planck equation, Appl. Math. Model. 33 (2009) 256–273. [29] B.I. Henry, T.A.M. Langlands, P. Straka, Fractional Fokker–Planck equations for subdiffusion with space- and time-dependent forces, Phys. Rev. Lett. 105 (2010) 170602. [30] Z.Z. Sun, X.N. Wu, A fully discrete difference scheme for a diffusion-wave system, Appl. Numer. Math. 56 (2006) 193–209.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Combined compact difference scheme for the time fractional convection–diffusion equation with variable coefficients q Mingrong Cui School of Mathematics, Shandong University, Jinan, Shandong 250100, China
a r t i c l e
i n f o
Keywords: Compact Combined Fractional differential equation Convection–diffusion Mixed Difference High-order
a b s t r a c t Fourth-order combined compact finite difference scheme is given for solving the time fractional convection–diffusion–reaction equation with variable coefficients. We introduce the flux as a new variable and transform the original equation into a system of two equations. Compact difference is used as a high-order approximation for spatial derivatives of integer order in the coupled partial differential equations. The Caputo fractional derivative is discretized by a high-order approximation. Both Dirichlet and Robin boundary conditions are discussed. Convergence analysis is given for the problem of integer order with constant coefficients under some assumption. Numerical results are provided to verify the accuracy and efficiency of the proposed algorithm. Ó 2014 Elsevier Inc. All rights reserved.
1. Introduction Fractional differential equations(FDEs) can be used to in mathematical models for physical, biological, geological and financial systems, and the review article [1] and monograph [2] have given detailed discussions on fractional differential equations. In the recent years, many numerical methods have been proposed and studied for solving FDEs up to now, e.g., [3–5]. Compact finite difference schemes have high-order of accuracy and the desirable tridiagonal nature of the finite-difference equations (see [6,7]), and one-dimensional fractional sub-diffusion equation was recently solved by the compact finite dif4 ference scheme with convergence order Oðs þ h Þ in [8,9], a higher order one for the temporal variable in [10], and for two dimensional problems, ADI schemes with error analysis were given in papers [11,12], with the compact schemes for the time fractional convection–diffusion equations with constant coefficients considered in [13]. While high order methods for partial different equations with constant coefficients have been discussed by many authors, it is not so easy to give corresponding results for variable coefficients. One way to give high order difference schemes is to start from lower order scheme, then follow a substitution of the lower order truncation error by replacing the derivatives using the successive derivation of the original continuous problem. Thus high order schemes for variable coefficients must be examined carefully as many derivatives appear, and our aim in this paper is to give some easily implemented schemes for problems with variable coefficients. We must admit for problems with variable coefficients, following this way would be a tedious task. As pointed in [14], the idea of posing a boundary-value problem in the form of a first-order system as a governing conservation equation and an associated constitutive equation, and then using a mixed approximation scheme, has been
q This work was partially supported by Natural Science Foundation of Shandong Province, Grant No. BS2010HZ012, and National Natural Science Foundation of China, Grant No. 11171189. E-mail address: [email protected]
http://dx.doi.org/10.1016/j.amc.2014.08.025 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.
M.R. Cui / Applied Mathematics and Computation 246 (2014) 464–473
465
receiving increasing attention. For the steady and unsteady convection–diffusion equations of integer order, many research papers, e.g., [15,16], and recently a generalization work in [17], discussed solving them by using the compact difference scheme combined with mixed method. For two point boundary problem, this method was also discussed in paper [18]. After introducing the first derivative of the unknown solution as a new variable, then we can establish an approximation system between the derivative and the original solution. These schemes are also called combined compact difference(CCD) schemes, see papers [19,20]. The CCD schemes have significantly better spectral representation [21], and they have been proved very useful solving the shallow water equations in spherical geometry [22], the one dimensional convection diffusion equation [23], high-order compact upwind difference schemes has been developed in [24], and the two-dimensional lid-driven cavity problem [25,26], and Navier–Stokes equations [27]. In this paper we solve the time fractional convection–diffusion–reaction problems using the CCD scheme. The problems considered here are convection–diffusion–reaction problems with variable coefficients, thus they can be used for solving the fractional Fokker–Planck equation [28,29]. We give the CCD schemes and the variable coefficients can be handled easily. Besides, both Dirichlet and Robin conditions are discussed, and as we use the flux as a new variable, the solution and the flux can be obtained simultaneously without losing accuracy. We consider the following model equation C c 0 Dt p
þ bðx; tÞ
@p @ @p þ cðx; tÞp ¼ f ðx; tÞ; aðx; tÞ @x @x @x
L1 < x < L2 ; 0 < t 6 T
ð1Þ
where the Caputo fractional derivative C0 Dct v ð0 < c < 1Þ of the function v ðx; tÞ is defined by [2] C c 0 Dt
1 v ðx; tÞ ¼ Cð1 cÞ
Rt 0
ðt sÞc @@vs ðx; sÞds:
The initial condition for (1) is
pðx; 0Þ ¼ wðxÞ;
L1 < x < L2
ð2Þ
with the Dirichlet boundary condition given by
pðL1 ; tÞ ¼ u1 ðtÞ;
pðL2 ; tÞ ¼ u2 ðtÞ; 0 < t 6 T
ð3Þ
or with the Robin boundary condition
@p þ a1 ðx; tÞp jx¼L1 ¼ b1 ðtÞ; @x @p þ a2 ðx; tÞp jx¼L2 ¼ b2 ðtÞ; 0 < t 6 T: @x
ð4Þ
The paper is organized as follows. In Section 2, after introducing the flux as a new variable and then changing the original problem into two coupled equations, we use the high-order compact difference to approximate the spatial derivative and then adopt the high order discretization for the time fractional derivative. Then we get an implicit CCD scheme for the time fractional convection–diffusion equation with the local truncation error discussed. Treatment of the Robin boundary is discussed in Section 3. We give the convergence analysis for the integer order problem with constant coefficients under some assumption in Section 4. Finally, some numerical examples are given in Section 5 to verify the order of convergence of the schemes. This paper closes with a summary in Section 6. Throughout this paper, the symbol C is a generic positive constant, it may take different values at different places. We use P the ‘‘empty sum’’ convention ql¼p Sl ¼ 0 for q < p. 2. CCD scheme for Dirichlet boundary condition 2.1. Partition of the domain For the numerical solution of Eq. (1) we introduce a uniform grid of mesh points ðxj ; t n Þ, with xj ¼ L1 þ jh; j ¼ 0; 1; . . . ; N x þ 1; and t n ¼ ns; n ¼ 0; 1; . . . ; N. Here both N x and N are positive integers, h ¼ ðL2 L1 Þ=ðN x þ 1Þ is the mesh-width in x, and s ¼ T=N is the time step. For any function v ðx; tÞ, we let v nj ¼ v ðxj ; tn Þ, e.g., the theoretical solution u (i.e., the flux, see below) at the mesh point ðxj ; t n Þ is denoted by unj , and U nj stands for the numerical solution of the flux at the same mesh point. 2.2. Derivation of the CCD scheme In order to give the CCD scheme for the model equation (1), we introduce the flux uðx; tÞ ¼ aðx; tÞ @p as a new variable. @x The unknowns p and u may correspond to pressure and velocity in the flow equations in the mathematical modeling of compressible/incompressible miscible displacement problems, respectively. Then Eq. (1) becomes
(
C c 0 Dt p
þ @u bðx;tÞ u þ cðx; tÞp ¼ f ðx; tÞ; @x aðx;tÞ
¼ 0; u þ aðx; tÞ @p @x
L1 < x < L2 ; 0 < t 6 T:
We use high-order compact approximations for the first order derivatives in (5), i. e.,
ð5Þ
466
M.R. Cui / Applied Mathematics and Computation 246 (2014) 464–473
8 c > < C0 Dt pnj þ > : a1n unj þ j
bn
dx
2
1þh6 d2x
n
j
4
dx
2
1þh6 d2x
4
unj anj unj þ cnj pnj ¼ f j þ Oðh Þ;
pnj ¼ Oðh Þ;
ð6Þ
1 6 j 6 Nx ; 1 < n 6 N
where dx ; d2x denote the approximations of the first and second order derivatives, respectively. That is, for a mesh function v nj ,
dx v nj ¼
v njþ1 v nj1 2h
;
d2x v nj ¼
v nj1 2v nj þv njþ1 h2
:
We need to notice that the subscript j runs from 1 to N x , that is, we try to find numerical solutions for the unknowns pnj ð1 6 j 6 N x Þ and unj ð0 6 j 6 N x þ 1Þ. 2 Multiply the operator 1 þ h6 d2x on both sides of Eq. (6), we obtain
8 C n 2 2 2 > > 1 þ h6 d2x Dct pnj þ dx unj 1 þ h6 d2x ba u j þ 1 þ h6 d2x ½cpnj > > > 0 < 2 n 4 ¼ 1 þ h6 d2x f j þ Oðh Þ; 1 6 j 6 Nx ; 1 < n 6 N; > > > > > : dx pn þ 1 þ h2 d2x 1 u n ¼ Oðh4 Þ; 1 6 j 6 Nx ; 1 < n 6 N: j 6 a j
ð7Þ
The Caputo fractional derivative C0 Dct p appeared in (7) needs to be discretized now, and a high-order approximation can be obtained by the following lemma. Lemma 2.1 [30]. Suppose f ðtÞ 2 C 2 ½0; t k . Let
1 Rðf ðt k ÞÞ ¼ Cð1 cÞ
Z
tk
0
" # 0 k1 X f ðsÞ sc ds a f ðt Þ ðdkj1 dkj Þf ðt j Þ dk1 f ðt 0 Þ ; Cð2 cÞ 0 k ðt k sÞc j¼1
then
jRðf ðt k ÞÞj 6
" # 1 1 c 22c 00 þ ð1 þ 2c Þ max jf ðtÞjs2c ; 06t6t k Cð2 cÞ 12 2c 1c
1c
where c 2 ð0; 1Þ and aj ¼ ðj þ 1Þ j . c Using Lemma 2.1 replace the fractional derivative C0 Dt p by the following equality, C c 0 Dt pðxj ; t n Þ
¼
s c Cð2 cÞ
"
a0 pnj
n1 X
#
ðdnl1
dnl Þplj
dn1 p0j
þ Oðs2c Þ:
ð8Þ
l¼1
Using (8) and multiply both sides of (7) by
sc Cð2 cÞ, then an approximation of (1)–(3) can be given by
8 n1 X > 2 2 > > 1 þ h2 d2x Pnj ðdnl1 dnl Þ 1 þ h6 d2x Plj dn1 1 þ h6 d2x P0j > 6 > > > l¼1 > > > > c Cð2 cÞfd U n 1 þ h2 d2 b U n þ 1 þ h2 d2 ½cPn g > > þ s x < j j 6 x 6 x a j n c h2 2 > ¼ s Cð2 cÞ 1 þ 6 dx f j ; 1 6 j 6 Nx ; 1 < n 6 N; > > > > > 2 n > > dx Pnj 1 þ h6 d2x 1a U j ¼ 0; 1 6 j 6 Nx ; 1 < n 6 N; > > > > > : 0 Pj ¼ wðxj Þ; 0 6 j 6 Nx þ 1; Pn0 ¼ u1 ðtn Þ; PnNx þ1 ¼ u2 ðt n Þ; 1 6 n 6 N: With the definitions of the discrete operators, we can write (9) on the nodal points explicitly.
8 h n b n b n i n n 3 c b c > > > h s Cð2 cÞðU jþ1 U j1 Þ s Cð2 cÞ a U j1 þ 4 a U j þ a U jþ1 > > h i > > > þsc Cð2 cÞ ðcPÞn þ 4ðcPÞn þ ðcPÞn þ ðPn þ 4Pn þ Pn Þ > > j1 j jþ1 j1 j jþ1 > > > > n1 > X > < ðdnl1 dnl ÞðPlj1 þ 4Plj þ Pljþ1 Þ dn1 ðP0j1 þ 4P 0j þ P0jþ1 Þ l¼1 > > > n n n > > ¼ sc Cð2 cÞðf j1 þ 4f j þ f jþ1 Þ; 1 6 j 6 Nx ; 1 < n 6 N; > > > h i > 1 n n n >3 n n 1 1 > ðP P Þ þ U þ 4 U þ ð UÞ ¼ 0; 1 6 j 6 Nx ; 1 < n 6 N; > jþ1 j1 > h a a a j1 j jþ1 > > > : 0 n n Pj ¼ wðxj Þ; 0 6 j 6 Nx þ 1; P0 ¼ u1 ðtn Þ; PNx þ1 ¼ u2 ðt n Þ; 1 6 n 6 N:
ð9Þ
M.R. Cui / Applied Mathematics and Computation 246 (2014) 464–473
467
For every time level t n , we have obtained 2N x equations, but the total number of unknowns is 2N x þ 2, that is, we need to seek the values of the variables Pnj ð1 6 j 6 N x Þ and U nj ð0 6 j 6 N x þ 1Þ. Hence we establish two additional boundary relations for the flux. These equations can be given as in paper [7], but we shift the subscript to obtain a tridiagonal coefficient matrix.
n n 1a 0 U n0 3 1a 1 U n1 ¼ 1h 17 Pn þ 32 Pn1 þ 32 Pn2 16 P n3 ; 6 0 n n n P 32 PnNx 32 PnNx 1 þ 16 P nNx 2 : 1a Nx þ1 U nNx þ1 3 1a Nx U nNx ¼ 1h 17 6 N x þ1
ð10Þ
Therefore, we can give the CCD scheme for (1)–(3) as follows. Step 1: As we have fP0 g, we first solve for fU0 g.
8 0 0 0 1 > U þ 3 1a 1 U 01 ¼ 1h ð17 P0 32 P01 32 P02 þ 16 P03 Þ; > a 0 0 6 0 > > > > < 1a 0j1 U 0j1 þ 4 1a 0j U 0j þ 1a 0jþ1 U 0jþ1 ¼ 3h ðP0j1 P 0jþ1 Þ; 1 6 j 6 Nx ; > 310 U 0 þ 10 U 0 > ¼ 1h 17 P0 þ 32 P0Nx þ 32 P0Nx 1 16 P0Nx 2 ; > a Nx Nx a Nx þ1 N x þ1 6 N x þ1 > > > : 0 Pj ¼ wðxj Þ; 0 6 j 6 Nx þ 1:
ð11Þ
Step 2: After we have obtained fPn1 ; Un1 gðn P 1Þ, we get the numerical solutions fPn ; Un gð1 6 n 6 NÞ by solving the following equations.
n 8 1n n U þ 3 1a 1 U n1 þ 1h ð32 Pn1 þ 32 Pn2 16 P n3 Þ ¼ 17 Pn ; > a 0 0 6h 0 > > > n n n > 1 > U n þ 4 1a j U nj þ 1a jþ1 U njþ1 þ 3h ðPnjþ1 Pnj1 Þ ¼ 0; 1 6 j 6 Nx ; > a j1 j1 > > > 1n n n n > n n n > Pn ; > 3 a Nx U Nx þ 1a Nx þ1 U Nx þ1 þ 1h 16 PNx 2 32 PNx 1 32 PNx ¼ 17 > 6h N x þ1 > h i > n n > n > c > s Cð2 cÞ ba j1 þ 3h U nj1 þ 4ðbaÞj U nj þ ba jþ1 3h U njþ1 > > > > < þð1 þ sc Cð2 cÞcnj1 ÞPnj1 þ 4ð1 þ sc Cð2 cÞcnj ÞPnj > > þð1 þ sc Cð2 cÞcn ÞPnjþ1 > jþ1 > > > > X n1 > > > > ðdnl1 dnl ÞðPlj1 þ 4Plj þ Pljþ1 Þ þ dn1 ðP0j1 þ 4P 0j þ P0jþ1 Þ >¼ > > > l¼1 > > > n n n > > þsc Cð2 cÞðf j1 þ 4f j þ f jþ1 Þ; 1 6 j 6 Nx ; > > > : 0 Pj ¼ wðxj Þ; 0 6 j 6 Nx þ 1; P n0 ¼ u1 ðt n Þ; PnNx þ1 ¼ u2 ðtn Þ; 1 6 n 6 N:
ð12Þ
On each time level t n we denote the exact solution vector of dimension N x by pn ¼ pðtn Þ ¼ ðpn1 ; pn2 ; . . . ; pnNx ÞT and the approxT imate solution vector Pn ¼ Pðt n Þ ¼ ðPn1 ; Pn2 ; . . . ; PnNx Þ . Similarly we have un ¼ uðt n Þ ¼ ðun0 ; un1 ; ; unNx þ1 ÞT and T Un ¼ Uðt n Þ ¼ ðU n0 ; U n1 ; . . . ; U nNx þ1 Þ , respectively. Then we can write (11) and (12) in the matrix form,
A1 U0 ¼ F0
ð13Þ
and
A1 B2
B1 A2
Un
Pn
! n1 X O O D11 Ul þ ¼ l D21 O C P l l¼0
D12 D22
Fn n
G
:
ð14Þ
Let sC ¼ sc Cð2 cÞ, the matrices A1 2 RðNx þ2ÞðNx þ2Þ ; A2 2 RNx Nx ; B1 2 RðNx þ2ÞNx ; B2 2 RNx ðNx þ2Þ , and C l 2 RNx Nx ð1 6 l 6 n 1Þ in (14) are given by
0
n
ð1aÞ0 B B 1n Bð Þ B a0 B B B 0 A1 ¼ B B . B . B . B B .. B . @ 0
n
0
n
ð1aÞ2
3ð1aÞ1 4ð1aÞ1 n
n
n
.. . n
ð1aÞ3 .. .
.. . .. . .. .
ð1aÞ1 .. . .. .
4ð1aÞ2 .. . .. .
ð1aÞNx 1
n
4ð1aÞNx
0
3ð1aÞNx
n n
0 .. . .. .
1
C C C C C C C C; C C 0 C C C n ð1aÞNx þ1 C A n ð1aÞNx þ1
468
M.R. Cui / Applied Mathematics and Computation 246 (2014) 464–473
0
3 2
3 2
16
0
0
0
0
0 C C C 0 C C C 0 C C C 0 C C; C 0 C C C C 3 C C C 0 A
B 0 3 0 0 B B B 3 0 3 0 0 B B .. B 0 3 0 3 0 . B .. 1B B 0 0 3 0 3 . B1 ¼ B hB . .. .. .. .. .. B .. . . . . . B B . B . 0 3 0 B . B . B . @ . 0 0 3
0
b1;1
B B 0 B B2 ¼ sC B B .. @ . 0 with b1;j ¼ 3h
0
bn a j1
1 6
b2;1
b3;1
0
b1;2 .. .
b2;2 .. . 0
b3;2 .. .
32
b1;Nx 1
ð1 6 j 6 N x 1Þ; b2;j
4ð1 þ sC cn1 Þ
B B B 1 þ sC cn1 B B A2 ¼ B 0 B B .. B . @ 0
1 þ sC cn2
1
32
.. . .. . b2;Nx
0 0 0
C C C C; C A
b3;Nx þ1
n n ¼ 4 ba j ð1 6 j 6 N x Þ; b3;j ¼ 3h ba jþ1 ðfor 1 6 j 6 N x þ 1Þ.
0
4ð1 þ sC cn2 Þ 1 þ sC cn3 .. .. . . .. .. . . 0
1
.. . .. . .. . 1þs
n C c N x 1
0 .. . 0 1 þ sC cnNx 4ð1 þ sC cnNx Þ
1 C C C C C C; C C C A
and
1 0 B .. .C B1 4 1 . .. C C B C B .. .. .. C l ¼ ðdnl1 dnl ÞB . . . 0C C: B0 C B. . C B. . . 1 4 1A @. 0 0 1 4 0
4
1
0
Remark 2.1. We can use other choices of fourth order approximations for the boundary points to replace (10), such as those two equations given in paper [18].
n n 1 1 1 11 4 1 U n1 U n2 ¼ P n0 4Pn1 þ 6Pn2 Pn3 þ Pn4 ; a 1 a 2 h 12 3 4 n n 1 1 1 1 4 11 n U nNx 1 4 U nNx ¼ PnNx 3 þ PnNx 2 6PnNx 1 þ 4PnNx þ PNx þ1 : a Nx 1 a Nx h 4 3 12
4
Remark 2.2. Using (7), let Rn1 and Rn2 be defined by
" # n1 n X 2 2 2 n sc h2 2 n l 0 þ dx unj 1 þ h6 d2x ba u j þ 1 þ h6 d2x ½cpnj 1 þ h6 d2x f j ; p Rn1 ¼ Cð2 1 þ d ðd d Þp d p n1 j nl1 nl j j cÞ 6 x l¼1
2 n Rn2 ¼ dx pnj þ 1 þ h6 d2x 1a u j ; 4
1 6 j 6 Nx ; 1 < n 6 N: 4
4
We can see that jRn1 j 6 Cðs2c þ h Þ, jRn2 j 6 Ch , hence the local truncation error of the scheme (11) and (12) is Oðs2c Þ þ Oðh Þ.
M.R. Cui / Applied Mathematics and Computation 246 (2014) 464–473
469
3. CCD scheme for Robin boundary condition We need to add P n0 and PnNx þ1 as unknowns now. That is, we seek solution vectors Pn ; Un 2 RNx þ2 now. Using (4), we have
n 1 U n þ an1;0 Pn0 ¼ bn1;0 ; a 0 0
n 1 Un þ an2;Nx þ1 PnNx þ1 ¼ bn2;Nx þ1 : a Nx þ1 Nx þ1
The CCD scheme now becomes the following two steps. Step 1 : We first solve for U 0j ð0 6 j 6 N x þ 1Þ and P00 ; P0Nx þ1 .
8 1n 0 U þ a01;0 P 00 ¼ b01;0 ; > > > a 0 0 > 10 0 17 0 1 3 0 3 0 1 0 > 1 0 0 > > > a 0 U 0 þ 3 a 1 U 1 6h P0 ¼ h 2 P 1 2 P2 þ 6 P3 ; > > > > > 10 U 0 þ 410 U 0 þ 10 U 0 3 P0 ¼ 3 P 0 ; > > a 1 1 a 2 2 h 0 h 2 > < a 0 0 10 0 10 0 0 1 0 3 U þ 4 a j U j þ a jþ1 U jþ1 ¼ h ðP0j1 P 0jþ1 Þ; 2 6 j 6 Nx 1; a j1 j1 > > 10 0 10 > 0 0 > 1 0 3 0 3 0 > > > a Nx 1 U Nx 1 þ 4 a Nx U Nx þ a Nx þ1 U Nx þ1þ h PNx þ1 ¼ h PNx 1 ; > > > 0 0 0 > > P0 ¼ 1h 32 P0Nx þ 32 P0Nx 1 16 P0Nx 2 ; 1a Nx þ1 U 0Nx þ1 þ a02;Nx þ1 PnNx þ1 ¼ b02;Nx þ1 ; 3 1a Nx U 0Nx þ 1a Nx þ1 U 0Nx þ1 þ 17 > 6h Nx þ1 > > > : 0 Pj ¼ wðxj Þ; 1 6 j 6 Nx :
ð15Þ
Step 2 : We get the numerical solutions fPn ; Un gð1 6 n 6 NÞ by solving the following equations.
8 1n n U þ an1;0 P n0 ¼ bn1;0 ; > > > a0n 0 1n n 1 17 n 3 n 3 n 1 n > > n 1 > > > a 0 U 0 þ 3 a 1 U 1 þ h 6 P0 þ 2 P1 þ 2 P2 6 P3 ¼ 0; > n n n n n n > n n > > > 1a j1 U j1 þ 4 1a j U j þ 1a jþ1 U jþ1 þ 3h ðPjþ1 Pj1 Þ ¼ 0; > > h i > > > sc Cð2 cÞ bn þ 3 U n þ 4bn U n þ bn 3 U n > > j1 j jþ1 a h a a h j1 j jþ1 > > > > > > þð1 þ sc Cð2 cÞcnj1 ÞPnj1 þ 4ð1 þ sc Cð2 cÞcnj ÞPnj > > > > n c n > > < þð1 þ s Cð2 cÞcjþ1 ÞPjþ1 n1 X > > ¼ ðdnl1 dnl ÞðPlj1 þ 4Plj þ Pljþ1 Þ þ dn1 ðP0j1 þ 4P 0j þ P0jþ1 Þ > > > > l¼1 > > > c n n n > þs Cð2 cÞðf j1 þ 4f j þ f jþ1 Þ; > > > > n n n > > > 3 1 U þ 1a Nx þ1 U nNx þ1 > > a Nx Nx > > 1 1 n 3 n 3 n 17 n > > > þ h 6 PNx 2 2 PNx 1 2 PNx þ 6 PNx þ1 ¼ 0; > > n > n n n > > 1a Nx þ1 U Nx þ1 þ an2;Nx þ1 PNx þ1 ¼ b2;Nx þ1 ; > > > : 0 Pj ¼ wðxj Þ; for 1 6 j 6 Nx and 1 6 n 6 N:
ð16Þ
Remark 3.1. As we substitute the derivatives appeared on the boundary conditions for the corresponding fluxes directly, and the equations on the interior points are unchanged, thus the local truncation error of the scheme (15) and (16) is also 4 Oðs2c Þ þ Oðh Þ. 4. Convergence analysis for the integer order problem with constant coefficients under some assumption In this section, we give the convergence analysis for the problem of integer order with constant coefficients under the assumption that both the errors of the primitive unknown and its derivative are zero on the boundary. To give the convergence analysis, we need to introduce the inner product and the induced norm. Let
ðf ; gÞ ¼ h
Nx X
f j gj ;
1=2
kf k ¼ ðf ; gÞ
;
j¼1
and we have the following theorem. Theorem 4.1. We assume that c ¼ 1 and the solution of equation (1) is sufficiently smooth, the coefficients appeared in equation (1) are constants, and the errors of both the primitive variable and its derivative are zero on the boundary. Then the convergence order (1) of the compact difference scheme defined by (11) and (12) is given by 4
kpn Pn k þ kun Un k 6 Cðs þ h Þ;
n ¼ 1; 2; . . . ; N:
ð17Þ
470
M.R. Cui / Applied Mathematics and Computation 246 (2014) 464–473
Proof. Denote the errors by nn ¼ pn Pn and fn ¼ un Un , then from (1) and (9) (with c ¼ 1 now), we get
! n ! n ! n1 2 2 2 h 2 nj nj h 2 b h 2 n 1 þ dx þ dx fj 1 þ dx f þ 1 þ dx ½cnnj ¼ Rn1;j ; a j 6 s 6 6 2
dx nnj
h þ 1 þ d2x 6
!
n 1 f ¼ Rn2;j ; a j
4
ð18Þ
1 6 j 6 Nx ; 1 < n 6 N;
ð19Þ
4
where jRn1;j j 6 Cðs þ h Þ; jRn2;j j 6 Ch , for 1 6 j 6 N x ; 1 < n 6 N.
For any real vector functions f ; g defined on the mesh points xj ð1 6 j 6 N x Þ, let A be the matrix so that 2 T 1 þ h6 d2x f ; g ¼ f Ag holds, it follows that
0
4
B B1 B B B 1B B0 A¼ B. 6B . B. B B. B .. @
1
0
4
1
1 .. . .. .
4 .. . .. .
1 0 .. .. .C . . .. C C C .. .. C . .C 1 C C; .. .. C . . 0C C C 1 4 1C A
0
0
1
4
which is symmetric, positive definite, and in fact, we have
! ! 2 h 2 1 þ dx f ; f 6 kf k2 : 6
1 kf k2 6 3
ð20Þ
Hence we can introduce a weighted norm kf kA ¼ T
1=2
ðf Af Þ
ðg T AgÞ
1þ
1=2
6 12 ðkf k2A þ kgk2A Þ. Now we can multiply (18) by hnnj , (19) by hfnj , and sum for 1 6 j 6 N x to get
! ! 2 h 2 nn nn1 n dx ; n þ ðdx fn ; nn Þ 6 s 2
n
1=2 2 2 T 1 þ h6 d2x f ; f . Besides, we have 1 þ h6 d2x f ; g ¼ jf Agj 6
n
ðdx n ; f Þ þ
h 1 þ d2x 6
!
2
1þ
n ! 1 f ; fn ¼ ðRn2 ; fn Þ; a
h 2 d 6 x
! n ! b f ; nn þ a
1þ
! ! 2 h 2 dx ½cnn ; nn ¼ ðRn1 ; nn Þ; 6
1 < n 6 N:
ð21Þ
ð22Þ
As we assume that nn0 ¼ nnNx þ1 ¼ fn0 ¼ fnNx þ1 ¼ 0 holds true, therefore,
ðdx nn ; fn Þ ¼ h
Nx X nnjþ1 nnj1
2h
j¼1
fnj ¼
Nx x þ1 x 1 1 NX 1 NX 1X 1 nnj fnj1 nnj fnjþ1 ¼ nn ðfn fnjþ1 Þ þ ðnnNx þ1 fnNx þ nnNx fnNx þ1 nn0 fn1 nn1 fn0 Þ 2 j¼2 2 j¼0 2 j¼1 j j1 2
n
¼ ðdx f ; nn Þ: Hence, we can add (21) and (22) together to obtain
1þ
! ! 2 h 2 nn nn1 n dx ;n þ 6 s
¼ ðRn1 ; nn Þ þ ðRn2 ; fn Þ; thus
2
1þ
h 2 d 6 x
! ! n 1 f ; fn þ a
1þ
! ! 2 h 2 dx ½cnn ; nn 6
2
1þ
h 2 d 6 x
!
n ! b f ; nn a
1 < n 6 N;
! n ! 2 1 1 h 2 b 1 1 1 2 2 n 2 n1 2 ðkn kA kn kA Þ þ kfkA þ cknkA 1 þ dx f ; nn 6 ðknn k2 þ kfn k2 Þ þ 3 kRn1 k2 þ akRn2 k2 2s a a 3 4a 4 6 1 1 6 knn k2A þ kfn k2A þ 3 kRn1 k2 þ akRn2 k2 ; 1 < n 6 N: 4a 4
Note that
2
h 1 þ d2x 6
! ! n b b n f ; n ¼ jj a a
! ! 2 2 h 2 n n 1 n 2 b 1 þ dx f ; n 6 kf kA þ knn k2A ; 4a 6 a
471
M.R. Cui / Applied Mathematics and Computation 246 (2014) 464–473
consequently, we get 2
1þ2 c1
b a
! !
1 a
s knn k2A þ skfk2A 6 knn1 k2A þ 6s
1 n 2 kR1 k þ akRn2 k2 : 4
ð23Þ
2
2
Based on the sign of c 1 ba , we divide the discussion of (23) in the following two cases. If c P 1 þ ba , then we have
knn k2A þ
n X 1 1 1 l 2 4 2 skfk2A 6 knn1 k2A þ 6s kRn1 k2 þ akRn2 k2 6 6s kR1 k þ akRl2 k2 6 Cðs þ h Þ : a 4 4 l¼0 2
Otherwise, let c0 ¼ j2ðc 1 ba Þj, for
16
1 2
1 þ 2ðc 1 ba Þs
¼
2
s small enough such that both 1 þ 2ðc 1 ba Þs > 0 and
1 6 1 þ 2c0 s 6 2 1 c0 s
hold, note that n0j ¼ 0ð0 6 j 6 N x þ 1Þ, we get
knn k2A þ
n X 2 1 1 1 skfk2A 6 ð1 þ 2c0 sÞknn1 k2A þ 12s kRn1 k2 þ akRn2 k2 6 12s ð1 þ 2c0 sÞnl kRl1 k2 þ akRl2 k2 6 Cðs þ h4 Þ : a 4 4 l¼1
Thus we can obtain
knn k2A þ
2 1 skfk2A 6 Cðs þ h4 Þ a
ð24Þ
again, hence we can prove (17) by using (20) and (24). The proof of the theorem is completed.
h
Remark 4.1. In the proof we have used the assumption nn0 ¼ nnNx þ1 ¼ fn0 ¼ fnNx þ1 ¼ 0. For problems with Dirichlet boundary conditions, we have nn0 ¼ nnNx þ1 ¼ 0, likewise we have fn0 ¼ fnNx þ1 ¼ 0 for Neumann boundary problems. However, we cannot give a condition to make all these conditions hold true at present. Remark 4.2. For the partial differential equations of fractional order c 2 ð0; 1Þ, in the estimate we will have a factor ð1 C sc Þn multiplied on the righthand side, but this term tends to infinity when c < 1, so we only discuss c ¼ 1 here. 5. Numerical experiments We give some computational results in this section, to test the convergence order of our scheme. In the numerical simulations, we use the discrete l2 -error norm which include errors on the boundary points for both test Problems 1 and 2. Table 1 Results for test Problem 1 with c ¼ 0:2; 0:4; 0:6; 0:8 and 1:0.
c
h
s
N
keN pk
order
keN uk
order
0.2 0.2 0.2 0.2
0.1 0.05 0.025 0.0125
0.1 0.0214 0.0046 9.8431e-4
10 47 218 1016
3.6451e-5 2.3881e-6 1.5520e-7 9.9066e-9
– 3.9320 3.9437 3.9696
2.2814e-4 1.3622e-5 8.4283e-7 5.2457e-8
– 4.0659 4.0146 4.0060
0.4 0.4 0.4 0.4
0.1 0.05 0.025 0.0125
0.1 0.0177 0.0031 5.5243e-4
10 57 320 1810
2.2911e-4 1.4357e-5 9.0508e-7 5.6675e-8
3.9962 3.9876 3.9973
0.0014 8.1324e-5 4.8879e-6 2.9865e-7
– 4.1056 4.0564 4.0327
0.6 0.6 0.6 0.6
0.1 0.05 0.025 0.0125
0.1 0.0138 0.0019 2.6287e-4
10 72 525 3804
9.0349e-4 5.6485e-5 3.5221e-6 2.2017e-7
– 3.9996 4.0034 3.9997
0.0056 3.1955e-4 1.9003e-5 1.1593e-6
– 4.1313 4.0717 4.0349
0.8 0.8 0.8 0.8
0.1 0.05 0.025 0.0125
0.1 0.0099 9.8431e-4 9.7656e-5
10 101 1016 10240
0.0029 1.8331e-4 1.1459e-5 7.1620e-7
– 3.9837 3.9997 4.0000
0.0182 0.0010 6.1813e-5 3.7703e-6
– 4.1859 4.0159 4.0352
1.0 1.0 1.0 1.0
0.1 0.05 0.025 0.0125
0.1 0.0063 3.9063e-4 2.4414e-5
10 160 2560 40960
0.0085 5.3033e-4 3.3146e-5 2.0716e-6
– 4.0025 4.0000 4.0000
0.0524 0.0030 1.7877e-4 1.0904e-5
– 4.1265 4.0688 4.0352
472
M.R. Cui / Applied Mathematics and Computation 246 (2014) 464–473
Table 2 Results for test Problem 2 with c ¼ 0:2; 0:4; 0:6; 0:8 and 1:0.
c
h
s
N
keN pk
order
keN uk
order
0.2 0.2 0.2 0.2
0.1 0.05 0.025 0.0125
0.1 0.0214 0.0046 9.8431e-4
10 47 218 1016
0.0082 5.7642e-4 3.9408e-5 2.6200e-6
– 3.8304 3.8706 3.9108
0.0053 3.7669e-4 2.5776e-5 1.7144e-6
– 3.8145 3.8693 3.9103
0.4 0.4 0.4 0.4
0.1 0.05 0.025 0.0125
0.1 0.0177 0.0031 5.5243e-4
10 57 320 1810
0.0554 0.0037 2.4204e-4 1.5378e-5
– 3.9043 3.9342 3.9763
0.0361 0.0024 1.5831e-4 1.0063e-5
– 3.9109 3.9222 3.9756
0.6 0.6 0.6 0.6
0.1 0.05 0.025 0.0125
0.1 0.0138 0.0019 2.6287e-4
10 72 525 3804
0.2006 0.0131 8.2802e-4 5.1725e-5
– 3.9367 3.9838 4.0007
0.1309 0.0086 5.4158e-4 3.3847e-5
– 3.9280 3.9891 4.0001
0.8 0.8 0.8 0.8
0.1 0.05 0.025 0.0125
0.1 0.0099 9.8431e-4 9.7656e-5
10 101 1016 10240
0.5495 0.0351 0.0022 1.3538e-4
– 3.9686 3.9959 4.0224
0.3584 0.0229 0.0014 8.8586e-5
– 3.9682 4.0318 3.9822
1.0 1.0 1.0 1.0
0.1 0.05 0.025 0.0125
0.1 0.0063 3.9063e-4 2.4414e-5
10 160 2560 40960
1.2862 0.0784 0.0048 3.0034e-4
– 4.0361 4.0297 3.9984
0.8390 0.0512 0.0032 1.9653e-4
– 4.0345 4.0000 4.0253
" keNp k ¼ kpN PN k ¼
N x þ1 X
" keNu k
N
N
¼ ku U k ¼
2
#1=2
ðpNj PNj Þ h
j¼0
N x þ1 X
; #1=2
ðuNj
2 U Nj Þ h
:
j¼0 4
As in paper [12], since the local truncation error of the schemes are proved to be Oðs2c Þ þ Oðh Þ, then if we decrease the mesh size of h to h=2 and s to s=24=ð2cÞ , we expect to get a factor 1/16 for the error kei s; h; tN k with index i indicates p or u. If we decrease h to h=4, then we multiply s by s=28=ð2cÞ . The numbers N in Table 1 and Table 2 are calculated by 1s , that is, we use the round number of 1s, hence N s – 1 but N s 1. The experimental order of convergence r i ðs; hÞ is computed by the formula
1 kei 24=ð2cÞ s; 2h; tN k A; ri ðs; hÞ ¼ log2 @ kei s; h; t N k 0
i ¼ p or u:
5.1. Test problem 1 Consider the Dirichlet problem (1)–(3) with the coefficients aðxÞ ¼ 1 þ x; bðxÞ ¼ 1 þ 2x; cðxÞ ¼ 1 þ 5x, the theoretical solution is pðx; tÞ ¼ ex t1þc , and the right hand side function f ðxÞ ¼ tex ðCð2 þ cÞ þ 6xt c Þ. The errors for c ¼ 0:2; 0:4; 0:5; 0:6; 0:8 and 1.0, respectively, together with the experimental convergence order are shown in Table 1. We see that the order is approximately 4, that is, the performance meets our expectations. 5.2. Test problem 2 Inspired by the example in paper [28], we use the Caputo derivative now and consider the Robin problem (1), (2) and (4) for the following equation C c 0 Dt p
¼
@ 1 @2p p þ 2 þ f ðx; tÞ; @x 1 þ x @x
0 < x < 1; 0 < t 6 1;
1þc 1 1 with bðx; tÞ ¼ 1þx ; cðx; tÞ ¼ ð1þxÞ . 2 and f ðx; tÞ ¼ 8ðx þ 1ÞCð2 þ cÞt now, and the exact solution is pðx; tÞ ¼ 8ðx þ 1Þt
For this test problem, the errors for p and u with c ¼ 0:2; 0:4; . . . ; 1:0, together with the experimental convergence order are shown in Table 2. The results show that the order is also approximately 4, that is, the l2 norm of both keNp k and keNu k is 4
Oðs2c Þ þ Oðh Þ.
M.R. Cui / Applied Mathematics and Computation 246 (2014) 464–473
473
6. Conclusion In this paper we propose a CCD scheme for solving the one dimensional time fractional convection–diffusion–reaction equation with variable coefficients. Both Dirichlet and Robin boundary conditions are considered. We solve both the unknown and its gradient simultaneously. The proposed algorithm has high accuracy and good efficiency. The local trunca4 tion error of the CCD scheme is Oðs2c Þ þ Oðh Þ. We give the convergence analysis for the problem of integer order with constant coefficients under some assumption.The conclusion is verified by some numerical experiments. Acknowledgement The author thanks the referees for their valuable comments and suggestions, which help to improve this manuscript significantly. References [1] R. Metzler, J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Phys. Rep. 339 (2000) 1–77. [2] I. Podlubny, Fractional Differential Equations, An Introduction to Fractional Derivatives, Fractional Differential Equations, Some Methods of Their Solution and Some of their Applications, Academic Press, San Diego, 1999. [3] C.-M. Chen, F. Liu, I. Turner, V. Anh, A Fourier method for the fractional diffusion equation describing sub-diffusion, J. Comput. phys. 227 (2007) 886– 897. [4] P. Zhuang, F. Liu, V. Anh, I. Turner, New solution and analytical techniques of the implicit numerical method for the anomalous subdiffusion equation, SIAM J. Numer. Anal. 46 (2008) 1079–1095. [5] P. Zhuang, F. Liu, V. Anh, I. Turner, Numerical methods for the variable-order fractional advection-diffusion equation with a nonlinear source term, SIAM J. Numer. Anal. 47 (2009) 1760–1781. [6] R.S. Hirsch, Higher order accurate difference solutions of fluid mechanics problems by a compact difference technique, J. Comput. Phys. 24 (1975) 90– 109. [7] S.K. Lele, Compact finite difference schemes with spectral-like resolution, J. Comput. Phys. 103 (1992) 16–42. [8] M.R. Cui, Compact finite difference method for the fractional diffusion equation, J. Comput. Phys. 228 (2009) 7792–7804. [9] C.-M. Chen, F. Liu, V. Anh, I. Turner, Numerical schemes with high spatial accuracy for a variable-order anomalous subdiffusion equation, SIAM J. Sci. Comput. 32 (2010) 1740–1760. [10] G.-H. Gao, Z.-Z. Sun, A compact finite difference scheme for the fractional sub-diffusion equations, J. Comput. Phys. 230 (2011) 586–595. [11] M.R. Cui, Compact alternating direction implicit method for two-dimensional time fractional diffusion equation, J. Comput. Phys. 231 (2012) 2621– 2633. [12] M.R. Cui, Convergence analysis of high-order compact alternating direction implicit schemes for the two-dimensional time fractional diffusion equation, Numer. Algor. 62 (2013) 383–409. [13] M.R. Cui, A high-order compact exponential scheme for the fractional convection-diffusion equation, J. Comput. Appl. Math. 255 (2014) 404–416. [14] G.F. Carey, W.F. Spotz, Higher-order compact mixed methods, Commun. Numer. Methods Eng. 13 (1997) 553–564. [15] W.F. Spotz, G.F. Carey, High-order compact scheme for the steady stream-function vorticity equations, Int. J. Numer. Methods Eng. 38 (1995) 3497– 3512. [16] W.F. Spotz, G.F. Carey, Extension of high-order compact schemes to time-dependent problems, Numer. Methods Partial Differ. Eqs. 17 (2001) 657–672. [17] S. Abide, X. Chesneau, B. Zeghmati, Compact mixed methods for convection/diffusion problems, Appl. Math. Comput. 218 (2012) 5867–5876. [18] J. Zhao, Highly accurate compact mixed methods for two point boundary value problems, Appl. Math. Comput. 188 (2007) 1402–1408. [19] F.C. Chu, C. Fan, A three-point combined compact difference scheme, J. Comput. Phys. 140 (1998) 370–399. [20] F.C. Chu, C. Fan, A three-point sixth-order nonuniform combined compact difference scheme, J. Comput. Phys. 148 (1999) 663–674. [21] K. Mahesh, A family of high order finite difference schemes with good spectral resolution, J. Comput. Phys. 145 (1998) 332–358. [22] T. Nihei, K. Ishii, A fast solver of the shallow water equations on a sphere using a combined compact difference scheme, J. Comput. Phys. 187 (2003) 639–659. [23] J. Zhang, J.J. Zhao, Truncation error and oscillation property of the combined compact difference scheme, Appl. Math. Comput. 161 (2005) 241–251. [24] Q. Zhou, Z.H. Yao, F. He, M.Y. Shen, A new family of high-order compact upwind difference schemes with good spectral resolution, J. Comput. Phys. 227 (2007) 1306–1339. [25] T.K. Sengupta, V. Lakshmanan, V.V.S.N. Vijay, A new combined stable and dispersion relation preserving compact scheme for non-periodic problems, J. Comput. Phys. 228 (2009) 3048–3071. [26] T.K. Sengupta, V.V.S.N. Vijay, S. Bhaumik, Further improvement and analysis of CCD scheme: dissipation discretization and de-aliasing properties, J. Comput. Phys. 228 (2009) 6150–6168. [27] W.J. Chen, J.C. Chen, Combined compact difference method for solving the incompressible Navier–Stokes equations, Int. J. Numer. Methods Fluids 68 (2012) 1234–1256. [28] S. Chen, F. Liu, P. Zhuang, V. Anh, Finite difference approximations for the fractional Fokker–Planck equation, Appl. Math. Model. 33 (2009) 256–273. [29] B.I. Henry, T.A.M. Langlands, P. Straka, Fractional Fokker–Planck equations for subdiffusion with space- and time-dependent forces, Phys. Rev. Lett. 105 (2010) 170602. [30] Z.Z. Sun, X.N. Wu, A fully discrete difference scheme for a diffusion-wave system, Appl. Numer. Math. 56 (2006) 193–209.
