An efficient high‐order two‐level explicit/implicit numerical scheme for two‐dimensional time fractional mobile/immobile advection‐dispersion model

This article constructs a new two‐level explicit/implicit numerical scheme in an approximate solution for the two‐dimensional time fractional mobile/immobile advection‐dispersion problem. The stability and error estimates of the proposed technique are deeply analyzed in the L∞(0,T;L2)$$ {L}^{\infty}\left(0,T;{L}^2\right) $$ ‐norm. The developed approach is less time consuming, fourth‐order in space and temporal accurate of order O(k2−λ2)$$ O\left({k}^{2-\frac{\lambda }{2}}\right) $$ , where k$$ k $$ is the time step and λ$$ \lambda $$ denotes a positive parameter less than 1. This result shows that the two‐level explicit/implicit formulation is faster and more efficient than a large class of numerical schemes widely discussed in the literature for the considered problem. Numerical experiments are performed to verify the theoretical studies and to demonstrate the efficiency of the new numerical method.


INTRODUCTION
In the last decades, a great deal of works has been devoted to fractional calculus and its applications in various fields of finance, biology, economy, physics, and chemistry are witnessing remarkable progress. 1-5Fractional mobile/immobile advection-diffusion models describe a special type of anomalous convection-diffusion-reaction equation and lie in an important class of fractional partial differential equations (FPDEs).The complexity of such equations comes from the adding of time fractional derivative that is nonlocal and has the history dependence and universal mutuality character.Specifically, find efficient solutions to these equations requires substantial amounts of calculations.This is equivalent to solve a certain class of complex ordinary/partial differential equations (e.g., see .Thus, numerical methods should use all previous solutions to compute the solution on the current time-level.This suggests that the computational cost becomes very expensive if implicit schemes are applied, especially for multidimensional problem.17  High-order explicit numerical methods are fast and computational efficient, but they suffer a severe time step restriction for stability.To overcome this drawback, suitable techniques such as high-order explicit/implicit methods in the temporal direction, so called high-order predictor-corrector schemes should be developed and deeply analyzed.Such approaches are not time consuming and they do not suffer a severely restricted time step for stability condition.18-25 In the literature 17,26-31, a wide set of explicit schemes, implicit methods, spectral approximations, collocation techniques and 1306 NGONDIEP compact difference schemes of order O(k + h 2 ), O(k 2− + h 2 ), O(k + h 4 ), O(k 2− + h 4 ) and O(k 2−  2 + h 4 ), where h is the grid space, have been discussed in a computed solution of time FPDEs.In this article, we construct a two-level high-order explicit-implicit approach in a numerical solution of the following two-dimensional time fractional mobile/immobile advection-dispersion equation with initial condition w(x, y, 0) = w 0 (x, y), on Ω ∪ Ω, (2) and boundary condition w(x, y, t) = g(x, y, t) on Ω × [0, T], (3) where Ω denotes the boundary of Ω, ∇ and △ are the gradient and Laplacian operators, respectively, w represents the solute concentration in the total (mobile+immobile), d means the fractional capacity coefficients whereas a and b designate the dispersion and velocity coefficients for the mobile phase, respectively.In Reference 32 the authors added the time derivative w t to conveniently describe the motion time and thus to distinguish the status of particles. is the reaction coefficient while f = f (x, y, t) denotes the source term, w 0 is the initial condition, g represents the boundary conditions whereas cD  0t w, for 0 <  < 1, is the Caputo fractional derivative of order  defined in Reference 33 by Under a suitable time step restriction, the proposed algorithm is stable, convergent with order O(k 2−  2 + h 4 ), more faster and efficient than a broad range of numerical schemes largely studied in the literature for the considered initial-boundary value problem (1)- (3).For more details, we refer the readers to References 34-40 and references therein.
The highlights of this work are the following: (i) construction of the new two-level explicit/implicit technique for solving the time fractional mobile/immobile advection-dispersion equation (1) subjected to initial condition (2) and boundary condition (3), (ii) analysis of the stability together with the convergence rate of the proposed numerical procedure, (iii) numerical experiments that validate the theoretical studies and demonstrate the efficiency and utility of the developed approach.
The article is organized as follows.In Section 2, we develop the two-level explicit/implicit method in an approximate solution of the initial-boundary value problem (1)- (3).The stability and error estimates of the new algorithm are deeply analyzed in Section 3. Section 4 considers a wide set of numerical examples whereas the general conclusions and future works are provided in Section 5.

CONSTRUCTION OF THE TWO-LEVEL EXPLICIT/IMPLICIT NUMERICAL SCHEME
In this section, we develop a high-order two-level explicit/implicit numerical method for solving the two-dimensional time fractional mobile/immobile advection-dispersion equation (1) subjected to initial and boundary conditions (2) and (3), respectively.For the sake of discretization, we set Ω = ( 1 ,  2 ) × ( 3 ,  4 ).Let M and N be positive integers.Set k = T N , h x =  2 − 1 M and h y =  4 − 3 M be the time step and the space steps in the x-direction and y-direction, respectively.For the temporal and spatial discretizations, we introduce uniform grid of mesh points (x i , y j , t n ), where x i =  1 + ih x , y j =  3 + jh y , and t n = nk, for i, j = 0, 1, 2, … , M, and n = 0, 1, 2, … , N. In addition, a regular partition of the domain Ω × I and the space of mesh functions are defined as Ω h x h y k = {(x i , y j , t n ), 0 ≤ i, j ≤ M; n = 0, 1, 2, … , N} and  h x h y k = {w(x i , y j , t n ), 0 ≤ i, j ≤ M; 0 ≤ n ≤ N}.For the sake of readability, we let w(x i , y j , t n ) = w n ij .Furthermore, the analytical solution of the initial-boundary value problem (1)-( 3) and the numerical one provided by the constructed algorithm at the discrete point (x i , y j , t n ) are denoted by w n ij and W n ij , respectively.The values of the source term f at the grid point (x i , y j , t n ) is given by f n ij .Finally, we define the following operators Let 0 <  < 2 −1 , be a positive parameter.Expanding the Taylor series for the function w at the grid point (x i , y j , t n ) using forward and backward difference formulations, this yields where w z and w 2z denote w z and  2 w z 2 , respectively.In the following, we will adopt these notations.Subtracting Equation ( 7) from ( 6) to get which is equivalent to Combining Equations ( 8) and (1), simple manipulations give Applying the Taylor series with time step k utilizing forward difference, it is easy to see that for z ∈ {0, x, y}, m = 1, 2 and s ∈ {0, 1 2 , 1}, where we set w m0 = w.Substituting (10) into Equation ( 9) and rearranging terms result in where w zt and w 2zt designate  2 w zt and  3 w z 2 t , respectively.Furthermore, the application of the Taylor series expansion for the function w about the mesh point (x i , y j , t n+ 1

2
) utilizing both forward and backward difference representations provides Plugging both equations, we obtain In addition, using the Taylor expansion, it holds Subtracting ( 13) from ( 14), utilizing (15) and rearranging terms to obtain Plugging this equation, ( 15) and ( 12), it is not hard to see that Applying Equation ( 1) at the grid points (x i , y j , t n+ 1 2 + ) and (x i , y j , t n+1+ ) to get Combining Equations ( 16)-( 18) and performing direct calculations yield which is equivalent to Substituting the equations given by relation (10) into (19), this provides Setting where "*" denotes the usual multiplication in the complex field C, Equations ( 11) and ( 20) become ) for n = 0, 1, 2, … , N − 1, and i, j = 0, 1, 2, … , M. In Reference 30 (p. 5, 6, 8, 11, 12, and 17), it has been established that the terms c D  0t w n+s+ ij , for s = 0, 1 2 , 1, can be approximated as where is the error term and the summation index "l" varies in the range: Furthermore, the sequences of coefficients {a  n+s,l } l≤n+s , are defined by: for n ≥ 1 For q = 0, 1, 2, … , n − 1, the terms f and d , are given by for n ≥ 1 and q = 0, 1, 2, … , n − 1.Furthermore, For n ≥ 1 For n ≥ 0 and q = 0, 1, 2, … , n, the terms f  n+1,q and d n+1,q are defined as follows: In addition, the error terms I ,n+s+ ij are given by and where ) is the error associated to the second-order polynomial approximating the function w ij (⋅) at the mesh points (w ), (w l+1 ij , t l+1 ) and (w ) (resp., (w l ij , t l ), (w ) and (w l+1 ij , t l+1 )).We should provide the spatial fourth-order approximations of the partial derivatives w n+s mz,ij , for m = 1, 2, z = x, y, and s{0, 1  2 , 1}.Utilizing the Taylor series expansion for w about the discrete point (x i , y j , t n+s ) with step size h z , for z = x, y, using both forward and backward representations, one easily proves as for the one-dimensional case discussed in Reference 30 that where Now, replacing in Equations ( 35)-( 38), x by y and rearranging indices to obtain where For s = 0 and n ≥ 1, combining Equations ( 23), ( 25), ( 26), ( 35), ( 36), (39), and (40), simple calculations yield where Furthermore, For s = 1, 1 2 and n ≥ 0, plugging Equations ( 24)-( 26), ( 35), ( 36), (39), and (40), direct computations give where We remind that We introduce the following linear operator: where  4 x ,  4 2x ,  4 y and  4 2y are defined by Equations ( 37), ( 38), (41), and (42), respectively.Tracking the error terms  n 1,ij and  n+1 2,ij , replacing the exact solution w n+s ij with the computed one W n+s ij at the grid point (x i , y j , t n+s ), for s = 0, 1 2 , 1, and rearranging terms, relations ( 43) and ( 45) can be approximated as for n = 0, 1, 2, … , N − 1, and i, j , Equation (49) can be rewritten as For n = 0, Equation ( 6) can be rewritten as: Plugging this together ( 1) and ( 47) to get This is equivalent to In Reference 29, replacing in Equations ( 16) and ( 17),   j and u j by  and w ij , respectively, we obtain where where 1,ij is the interpolating polynomial of degree one approximating w ij (⋅) at the grid points (w 0 ij , t 0 ) and (w  ij , t  ).Utilizing the Taylor series expansion, it is easy to show that Substituting ( 54) and ( 52) into (51) and using (53), direct calculations give Truncating the error term: , replacing the analytical solution w with the approximate one W, and the operator L by L h x h y , Equation (55) implies where A combination of Equations ( 48), (50), and (56) provides the new two-level explicit/implicit numerical approach in a numerical solution of the two-dimensional mobile/immobile advection-dispersion equation (1) subjected to initial and boundary conditions ( 2) and (3), respectively.That is, for i, j = 2, 3, … , M − 2, with initial and boundary conditions To start the algorithm, we set be the vector of computed solution at time t s .Let and be a block pentadiagonal matrix of size (M − 2) 2 , where I is the identity matrix of size M − 2, B h is the pentadiagonal matrix of size M − 2, given by .
Furthermore, we introduce the block pentadiagonal matrices M(0) kh , M (n) kh and M(n) kh , of size (M − 2) 2 , defined as where where where Utilizing ( 64)-( 66), ( 68), (70), and (72), it is not difficult to observe that Equations (58), (59), and (61) can be written in the matrix form as: ] , (75) ] , (76) where the summation index "l" varies in the range: l = 0, 1 2 , 1, … , n − 1 2 , n.The sequences of matrices {M (n)  kh } n and { M(n) kh } n can be observed as the g-Toeplitz matrix sequences generated by Lebesgue integrable functions u and v, respectively, defined over the domain (−, ).For more details about the Toeplitz matrices, g-Toeplitz and g-circulant structures, the readers can consult the works discussed in References 41-43.Furthermore, after discretization each system of linear (or nonlinear) equations can be solved efficiently using the notion of symbol and the related spectral analysis using either preconditioned Krylov methods or multigrid solvers. 44-46

STABILITY ANALYSIS AND ERROR ESTIMATES OF THE PROPOSED TWO-STEP APPROACH (58)-(62)
In this section, we analyze the stability and error estimates of the proposed numerical technique (58)-(62) for solving the two-dimensional time-fractional mobile/immobile advection-dispersion equation ( 1) subjected to initial-boundary conditions ( 2) and ( 3), under the time step restriction where: 0 <  < 1, C p > 0 denotes the Poincaré-Friedrich constant, a 1 = 2 2− − 1 2 > 0 and h is the space step in the x, y-directions.When  tends to 1, it is easy to see that estimate (77) provides The left side of this inequality depends on k  h 2 , so that in the usual terminology of explicit/implicit numerical scheme, the developed two-level procedure is stable under a suitable time step constraint.Furthermore, when the parameter  approaches 1 and  = 0, requirement (77) is more attractive and it is equivalent to restriction (78) which is well known in the literature as the Courant Friedrich-Levy condition ([CFL] condition).For the convenient of writing, we set in the following h = h x = h y .
We introduce the discrete norms In addition, the Hilbert space L 2 (Ω) is equipped with the scalar product (⋅, ⋅) 2 defined as Analogously, the quantities (  2 y w n ,  2 y v n ) 2 satisfy similar relations.Furthermore, where , respectively.Finally, we assume that the analytical solution w satisfies the regularity condition: w ∈  12,3  D .Thus, there is a positive constant  0 , independent of the mesh grid h and time step k, so that Theorem 1 (Stability and error estimates).Let w be the analytical solution of the initial-boundary value problem ( 1)-( 3) and W be the approximate one provided by the new algorithm ( 58)-( 62).Suppose  ∈ (0, 2 3 ) be a parameter and let 0 <  < 1  2 .Consider the generalized sequences (a  n+s,l ) 1 2 ≤l≤n+s , for s ∈ {0, 1 2 , 1}, defined by relations ( 27), ( 28), (30), and (31).Under the time step restriction (77), it holds where e n = w n − W n is the error term, Ĉ > 0, is a constant independent of the space step h and time step k.
It is worth mentioning that relation (83) shows that the constructed numerical scheme (58)-( 62) is stable whereas estimate (84) indicates that the developed approach is fourth-order spatial convergent and temporal accurate of order O(k 2−  2 ).The proof of Theorem 1 requires some intermediate results (namely, Lemmas 1-5).
Proof.The proof of Lemma 2 is similar to the proof of Lemma 1 established in Reference 30 (p. 10), by replacing the domain Consider w, v ∈  0,0 D , be functions defined over the domain D = Ω × [0, T].Suppose w mj = v mj = 0 and w im = v im = 0, for any m = {0, 1, M − 1, M}, and i, j = 0, 1, 2, … , M. Thus, the linear operator L h defined in (47) satisfies where C p is a positive constant that does not dependent on the time step k and the space step h.
Proof.Expanding the left side and rearranging terms gives the result.▪
(103) But it comes from the initial condition ( 62) that e 0 ij = 0, for i, j = 0, 1, 2, … , M. Furthermore, Utilizing this, it is not difficult to observe that ( 102) and ( 103) can be rewritten as (105) It follows from relations ( 27)-( 32) that Thus, the sequences {a  n+s,n+s } n+s and {a  n+s, 1   2   } n+s are bounded for any value of the parameter  ∈ (0, and taking the limit in both sides of estimates ( 104) and (105) when  tends to zero, this results in (108) The application of the Cauchy-Schwarz inequality gives In addition, it follows from the second estimate of relation ( 90) and the Cauchy-Schwarz inequality that Plugging estimates (107) and (109) (respectively, estimates (108), ( 110)-( 112)) and rearranging terms, to obtain (114) A combination of (113) and ( 114) provides Using the boundary condition (62) and the definition of the discrete norm and scalar product given by ( 79) and (80), it is not difficult to show that Substituting this into (115) and rearranging terms, this yields and  ∈ (0, 2 3 ), utilizing relations (106) and ( 27)- (32), it holds (117) Since for k small and satisfying the time step restriction (77), it holds: ) ≥ 0. In addition, utilizing the Taylor series expansion together with (117), it is not difficult to show that where "*" denotes the usual multiplication in the set of real numbers, t ] .
This estimate completes the proof of Theorem 1. ▪

NUMERICAL EXPERIMENTS AND CONVERGENCE RATE
In this section, we present some computational results to verify stability and the convergence of the proposed approach (58)-(62) for solving the initial-boundary value problem (1)-(3).To demonstrate the efficiency of the new algorithm, we set k = h In each example, the computations are performed with different values of the parameter :  ∈ {0.01, 0.25, 0.9}.Furthermore, the convergence order R(k, h) of the developed technique is estimate using the formula where we set k = h 8 4− .Finally, the numerical calculations are carried out by the use of MATLAB R2013b.Figures 1-5 suggest that the developed two-step explicit/implicit approach (58)-( 62) is stable whereas Tables 1-7 indicate that the new algorithm is spatial fourth-order accurate and temporal convergent with order O(k 2−  2 ).These numerical results confirm the theory (see Section 3, Theorem 1).

F I G U R E 5
Exact solution (w: in green), computed solution (W: in blue) and error (E: in red) for Example 4. Stability analysis and convergence order of the proposed approach for time fractional mobile-immobile with  = 0.01,  = 0.25 and k = h  Note: We take  = 0.25,  = 0.9 and k = h
The initial and boundary conditions are directly obtained from the exact solution.
Example 3. Suppose D be the bounded domain (0, 1) 2 × (0, 1].We assume that the parameters  = 0.01 and  = 0.9.We consider the time-fractional mobile/immobile solute transport model defined as .The initial and boundary conditions are directly obtained from the exact solution w which is given by w(x, y, t) = e −(x− 1 2 ) 2 −(y−

GENERAL CONCLUSIONS AND FUTURE WORKS
In this article, we have developed a two-level fourth-order explicit/implicit approach for solving the two-dimensional time fractional mobile/immobile advection-dispersion equation ( 1) with initial and boundary conditions ( 2) and (3), respectively.Under the suitable time step requirement (77), the stability and the error estimates of the proposed two-level explicit/implicit numerical scheme have been analyzed in L ∞ -norm.Both theoretical and computational results have suggested that the new algorithm stable, spatial fourth-order convergent and temporal accurate with order O(k 2−  2 ).This indicates that the constructed numerical technique (58)-( 62) is more faster and efficient than a large set of methods discussed in the literature 35-40 for solving the time fraction equation (1) subjected to initial and boundary conditions (2) and (3), respectively.Furthermore, the new formulation (58)-(62) can be observed as a robust numerical approach for solving general time FPDEs.Finally, it's well known in the literature 14-16,18 that the Von Neumann stability approach, based on a Fourier analysis in the space domain is a powerful tool to study the global amplification factor of numerical schemes applied to linear time-dependent PDEs while this technique is not suitable for nonlinear unsteady PDEs.Regarding the time fractional PDEs as the one considered in this work, for time steps sufficiently small, Equations (66)-(73) indicate that both matrices M (n)  kh and M(n) kh are good approximations of the scalar matrix 12h 2 I (M−2) 2 , where I (M−2) 2 denotes the identity matrix of size (M − 2) 2 and h represents the space step.This suggests that the product (M (n) kh ) −1 M(n) kh is a good approximation of I (M−2) 2 , for any values of the mesh grid h.So the numerical amplification factor of the developed approach (74)-( 76) is close to one.In the case where the time step and the order of time fractional derivative are not small, neither M (n)  kh nor M(n) kh are approximations of 12h 2 I (M−2) 2 whereas their coefficients depend on the order of time fractional derivative.Additionally, the coefficients of both matrices show that M (n)  kh cannot be considered as a good preconditioner for M(n) kh .Thus, the amplification factor of the proposed algorithm which equals the spectral radius of (M (n) kh ) −1 M(n) kh strongly depends on the order of the time fractional derivative appearing in the governing PDE.Our future works will be focused on the development of a two-level explicit/implicit scheme with convergence order O(k 2 + h 4 ) in a computed solution of multidimensional time FPDEs.

APPENDICES
In this appendix, provide the proof of some intermediate results (Lemmas 3 and 4).

APPENDIX A. PROOF OF LEMMA 3
This subsection establishes the Proof of Lemma 3.

) 2 .F I G U R E 2
Exact solution (w: in green), numerical solution (W: in blue) and error (E: in red) for Example 1. Stability and convergence of the developed two-level explicit/implicit approach for time fractional mobile-immobile with  = 0.01,  = 0.25 and k = h 8 4− .[Colour figure can be viewed at wileyonlinelibrary.com]

) 1 2 ,F I G U R E 4
Exact solution (w: in green), computed solution (W: in blue) and error (E: in red) for Example 3. Stability and convergence rate of the new algorithm for time fractional mobile-immobile with  = 0.01,  = 0.9 and h = k 4− 8 .[Colour figure can be viewed at wileyonlinelibrary.com] and |||e
Exact solution (w: in green), approximate solution (W: in blue) and error (E: in red) for Example 1. Analysis of stability and convergence of the proposed two-level explicit/implicit scheme for time fractional mobile-immobile with  =  = 0.01 and k = h
TA B L E 1Note: We take  = 0.01,  = 0.25 and k = h TA B L E 3 Analysis of stability and convergence rate R(k, h) of the proposed two-step explicit/implicit scheme with varying space step h.Note: We take  = 0.25,  = 0.9 and k = h TA B L E 4 Analysis of stability and convergence rate R(k, h) of the two-step fourth-order algorithm with varying space step h.
Stability analysis and convergence order R(k, h) of the proposed technique with varying time step k.Analysis of stability and convergence rate R(k, h) of the two-step fourth-order algorithm with varying mesh size h.