Sinc-Galerkin method for solving the time fractional convection–diffusion equation with variable coefficients

In this paper, a new numerical algorithm for solving the time fractional convection–diffusion equation with variable coefficients is proposed. The time fractional derivative is estimated using the L1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$L_{1}$\end{document} formula, and the spatial derivative is discretized by the sinc-Galerkin method. The convergence analysis of this method is investigated in detail. The numerical solution is 2−α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$2-\alpha$\end{document} order accuracy in time and exponential rate of convergence in space. Finally, some numerical examples are given to show the effectiveness of the numerical scheme.


Introduction
In the last few decades, fractional differential equations have been widely applied in various fields of science and engineering to model many phenomena [1][2][3][4][5][6][7][8][9][10][11]. The current applications of fractional differential equations make it important to find efficient methods to solve them. Though many analytic methods have been proposed, the solutions of most fractional differential equations cannot be easily obtained due to the nonlocality and complexity of fractional differential operators. Hence, numerical algorithms for fractional differential equations have been studied extensively by many researchers, which mainly cover finite difference methods [12][13][14], finite element methods [15], spectral methods [16], local fractional series expansion methods [17], and other methods [18][19][20]. In particular, operator-based approach [20] was proposed and developed to deal with nonlinear fractional order differential equations.
In this paper, we consider the following time fractional convection-diffusion equation with variable coefficients: ∂u(x, t) ∂x + p 2 (x)u(x, t) = f (x, t), where (·) is the gamma function.
The convection-diffusion equations widely occur in the mathematical modeling of various physical phenomena that arise in oil reservoir simulations, transport of mass and energy, global weather production, and dispersion of chemicals in reactors. Time fractional convection-diffusion equations can be utilized to simulate time-related abnormal diffusion processes. It is a generalization of the classical convection diffusion by replacing the integer order time derivative with a fractional order time derivative.
Sinc methods have been proposed and studied by Stenger [21]. The sinc method has been increasingly used for solving ordinary differential equations [22][23][24][25], partial differential equations [26][27][28][29], and integral equations [30][31][32][33]; it not only has exponential convergence rate, but also can deal with singular problems well. In recent years, the sinc method has been also frequently employed for the numerical solution of fractional partial differential equations. In [34], Nagy applied the sinc-Chebyshev collocation method for numerical investigation of the time fractional nonlinear Klein-Gordon equation. In [35], Saadatmandi et al. proposed the sinc-Legendre collocation method for a class of fractional convection-diffusion equations with variable coefficients. In [36], Mao et al. developed the sinc-Chebyshev collocation method to solve a class of fractional diffusion-wave equations. In [37], Pirkhedri et al. adopted the sinc-Haar collocation method for the time fractional diffusion equation. In [38], a new reliable algorithm based on the sinc function is employed for the time fractional diffusion equation. In [18], Sweilam et al. investigated the numerical solution of the time fractional order telegraph equation by means of the sinc-Legendre collocation method. In [39], Jalilian et al. adopted an algorithm based on sinc basis functions for the numerical solution of the nonlinear fractional integro-differential equation of pantograph type. In this paper, we apply the sinc-Galerkin method to solve the time fractional convection-diffusion equation with variable coefficients.
The remainder of this paper is organized as follows. In Sect. 2, some preliminaries are presented. In Sect. 3, the time fractional derivative is discretized and the sinc-Galerkin method is applied to find the approximate solution at each step. In Sect. 4, convergence analysis is presented for this method and an exponential convergence is obtained as well. In Sect. 5, illustrative examples are carried out to justify the theoretical results.

Preliminaries
In this section, we recall notations and definitions of the sinc function and state some known lemmas that are necessary for this paper. These are discussed thoroughly in [21].
For all x ∈ R, the sinc function is basically defined as For h > 0, the translated sinc basis functions are given as Let f be a function defined on R, then for h > 0 the series is called the Whittaker cardinal expansion of f whenever this series converges. The properties of Whittaker cardinal expansions have been extensively studied in [40]. The properties are derived in the infinite strip D d of the complex plane where d > 0 To extend the approximation on R to the finite interval (0, 1), we consider the conformal map which maps the eye-shaped domain The basis functions on (0, 1) for z ∈ D E are taken to be the composite translated sinc functions The inverse map of w = φ(z) is Let ψ denote the inverse map of φ, so we define the range of φ -1 on R as For the uniform grid {kh} ∞ k=-∞ on R, the sinc points which correspond to these nodes are denoted by Definition 1 Let B(D E ) be the class of functions f which are analytic in D E and satisfy where ∂D E represents the boundary of D E .
, for some positive constants c and β, and if we select where C 1 depends only on f , d, and β.
and φ be a conformal map with constants β and C 2 so that by selecting h = √ πd/βN , then the sinc trapezoidal quadrature rule is

Lemma 3 ([41]) Let φ be the conformal one-to-one mapping of the simply connected domain D E onto D d .
Then In relations (5)- (7), h is the step size and x k is the sinc grid given by (3). For the assembly of the discrete system, it is convenient to define the following matrices: where δ (i) jk denotes the (j, k) the element of the matrix I (i) . The matrix I (0) is the m × m identity matrix. The matrix I (1) is the skew symmetric Toeplitz matrix and I (2) is the symmetric Toeplitz matrix.

Temporal discretization
In order to discretize equation (1) in time direction, let t n = nτ , n = 0, 1, . . . , where τ is the time step size. Denote u n (x) = u(x, t n ), u n . The time Caputo derivative is replaced by the L 1 -approximation, and the approximation order can be given in the following lemma.
Omitting the small error term R n,1 , we reach the following semi-discrete scheme of equation (1): subjected to the initial and boundary conditions

Space discretization
Next we consider space discretization to (12) by the sinc-Galerkin method. The approximation solution of (12) can be given by where S j (x) is the function S(j, h) • φ(x) for some fixed step size h. The unknown coefficients c n j in (13) are determined by orthogonalizing the residual Lu ng(x) with respect to the basis function {S k (x)} N k=-N . Taking the inner product on both sides of (12), we have where the inner product of functions f and η is defined by Equation (14) can be rewritten as Integrating by parts for the first two integral terms on the left-hand side of (15), we have where We can choose the weight function ω(x) = 1 φ (x) such that B 1 = B 2 = 0, then we apply the sinc quadrature rule (4) in Lemma 1 to (15) and obtain the following approximations: where Substituting (17)- (20) into (15) and replacing u n (x j ) with c n j , we obtain the discrete sinc-Galerkin system for determination of the unknown coefficients {c n j } N j=-N as follows: for k = -N, -N + 1, . . . , N .
We define the diagonal matrix as follows: Based on (8) and (22), system (21) can be represented by the following matrix form: where C and R are 2N + 1-vectors and M is (2N + 1) × (2N + 1) matrix as By solving the system of linear algebraic equations (23), we can obtain an approximate solution u n M of equation (12) from (13).

Convergence analysis
In this section, we show that the approximate solution u n M (x) converges to the exact solution u n (x) of (12) at an exponential rate. In order to establish a bound of |u n (x)u n M (x)|, we first need to get a bound of MU n -R 2 where where components u n (x j ), j = -N, . . . , N , are the values of the exact solution of (12) at the sinc points.

Lemma 5 Assume that (12) has a unique solution u n
then there exists a constant C 3 independent of N such that Proof For simplicity, we denote r k = (MU n -R) k for k = -N, . . . , N . By orthogonalizing the residual Lu ng(x) with respect to the basis function {S k (x)} N k=-N and using Theorem 4.4 of Lund et al. [41], we get where C 4 is a constant independent of N , then Utilizing Euclidean norm, we obtain Theorem 1 Let u n (x) be the exact solution of (12) and u n M (x) be its sinc approximation defined by (13), then under the assumptions of Lemma 2 and Lemma 5, there exists a constant C independent of N such that Proof Suppose that exact solution of (12) at sinc points x j (j = -N, . . . , N ) is denoted by Then, by making use of the triangular inequality, we have According to Lemma 1, there exists a constant C 5 independent of N such that sup x∈(0,1) The second term on the right-hand side of (25) satisfies We know that ( N j=-N |S j (x)| 2 ) 1 2 ≤ C 6 , where C 6 is a constant, then by using (22) and (24) Finally, by applying relations (25)-(28), we conclude sup x∈(0,1) where C = max{C 5 , C 2 C 6 M -1 2 }.

Numerical experiments
In this section, we give some numerical examples to confirm our analysis. In all the examples considered in this paper, we set parameters β = 1 and d = π Example 1 Consider the following fractional convection-diffusion equation: with the initial condition u(x, 0) = 0, 0 < x < 1, and the boundary conditions The exact solution of (29) is given by To illustrate the accuracy of our method, the maximum absolute error at sinc grid points is taken as 1+e ih . Furthermore, the temporal convergence order of presented method is denoted by for sufficiently small h. Table 1 shows the maximum absolute errors and temporal convergence order for α = 0.2, 0.4, 0.6, 0.8 and N = 64 with different values of time step size. Convergence curves of our method and the finite difference method with upwind strategy for the convection term are plotted in Fig. 1. It is clearly observed that the spatial errors appear in an exponential decay for our method. Figure 2 presents the graph of the numerical solution and the exact solution with α = 0.5, N = 128, and τ = 1 1000 . From these diagrams, it can be seen that our scheme gives a good approximation to the exact solution at mesh points.
Example 2 Consider the following fractional convection-diffusion equation: with the initial condition u(x, 0) = 0, 0 < x < 1, As the exact solution u(x, t) of (30) is unknown, therefore for each the maximum absolute error at the sinc grid points is taken as 1+e ih . Table 2 shows the maximum absolute errors and temporal convergence order for α = 0.8 and N = 128 with different values of and time step size. It is observed that the scheme generates temporal convergence order, which is consistent with our theoretical analysis. Convergence curves of our method and the finite difference method with upwind strategy for the convection term are plotted in Figs. 3 and 4. The two figures clearly show that the proposed method converges at the exponential rate with increasing N in space. Numerical solution profiles for = 2 -8 and = 2 -12 are plotted in Figs. 5 and 6. The two figures show that the boundary layer is located on the right-hand side of the rectangular domain.

Conclusion
In this paper, we have developed and analyzed an efficient numerical scheme for the time fractional convection-diffusion equations with variable coefficients. The problem is reduced to the solution of a system of linear algebraic equations at each time step. The convergence analysis of the proposed method is proven, and it is shown that the sinc-Galerkin method converges to the solution at the exponential rate in space. Two examples are tested and the obtained numerical results confirm the theoretical analysis.