Skip to main content

Theory and Modern Applications

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

Abstract

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 \(L_{1}\) 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-\alpha\) 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.

1 Introduction

In the last few decades, fractional differential equations have been widely applied in various fields of science and engineering to model many phenomena [111]. 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 [1214], finite element methods [15], spectral methods[16], local fractional series expansion methods[17], and other methods[1820]. 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:

$$ \begin{aligned} &\frac{\partial ^{\alpha }u(x,t)}{\partial t^{\alpha }}-\epsilon \frac{\partial ^{2}u(x,t)}{\partial x^{2}}+p_{1}(x) \frac{\partial u(x,t)}{\partial x}+p_{2}(x)u(x,t)=f(x,t), \\ & 0< x < 1, \qquad 0< t \leq T, \end{aligned} $$
(1)

with the initial condition

$$ u(x,0)=h(x), \quad 0< x< 1, $$

and the boundary conditions

$$ u(a,t)=u(b,t)=0,\quad 0< t\leq T, $$

where \(0<\alpha < 1 \), ϵ is a known positive constant, \(p_{1}(x)\in C^{1}([0,1])\), \(p_{2}(x),h(x)\in C([0,1])\), and \(f(x,t)\in C(\Omega )\), \(\Omega =(0,1)\times (0,T]\). The time fractional derivative is defined in the Caputo sense as follows:

$$ \frac{\partial ^{\alpha }u(x,t)}{\partial t^{\alpha }}= \frac{1}{\Gamma (1-\alpha )} \int _{0}^{t} \frac{\partial u(x,\xi )}{\partial \xi } \frac{d\xi }{(t-\xi )^{\alpha }}, $$

where \(\Gamma (\cdot )\) 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[2225], partial differential equations[2629], and integral equations[3033]; 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.

2 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\in \mathbb{R}\), the sinc function is basically defined as

$$ \operatorname{sinc}(x)=\textstyle\begin{cases} \frac{\sin (\pi x)}{\pi x}, & x\neq 0, \\ 1, & x = 0. \end{cases} $$

For \(h>0\), the translated sinc basis functions are given as

$$ S(j,h) (x)=\operatorname{sinc} \biggl(\frac{x-jh}{h} \biggr), \quad j=0,\pm 1,\pm 2,\ldots $$

Let f be a function defined on \(\mathbb{R}\), then for \(h>0\) the series

$$ C(f,h) (x)=\sum_{j=-\infty }^{\infty } f(jh) S(j,h) (x) $$

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\)

$$ D_{d} = \biggl\{ w=u+iv: \vert v \vert < d \leq \frac{\pi }{2} \biggr\} . $$

To extend the approximation on \(\mathbb{R}\) to the finite interval \((0,1)\), we consider the conformal map

$$ \phi (z)=\log \biggl(\frac{z}{1-z} \biggr), $$

which maps the eye-shaped domain

$$ D_{E}= \biggl\{ z=x+iy: \biggl\vert \arg \biggl( \frac{z}{1-z} \biggr) \biggr\vert < d \leq \frac{\pi }{2} \biggr\} $$

onto the infinite strip \(D_{d}\).

The basis functions on \((0,1)\) for \(z\in D_{E}\) are taken to be the composite translated sinc functions

$$ S_{j}(z)=S(j,h)\circ \phi (z)=\operatorname{sinc} \biggl(\frac{\phi (z)-jh}{h} \biggr), \quad j=0, \pm 1, \pm 2,\ldots . $$
(2)

The inverse map of \(w=\phi (z)\) is

$$ z=\phi ^{-1}(w)=\frac{e^{w}}{1+e^{w}}. $$

Let ψ denote the inverse map of ϕ, so we define the range of \(\phi ^{-1}\) on \(\mathbb{R}\) as

$$ (0,1)=\bigl\{ \psi (u)=\phi ^{-1}(u)\in D_{E}: -\infty < u < \infty \bigr\} . $$

For the uniform grid \(\{kh\}^{\infty }_{k=-\infty }\) on \(\mathbb{R}\), the sinc points which correspond to these nodes are denoted by

$$ x_{k}=\psi (kh)=\frac{e^{kh}}{1+e^{kh}},\quad k=0,\pm 1, \pm 2, \ldots . $$
(3)

Definition 1

Let \(B(D_{E})\) be the class of functions f which are analytic in \(D_{E}\) and satisfy

$$ \int _{\psi (t+\Sigma )} \bigl\vert f(z) \bigr\vert \,\mathrm{d}z \rightarrow 0, \quad \mbox{as } t \rightarrow \pm \infty , $$

where \(\Sigma = \{\mathrm{i}\eta : |\eta |< d \leq \frac{\pi }{2} \}\), and satisfy

$$ N(f)= \int _{\partial D_{E}} \bigl\vert f(z)\,\mathrm{d}z \bigr\vert < \infty , $$

where \(\partial D_{E}\) represents the boundary of \(D_{E}\).

Definition 2

Let \(L_{\beta }(D_{E})\) be the set of all analytic function f in \(D_{E}\), for which there exists a constant C such that

$$ \bigl\vert f(z) \bigr\vert \leq C \frac{ \vert \mathrm{e}^{\phi (z)} \vert ^{\beta }}{(1+ \vert \mathrm{e}^{\phi (z)} \vert )^{2\beta }}, \quad z \in D_{E}, 0< \beta \leq 1. $$

Lemma 1

([40])

If \(f \phi ' \in B(D_{E}) \), then for all \(x \in (0,1)\)

$$ \Biggl\vert f(x)-\sum_{j=-N}^{N}f(x_{j})S_{j}(x) \Biggr\vert \leq \frac{2N(f\phi ')}{\pi d}\mathrm{e}^{-\pi d/h}. $$

Moreover, if \(|f(x)|\leq ce^{-\beta |\phi (x)|} \), \(x\in (0,1)\), for some positive constants c and β, and if we select \(h=\sqrt{\pi d/\beta N}\), then

$$ \mathop{\sup }_{x\in (0,1)} \Biggl\vert f(x)- \biggl( \frac{\mathrm{d}}{\mathrm{d}x} \biggr)^{l}\sum_{j=-N}^{N}f(x_{j})S_{j}(x) \Biggr\vert \leq C_{1}N^{(l+1)/2}\mathrm{e}^{-\sqrt{\pi d\beta N}}, $$

where \(C_{1}\)depends only on f, d, and β.

Lemma 2

([41])

Let \(F\in B(D_{E})\)and ϕ be a conformal map with constants β and \(C_{2}\)so that

$$ \biggl\vert \frac{F(x)}{\phi '(x)} \biggr\vert \leq C_{2}e^{-\beta \vert \phi (x) \vert }, \quad x \in (0,1), $$

by selecting \(h=\sqrt{\pi d/\beta N}\), then the sinc trapezoidal quadrature rule is

$$ \int _{0}^{1}F(x)\,dx=h\sum _{j=-N}^{N} \frac{F(x_{j})}{\phi '(x_{j})}+ o \bigl(e^{-\sqrt{\pi d\beta N}} \bigr). $$
(4)

Lemma 3

([41])

Let ϕ be the conformal one-to-one mapping of the simply connected domain \(D_{E}\)onto \(D_{d}\). Then

$$\begin{aligned}& \delta _{jk}^{(0)}=\bigl[S(j,h)\circ \phi (x)\bigr]|_{x=x_{k}}=\textstyle\begin{cases} 1, & j=k, \\ 0, & j\neq k, \end{cases}\displaystyle \end{aligned}$$
(5)
$$\begin{aligned}& \delta _{jk}^{(1)}=h \frac{d}{d\phi }\bigl[S(j,h)\circ \phi (x)\bigr]|_{x=x_{k}}= \textstyle\begin{cases} 0, & j=k, \\ \frac{(-1)^{k-j}}{k-j}, & j\neq k, \end{cases}\displaystyle \end{aligned}$$
(6)
$$\begin{aligned}& \delta _{jk}^{(2)}=h^{2} \frac{d^{2}}{d\phi ^{2}}\bigl[S(j,h)\circ \phi (x)\bigr]|_{x=x_{k}}= \textstyle\begin{cases} \frac{-\pi ^{2}}{3}, & j=k, \\ \frac{-2(-1)^{k-j}}{(k-j)^{2}}, & j\neq k. \end{cases}\displaystyle \end{aligned}$$
(7)

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:

$$ I^{(i)}=\bigl[\delta _{jk}^{(i)} \bigr], \quad i=0,1,2, $$
(8)

where \(\delta _{jk}^{(i)}\) denotes the \((j,k)\) the element of the matrix \(I^{(i)}\). The matrix \(I^{(0)}\) is the \(m\times m\) identity matrix. The matrix \(I^{(1)}\) is the skew symmetric Toeplitz matrix and \(I^{(2)}\) is the symmetric Toeplitz matrix.

3 Derivation of numerical method

3.1 Temporal discretization

In order to discretize equation (1) in time direction, let \(t_{n}=n\tau \), \(n=0,1,\ldots \) , where τ is the time step size. Denote \(u^{n}(x)=u(x,t_{n})\), \(u^{n}_{x}(x)= \frac{\partial u(x,t_{n})}{\partial x}\), \(u^{n}_{xx}(x)= \frac{\partial ^{2} u(x,t_{n})}{\partial x^{2}}\), and \(f^{n}(x)=f(x,t_{n})\). The time Caputo derivative is replaced by the \(L_{1}\)-approximation, and the approximation order can be given in the following lemma.

Lemma 4

([42])

Suppose \(0< \alpha < 1\)and \(f(t)\in C^{2}[0,t_{k}]\), it holds that

$$ \begin{aligned} & \Biggl\vert \frac{1}{\Gamma (1-\alpha )} \int ^{t_{k}}_{0} \frac{f'(t)}{(t_{k}-t)^{\alpha }}\,dt \\ &\qquad {}- \frac{\tau ^{-\alpha }}{\Gamma (2-\alpha )} \Biggl[b_{0}f(t_{k})- \sum^{k-1}_{m=1}(b_{k-m-1}-b_{k-m})f(t_{m})-b_{k-1}f(t_{0}) \Biggr] \Biggr\vert \\ &\quad \leq \frac{1}{\Gamma (2-\alpha )} \biggl[\frac{1-\alpha }{12}+ \frac{2^{2-\alpha }}{2-\alpha }-\bigl(1+2^{-\alpha }\bigr) \biggr]\max _{0 \leq t\leq t_{k}} \bigl\vert f''(t) \bigr\vert \tau ^{2-\alpha }, \end{aligned} $$
(9)

where \(b_{m}=(m+1)^{1-\alpha }-m^{1-\alpha }\), \(m\geq 0\).

Based on Lemma 4, we can approximate the Caputo fractional derivative as follows:

$$\begin{aligned} \frac{\partial ^{\alpha }u^{n}(x)}{\partial t^{\alpha }}=\mu \Biggl[b_{0}u^{n}(x)- \sum^{n-1}_{m=1}(b_{n-m-1}-b_{n-m})u^{m}(x)-b_{n-1}u^{0}(x) \Biggr]+O\bigl(\tau ^{2-\alpha }\bigr), \end{aligned}$$
(10)

where \(\mu =\frac{\tau ^{-\alpha }}{\Gamma (2-\alpha )}\).

Substituting (10) into (1), we have

$$ \begin{aligned} &\mu \Biggl[b_{0}u^{n}(x)- \sum^{n-1}_{m=1}(b_{n-m-1}-b_{n-m})u^{m}(x)-b_{n-1}u^{0}(x) \Biggr] \\ &\quad =\varepsilon u^{n}_{xx}(x)-p_{1}(x)u^{n}_{x}(x)-p_{2}(x)u^{n}(x)+f^{n}(x)+R_{n,1}, \end{aligned} $$
(11)

where \(R_{n,1}=O(\tau ^{2-\alpha })\).

Omitting the small error term \(R_{n,1}\), we reach the following semi-discrete scheme of equation (1):

$$ Lu^{n}=-\epsilon u^{n}_{xx}(x)+p_{1}(x)u^{n}_{x}(x)+ \bigl(\mu +p_{2}(x)\bigr)u^{n}(x)=g(x), \quad 0< x< 1, $$
(12)

subjected to the initial and boundary conditions

$$\begin{aligned}& u^{0}(x)=h(x), \\& u^{n}(0)=u^{n}(1)=0, \end{aligned}$$

where

$$ g(x)=f^{n}(x)+\mu \Biggl[\sum^{n-1}_{m=1}(b_{n-m-1}-b_{n-m})u^{m}(x)+b_{n-1}u^{0}(x) \Biggr]. $$

3.2 Space discretization

Next we consider space discretization to (12) by the sinc-Galerkin method. The approximation solution of (12) can be given by

$$ u^{n}(x)\approx u_{M}^{n}= \sum_{j=-N}^{N}c_{j}^{n}S_{j}(x), \quad M=2N+1, $$
(13)

where \(S_{j}(x) \) is the function \(S(j,h)\circ \phi (x)\) for some fixed step size h. The unknown coefficients \(c_{j}^{n}\) in (13) are determined by orthogonalizing the residual \(Lu^{n}-g(x)\) with respect to the basis function \(\{S_{k}(x)\}_{k=-N}^{N}\). Taking the inner product on both sides of (12), we have

$$ \begin{aligned} & \bigl( -\epsilon u^{n}_{xx}(x),S_{k}(x) \bigr)+ \bigl(p_{1}(x)u^{n}_{x}(x),S_{k}(x) \bigr) + \bigl(\bigl(\mu +p_{2}(x)\bigr)u^{n}(x),S_{k}(x) \bigr) \\ &\quad = \bigl(g(x),S_{k}(x) \bigr), \end{aligned} $$
(14)

where the inner product of functions f and η is defined by

$$ ( f,\eta )= \int _{0}^{1}f(x)\eta (x)\omega (x)\,dx. $$

Equation (14) can be rewritten as

$$ \begin{aligned} & { \int _{0}^{1}p_{1}(x)u^{n}_{x}(x)S_{k}(x) \omega (x)\,dx}- { \int _{0}^{1}\epsilon u^{n}_{xx}(x)S_{k}(x) \omega (x)\,dx} \\ &\quad {}+ { \int _{0}^{1}\bigl(\mu +p_{2}(x) \bigr)u^{n}(x)S_{k}(x)\omega (x)\,dx}= \int _{0}^{1}g(x)S_{k}(x)\omega (x) \,dx. \end{aligned} $$
(15)

Integrating by parts for the first two integral terms on the left-hand side of (15), we have

$$ \begin{aligned} & { \int _{0}^{1}p_{1}(x)u^{n}_{x}(x)S_{k}(x) \omega (x)dx}=B_{1}- \int _{0}^{1}u^{n}(x) \bigl( p_{1}(x)S_{k}(x)\omega (x) \bigr)'dx, \\ & { \int _{0}^{1}\epsilon u^{n}_{xx}(x)S_{k}(x) \omega (x)dx}=B_{2}- \int _{0}^{1}u^{n}(x) \bigl( \epsilon S_{k}(x)\omega (x) \bigr)''dx, \end{aligned} $$
(16)

where

$$ \begin{aligned} &B_{1}=\bigl\{ u^{n}(x) p_{1}(x)S_{k}(x)\omega (x)\bigr\} |_{0}^{1}, \\ &B_{2}=\bigl\{ u_{x}^{n}(x)\epsilon S_{k}(x)\omega (x)-u^{n}(x) \bigl(\epsilon S_{k}(x) \omega (x)\bigr)'\bigr\} |_{0}^{1}. \end{aligned} $$

We can choose the weight function \(\omega (x)=\frac{1}{\phi '(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:

$$\begin{aligned}& { \int _{0}^{1}p_{1}(x)u^{n}_{x}(x)S_{k}(x) \omega (x)\,dx} \approx -h\sum_{j=-N}^{j=N} \sum_{l=0}^{1} \frac{u^{n}(x_{j})}{\phi '(x_{j})h^{l}} \delta _{kj}^{(l)}g_{1,l}(x_{j}), \end{aligned}$$
(17)
$$\begin{aligned}& { \int _{0}^{1}\epsilon u^{n}_{xx}(x)S_{k}(x) \omega (x)\,dx} \approx -h\sum_{j=-N}^{j=N} \sum_{l=0}^{2} \frac{u^{n}(x_{j})}{\phi '(x_{j})h^{l}} \delta _{kj}^{(l)}g_{2,l}(x_{j}), \end{aligned}$$
(18)
$$\begin{aligned}& { \int _{0}^{1}\bigl(\mu +p_{2}(x) \bigr)u^{n}(x)S_{k}(x)\omega (x)\,dx} \approx h \frac{(\mu +p_{2}(x_{k}))u^{n}(x_{k})\omega (x_{k})}{\phi '(x_{k})}, \end{aligned}$$
(19)
$$\begin{aligned}& { \int _{0}^{1}g(x)S_{k}(x)\omega (x) \,dx} \approx h \frac{g(x_{k})\omega (x_{k})}{\phi '(x_{k})}, \end{aligned}$$
(20)

where

$$ \begin{aligned} &g_{2,2}=-\epsilon \omega (x) \bigl(\phi '(x)\bigr)^{2}, \qquad g_{2,1}=- \epsilon \omega (x)\phi ''(x)-2\epsilon \omega '(x)\phi '(x), \qquad g_{2,0}=- \epsilon \omega ''(x), \\ &g_{1,1}=p_{1}(x)\omega (x)\phi '(x), \qquad g_{1,0}=\bigl(p_{1}(x)\omega (x)\bigr)'. \end{aligned} $$

Substituting (17)–(20) into (15) and replacing \(u^{n}(x_{j})\) with \(c_{j}^{n}\), we obtain the discrete sinc-Galerkin system for determination of the unknown coefficients \(\{c_{j}^{n}\}_{j=-N}^{N}\) as follows:

$$ \begin{aligned} &\sum _{j=-N}^{N}\sum_{l=0}^{2} \frac{\delta _{kj}^{(l)}g_{2,l}(x_{j})}{\phi '(x_{j})h^{l}}c_{j}^{n}- \sum _{j=-N}^{N}\sum_{l=0}^{1} \frac{\delta _{kj}^{(l)}g_{1,l}(x_{j})}{\phi '(x_{j})h^{l}}c_{j}^{n}+ \frac{(\mu +p_{2}(x_{k}))\omega (x_{k})}{\phi '(x_{k})}c_{k}^{n} \\ &\quad = \frac{f^{n}(x_{k})\omega (x_{k})}{\phi '(x_{k})}+\mu \Biggl[\sum^{n-1}_{m=1}(b_{n-m-1}-b_{n-m}) \frac{\omega (x_{k})}{\phi '(x_{k})}c_{k}^{m}+b_{n-1} \frac{\omega (x_{k})}{\phi '(x_{k})}c_{k}^{0} \Biggr] \end{aligned} $$
(21)

for \(k=-N, -N+1, \ldots , N\).

We define the diagonal matrix as follows:

$$ \mathbf{D}\bigl(g(x)\bigr)_{ij}= \textstyle\begin{cases} g(x_{i}), & i=j, \\ 0, & i\neq j. \end{cases} $$
(22)

Based on (8) and (22), system (21) can be represented by the following matrix form:

$$ \mathbf{M}\mathbf{C^{n}}=\mathbf{R}, $$
(23)

where C and R are \(2N+1\)-vectors and M is \((2N+1)\times (2N+1) \) matrix as

$$\begin{aligned}& \mathbf{C}^{n}= \bigl(c_{-N}^{n}, c_{-N+1}^{n}, \ldots , c_{N}^{n} \bigr)^{T}, \\& \mathbf{R}=\mathbf{D} \biggl(\frac{1}{(\phi '(x))^{2}} \biggr)\mathbf{F}^{n} \\& \hphantom{\mathbf{R}=}{}+ \mu \Biggl[\sum^{n-1}_{m=1}(b_{n-m-1}-b_{n-m}) \mathbf{D} \biggl( \frac{1}{(\phi '(x))^{2}} \biggr)\mathbf{C}^{m} +b_{n-1}\mathbf{D} \biggl(\frac{1}{(\phi '(x))^{2}} \biggr) \mathbf{C}^{0} \Biggr], \\& \mathbf{M}=\epsilon \mathbf{B}-\frac{1}{h}\mathbf{I}^{(1)} \mathbf{D} \biggl(\frac{p_{1}(x)}{\phi '(x)} \biggr) +\mathbf{D} \biggl( \frac{-1}{\phi '(x)} \biggl(\frac{p_{1}(x)}{\phi '(x)} \biggr)'+ \frac{d}{(\phi '(x))^{2}} \biggr), \end{aligned}$$

in which

$$\begin{aligned}& \mathbf{B}=\frac{-1}{h^{2}}\mathbf{I}^{(2)}+ \frac{1}{h}\mathbf{I}^{(1)} \mathbf{D} \biggl( \frac{\phi ''(x)}{(\phi '(x))^{2}} \biggr) +\mathbf{D} \biggl(\frac{-1}{\phi '(x)}\biggl( \frac{1}{\phi '(x)}\biggr)'' \biggr), \\& \mathbf{F}^{n}=\bigl(\mathbf{F}^{n}(x_{-N}), \mathbf{F}^{n}(x_{-N+1}), \ldots , \mathbf{F}^{n}(x_{N}) \bigr)^{T}. \end{aligned}$$

By solving the system of linear algebraic equations (23), we can obtain an approximate solution \(u_{M}^{n}\) of equation (12) from (13).

4 Convergence analysis

In this section, we show that the approximate solution \(u_{M}^{n}(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_{M}^{n}(x)|\), we first need to get a bound of \(\|MU^{n}-R\|_{2}\) where

$$ U^{n}= \bigl(u^{n}(x_{-N}),u^{n}(x_{-N+1}), \ldots ,u^{n}(x_{N}) \bigr)^{\mathrm{T}}, $$

where components \(u^{n}(x_{j})\), \(j=-N,\ldots ,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}(x)\in L_{\beta }(D_{E})\). If \(G \in L_{\beta }(D_{E})\)where

$$ G=\frac{g(x)}{\phi '(x)},\frac{\mu +p_{2}(x)}{\phi '(x)}, \frac{\phi ''(x)}{(\phi '(x))^{2}}, \frac{p_{1}(x)}{\phi '(x)}, - \frac{1}{\phi '(x)} \biggl(\frac{1}{\phi '(x)} \biggr)'', - \frac{1}{\phi '(x)} \biggl( \frac{p_{1}(x)}{\phi '(x)} \biggr)', $$

then there exists a constant \(C_{3}\)independent of N such that

$$ \bigl\Vert MU^{n}-R \bigr\Vert _{2} \leq C_{3}N^{\frac{3}{2}}\mathrm{e}^{-\sqrt{\pi d \beta N}}. $$
(24)

Proof

For simplicity, we denote \(r_{k}= (MU^{n}-R )_{k}\) for \(k=-N,\ldots ,N\). By orthogonalizing the residual \(Lu^{n}-g(x)\) with respect to the basis function \(\{S_{k}(x)\}_{k=-N}^{N}\) and using Theorem 4.4 of Lund et al. [41], we get

$$ 0=h \int _{0}^{1}\frac{(Lu^{n}(x)-g(x))S_{k}(x)}{\phi '(x)} \, \mathrm{d}x=r_{k}+C_{4}N\mathrm{e}^{-\sqrt{\pi d\beta N}}, $$

where \(C_{4}\) is a constant independent of N, then

$$ \vert r_{k} \vert \leq C_{4}N\mathrm{e}^{-\sqrt{\pi \beta dN}}. $$

Utilizing Euclidean norm, we obtain

$$ \bigl\Vert MU^{n}-R \bigr\Vert _{2}= \Biggl(\sum _{k=-N}^{k=N} \vert r_{k} \vert ^{2} \Biggr)^{ \frac{1}{2}} \leq C_{3}N^{\frac{3}{2}} \mathrm{e}^{-\sqrt{\pi d\beta N}}. $$

 □

Theorem 1

Let \(u^{n}(x)\)be the exact solution of (12) and \(u_{M}^{n}(x)\)be its sinc approximation defined by (13), then under the assumptions of Lemma 2and Lemma 5, there exists a constant C independent of N such that

$$ \sup_{x\in (0,1)} \bigl\vert u^{n}(x)-u^{n}_{M}(x) \bigr\vert \leq C N^{ \frac{3}{2}}\mathrm{e}^{-\sqrt{\pi d\beta N}}. $$

Proof

Suppose that exact solution of (12) at sinc points \(x_{j}\) (\(j=-N,\ldots ,N\)) is denoted by \(U_{N}^{n}(x)\) defined by

$$ U_{N}^{n}(x)=\sum_{j=-N}^{N}u^{n}(x_{j})S_{j}(x). $$

Then, by making use of the triangular inequality, we have

$$ \bigl\vert u^{n}(x)-u^{n}_{M}(x) \bigr\vert \leq \bigl\vert u^{n}(x)-U_{N}^{n}(x) \bigr\vert + \bigl\vert U_{N}^{n}(x)-u^{n}_{M}(x) \bigr\vert . $$
(25)

According to Lemma 1, there exists a constant \(C_{5}\) independent of N such that

$$ \sup_{x\in (0,1)} \bigl\vert u^{n}(x)-U_{N}^{n}(x) \bigr\vert \leq C_{5} N^{ \frac{1}{2}}\mathrm{e}^{-\sqrt{\pi d\beta N}}. $$
(26)

The second term on the right-hand side of (25) satisfies

$$ \begin{aligned} \bigl\vert U_{N}^{n}(x)-u^{n}_{M}(x) \bigr\vert &\leq \Biggl\vert {\sum_{j=-N}^{N}} \bigl(u^{n}(x_{j})-c_{j}^{n} \bigr)S_{j}(x) \Biggr\vert \\ & \leq \Biggl( {\sum_{j=-N}^{N}} \bigl\vert u^{n}(x_{j})-c_{j}^{n} \bigr\vert ^{2} \Biggr)^{\frac{1}{2}} \Biggl( {\sum _{j=-N}^{N}} \bigl\vert S_{j}(x) \bigr\vert ^{2} \Biggr)^{\frac{1}{2}}. \end{aligned} $$
(27)

We know that \(({\sum_{j=-N}^{N}}|S_{j}(x)|^{2} )^{\frac{1}{2}}\leq C_{6}\), where \(C_{6}\) is a constant, then by using (22) and (24) we get

$$ \begin{aligned} & \Biggl( {\sum _{j=-N}^{N}} \bigl\vert u^{n}(x_{j})-c_{j}^{n} \bigr\vert ^{2} \Biggr)^{\frac{1}{2}} \Biggl( {\sum _{j=-N}^{N}} \bigl\vert S_{j}(x) \bigr\vert ^{2} \Biggr)^{\frac{1}{2}} \\ &\quad \leq C_{6} \bigl\Vert U^{n}-C^{n} \bigr\Vert _{2} \\ &\quad =C_{6} \bigl\Vert M^{-1}\bigl(MU^{n}-R \bigr) \bigr\Vert _{2} \\ &\quad \leq C_{2}C_{6} \bigl\Vert M^{-1} \bigr\Vert _{2}N^{\frac{3}{2}}\mathrm{e}^{-\sqrt{\pi d \beta N}}. \end{aligned} $$
(28)

Finally, by applying relations (25)–(28), we conclude

$$ \sup_{x\in (0,1)} \bigl\vert u^{n}(x)-u^{n}_{M}(x) \bigr\vert \leq C N^{ \frac{3}{2}}\mathrm{e}^{-\sqrt{\pi d\beta N}}, $$

where \(C=\max \{C_{5},C_{2}C_{6}\|M^{-1}\|_{2} \}\). □

5 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 \(\beta =1\) and \(d=\frac{\pi }{2}\) which yield \(h=\frac{\pi }{\sqrt{2N}}\).

Example 1

Consider the following fractional convection–diffusion equation:

$$ \frac{\partial ^{\alpha }u(x,t)}{\partial t^{\alpha }}- \frac{\partial ^{2}u(x,t)}{\partial x^{2}}+x \frac{\partial u(x,t)}{\partial x}=f(x,t), \quad 0< x < 1, 0< t \leq 1, $$
(29)

with the initial condition

$$ u(x,0)=0, \quad 0< x< 1, $$

and the boundary conditions

$$ u(0,t)=u(1,t)=0, \quad 0< t\leq 1, $$

where

$$ f(x,t)=\frac{2t^{2-\alpha }}{\Gamma (3-\alpha )}\bigl(x^{2}-x^{3}\bigr)+ t^{2}\bigl(2x^{2}-2x^{3}+6x-2\bigr). $$

The exact solution of (29) is given by

$$ u(x,t)=t^{2}\bigl(x^{2}-x^{3}\bigr). $$

To illustrate the accuracy of our method, the maximum absolute error at sinc grid points is taken as

$$ e_{\infty }(h,\tau )=\max_{i,j} \bigl\vert u(x_{i},t_{j})-u_{M}(x_{i},t_{j}) \bigr\vert , $$

where \(x_{i}=\frac{e^{ih}}{1+e^{ih}}\).

Furthermore, the temporal convergence order of presented method is denoted by

$$ rate_{1}=\log _{2} \biggl( \frac{e_{\infty }(h,2\tau )}{e_{\infty }(h,\tau )} \biggr) $$

for sufficiently small h.

Table 1 shows the maximum absolute errors and temporal convergence order for \(\alpha =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 \(\alpha =0.5\), \(N=128\), and \(\tau =\frac{1}{1000}\). From these diagrams, it can be seen that our scheme gives a good approximation to the exact solution at mesh points.

Figure 1
figure 1

Convergence of the sinc-Galerkin and FDM-upwind methods for various values of N with \(\alpha =0.8\), \(\tau =1/4000\)

Figure 2
figure 2

Comparison of numerical solution and analytical solution for \(\alpha =0.5\), \(N=128\), \(\tau =1/1000\)

Table 1 The maximum absolute errors and temporal convergence orders with \(N=64\)

Example 2

Consider the following fractional convection–diffusion equation:

$$ \begin{aligned} &\frac{\partial ^{\alpha }u(x,t)}{\partial t^{\alpha }}-\epsilon \frac{\partial ^{2}u(x,t)}{\partial x^{2}}+ \bigl(2-x^{2}\bigr) \frac{\partial u(x,t)}{\partial x}+xu(x,t) \\ &\quad =f(x,t) \quad 0< x < 1, 0< t \leq 1, \end{aligned}$$
(30)

with the initial condition

$$ u(x,0)=0, \quad 0< x< 1, $$

and the boundary conditions

$$ u(0,t)=u(1,t)=0, \quad 0< t\leq 1, $$

where

$$ f(x,t)=10t^{2}\mathrm{e}^{-t}x(1-x). $$

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

$$ e_{\infty }(h,\tau )=\max_{i,j} \bigl\vert u_{M}(x_{i},t_{j})-u_{2M}(x_{i},t_{j}) \bigr\vert , $$

where \(x_{i}=\frac{e^{ih}}{1+e^{ih}}\).

Table 2 shows the maximum absolute errors and temporal convergence order for \(\alpha =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 \(\epsilon =2^{-8}\) and \(\epsilon =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.

Figure 3
figure 3

Convergence of the sinc-Galerkin and FDM-upwind methods for various values of N with \(\alpha =0.5\), \(\epsilon =2^{-5}\), \(\tau =1/500\)

Figure 4
figure 4

Convergence of the sinc-Galerkin and FDM-upwind methods for various values of N with \(\alpha =0.5\), \(\epsilon =2^{-12}\), \(\tau =1/500\)

Figure 5
figure 5

Numerical solutions with \(\alpha =0.2\), \(\epsilon =2^{-8}\), \(N=32\), \(\tau =1/100\)

Figure 6
figure 6

Numerical solutions with \(\alpha =0.2\), \(\epsilon =2^{-12}\), \(N=32\), \(\tau =1/100\)

Table 2 The maximum absolute errors and temporal convergence orders with \(\alpha =0.8\) and \(N=128\)

6 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.

References

  1. Podlubny, I.: Fractional Differential Equations. Academic Press, New York (1999)

    MATH  Google Scholar 

  2. Raberto, M., Scalas, E., Mainardi, F.: Waiting-times and returns in high-frequency financial data: an empirical study. Physica A 314, 749–755 (2002)

    Article  MATH  Google Scholar 

  3. Koeller, R.C.: Application of fractional calculus to the theory of viscoelasticity. J. Appl. Mech. 51, 229–307 (1984)

    Article  MathSciNet  MATH  Google Scholar 

  4. Meerschaert, M.M., Zhang, Y., Baeumer, B.: Particle tracking for fractional diffusion with two time scales. Comput. Math. Appl. 59, 1078–1086 (2010)

    Article  MathSciNet  MATH  Google Scholar 

  5. Jiang, X., Xu, M., Qi, H.: The fractional diffusion model with an absorption term and modified Fick’s law for non-local transport processes. Nonlinear Anal. 11, 262–269 (2011)

    Article  MathSciNet  MATH  Google Scholar 

  6. Yang, X., Tenreiro Machado, J.A.: A new fractal nonlinear Burgers equation arising in the acoustic signals propagation. Math. Methods Appl. Sci. 42, 7539–7544 (2019)

    Article  MathSciNet  MATH  Google Scholar 

  7. Yang, X., Feng, Y., Cattani, C., Inc, M.: Fundamental solutions of anomalous diffusion equations with the decay exponential kernel. Math. Methods Appl. Sci. 42, 4054–4060 (2019)

    Article  MathSciNet  MATH  Google Scholar 

  8. Liu, J., Yang, X., Feng, Y.: On integrability of the time fractional nonlinear heat conduction equation. J. Geom. Phys. 144, 190–198 (2019)

    Article  MathSciNet  MATH  Google Scholar 

  9. Liu, J., Yang, X., Feng, Y., Zhang, H.: Analysis of the time fractional nonlinear diffusion equation from diffusion process. J. Appl. Anal. Comput. 10, 1060–1072 (2020)

    MathSciNet  Google Scholar 

  10. Liu, J., Yang, X., Feng, Y., Iqbal, M.: Group analysis to the time fractional nonlinear wave equation. Int. J. Math. 31, 20500299 (2020)

    MathSciNet  Google Scholar 

  11. Yang, X.: General Fractional Derivatives: Theory, Methods and Applications. CRC Press, New York (2019)

    Book  MATH  Google Scholar 

  12. Mohebbi, A., Abbaszadeh, M.: Compact finite difference scheme for the solution of time fractional advection-dispersion equation. Numer. Algorithms 63, 431–452 (2013)

    Article  MathSciNet  MATH  Google Scholar 

  13. Sousa, E., Li, C.: A weighted finite difference method for the fractional diffusion equation based on the Riemann–Liouville derivative. Appl. Numer. Math. 90, 22–37 (2015)

    Article  MathSciNet  MATH  Google Scholar 

  14. Guo, B.L., Xu, Q., Yin, Z.: Implicit finite difference method for fractional percolation equation with Dirichlet and fractional boundary conditions. Appl. Math. Mech. 37, 403–416 (2016)

    Article  MathSciNet  MATH  Google Scholar 

  15. Jin, B., Lazarov, R., Pasciak, J., Zhou, Z.: Error analysis of semidiscrete finite element methods for inhomogeneous time-fractional diffusion. IMA J. Numer. Anal. 35, 561–587 (2015)

    Article  MathSciNet  MATH  Google Scholar 

  16. Bhrawy, A.H., Doha, E.H., Baleanud, D., Ezz-Eldien, S.S.: A spectral tau algorithm based on Jacobi operational matrix for numerical solution of time fractional diffusion-wave equations. J. Comput. Phys. 93, 142–156 (2015)

    Article  MathSciNet  MATH  Google Scholar 

  17. Yang, A., Yang, X., Li, Z.: Local fractional series expansion method for solving wave and diffusion equations on Cantor sets. Abstr. Appl. Anal. 2013, 351057 (2013)

    MathSciNet  MATH  Google Scholar 

  18. Sweilam, N.H., Nagy, A.M., El-Sayed, A.A.: Solving time-fractional order telegraph equation via sinc-Legendre collocation method. Mediterr. J. Math. 13, 5119–5133 (2016)

    Article  MathSciNet  MATH  Google Scholar 

  19. Yang, X., Tenreiro Machado, J.A., Srivastava, H.M.: A new numerical technique for solving the local fractional diffusion equation:two-dimensional extended differential transform approach. Appl. Math. Comput. 274, 143–151 (2016)

    MathSciNet  MATH  Google Scholar 

  20. Navickas, Z., Telksnys, T., Marcinkevicius, R., Ragulskis, M.: Operator-based approach for the construction of analytical soliton solutions to nonlinear fractional-order differential equations. Chaos Solitons Fractals 104, 625–634 (2017)

    Article  MathSciNet  MATH  Google Scholar 

  21. Stenger, F.: Numerical methods based on the Whittaker cardinal or sinc functions. SIAM Rev. 23, 165–224 (1981)

    Article  MathSciNet  MATH  Google Scholar 

  22. Parand, K., Dehghan, M., Pirkhedria, A.: Sinc-collocation methods for solving the Blasius equation. Phys. Lett. A 373, 4060–4065 (2009)

    Article  MathSciNet  MATH  Google Scholar 

  23. Saadatmandi, A., Razzaghi, M., Dehghan, M.: Sinc-Galerkin solution for nonlinear two point boundary value problems with applications to chemical reactor theory. Math. Comput. Model. 42, 1237–1244 (2005)

    Article  MathSciNet  MATH  Google Scholar 

  24. Okayama, T.: Theorectical analysis of sinc-collocation methods and sinc-Nyström methods for systems of initial value problems. BIT Numer. Math. 58, 199–220 (2018)

    Article  MATH  Google Scholar 

  25. Rashidinia, J., Nabati, M., Barati, A.: Sinc-Galerkin method for solving nonlinear weakly singular two point boundary value problems. Int. J. Comput. Math. 94, 79–94 (2017)

    Article  MathSciNet  MATH  Google Scholar 

  26. Mueller, J.L., Shores, T.S.: A new sinc-Galerkin method for convection-diffusion equations with mixed boundary conditions. Comput. Math. Appl. 47, 803–822 (2004)

    Article  MathSciNet  MATH  Google Scholar 

  27. Rashidinia, J., Barati, A., Babati, M.: Application of sinc-Galerkin method to singularly perturbed parabolic convection-diffusion problems. Numer. Algorithms 66, 643–662 (2014)

    Article  MathSciNet  MATH  Google Scholar 

  28. Dehghan, M., Emami-Naeini, F.: The sinc-collocation and sinc-Galerkin methods for solving the two-dimensional Schrödinger equation with nonhomogeneous boundary conditions. Appl. Math. Model. 37, 9379–9397 (2013)

    Article  MathSciNet  MATH  Google Scholar 

  29. Abdrabou, A., El-Gamel, M.: On the sinc-Galerkin method for triharmonic boundary-value problems. Comput. Math. Appl. 76, 520–533 (2018)

    Article  MathSciNet  MATH  Google Scholar 

  30. Maleknejad, K., Ostadi, A.: Numerical solution of system of Volterra integral equations with weakly singular kernels and its convergence analysis. Appl. Numer. Math. 115, 82–98 (2017)

    Article  MathSciNet  MATH  Google Scholar 

  31. Okayama, T.: Theorectical analysis of a sinc-Nyström method for Volterra integro-differential equations and its improvement. Appl. Math. Comput. 324, 1–15 (2018)

    MathSciNet  MATH  Google Scholar 

  32. Maleknejad, K., Ostadi, A.: Using sinc-collocation method for solving weakly singular Fredholm integral equations of the first kind. Appl. Anal. 96, 702–713 (2016)

    Article  MathSciNet  MATH  Google Scholar 

  33. Ma, Y., Huang, J., Wang, C., Li, H.: sinc Nyström method for a class of nonlinear Volterra integral equations of the first kind. Adv. Differ. Equ. 2016, 151 (2016)

    Article  MATH  Google Scholar 

  34. Nagy, A.M.: Numerical solution of time fractional nonlinear Klein–Gordon equation using sinc-Chebyshev collocation method. Appl. Math. Comput. 310, 139–148 (2017)

    MathSciNet  MATH  Google Scholar 

  35. Saadatmandi, A., Dehghan, M., Azizi, M.: The sinc-Legendre collocation method for a class of fractional convection-diffusion equations with variable coefficients. Commun. Nonlinear Sci. Numer. Simul. 17, 4125–4136 (2012)

    Article  MathSciNet  MATH  Google Scholar 

  36. Mao, Z., Xiao, A., Yu, Z., Shi, L.: Sinc-Chebyshev collocation method for a class of fractional diffusion-wave equation. Sci. World J. 2014, 143983 (2014)

    Google Scholar 

  37. Pirkhedri, A., Javadi, H.H.S.: Solving the time-fractional diffusion equation via sinc-Haar collocation method. Appl. Math. Comput. 257, 317–326 (2015)

    MathSciNet  MATH  Google Scholar 

  38. Hesameddini, E., Asadollahifard, E.: A new reliable algorithm based on the sinc function for the time fractional diffusion equation. Numer. Algorithms 72, 893–913 (2016)

    Article  MathSciNet  MATH  Google Scholar 

  39. Jalilian, Y., Ghasemi, M.: On the solutions of a nonlinear fractional integro-differential equation of pantograph type. Mediterr. J. Math. 14, 194 (2017)

    Article  MathSciNet  MATH  Google Scholar 

  40. Stenger, F.: Numerical Method Based on Sinc and Analytic Functions. Springer, New York (1993)

    Book  MATH  Google Scholar 

  41. Lund, J., Bowers, K.: Sinc Method for Quadrature and Differential Equations. SIAM, Philadelphia (1992)

    Book  MATH  Google Scholar 

  42. Sun, Z.Z., Wu, X.N.: A fully discrete difference scheme for a diffusion-wave system. Appl. Numer. Math. 56, 193–209 (2006)

    Article  MathSciNet  MATH  Google Scholar 

Download references

Acknowledgements

The authors are very grateful to the editors and the referees for carefully reading and comments on this paper.

Availability of data and materials

Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.

Funding

The work is supported by the Natural Science Foundation of China (11271101).

Author information

Authors and Affiliations

Authors

Contributions

All authors contributed equally to the writing of this paper. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Qiang Xu.

Ethics declarations

Competing interests

The authors declare that they have no competing interests.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Chen, L.J., Li, M. & Xu, Q. Sinc-Galerkin method for solving the time fractional convection–diffusion equation with variable coefficients. Adv Differ Equ 2020, 504 (2020). https://doi.org/10.1186/s13662-020-02959-5

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13662-020-02959-5

Keywords