1. Introduction
This work includes the numerical solution of the following generalized fractional advection–diffusion equation,
where
is a bounded domain with boundary
and the notation
denotes the GCFD (defined in [
1] and related references therein) with respect to
t of order
:
where
, parameter
is the diffusivity,
is the advection constant and
U is the solute concentration;
g,
,
and
are continuous functions on their respective domains with
and
. Here scale and weight are sufficiently regular functions and our model (
1) reduces to the diffusion problem when
. The advection–diffusion equation is basically a transport problem that transports a passive scalar quantity in a fluid flow. Due to diffusion and advection, this model represents the physical phenomenon of species concentration for mass transfer and temperature in heat transfer; for more details, refer to [
2,
3,
4,
5,
6], and [
7,
8,
9,
10,
11] for further history and significance of the advection–diffusion equation in physics, chemistry and biology.
The most used fractional derivatives in the problem formulation are the Riemann–Liouville and the Caputo derivatives [
12]. In the year 2012, the generalizations of fractional integrals and derivatives were discussed by Agrawal [
1]. Two functions, scale
and weight
in one parameter, appear in the definition of the generalized fractional derivative of a function
. If
and
then the generalized fractional derivative reduces to the Riemann–Liouville (R-L) and the Caputo derivative; whereas if
,
, and
,
then it will convert to Hadamard [
13], and modified Erdélyi–Kober fractional derivatives, respectively. The authors of [
14] studied the generalized form of R-L and the Hadamard fractional integrals, which is a special case of the Erdélyi–Kober generalized fractional derivative, and some properties of this operator. Atangana and Baleanu [
15] discussed a new fractional derivative with a non-singular kernel and used this derivative in the formation of a fractional heat transfer model. Therefore, we obtain different types of fractional derivatives for different choices of weight and scale functions. In the generalized derivative, the scale function
manages the considered time domain; it can stretch or contract accordingly to capture the phenomena accurately over the desired time range. The weight function
allows the events to be estimated differently at different times.
Over the last decade, many numerical methods were investigated to approximate the Caputo fractional derivative. For example, Mustapha [
16] presented an
approximation formula to solve a fractional reaction–diffusion equation and second-order error bound discussed on non-uniform time meshes. Alikhanov [
17] constructed an
formula to approximate the Caputo fractional time derivative and then used this derived scheme in solving the time fractional diffusion equation with variable coefficients. Abu Arqub [
18] considered reproducing a kernel algorithm for approximate solution of the nonlinear time-fractional PDEs with initial and Robin boundary conditions. Li and Yan [
19] discussed the idea of [
20] (i.e.,
approximation formula for time discretization), and also derived a new time discretization method with accuracy order
and finite element method for spatial discretization. Cao et al. [
21] presented a high-order approximation formula based on the cubic interpolation to approximate the Caputo derivative for the time fractional advection–diffusion equation. Xu and Agrawal [
22] used the finite difference method (FDM) to approximate the GCFD for solving the generalized fractional Burgers equation. Kumar et al. [
23] presented
and
methods to approximate the generalized time fractional derivative which are defined as follows, respectively
where,
where the domain
was discretized into
n equal subintervals, i.e.,
with step-size
, and errors
were shown in [
23]. In this work, we discuss the numerical scheme for GCFD with convergence rate
; for this accuracy we have to assume that
,
,
. This idea is discussed in [
21] for Caputo derivative approximation, but the error bound was discussed directly.
Due to the non-local property of the fractional derivatives, the numerical solution of the fractional partial differential equations is a very difficult task [
24]. Several authors have presented some precise and efficient numerical methods for the fractional advection–diffusion equation. For examples, Zheng et al. [
25] used the finite element method (FEM) for space fractional advection–diffusion equation. Mardani et al. [
26] discussed meshless moving least square method for solving the time-fractional advection–diffusion equation with variable coefficients. Cao et al. [
21] proposed the higher-order approximation of the Caputo derivatives and further applied it in solving the fractional advection–diffusion equation. They used the Lagrange interpolation method to discretize the time derivative and second-order central difference for the spatial derivatives. Li and Cai [
27] considered a three-step process for the Caputo fractional derivative approximation; the first two steps include the shifted Lubich formula derivation for infinite interval then for finite interval, and after that it is generalized to the Caputo derivative. Yadav et al. [
28] discussed the Taylor expansion for the approximation of the generalized time-fractional derivative to solve the generalized fractional advection–diffusion equation. Tian et al. [
29] presented a polynomial spectral collocation method for the space fractional advection–diffusion equation. In [
30], the authors developed explicit and implicit Euler approximations to solve a variable-order fractional advection–diffusion equation on a finite domain. Singh et al. [
31] investigated the numerical approximation of the Caputo–Prabhakar derivative and then used this approximation to solve the fractional advection–diffusion equation. Kannan et al. [
32] used a variant of the local discontinuous Galerkin (LDG2) flux formulation to discuss the high accuracy of the 1D and 2D diffusion equations. In [
33], the authors discussed three approaches, penalty approach, BR2 (Bassi and Rebay) method and LDG method, to study the 2D diffusion equation.
1.1. Motivation and Literature Review for Approximation of the GCFD
Up to now, there are only limited works available in the literature to approximate the GCFD. These works include the finite difference method, -method, -method and collocation method used to approximate GCFD. In this article, we present the approximation of the GCFD based on cubic interpolation polynomials (say -method).
Xu and Agrawal [
22] considered the FDM for approximation of the GCFD for the generalized fractional Burgers equation.
Kumar et al. [
34] presented a numerical scheme for the generalized fractional telegraph equation in time.
Cao et al. [
35] worked on the generalized time-fractional Kdv equation. They used the collocation method which is constructed using Jacobi–Gauss–Lobatto (JGL) nodes.
Xu et al. [
36] considered the FDM to approximate the GCFD of order
and discussed the analytical and numerical solutions of the generalized fractional diffusion equation.
In [
37], the authors used the generalized weighted and shifted Grünwald–Letnikov difference operator to approximate the GCFD and then adopt it to solve the generalized fractional diffusion equation.
Yadav et al. [
28] discussed the Taylor expansion and finite difference approach to approximate the GCFD and then presented the numerical solution of the generalized fractional advection–diffusion equation.
Sultana et al. [
38] presented a numerical scheme based on non-uniform meshes for solution of the generalized time-fractional telegraph equation using quadratic interpolation and compact finite difference approximations for temporal and spatial directions respectively.
Motivated by all these works, the main focus of this paper is to present a much higher-order numerical scheme to approximate the GCFD, and also establish the error analysis in both time and space discretization. To the best of our knowledge, no work has been done yet for a third-order error bound of cubic interpolation formula to approximate the GCFD.
1.2. The Main Contributions of This Work Are as Follows
(1) We extend the approximation method of Cao et al. [
21] for approximating the GCFD and obtain the convergence order
. Further, we show that the obtained scheme reduces to the approximation scheme discussed by Cao et al. [
21] for choice of the scale and the weight functions as
and
.
(2) We establish the full error analysis of the presented higher-order numerical scheme for the generalized time fractional derivative by using the Lagrange interpolation formula.
(3) We introduce some numerical results for different choices of scale and weight functions for the high-order time discretization scheme with convergence order
for all
which achieves higher accuracy than the numerical methods developed in Gao et al. [
39] for
and
.
The remaining sections of the paper are arranged as follows: In
Section 2, we discuss the
-th order scheme to approximate the GCFD of order
of the function
. In
Section 3, a higher-order difference scheme to solve the generalized fractional advection–diffusion equation is presented. Stability and convergence analysis are also discussed in this section. We present three numerical examples which illustrate the error and convergence order of our established numerical scheme in
Section 4. Finally,
Section 5 concludes with some remarks.
2. Numerical Scheme for the Generalized Caputo Fractional Derivative
Motivated by the research carried out in [
21,
23], this section is devoted to presenting a high-order approximation formula for the generalized Caputo-type fractional derivative using cubic interpolation polynomials.
Suppose that
and grid points
with step length
for
. For simplicity, we use
,
and
. The generalized Caputo fractional derivative of order
of the function
at grid point
is given by,
On the first interval
of the domain, we use continuous linear polynomial
to approximate the function
. Let
and the difference operator
for
. Then, we have
Thus, Equation (
5) yields,
where
is the truncation error on the first interval and coefficients for this approximation are
Here, we denote notation
. Then
Remark 1. If the scale function ζ is a positive strictly increasing function on the domain , then the following inequality holdsSince then it implies for also . Remark 2. To estimate , suppose thattherefore, On the second interval
, we use the continuous quadratic polynomial
to approximate the function
, then we get
Thus from Equation (
5), we get,
Here, the truncation error on the second interval is
and
On the other subdomains , we use the cubic interpolation polynomial to approximate the function using four points .
As we know the cubic interpolation polynomial is defined as:
then we get,
Therefore, from Equation (
5), we have,
where
is the truncation error, and
where
. After simplifying Equation (
11), we obtain the following form
Here, we introduce another coefficient
which is defined as
and coefficients
are defined in Equations (
7) and (
10), respectively, and Equation (
16) gives a more compact form of Equation (
11). Such forms of coefficients were missing in [
21]. From this we can easily discuss properties of coefficients.
Motivated by [
21] (developed for Caputo derivative), a new numerical scheme for the generalized Caputo-type fractional derivative of order
of the function
at grid point
, with the help of Equations (
6), (
9) and (
11), is defined by
with
.
Lemma 1. For any and thenwhere, is the approximation of GCFD and . For distinct values of
n, the coefficients in (
18) can be expressed as below.
Lemma 2. If the scale function ζ fulfills (8), and the weight function ω is non-negative and non-decreasing on uniform time grids, then - (i)
Ref. [35] The linear approximation coefficient satisfies, - (ii)
The quadratic approximation coefficient satisfies,
Proof. (i) If scale function
is continuous on its respective domain then, by using mean-value theorem, there exist
such that
Since is a monotone increasing function, we easily get our required result.
(ii) Let
; suppose scale function
is sufficiently smooth on domain
then mean-value theorem yields
Using (
21), we can easily get that
for positive strictly increasing weight function
. Since
is a monotone increasing function on temporal domain
, so we get the desired result. □
Lemma 3. Suppose that the scale function ζ is positive and strictly increasing, and the weight function ω is non-negative and non-decreasing, then, for each , the following conditions hold for coefficients in (19) - (1)
, ,
- (2)
.
Proof. (1) If
, then
as the scale function is strictly increasing, therefore
, for
. This implies that
.
If , then
since
and we have already shown that
, therefore
.
If
then, for
,
(2) If , then for . This implies that .
If
, then there exists an
, by numerical analysis
□
Lemma 4. If scale function ζ is a Lipschitz function on interval with Lipschitz constant L, then Truncation Error for Generalized Caputo Derivative Term
For truncation error of approximation of the generalized Caputo derivative defined in (
18), for simplicity, suppose that function
such that
and
; this implies
. Therefore,
.
Theorem 1. A triangle inequality gives the boundwith Proof. Ref. [
23] For
the scheme (
18) is a linear approximation of the GCFD and the convergence order for this approximation formula is
.
Ref. [
23] For
the scheme (
18) is a quadratic approximation of GCFD and here the convergence order of the approximation formula is
.
For
by the Lagrange interpolation remainder theorem, we use quadratic interpolation function
to interpolate
using node points
,
,
on interval
and cubic interpolation function
depends on
,
,
,
to interpolate
on
as follows
where
.
Now,
Since, from (
22) and (
23),
Consider the first integration of Equation (
24)
Using Lemma (4), we have the following inequality
Consider the second integration of Equation (
24)
Consider the first part of the RHS of Equation (
26)
Here,
is the
. Using Lemma 4, we get the following inequality
Remark 3. If scale function ζ fulfills the condition of the Lipschitz function such that , then is obtained at , therefore Consider the second part of the RHS of Equation (
26)
Now combining (
25), (
27) and (
28), then the error bound is
□
3. Numerical Scheme for the Generalized Fractional Advection–Diffusion Equation
In this section, we study the numerical scheme for solving the generalized fractional advection–diffusion equation defined by (
1).
Let
. We rewrite Equation (
1) to a similar equation using a unity weight function and the same scale function
. Then, Equation (
1) converts into the following form:
where
,
, and
.
For the uniform spatial mesh, let of the interval with step size where M denotes number of subintervals, the grid points and be the step size in the temporal direction with grids .
Now, we discretize our problem (
30) at
, then we get
In Equation (
31), for fixed
and
, the first- and second-order spatial derivatives are discretized by using the following central difference approximations:
With the help of Equation (
18), we get an approximation of the generalized Caputo-type fractional derivative term in (
31) as follows:
where
is defined in Equation (
19). Next, we use Equations (
32)–(
34) in discretized Equation (
31), yielding
where
for some constant
.
Now, we use a numerical approximation
of
to neglect the truncation error term in Equation (
35), then we determine the following finite difference scheme:
where
,
. That is,
We rewrite the matrix form of the above equation as follows:
where the matrix as well as vectors of Equation (
38) are defined as follows:
.
Remark 4. Since the coefficient matrix A is a tridiagonal and strictly diagonally dominant then (Levy–Desplanques theorem). Therefore, at each time level the proposed scheme (37) has a unique solution for . Theorem 2. The local truncation error of difference scheme (36) at , iswhere C is the positive constant independent of the time and space step sizes. Proof. From Equations (4.7) the LTE of difference scheme (
36) is
□
In the next theorem, we use -space with the norm and the inner product .
Theorem 3. (see [22]) If the tridiagonal matrix elements satisfy the inequalitythen the finite difference scheme is stable. Proof. Equation (
38) can be rewritten as follows, for
let,
.
Since matrix
K is invertible then the above equation can be rewritten as
Using the recurrence relation, we get the following equation
Let
be the approximate solution of Equation (
38), then we define the error at grid point
From the definition of the compatible matrix norm
According to the Ostrowski therorem ([
40], Theorem 3.1), let
, then the determinant of
K satisfies
Using assumed inequality (
40) into (
42), we get
Therefore, .
Thus, the numerical scheme is stable. □
Theorem 4. The solution of the difference scheme (36) satisfiesfor some constant . Proof. The truncation error for the difference scheme at
is
by Theorem 2. Now, using the theorem of stability, it implies
We obtain the desired result easily after using the truncation error as discussed in Theorem 2. □
4. Numerical Results
In this section, we will check the numerical accuracy of the proposed schemes (
18) and difference scheme (
37), and also verify the theoretical convergence order discussed in Theorem 4. Here, we provide three examples to numerically support our theory; in the first example we check the convergence order and absolute error of approximation for the GCFD, while the last two problems get the form (
30) to describe the accuracy and maximum absolute error of the difference scheme. All numerical results are implemented in MATLAB R2018b.
To calculate the maximum absolute error
and error
corresponding to the
-norm, we use the following formulas [
19,
21,
39], respectively, due to the fact that the exact solutions of the considered test problems are known.
where
is the exact solution of advection diffusion equation and
is the approximate solution at the point
.
Moreover, the convergence order in the space and time directions for the described difference scheme corresponding to
-norm can be evaluated using the following formulas.
is the order of convergence on the space side and
on the temporal side.
and
Example 1. Ref. [39] Take function , for . Determine the α-th order GCFD for at numerically. The maximum absolute error and rate of convergence for the scheme (
18) to approximate the GCFD of function
for
with uniform time steps
,
,
,
,
are calculated and shown in the following tables.
Table 1 shows the errors with respect to
-norm and convergence rates in time, and these data are found after calculating the classical Caputo derivative (i.e., taking
) of function
with the help of scheme (
18). From this table, we can see that the errors of our scheme (
18) obtained from the GCFD approximation are lesser compared to the scheme developed in Gao et al. [
39] for approximation of the Caputo derivative and the convergence rate of our scheme is
while accuracy for the time derivative in [
39] is
. In
Table 2, to compute the maximum errors
and order of convergence in the time direction, we take
while the scale function is fixed with
t. From
Table 3, we validate the convergence for
and the CPU time in seconds for
are discussed.
Table 4 shows maximum errors and order of convergence for different choices of weight
,
; scale is
and
. In all cases, accuracy in time is obtained as
for scheme (
18), which is higher than [
28,
40].
Figure 1a,b show that the comparison of absolute error for the scheme defined in [
39] and the current scheme (
18) for
and
, respectively.
Example 2. Ref. [21] We take the following generalized fractional advection–diffusion equation:where . When and , then is the exact solution. To solve this example, we use the numerical scheme defined in (
37). The maximum errors at time
for different values of
with different step sizes, and rate of convergence in time direction and space direction are displayed in
Table 5 and
Table 6. In
Table 5, we set
and describe the numerical errors and convergence rates in time
for different values of
N. In
Table 6, we fix
and present the numerical errors and spatial convergence rates
for different values of
M. It is shown that our scheme (
37) gives
-order convergence in the temporal direction and second-order convergence in the spatial direction.
Figure 2 represents the exact and numerical solution of the Example 2 for different values of
.
It is clearly visible from the above
Figure 3a,b that scale function
can stretch or contract the domain.
In
Table 7, we compare our scheme (
37) with Gao et al. ([
39], Example 4.1) in temporal direction to fix
for particular choices of scale
and weight function
.
Example 3. Take the following fractional advection–diffusion equation:where . When and , the exact solution is . To solve this PDE, we use our difference scheme (
37). Here, two
Table 8 and
Table 9 are given in support of the numerical results. In
Table 8, we take fixed space step size
and display maximum absolute errors and errors
with respect to norm
, and also the rate of convergence for the temporal direction with different time steps
. In
Table 9, we express maximum errors and convergence order in the space direction to set time steps fixed with
and taking different
.
In
Table 10, we validate the proposed scheme with [
21] (Example 5.1). Firstly, we present the max error, the convergence order with fixed
and varying
N with 8, 16, 32, 64, 128 for
. It is noted that the convergence order of our scheme (
37) for the temporal direction is almost the same as [
21]. After that, we express maximum error and convergence rate for the spatial dimension to set
and changing
for the same value of
; we observe that the spatial convergence order of our numerical scheme is same as [
21]. Thus, the scheme presented in [
21] becomes a particular case of the proposed scheme (
37) for
and
. Furthermore, we can compute numerical results for different suitable choices of scale and weight functions.