A Finite Difference Method on Non-Uniform Meshes for Time-Fractional Advection–Diffusion Equations with a Source Term

Featured Application

In this paper, the authors develop an unconditionally-stable, finite difference method on non-uniform grids for solving time-fractional advection–diffusion equations with a source term involving the Caputo fractional derivative. In previous years, the interest in these equations has come from their mathematical structures and from their applications. Fractional partial differential equations (FPDEs) are applied in various areas of engineering, science, finance, applied mathematics, bio-engineering and so on. In regard to its simplicity, the method is applicable to a wide class of fractional advection–diffusion equations in many fields of the applied sciences.

Abstract

The present paper deals with the numerical solution of time-fractional advection–diffusion equations involving the Caputo derivative with a source term by means of an unconditionally-stable, implicit, finite difference method on non-uniform grids. We use a special non-uniform mesh in order to improve the numerical accuracy of the classical discrete fractional formula for the Caputo derivative. The stability and the convergence of the method are discussed. The error estimates established for a non-uniform grid and a uniform one are reported, to support the theoretical results. Numerical experiments are carried out to demonstrate the effectiveness of the method.

1. Introduction

The fractional partial differential equations (FPDEs) have become increasingly popular in recent years. The interest in these equations comes from their mathematical structure and from their applications. FPDEs are applied in various areas of engineering, science, finance, applied mathematics, bio-engineering and so on. Several ways to solve them theoretically have been proposed [1,2,3], including the Green function method and Laplace and Fourier transform methods [4,5], the Adomian decomposition [6,7], and the Homotopy Perturbation methods [8,9].

Generally, numerical solution techniques are preferred when dealing with fractional models since the analytical solutions are available for a few simple cases, but also the memory effect of the problem under consideration may lead to arbitrarily wild solutions, as shown in [10]. In recent years, efficient numerical methods have been developed to solve fractional differential equations, including finite difference methods [11,12,13,14,15], the finite volume method [16], the finite element method [17,18,19], and the spectral method [20,21].

Generally, Riemann–Liouville derivatives and Caputo derivatives are used for the formulation of fractional differential problems, and often, the Riemann–Liouville formula is approximated by the Caputo one. As for the non-fractional differential equations involving the Caputo derivative, finite difference methods are one of the most important classes of numerical methods for solving FPDEs. Zhuang and Liu [22] obtained an implicit difference approximation to solve time-fractional diffusion equations. Lin and Xu [21,23] proposed the numerical solution by finite/difference approximations for a time-fractional diffusion equation. Liu et al. [24] developed an explicit difference method and an implicit difference method for solving a space–time fractional advection dispersion equation on a finite domain. In [25,26,27], high order numerical difference schemes were constructed in order to solve the Caputo-type advection–diffusion equations. Zhang et al. [28] obtained a finite difference method for FPDEs involving the Caputo derivative on a non-uniform mesh. Recently Jannelli et al. [29,30,31,32] determined exact and numerical solutions for the time and space fractional advection–diffusion differential equations involving the Riemann–Liouville derivative with a non-linear source term by means of the Lie symmetries. The authors transformed the fractional partial differential equations into a fractional ordinary differential equations, written in terms of Caputo derivatives, which were then solved using implicit finite difference methods.

The main goal of this paper was to construct an unconditionally-stable, implicit finite difference method for solving the time-fractional advection–diffusion equations (TFADEs) with a non-homogeneous source term involving the Caputo fractional derivative on non-uniform grids. We chose to use a special non-uniform mesh in order to improve the numerical accuracy of the classical discrete fractional formula for the Caputo derivative, since the fractional derivatives are integrals with weakly singular kernels, and the discretization on the uniform mesh may lead to poor accuracy. The consistency, stability, and convergence of the proposed difference method were investigated. Three numerical examples are given to show the reliability and efficiency of the derived difference method.

2. The Mathematical Model

We considered the following linear time-fractional advection–diffusion equation

$$\begin{array}{c}\hfill \frac{{\partial }^{\alpha }}{\partial t^{\alpha }}u(x,t)+K_1\frac{\partial }{\partial x}u(x,t)-K_2\frac{{\partial }^2}{\partial x^2}u(x,t)=f(x,t),a<x<b,0<t\le T,\end{array}$$

(1)

with the initial and boundary conditions given by

$$\begin{array}{c}u(x,0)=\varphi \left(x\right),a\le x\le b,\hfill \\ u(a,t)=\phi \left(t\right),u(b,t)=\psi \left(t\right),0<t\le T,\hfill \end{array}$$

where u is the field variable that can represent, for example, the solute concentration, and $K_1$ and $K_2$ are the constant fluid velocity and the dispersion coefficient, respectively. The time fractional derivative, ${\displaystyle \frac{{\partial }^{\alpha }}{\partial t^{\alpha }}u(x,t)}$ is the $\alpha$ order Caputo fractional derivative defined by

$$\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}u(x,t)=\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^t\frac{\partial }{\partial s}u(x,s){(t-s)}^{-\alpha }ds0<\alpha <1.$$

(2)

The function $f(x,t)$ can be used to represent sources and sinks. $\varphi \left(x\right)$, $\phi \left(t\right)$ and $\psi \left(t\right)$ are known smooth functions. We took $K_1\ne 0$ and $K_2>0$ and we assumed problem (1) had a unique and sufficiently smooth solution under the above initial and boundary conditions.

The fractional Equation (1) represents a suitable mathematical model for describing anomalous physical processes that exhibit fractional order behavior that varies with time or space and that cannot be modeled accurately by normal integer order equations. In fact, it has been treated by a number of authors as a useful approach for the description of transport dynamics in complex systems which are governed by anomalous diffusion and non-exponential relaxation patterns [33]. The TFADE is also used in groundwater hydrology research to model the transport of passive tracers carried by fluid flow in porous mediums [34] and in neurology [35]. It is notable that, when $\alpha =1$, the model (1) reduces to the classical advection–diffusion–reaction equation

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}u(x,t)+K_1\frac{\partial }{\partial x}u(x,t)-K_2\frac{{\partial }^2}{\partial x^2}u(x,t)=f(x,t),a\le x\le b,0<t\le T,\end{array}$$

used in order to describe several phenomena of relevant interest in many fields of applied sciences. It is a mathematical model that describes how the concentration of the substance distributed in a medium changes under the influence of three processes: advection, diffusion and reaction. In [36,37,38], a fractional step approach with variable time step was used in order to solve numerically mathematical models that describe evolution problems on computational domains of one or three space dimensions.

3. Discretization in Time on a Non-Uniform Mesh

The main goal of this work was to construct an unconditionally-stable implicit finite difference method defined on a non-uniform mesh. In general, the existence of a weakly singular kernel ${(t-s)}^{-\alpha }$, $0<\alpha <1$, in fractional derivatives makes it more difficult to obtain a higher-order scheme. In particular, when the solutions are not suitably smooth, numerical methods on uniform meshes seem to have a poor convergent rate. For these reasons, numerical schemes on non-uniform meshes have been developed in the last years.

In this section, first, we construct the non-uniform mesh and then, in order to approximate the Caputo derivative on the non-uniform grid, we define a suitable discrete fractional derivative formula.

We divide the interval, $[0,T]$, into N subintervals, $[t^{n-1},t^n]$, for $n=1,\cdots ,N$ and with $0=t^0<t^1<\cdots t^N=T$. We denote the time step $\Delta t^n$ as

$$\begin{array}{c}\hfill \Delta t^n=t^n-t^{n-1},1\le n\le N,\end{array}$$

and let

$$\begin{array}{c}\hfill \Delta t^{max}=\underset{1\le n\le N}{max}\Delta t^n\Delta t^{min}=\underset{1\le n\le N}{min}\Delta t^n.\end{array}$$

A sequence of mesh is described as quasi-uniform if there a finite constant, C, exists such that

$$\begin{array}{c}\hfill \Delta t^{max}/\Delta t^{min}\le C.\end{array}$$

In this case, it holds that $\Delta t^{max}\le \beta T{N}^{-1}$. When $C=1$, it holds that $\Delta t^n=T{N}^{-1}$ for all $n=1,\cdots ,N$, and the mesh is reduced to a uniform mesh.

A sequence of meshes is not quasi-uniform if

$$\begin{array}{c}\hfill \Delta t^{max}/\Delta t^{min}\to +\infty \mathrm{as}N\to +\infty .\end{array}$$

The non-uniform mesh and quasi-uniform mesh methods have been used for solving differential equations by several authors. In [39,40,41,42], quasi-uniform meshes were used for solving a class of boundary values problems on infinite domains.

In this work, we are interested on the non-uniform mesh [28], which is defined as follows

$$\begin{array}{c}\hfill \Delta t^n=(N+1-n)\sigma ,1\le n\le N,\end{array}$$

(3)

where $\sigma ={\displaystyle \frac{2T}{N(N+1)}}$. Note that the time steps, ${\{\Delta t^n\}}_{n=1}^N$, are a monotonically decreasing sequence, with $\Delta t^1=O\left(N^{-1}\right)$ and $\Delta t^N=O\left(N^{-2}\right)$. Figure 1 shows a sample of the non-uniform grid (3) obtained for $N=10$. The fractional derivatives are integrals with weakly singular kernels, and it is well known that the discretization on the uniform mesh may lead to poor accuracy. We decided to use a non-uniform mesh in order to improve the accuracy of the numerical solution. Another interesting choice could be to construct non-uniform meshes that refine the computational domain according to the fractional order of $\alpha$ or to use adaptive procedures to dynamically choose the size of the time-steps according to the local behavior of the solution.

To motivate the choice of the non-uniform mesh (3), we discretized the Caputo derivative (2) of a function $v\left(t\right)$ by means of the following classical approximate formula

$$\begin{array}{ccc}\hfill \frac{d^{\alpha }}{d{t}^{\alpha }}v\left(t^n\right)& =& \frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^{t^n}{v}^{\prime }\left(s\right){(t^n-s)}^{-\alpha }ds\hfill \\ & =& \frac{1}{\mathsf{\Gamma }(1-\alpha )}\sum _{k=1}^n\frac{v\left(t^k\right)-v\left(t^{k-1}\right)}{t^k-t^{k-1}}{\int }_{t^{k-1}}^{t^k}{(t^n-s)}^{-\alpha }ds+R^n,\hfill \end{array}$$

(4)

which is the well-known, so-called L1 formula defined in [2] where $R^n$ is the local truncation error. The L1 formula has been used to solve the fractional differential equations with Caputo derivatives (see [22,43]). Moreover, using the relationship between the Caputo derivative and the Riemann–Liouville fractional derivative, the L1 formula was also applied to the time fractional diffusion equation with the Riemann–Liouville fractional derivative (see [44,45]). High-order approximations, such as the compact difference scheme [43,45,46,47] and spectral method [21,23,48] were applied to improve the spatial accuracy of fractional diffusion equations. It is important to note that it is rather difficult to get a high-order time approximation due to the singularity of fractional derivatives.

For any temporal mesh, for $0<\alpha <1$ and $v\left(t\right)\in C^2[0,T]$, it can be verified that [28]

$$\begin{array}{c}\hfill {\int }_0^{t^n}{v}^{\prime }\left(s\right){(t^n-s)}^{-\alpha }ds=\sum _{k=1}^n\frac{v\left(t^k\right)-v\left(t^{k-1}\right)}{t^k-t^{k-1}}{\int }_{t^{k-1}}^{t^k}{(t^n-s)}^{-\alpha }ds+r^n,\end{array}$$

where

$$\begin{array}{c}\hfill |r^n|\le \left(\frac{{(\Delta t^n)}^2}{2(1-\alpha )}+\frac{{(\Delta t_{max})}^2}{8}\right){(\Delta t^n)}^{-\alpha }\underset{0\le t\le t^n}{max}|v^{″}\left(t\right)|1\le n\le N.\end{array}$$

Since

$$\begin{array}{ccc}\hfill R^n& =& \frac{1}{\mathsf{\Gamma }(1-\alpha )}\left[{\int }_0^{t^n}{v}^{\prime }\left(s\right){(t^n-s)}^{-\alpha }ds-\sum _{k=1}^n\frac{v\left(t^k\right)-v\left(t^{k-1}\right)}{t^k-t^{k-1}}{\int }_{t^{k-1}}^{t^k}{(t^n-s)}^{-\alpha }ds\right]\hfill \\ & =& \frac{1}{\mathsf{\Gamma }(1-\alpha )}{r}^n,\hfill \end{array}$$

we obtain

$$\begin{array}{c}\hfill |R^n|\le \frac{1}{\mathsf{\Gamma }(1-\alpha )}\left(\frac{{(\Delta t^n)}^2}{2(1-\alpha )}+\frac{{(\Delta t_{max})}^2}{8}\right){(\Delta t^n)}^{-\alpha }\underset{0\le t\le t^n}{max}|v^{″}\left(t\right)|1\le n\le N.\end{array}$$

For the uniform mesh, i.e., $\Delta t^n=\Delta t$ for all $n=1,2,\cdots ,N$, then $R^n=O(\Delta t^{2-\alpha })$ (see [21,49]). From the truncation error estimate of the L1 formula, it is clear that the accuracy is dependent on the fractional order, $\alpha$. This is justified since a weakly singular kernel, ${(t-s)}^{-\alpha }$, is contained in the integral.

In order to improve the accuracy of the L1 numerical approximation of the derivative of fractional order, it is possible to consider non-uniform meshes. Numerical methods developed with the non-uniform meshes have been developed to solve weakly singular integro-differential equations. Mustapha [50], Yuste and Quintana-Murillo [51,52] proposed an implicit finite-difference time-stepping method for the discretization of the time diffusion equation.

Zhang et al. [28] obtained a numerical integration formula for any $\alpha \in (0,1)$ by employing the special non-uniform grid (3) and obtained the following results for sufficiently smooth solutions: for the non-uniform mesh (3), for $0<\alpha <1$ and $v\left(t\right)\in C^2[0,T]$, it held that

$$\begin{array}{c}\hfill {\int }_0^{t^n}{v}^{\prime }\left(s\right){(t^n-s)}^{-\alpha }ds=\sum _{k=1}^n\frac{v\left(t^k\right)-v\left(t^{k-1}\right)}{t^k-t^{k-1}}{\int }_{t^{k-1}}^{t^k}{(t^n-s)}^{-\alpha }ds+r^{n}1\le n\le N,\end{array}$$

where

$$\begin{array}{c}|r^n|\le \left(1+\alpha +\frac{2^{1-\alpha }}{1-\alpha }\right)\underset{0\le t\le t^n}{max}|v^{″}\left(t\right)|T^{2-\alpha }{(N+1)}^{\alpha -2},1\le n\le N-1,\hfill \end{array}$$

and

$$\begin{array}{c}\hfill {\int }_0^{t^N}{v}^{\prime }\left(s\right){(t^N-s)}^{-\alpha }ds=\sum _{k=1}^N\frac{v\left(t^k\right)-v\left(t^{k-1}\right)}{t^k-t^{k-1}}{\int }_{t^{k-1}}^{t^k}{(t^N-s)}^{-\alpha }ds+r^N,\end{array}$$

where

$$\begin{array}{c}|r^N|\le \frac{1+\alpha }{1-\alpha }{2}^{1-\alpha }\underset{0\le t\le T}{max}|v^{″}\left(t\right)|T^{2-\alpha }{N}^{-2}.\hfill \end{array}$$

Then, we obtain

$$\begin{array}{c}|R^n|\le \frac{1}{\mathsf{\Gamma }(1-\alpha )}\left(1+\alpha +\frac{2^{1-\alpha }}{1-\alpha }\right)\underset{0\le t\le t^n}{max}|v^{″}\left(t\right)|T^{2-\alpha }{(N+1)}^{\alpha -2},1\le n\le N-1,\hfill \\ |R^N|\le \frac{1}{\mathsf{\Gamma }(1-\alpha )}\frac{1+\alpha }{1-\alpha }{2}^{1-\alpha }\underset{0\le t\le T}{max}|v^{″}\left(t\right)|T^{2-\alpha }{N}^{-2}.\hfill \end{array}$$

(5)

We use these estimates to study the consistency, stability and convergence of the method presented in the following sections.

4. An Implicit Finite Difference Method

For the derivation of the implicit difference method, first, we constructed a computational uniform grid in the x direction, that is, the spatial size of the mesh, $\Delta x=x_j-x_{j-1}$, is constant, for $1\le j\le J$, and it is non-uniform in the time direction. We defined the mesh points $(x_j,t^n)$ with $x_j=a+j\Delta x$, $j=0,\cdots ,J$ and $t^n=t^{n-1}+\Delta t^n$, for $n=1,\cdots ,N$, with $\Delta t^n$ defined by (3). J and N are positive integers. We denoted the numerical approximation provided by the difference method of the exact solution $u(x_j,t^n)$ by $U_j^n$ at the mesh points $(x_j,t^n)$, for $j=0,\cdots ,J$ and $n=0,\cdots ,N$.

As usual, we supposed that the solution was sufficiently smooth and discretized its first $\partial /\partial x$ and second order ${\partial }^2/\partial x^2$ spatial derivatives by means of the second order, three-point central difference formula so that

$$\begin{array}{c}\frac{\partial }{\partial x}u(x_j,t^n)=\frac{u(x_{j+1},t^n)-u(x_{j-1},t^n)}{2\Delta x}+O(\Delta x^2)\hfill \\ \hfill \end{array}$$

(6)

$$\begin{array}{c}\frac{{\partial }^2}{\partial x^2}u(x_j,t^n)=\frac{u(x_{j+1},t^n)-2u(x_j,t^n)+u(x_{j-1},t^n)}{\Delta x^2}+O(\Delta x^2).\hfill \end{array}$$

(7)

According to the discretization for the Caputo derivative (4), we approximated the time fractional derivative in (1) as follows

$$\begin{array}{c}\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}u(x_j,t^n)=\hfill \\ \frac{1}{\mathsf{\Gamma }(2-\alpha )}\sum _{k=1}^n\frac{u(x_j,t^k)-u(x_j,t^{k-1})}{t^k-t^{k-1}}\left[{(t^n-t^{k-1})}^{1-\alpha }-{(t^n-t^k)}^{1-\alpha }\right]+R_j^n,\hfill \end{array}$$

(8)

where

$$\begin{array}{c}\hfill \frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_{t^{k-1}}^{t^k}{(t^n-s)}^{-\alpha }ds=\frac{1}{\mathsf{\Gamma }(2-\alpha )}[{(t^n-t^{k-1})}^{1-\alpha }-{(t^n-t^k)}^{1-\alpha }].\end{array}$$

Replacing $u(x_j,t^n)$ with its numerical approximation, $U_j^n$, and neglecting the local truncation errors, the time-fractional advection–diffusion Equation (1) was discretized as follows

$$\begin{array}{c}\hfill \frac{1}{\mathsf{\Gamma }(2-\alpha )}\sum _{k=1}^n{T}_{n,k}(U_j^k-U_j^{k-1})+K_1\frac{U_{j+1}^n-U_{j-1}^n}{2\Delta x}-K_2\frac{U_{j+1}^n-2U_j^n+U_{j-1}^n}{\Delta x^2}=f_j^n,\end{array}$$

(9)

for $1\le n\le N$ and $1\le j\le J-1$, where $f_j^n=f(x_j,t^n)$ is the source term and

$$\begin{array}{c}\hfill T_{n,k}=\frac{{(t^n-t^{k-1})}^{1-\alpha }-{(t^n-t^k)}^{1-\alpha }}{t^k-t^{k-1}}1\le k\le n,1\le n\le N.\end{array}$$

The initial and boundary conditions can be rewritten as

$$\begin{array}{c}{U}_j^0=\varphi \left(x_j\right)0\le j\le J\hfill \\ U_0^n=\phi \left(t^n\right)U_J^n=\psi \left(t^n\right)0\le n\le N.\hfill \end{array}$$

(10)

For any temporal meshes on $[0,T]$ and for any $n=1,2,\cdots ,N$, it holds that [22]

$$\begin{array}{c}{T}_{n,k}>0T_{n,k}>T_{n,k-1},1\le k\le n.\hfill \end{array}$$

Taking into account that $T_{n,n}={(\Delta t^n)}^{-\alpha }$, we can write

$$\begin{array}{c}\sum _{k=1}^n{T}_{n,k}(U_j^k-U_j^{k-1})=\sum _{k=1}^{n-1}{T}_{n,k}(U_j^k-U_j^{k-1})+{(\Delta t^n)}^{-\alpha }(U_j^n-U_j^{n-1}).\hfill \end{array}$$

Then, the following implicit finite difference scheme was obtained:

$$\begin{array}{c}(-K1-K2)U_{j-1}^n+(1+2K2)U_j^n+(K1-K2)U_{j+1}^n\hfill \\ =U_j^{n-1}-{(\Delta t^n)}^{\alpha }\sum _{k=1}^{n-1}{T}_{n,k}(U_j^k-U_j^{k-1})+F_j^n,1\le n\le N,1\le j\le J-1,\hfill \end{array}$$

(11)

where we set

$$\begin{array}{c}\hfill K1=\frac{K_1{(\Delta t^n)}^{\alpha }\mathsf{\Gamma }(2-\alpha )}{2\Delta x},K2=\frac{K_2{(\Delta t^n)}^{\alpha }\mathsf{\Gamma }(2-\alpha )}{\Delta x^2},F_j^n={(\Delta t^n)}^{\alpha }\mathsf{\Gamma }(2-\alpha )f_j^n.\end{array}$$

The Equation (11) can be written in vectorial form as

$$\begin{array}{c}K{U}^n=\mathcal{L}{U}^{n-1}+F^n,\hfill \end{array}$$

(12)

where K is a tridiagonal matrix and where $\mathcal{L}$ denotes the following difference operator

$$\begin{array}{c}\mathcal{L}{U}^{n-1}=U_j^{n-1}-{(\Delta t^n)}^{\alpha }\sum _{k=1}^{n-1}{T}_{n,k}(U_j^k-U_j^{k-1}).\hfill \end{array}$$

(13)

Here, and in the following text, we assume the convention that the summation is equal to zero if the lower bound is larger that the upper bound. The obtained method (12) is implicit. In order to compute the numerical solution, $U^n$, a system with the tridiagonal coefficients matrix, K, has to be solved.

It is interesting to note that the operator, $\mathcal{L}$, is a kind operator with memory, due to the non-local character of the fractional derivative. This means that the effect on U at time $t^n$, $U^n$ depends on all the previous values, $U^0,U^1,\cdots ,U^{n-1}$, evaluated at all the previous time, $t^0,t^1,\cdots ,t^{n-1}$.

The main difference compared to the non-fractional case is that, in order to evaluate $\mathcal{L}$, the numerical solutions for all the n previous time values $t^0,t^1,\cdots ,t^{n-1}$ are required, while for non-fractional equations, only the solution to the previous value $t^{n-1}$ is used. The computational cost to compute the solution at the time $t^n$ from the solution at the time $t^{n-1}$ grows as n. That is, it grows as the number of terms in the summation that compares in the second term of the (13). This implies that the computational cost goes from $t^0$ to $t^n$, thus growing by $n^2$.

5. Consistency, Stability, and Convergence

In this section, we discuss the consistency, the stability, and the convergence of the implicit finite difference scheme.

Consistency. According to Equations (6)–(8), the local truncation error of the difference scheme (9) is

$$\begin{array}{ccc}\hfill R_j^n& =& \frac{1}{\mathsf{\Gamma }(2-\alpha )}\sum _{k=1}^n{T}_{n,k}(u(x_j,t^k)-u(x_j,t^{k-1}))\hfill \\ & & +K_1\frac{u(x_{j+1},t^n)-u(x_{j-1},t^n)}{2\Delta x}-K_2\frac{u(x_{j+1},t^n)-2u(x_j,t^n)+u(x_{j-1},t^n)}{\Delta x^2}-f(x_j,t^n)\hfill \\ & =& \frac{1}{\mathsf{\Gamma }(2-\alpha )}\sum _{k=1}^n{T}_{n,k}(u(x_j,t^k)-u(x_j,t^{k-1}))-\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}u(x_j,t^n)\hfill \\ & & +K_1\left[\frac{u(x_{j+1},t^n)-u(x_{j-1},t^n)}{2\Delta x}-\frac{\partial }{\partial x}u(x_j,t^n)\right]\hfill \\ & & -K_2\left[\frac{u(x_{j+1},t^n)-2u(x_j,t^n)+u(x_{j-1},t^n)}{\Delta x^2}-\frac{{\partial }^2}{\partial x^2}u(x_j,t^n)\right]\hfill \\ & =& O\left(N^{\alpha -2}\right)+K_{1}O(\Delta x^2)+K_{2}O(\Delta x^2)=O(N^{\alpha -2}+\Delta x^2).\hfill \end{array}$$

(14)

The implicit finite difference scheme defined by (9) or (11) is consistent with model (1) of order $O(N^{\alpha -2}+\Delta x^2)$.

Stability. For the stability analysis of the implicit finite difference scheme Equation (11) was rewritten as

$$\begin{array}{c}(-K1-K2)U_{j-1}^n+(1+2K2)U_j^n+(K1-K2)U_{j+1}^n\hfill \\ ={(\Delta t^n)}^{\alpha }\sum _{k=0}^{n-1}(T_{n,k+1}-T_{n,k})U_j^k+F_j^n,1\le n\le N,1\le j\le J-1.\hfill \end{array}$$

(15)

Let ${\overline{U}}_j^n$ be another approximate solution of the difference scheme (11), and let

$$\begin{array}{c}\hfill {\rho }_j^n=U_j^n-{\overline{U}}_j^n,0\le j\le J,1\le n\le N,\end{array}$$

be the corresponding round-off error. We let

$$\begin{array}{c}\hfill {\rho }^n={({\rho }_0^n,{\rho }_1^n,\cdots ,{\rho }_J^n)}^T,\end{array}$$

and we considered the infinity norm

$$\begin{array}{c}\hfill ||{\rho }^n{||}_{\infty }=\underset{0\le j\le J}{max}|{\rho }_j^n|=|{\rho }_i^n|.\end{array}$$

The round-off error satisfied the following round-off equations

$$\begin{array}{c}(-K1-K2){\rho }_{j-1}^n+(1+2K2){\rho }_j^n+(K1-K2){\rho }_{j+1}^n\hfill \\ ={(\Delta t^n)}^{\alpha }\sum _{k=0}^{n-1}(T_{n,k+1}-T_{n,k}){\rho }_j^k.\hfill \end{array}$$

(16)

In order to check whether the finite difference scheme was stable, we studied how the size of the round-off error, ${\rho }^n$, evolved over time.

We defined

$$\begin{array}{c}{L}_1{\rho }_j^n=(-K1-K2){\rho }_{j-1}^n+(1+2K2){\rho }_j^n+(K1-K2){\rho }_{j+1}^n,\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill & L_2{\rho }_j^{n-1}={(\Delta t^n)}^{\alpha }\sum _{k=0}^{n-1}(T_{n,k+1}-T_{n,k}){\rho }_j^k.\hfill \end{array}$$

(17)

Equation (16) can be written as

$$\begin{array}{c}{L}_1{\rho }_j^n=L_2{\rho }_j^{n-1}.\hfill \end{array}$$

From (17), and taking into account that $T_{n,k+1}-T_{n,k}>0$, we obtained

$$\begin{array}{c}\hfill |L_2{\rho }_j^{n-1}|=|{(\Delta t^n)}^{\alpha }\sum _{k=0}^{n-1}(T_{n,k+1}-T_{n,k}){\rho }_j^k|\le |{\rho }_j^{n-1}|{(\Delta t^n)}^{\alpha }\sum _{k=0}^{n-1}(T_{n,k+1}-T_{n,k}),\end{array}$$

(18)

where we defined

$$\begin{array}{c}\hfill |{\rho }_j^{n-1}|=\underset{0\le k\le n-1}{max}|{\rho }_j^k|.\end{array}$$

But, since ${\sum }_{k=0}^{n-1}(T_{n,k+1}-T_{n,k})={(\Delta t^n)}^{-\alpha }$ and recalling that $T_{n,n}={(\Delta t^n)}^{-\alpha }$ and $T_{n,0}=0$,

$$\begin{array}{c}\hfill |L_2{\rho }_j^{n-1}|=|{(\Delta t^n)}^{\alpha }\sum _{k=0}^{n-1}(T_{n,k+1}-T_{n,k}){\rho }_j^k|\le |{\rho }_j^{n-1}|,\end{array}$$

Thus, we concluded that

$$\begin{array}{ccc}\hfill ||{\rho }^n{||}_{\infty }& =& |{\rho }_i^n|=|(-K1-K2){\rho }_i^n+(1+2K2){\rho }_i^n+(K1-K2){\rho }_i^n|\hfill \\ & =& |L_1{\rho }_i^n|=|L_2{\rho }_i^{n-1}|\le |{\rho }_i^{n-1}|=||{\rho }^{n-1}{||}_{\infty },\hfill \end{array}$$

i.e.,

$$||{\rho }^n{||}_{\infty }\le ||{\rho }^0{||}_{\infty },$$

for $1\le n\le N$. This means that the present method is unconditionally-stable.

Convergence. We let $u(x_j,t_n)$ be the exact solution of Equation (1) at mesh point $(x_j,t^n)$ for $j=0,1,\cdots ,J$ and $n=0,1,\cdots ,N$. Denoting ${ϵ}_j^n=u(x_j,t^n)-U_j^n$, we obtained the error equations

$$\begin{array}{c}{L}_1{ϵ}_j^n=L_2{ϵ}_j^{n-1}+R_j^n,j=0,1,\cdots ,J,n=0,1,\cdots ,N,\hfill \end{array}$$

(19)

with ${ϵ}_j^0=0$, for $j=0,1,\cdots ,J$ and ${ϵ}_J^n=0$, for $n=1,2,\cdots ,N$. $R_j^n$ is the local truncation error. We introduced the following norm

$$\begin{array}{c}\hfill ||{ϵ}^n{||}_{\infty }=\underset{0\le j\le J}{max}|{ϵ}_j^n|=|{ϵ}_i^n|,\end{array}$$

then, we obtained

$$\begin{array}{ccc}\hfill ||{ϵ}^n{||}_{\infty }& =& |{ϵ}_i^n|=|L_1{ϵ}_i^n|=|L_2{ϵ}_i^{n-1}+R_i^n|\hfill \\ & \le & |L_2{ϵ}_i^{n-1}|+|R_i^n|\le |{ϵ}_i^{n-1}|+|R_i^n|\le ||{ϵ}^{n-1}{||}_{\infty }+R_{max},\hfill \end{array}$$

where $R_{max}={max}_{n,i}|R_i^n|$. For $n=1$, we obtained

$$\begin{array}{c}\hfill ||{ϵ}^1{||}_{\infty }\le ||{ϵ}^0{||}_{\infty }+R_{max}=R_{max}.\end{array}$$

Then,

$$\begin{array}{c}\hfill ||{ϵ}^n{||}_{\infty }\le R_{max},1\le n\le N.\end{array}$$

To obtain the error estimates of the numerical solutions, we needed the uniform error bounds on all time levels. Then, applying the first of the estimates (5), we were able to consider

$$\begin{array}{c}{R}_{max}\le C_R(N^{\alpha -2}+\Delta x^2),1\le n\le N,\hfill \end{array}$$

(20)

where $C_R$ is a positive constant that is dependent on T, $\alpha$ and the exact solution $u(x,t)$, but is independent of N and $\Delta x$. Thus, we proved that the solution of the difference method (11), with initial and boundary conditions given by (10), is convergent.

6. Numerical Experiments

In this section, we report some numerical examples of the TFPDEs to demonstrate the accuracy and efficiency of the numerical method. At first, we wanted to show that the approach proposed in this paper properly works, so we chose two examples in such a way that the exact solutions of FPDE could be evaluated analytically. This allowed us to check the accuracy and the order of convergence of the numerical solution. In both the examples, we assumed suitable smooth solutions, and we used the proposed method on the non-uniform grid and on a uniform grid. We compared the numerical results, observing that the finite difference scheme generated more accurate numerical solutions on the non-uniform grid than on the uniform one. In the third test, we solved a TFADE of physical interest with the source term chosen as a linear function of the field variable in order to show that the method is applicable to a wide class of TFADEs. The presented problem is one of the most used mathematical models in the applied sciences. It describes how the concentration of one or more substances, chemical or biological species, distributed in a medium changes under the influence of three processes, namely, advection, diffusion and reaction. With this last test, we illustrate how the changes in the solution behavior arise when the fractional order is varied.

Example 1.

We consider the following TFADE

$$\begin{array}{c}\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}u(x,t)+\frac{\partial }{\partial x}u(x,t)-\frac{{\partial }^2}{\partial x^2}u(x,t)=f(x,t),0<x<1,0<t\le T,\hfill \\ u(x,0)=0,0\le x\le 1,\hfill \\ u(0,t)=t^{\beta },u(1,t)=e{t}^{\beta },0<t\le T,\hfill \end{array}$$

(21)

where the source term is given by

$$\begin{array}{c}f(x,t)=\frac{\mathsf{\Gamma }(\beta +1)}{\mathsf{\Gamma }(\beta +1-\alpha )}{e}^x{t}^{\beta -\alpha }.\hfill \end{array}$$

The exact solution is

$$\begin{array}{c}u(x,t)=e^x{t}^{\beta }.\hfill \end{array}$$

In this example, we took $K_1=K_2=1$ with $\alpha =0.5$ and $\beta =5$. Figure 2 shows the comparison between the numerical solution, $U_j^n$, and the exact solution, $u(x_j,t^n)$, at different times computed with $N=J=20$. From Figure 2, it can be seen that the numerical solution, $U_j^n$, is in good agreement with the exact solution $u(x_j,t^n)$. The exact solution is reported with the solid line. We report only the numerical results obtained for the value parameter $\alpha =0.5$. Analogous results were obtained for $0.1\le \alpha \le 0.9$.

In order to investigate the temporal error and the convergence order of the numerical difference method, we defined the maximum error between the exact solution, $u(x_j,t^N)$, and the numerical solution, $U_j^N$, at the final time, $t^N$

$$\begin{array}{ccc}\hfill e_{\infty }(N,J)& =& \underset{1\le j\le J}{max}|u(x_j,t^N)-U_j^N|,\hfill \end{array}$$

(22)

and the convergence order, as follows

$$\begin{array}{c}\hfill Order={log}_2\left(\frac{e_{\infty }(N,J)}{e_{\infty }(2N,J)}\right).\end{array}$$

(23)

In order to show the efficiency of the method, we defined ${\overline{U}}_j^n$ as the numerical solution obtained by the implicit difference method on a uniform grid defined by $\Delta t^n=T/N$, for $n=1,\cdots ,N$. We used notations similar to the (22) and (23) to define the maximum error, ${\overline{e}}_{\infty }(N,J)$, between the exact solution, $u(x_j,t^N)$, and the numerical solution, ${\overline{U}}_j^N$, at the final time, $t^N$, and the convergence $Orde{r}_{\overline{U}}$. In this test, we fixed $J=100$, a value large enough such that the spatial error is negligible as compared with the temporal error.

Table 1. $e_{\infty }$, ${\overline{e}}_{\infty }$ and convergence orders, related to the numerical solutions $U_j^n$ and ${\overline{U}}_j^n$ respectively, for different values of N and $\alpha$, with $J=100$.

$\mathit{\alpha }$	N	${\mathit{e}}_{\mathbf{\infty }}$	$\mathit{Order}$	${\overline{\mathit{e}}}_{\mathbf{\infty }}$	${\mathit{Order}}_{\overline{\mathit{U}}}$
$0.1$	10	$3.6363\times {10}^{-4}$		$1.3741\times {10}^{-3}$
	20	$9.1021\times {10}^{-5}$	$1.9982$	$4.3575\times {10}^{-4}$	$1.6570$
	40	$2.2054\times {10}^{-5}$	$2.0452$	$1.3222\times {10}^{-4}$	$1.7206$
	80	$4.4649\times {10}^{-6}$	$2.3043$	$3.8355\times {10}^{-5}$	$1.7854$
$0.5$	10	$5.5793\times {10}^{-3}$		$1.9875\times {10}^{-2}$
	20	$1.7121\times {10}^{-3}$	$1.7044$	$7.7106\times {10}^{-3}$	$1.3660$
	40	$5.4236\times {10}^{-4}$	$1.6584$	$2.8895\times {10}^{-3}$	$1.4160$
	80	$1.7544\times {10}^{-4}$	$1.6283$	$1.0602\times {10}^{-3}$	$1.4465$
$0.9$	10	$4.6556\times {10}^{-2}$		$1.0115\times {10}^{-1}$
	20	$2.1358\times {10}^{-2}$	$1.1242$	$4.9512\times {10}^{-2}$	$1.0307$
	40	$9.8735\times {10}^{-3}$	$1.1131$	$2.3691\times {10}^{-2}$	$1.0634$
	80	$4.5790\times {10}^{-3}$	$1.1085$	$1.1201\times {10}^{-2}$	$1.0807$

shows the values of $e_{\infty }$ and ${\overline{e}}_{\infty }$ and the corresponding numerical convergence orders for $\alpha =0.1,0.5$ and $0.9$. It can be seen that the method is stable and convergent when solving the problem (21) on both the computational grids. By using the L1 formula on non-uniform mesh, we improvde the order of accuracy in time of the proposed method. In fact, we observe that the solutions are more accurate on the non-uniform mesh and the convergence order on the non-uniform mesh is greater than the convergence order obtained with the uniform mesh. The numerical results agree well with the theoretical results.

Example 2.

We consider the following TFADE

$$\begin{array}{c}\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}u(x,t)+\frac{\partial }{\partial x}u(x,t)-\frac{{\partial }^2}{\partial x^2}u(x,t)=f(x,t),0<x<1,0<t\le T,\hfill \\ u(x,0)=0,0\le x\le 1,\hfill \\ u(0,t)=0,u(1,t)=t^3,0<t\le T,\hfill \end{array}$$

(24)

where the source term is given by

$$\begin{array}{c}f(x,t)=\frac{6}{\mathsf{\Gamma }(4-\alpha )}{x}^2{t}^{3-\alpha }+2t^3(x-1)\hfill \end{array}$$

and the exact solution is

$$\begin{array}{c}u(x,t)=x^2{t}^3.\hfill \end{array}$$

We took $K_1=K_2=1$ and set $\alpha =0.1$. Figure 3 shows the comparison between the numerical solution, $U_j^n$, and the exact solution, $u(x_j,t^n)$, at different times computed with $N=J=20$. The exact solution is reported with the solid line. The numerical solution $U_j^n$ is in good agreement with the exact solution $u(x_j,t^n)$. We report only the numerical results obtained for the value of the parameter, $\alpha =0.1$. Analogous results were obtained for $0.1<\alpha \le 0.9$.

The values of $e_{\infty }$ and ${\overline{e}}_{\infty }$ and the corresponding numerical convergence orders obtained for $\alpha =0.1,0.5$ and $0.9$ are reported in

Table 2. $e_{\infty }$, ${\overline{e}}_{\infty }$ and convergence orders, related to the numerical solutions, $U_j^n$ and ${\overline{U}}_j^n$, respectively, for different values of N and $\alpha$, with $J=100$.

$\mathit{\alpha }$	N	${\mathit{e}}_{\mathbf{\infty }}$	$\mathit{Order}$	${\overline{\mathit{e}}}_{\mathbf{\infty }}$	${\mathit{Order}}_{\overline{\mathit{U}}}$
$0.1$	10	$3.7836\times {10}^{-5}$		$9.4723\times {10}^{-5}$
	20	$9.6542\times {10}^{-6}$	$1.9705$	$2.8745\times {10}^{-5}$	$1.7204$
	40	$2.4674\times {10}^{-6}$	$1.9682$	$8.5307\times {10}^{-6}$	$1.7526$
	80	$6.2843\times {10}^{-7}$	$1.9732$	$2.4911\times {10}^{-6}$	$1.7759$
$0.5$	10	$4.4597\times {10}^{-4}$		$1.3182\times {10}^{-3}$
	20	$1.3436\times {10}^{-4}$	$1.7309$	$4.8998\times {10}^{-4}$	$1.4278$
	40	$4.1920\times {10}^{-5}$	$1.6803$	$1.7898\times {10}^{-4}$	$1.4529$
	80	$1.3435\times {10}^{-5}$	$1.6416$	$6.4671\times {10}^{-5}$	$1.4686$
$0.9$	10	$3.3229\times {10}^{-3}$		$6.9875\times {10}^{-3}$
	20	$1.4844\times {10}^{-3}$	$1.1625$	$3.3290\times {10}^{-3}$	$1.0682$
	40	$6.7617\times {10}^{-4}$	$1.1345$	$1.5524\times {10}^{-3}$	$1.0837$
	80	$3.1119\times {10}^{-4}$	$1.1196$	$7.2847\times {10}^{-4}$	$1.0915$

. The results confirm that the method is stable and convergent when solving problems (24) on both the computational grids. The numerical solutions computed on non-uniform mesh were more accurate and the convergence order on the non-uniform mesh was greater than the convergence order obtained with the uniform mesh. The errors $e_{\infty }$ and ${\overline{e}}_{\infty }$ satisfied the relationships (22). The convergence orders, $Order$ and $Orde{r}_{\overline{U}}$, satisfied the relationships (23).

Example 3.

In this test, in order to show the efficiency of the method, we solve the following TFADE

$$\begin{array}{c}\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}u(x,t)+K_1\frac{\partial }{\partial x}u(x,t)-K_2\frac{{\partial }^2}{\partial x^2}u(x,t)=f\left(u(x,t)\right),\hfill \\ u(x,0)=x^2{(5-x)}^2,a\le x\le b,\hfill \\ u(a,t)=u(b,t)=0,0<t\le T,\hfill \end{array}$$

(25)

where the source term is chosen as a linear function of the field variable

$$\begin{array}{c}f\left(u\right)=\beta u(x,t).\hfill \end{array}$$

The model describes the one-dimensional transport problem of a concentration, $u(x,t)$, of a chemical or biological species in a flowing medium, such as air or water. The species concentration is assumed to be horizontally and vertically well mixed, such that it varies only in the longitudinal or downstream direction. Moreover, a steady and uniform flow field is imposed and the effects of the dispersion are constant in time and space. A reaction where the transformation rate $\beta$ is proportional to the species concentration was considered; according to the sign of rate, decay effects may or may not have occurred. This model has been used by several authors, because it is able to model a wide variety of physical and biological phenomena, see, for example, [53,54,55].

In this example, we used $a=0$, $b=5$, $K_1=1$, $K_2=1$ and $\beta =0.2$. Figure 4 shows the solution behavior at different times obtained for $\alpha =0.1$ and $\alpha =0.9$ and with $N=J=100$. The height of the solution profile decreased as the time increased. Analogous results were obtained for other values of $\alpha$. Lower profiles were obtained for $\alpha =0.1$, and higher ones were obtained for $\alpha =0.9$ at the final time, $t^N=1$, that the behavior of solutions were comparable. Figure 5 shows the solution behavior with different values of $\alpha$ between 0 and 1 at the final time, $t^N=1$. Figure 5 also shows that the solution exhibited an anomalous diffusion behavior. The height of the solution profile decreased when $0.1\le \alpha \le 0.5$ and increased when $0.6\le \alpha \le 0.9$. Thus, the solution continuously depends on the $\alpha$-order of the time-fractional derivative.

7. Conclusions and Final Remarks

In this paper, we developed an unconditionally-stable, finite difference method on a non-uniform grid for solving time-fractional advection–diffusion equations involving the Caputo fractional derivative. The Caputo time derivative was discretized by means of a direct generalization of the well-known fractional L1 formula (4) to the case of non-uniform meshes. The L1 formula takes into account the non-local character of the time-fractional operator and allows the order of accuracy in time of the proposed method to be improved by using of the non-uniform time discretization obtained by the mesh (3). We proved the stability and convergence of the proposed method. Numerical experiments were carried out to support the theoretical results. The reported numerical experiments pointed out that the difference method is more accurate on the non-uniform grid than on the uniform mesh, and the convergence order is greater on the non-uniform mesh than on the uniform mesh. Moreover, it is important to note that, in regard to its simplicity, the method is applicable to a wide class of TFADEs occurring in applied sciences.

Recent convergence results [56] have shown how the grading of the mesh and the regularity of the solution affect the order of convergence of the finite difference schemes for time-fractional diffusion equations. One of the future works in this area is to investigate the convergence results of the proposed method on different non-uniform grids, such as, for example, temporal meshes appropriately graded toward the initial time, $t=0$, or using meshes that refine the computational domain according to the fractional order of $\alpha$. Moreover, other future works include investigating the error analysis and accuracy of the method applied by using adaptive procedures to dynamically choose the size of the time-steps according to the local behavior of the solution.