Next Article in Journal
Numerical Study on the Orthogonality of the Fields Radiated by an Aperture
Next Article in Special Issue
On Solutions of Fractional Integrodifferential Systems Involving Ψ-Caputo Derivative and Ψ-Riemann–Liouville Fractional Integral
Previous Article in Journal
Finite-Time Bounded Tracking Control for a Class of Neutral Systems
Previous Article in Special Issue
On Stability Criteria Induced by the Resolvent Kernel for a Fractional Neutral Linear System with Distributed Delays
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

High-Order Approximation to Generalized Caputo Derivatives and Generalized Fractional Advection–Diffusion Equations

1
Department of Mathematical Sciences, Indian Institute of Technology (BHU) Varanasi, Varanasi 221005, Uttar Pradesh, India
2
Department of Mathematics, Texas A & M University-Kingsville, Kingsville, TX 77553, USA
*
Authors to whom correspondence should be addressed.
Mathematics 2023, 11(5), 1200; https://doi.org/10.3390/math11051200
Submission received: 21 December 2022 / Revised: 19 February 2023 / Accepted: 23 February 2023 / Published: 28 February 2023
(This article belongs to the Special Issue Stability Analysis of Fractional Systems-II)

Abstract

:
In this article, a high-order time-stepping scheme based on the cubic interpolation formula is considered to approximate the generalized Caputo fractional derivative (GCFD). Convergence order for this scheme is ( 4 α ) , where α ( 0 < α < 1 ) is the order of the GCFD. The local truncation error is also provided. Then, we adopt the developed scheme to establish a difference scheme for the solution of the generalized fractional advection–diffusion equation with Dirichlet boundary conditions. Furthermore, we discuss the stability and convergence of the difference scheme. Numerical examples are presented to examine the theoretical claims. The convergence order of the difference scheme is analyzed numerically, which is ( 4 α ) in time and second-order in space.

1. Introduction

This work includes the numerical solution of the following generalized fractional advection–diffusion equation,
0 C D t ; [ ζ ( t ) , ω ( t ) ] α U ( x , t ) = D 2 U ( x , t ) x 2 A U ( x , t ) x + g ( x , t ) , x Ω , t ( 0 , T ] , U ( x , 0 ) = U 0 ( x ) , x Ω ¯ = Ω Ω , U ( 0 , t ) = ρ 1 ( t ) , U ( a , t ) = ρ 2 ( t ) , t ( 0 , T ] ,
where Ω = ( 0 , a ) is a bounded domain with boundary Ω and the notation 0 C D t , [ ζ ( t ) , ω ( t ) ] α denotes the GCFD (defined in [1] and related references therein) with respect to t of order α :
0 C D t ; [ ζ ( t ) , ω ( t ) ] α U ( t ) = [ ω ( t ) ] 1 Γ ( 1 α ) 0 t [ ω ( s ) U ( s ) ] [ ζ ( t ) ζ ( s ) ] α d s , 0 < α < 1 ,
where s = s , parameter D > 0 is the diffusivity, A > 0 is the advection constant and U is the solute concentration; g, U 0 , ρ 1 and ρ 2 are continuous functions on their respective domains with U 0 ( 0 ) = ρ 1 ( 0 ) and U 0 ( a ) = ρ 2 ( 0 ) . Here scale and weight are sufficiently regular functions and our model (1) reduces to the diffusion problem when A = 0 . 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 ζ ( t ) and weight ω ( t ) in one parameter, appear in the definition of the generalized fractional derivative of a function V ( t ) . If ω ( t ) = 1 and ζ ( t ) = t then the generalized fractional derivative reduces to the Riemann–Liouville (R-L) and the Caputo derivative; whereas if ω ( t ) = 1 , ζ ( t ) = l n ( t ) , and ω ( t ) = t σ η , ζ ( t ) = t σ 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 ζ ( t ) manages the considered time domain; it can stretch or contract accordingly to capture the phenomena accurately over the desired time range. The weight function ω ( t ) 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 L 1 approximation formula to solve a fractional reaction–diffusion equation and second-order error bound discussed on non-uniform time meshes. Alikhanov [17] constructed an L 2 1 σ 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., L 2 approximation formula for time discretization), and also derived a new time discretization method with accuracy order O ( τ 3 α ) 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 L 1 and L 2 methods to approximate the generalized time fractional derivative which are defined as follows, respectively
0 C D t ; [ ζ ( t ) , ω ( t ) ] α V ( t ) | t = t n = [ ω ( t ) ] 1 Γ ( 1 α ) l = 1 n ω l V l ω l 1 V l 1 ζ l ζ l 1 t l 1 t l [ ζ ( t ) ζ ( s ) ] α ( s ζ ( s ) ) d s + r 1 n ,
0 C D t ; [ ζ ( t ) , ω ( t ) ] α V ( t ) | t = t n = [ ω ( t ) ] 1 Γ ( 1 α ) l = 1 n t l 1 t l [ ζ ( t ) ζ ( s ) ] α s ( Π l ω ( s ) V ( s ) ) d s + r 2 n ,
where,
Π l ( ω ( t ) V ( t ) ) = ζ ( t ) { [ 2 ζ ( t ) ζ l ζ l 1 ( ζ l 2 ζ l 1 ) ( ζ l 2 ζ l 1 ) ] ω l 2 V l 2 + 2 ζ ( t ) ζ l ζ l 2 ( ζ l 1 ζ l 2 ) ( ζ l 1 ζ l ) ω l 1 V l 1 + 2 ζ ( t ) ζ l 1 ζ l 2 ( ζ l ζ l 2 ) ( ζ l ζ l 1 ) ω l V l } ,
where the domain [ 0 , T ] was discretized into n equal subintervals, i.e., 0 = t 0 < t 1 < < t N = T with step-size τ = T N , and errors r 1 n = O ( τ 2 α ) , r 2 n = O ( τ 3 α ) were shown in [23]. In this work, we discuss the numerical scheme for GCFD with convergence rate ( 4 α ) ; for this accuracy we have to assume that V ( t 0 ) = 0 , V ( t 0 ) = 0 , V ( t 0 ) = 0 . 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, L 1 -method, L 1 2 -method and collocation method used to approximate GCFD. In this article, we present the approximation of the GCFD based on cubic interpolation polynomials (say L 1 2 3 -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 0 < α < 1 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 ( 4 α ) . 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 ζ ( t ) = t and ω ( t ) = 1 .
(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 O ( τ 4 α ) for all α ( 0 , 1 ) which achieves higher accuracy than the numerical methods developed in Gao et al. [39] for ζ ( t ) = t and ω ( t ) = 1 .
The remaining sections of the paper are arranged as follows: In Section 2, we discuss the ( 4 α ) -th order scheme to approximate the GCFD of order α of the function V . 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 V ( t ) C 4 [ 0 , T ] and grid points 0 = t 0 < t 1 < < t N = T with step length τ = t n t n 1 for 1 n N . For simplicity, we use g ( s ) = ω ( s ) V ( s ) , V ( t l ) = V l , ω ( t l ) = ω l and ζ ( t l ) = ζ l . The generalized Caputo fractional derivative of order α of the function V ( t ) at grid point t n is given by,
0 C D t ; [ ζ ( t ) , ω ( t ) ] α V n = [ ω n ] 1 Γ ( 1 α ) l = 1 n t l 1 t l g ( s ) [ ζ ( t n ) ζ ( s ) ] α d s .
On the first interval [ 0 , t 1 ] of the domain, we use continuous linear polynomial Π 1 g ( t ) to approximate the function g ( t ) ( = ω ( t ) V ( t ) ) . Let g ( t l ) = g l and the difference operator τ g l = g l g l 1 for l 1 . Then, we have
( Π 1 g ( t ) ) = ( g 1 g 0 ) τ = ( τ g 1 ) τ .
Thus, Equation (5) yields,
[ ω n ] 1 Γ ( 1 α ) 0 t 1 [ ζ ( t n ) ζ ( s ) ] α g ( s ) d s = [ ω n ] 1 Γ ( 1 α ) ( τ g 1 ) τ 0 t 1 [ ζ ( t n ) ζ ( s ) ] α d s + E τ 1 = a n 1 ( τ g 1 ) + E τ 1 ,
where E τ 1 is the truncation error on the first interval and coefficients for this approximation are
a n l = [ ω n ] 1 Γ ( 2 α ) [ ζ ( t n ) ζ ( t l 1 ) ] 1 α [ ζ ( t n ) ζ ( t l ) ] 1 α [ ζ ( t l ) ζ ( t l 1 ) ] .
Here, we denote notation α 0 ( n ) = [ ω n ] 1 Γ ( 2 α ) , 1 n N . Then
a n l = α 0 ( n ) [ ζ ( t n ) ζ ( t l 1 ) ] 1 α [ ζ ( t n ) ζ ( t l ) ] 1 α [ ζ ( t l ) ζ ( t l 1 ) ] , 1 l n .
Remark 1.
If the scale function ζ is a positive strictly increasing function on the domain [ 0 , T ] , then the following inequality holds
0 < [ ζ ( t n ) ζ ( t l ) ] 1 α [ ζ ( t n ) ζ ( t l 1 ) ] 1 α , 1 l n .
Since t l > t l 1 , then it implies ζ ( t l ) > ζ ( t l 1 ) for 1 l n , also ( 1 α ) > 0 .
Remark 2.
To estimate a n 1 , suppose that
Θ = [ ζ ( t n ) ζ ( s ) ] d Θ = ζ ( s ) d s = ζ ( t l ) ζ ( t l 1 ) τ d s , s ( t l 1 , t l ) ,
therefore,
[ ζ ( t n ) ζ ( s ) ] α d s = τ ζ ( t l ) ζ ( t l 1 ) [ ζ ( t n ) ζ ( s ) ] 1 α ( 1 α ) .
On the second interval [ t 1 , t 2 ] , we use the continuous quadratic polynomial Π 2 g ( t ) to approximate the function g ( t ) , then we get
( Π 2 g ( t ) ) = ( 2 t t 0 t 1 ) 2 τ 2 g 2 ( 2 t t 0 t 2 ) τ 2 g 1 + ( 2 t t 1 t 2 ) 2 τ 2 g 0 .
Thus from Equation (5), we get,
[ ω n ] 1 Γ ( 1 α ) t 1 t 2 [ ζ ( t n ) ζ ( s ) ] α g ( s ) d s = [ ω n ] 1 Γ ( 1 α ) t 1 t 2 [ ζ ( t n ) ζ ( s ) ] α ( Π 2 g ( s ) ) d s + E τ 2 = a n 2 ( g 2 g 1 ) + b n 2 ( g 2 2 g 1 + g 0 ) + E τ 2 = a n 2 τ g 2 + b n 2 τ 2 g 2 + E τ 2 .
Here, the truncation error on the second interval is E τ 2 , and
b n l = α 0 ( n ) { 1 ( 2 α ) [ ζ ( t n ) ζ ( t l 1 ) ] 2 α [ ζ ( t n ) ζ ( t l ) ] 2 α [ ζ ( t l ) ζ ( t l 1 ) ] 2 1 2 [ ζ ( t n ) ζ ( t l 1 ) ] 1 α + [ ζ ( t n ) ζ ( t l ) ] 1 α [ ζ ( t l ) ζ ( t l 1 ) ] } , 2 l n .
On the other subdomains ( l 3 ) , we use the cubic interpolation polynomial Π l g ( t ) to approximate the function g ( t ) using four points ( t l 3 , g l 3 ) , ( t l 2 , g l 2 ) , ( t l 1 , g l 1 ) , ( t l , g l ) .
As we know the cubic interpolation polynomial is defined as:
Π l g ( t ) = r = 0 3 g l r s = 0 , s r 3 t t l s t l r t l s ,
then we get,
( Π l g ( t ) ) = g l 3 ( t t l 2 ) ( t l + t l 1 2 t ) + ( t t l 1 ) ( t l t ) 6 τ 3 + g l 2 ( t t l 3 ) ( 2 t t l 1 t l ) + ( t t l ) ( t t l 1 ) 2 τ 3 + g l 1 ( t t l 2 ) ( t l 3 + t l 2 t ) + ( t t l 3 ) ( t l t ) 2 τ 3 + g l ( t t l 2 ) ( 2 t t l 3 t l 1 ) + ( t t l 1 ) ( t t l 3 ) 6 τ 3 .
Therefore, from Equation (5), we have,
[ ω n ] 1 Γ ( 1 α ) l = 3 n t l 1 t l [ ζ ( t n ) ζ ( s ) ] α g ( s ) d s = [ ω n ] 1 Γ ( 1 α ) l = 3 n t l 1 t l [ ζ ( t n ) ζ ( s ) ] α ( Π l g ( s ) ) d s + E τ n = l = 3 n A 1 , n l g l + A 2 , n l g l 1 + A 3 , n l g l 2 + A 4 , n l g l 3 + E τ n ,
where E τ n ( 3 n N ) is the truncation error, and
A 1 , n l = α 0 ( n ) { 1 6 2 [ ζ ( t n ) ζ ( t l 1 ) ] 1 α 11 [ ζ ( t n ) ζ ( t l ) ] 1 α [ ζ ( t l ) ζ ( t l 1 ) ] + 1 ( 2 α ) [ ζ ( t n ) ζ ( t l 1 ) ] 2 α 2 [ ζ ( t n ) ζ ( t l ) ] 2 α [ ζ ( t l ) ζ ( t l 1 ) ] 2 + 1 ( 2 α ) ( 3 α ) [ ζ ( t n ) ζ ( t l 1 ) ] 3 α [ ζ ( t n ) ζ ( t l ) ] 3 α [ ζ ( t l ) ζ ( t l 1 ) ] 3 } ,
A 2 , n l = α 0 ( n ) { 1 2 6 [ ζ ( t n ) ζ ( t l ) ] 1 α + [ ζ ( t n ) ζ ( t l 1 ) ] 1 α [ ζ ( t l ) ζ ( t l 1 ) ] + 1 ( 2 α ) 5 [ ζ ( t n ) ζ ( t l ) ] 2 α 2 [ ζ ( t n ) ζ ( t l 1 ) ] 2 α [ ζ ( t l ) ζ ( t l 1 ) ] 2 + 3 ( 2 α ) ( 3 α ) [ ζ ( t n ) ζ ( t l ) ] 3 α [ ζ ( t n ) ζ ( t l 1 ) ] 3 α [ ζ ( t l ) ζ ( t l 1 ) ] 3 } ,
A 3 , n l = α 0 ( n ) { 1 2 2 [ ζ ( t n ) ζ ( t l 1 ) ] 1 α + 3 [ ζ ( t n ) ζ ( t l ) ] 1 α [ ζ ( t l ) ζ ( t l 1 ) ] + 1 ( 2 α ) [ ζ ( t n ) ζ ( t l 1 ) ] 2 α 4 [ ζ ( t n ) ζ ( t l ) ] 2 α [ ζ ( t l ) ζ ( t l 1 ) ] 2 + 3 ( 2 α ) ( 3 α ) [ ζ ( t n ) ζ ( t l 1 ) ] 3 α [ ζ ( t n ) ζ ( t l ) ] 3 α [ ζ ( t l ) ζ ( t l 1 ) ] 3 } ,
A 4 , n l = α 0 ( n ) { 1 6 [ ζ ( t n ) ζ ( t l 1 ) ] 1 α + 2 [ ζ ( t n ) ζ ( t l ) ] 1 α [ ζ ( t l ) ζ ( t l 1 ) ] + 1 ( 2 α ) [ ζ ( t n ) ζ ( t l ) ] 2 α [ ζ ( t l ) ζ ( t l 1 ) ] 2 + 1 ( 2 α ) ( 3 α ) [ ζ ( t n ) ζ ( t l ) ] 3 α [ ζ ( t n ) ζ ( t l 1 ) ] 3 α [ ζ ( t l ) ζ ( t l 1 ) ] 3 } ,
where 3 l n . After simplifying Equation (11), we obtain the following form
[ ω n ] 1 Γ ( 1 α ) l = 3 n t l 1 t l [ ζ ( t n ) ζ ( s ) ] α g ( s ) d s = l = 3 n [ a n l ( g l g l 1 ) + b n l ( g l 2 g l 1 + g l 2 ) + c n l ( g l 3 g l 1 + 3 g l 2 g l 3 ) ] + E τ n = l = 3 n a n l ( τ g l ) + b n l ( τ 2 g l ) + c n l ( τ 3 g l ) + E τ n .
Here, we introduce another coefficient c n l , which is defined as
c n l = α 0 ( n ) { 1 ( 2 α ) ( 3 α ) [ ζ ( t n ) ζ ( t l 1 ) ] 3 α [ ζ ( t n ) ζ ( t l ) ] 3 α [ ζ ( t l ) ζ ( t l 1 ) ] 3 1 ( 2 α ) [ ζ ( t n ) ζ ( t l ) ] 2 α [ ζ ( t l ) ζ ( t l 1 ) ] 2 1 6 [ ζ ( t n ) ζ ( t l 1 ) ] 1 α + 2 [ ζ ( t n ) ζ ( t l ) ] 1 α [ ζ ( t l ) ζ ( t l 1 ) ] } , 3 l n ,
and coefficients a n l , b n l 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 V ( t ) at grid point t n , with the help of Equations (6), (9) and (11), is defined by
H D t ; [ ζ ( t ) , ω ( t ) ] α V ( t ) | t = t n = [ ω n ] 1 Γ ( 1 α ) 0 t n g ( s ) [ ζ ( t n ) ζ ( s ) ] α d s = [ ω n ] 1 Γ ( 1 α ) [ 0 t 1 g ( s ) [ ζ ( t n ) ζ ( s ) ] α d s + t 1 t 2 g ( s ) [ ζ ( t n ) ζ ( s ) ] α d s + l = 3 n l 1 t l g ( s ) [ ζ ( t n ) ζ ( s ) ] α d s ] = l = 0 n λ l ω n l V n l ,
with g ( s ) = ω ( s ) V ( s ) .
Lemma 1.
For any α ( 0 , 1 ) and V ( t ) C 4 [ 0 , T ] , then
0 C D t ; [ ζ ( t ) , ω ( t ) ] α V n = H D t ; [ ζ ( t ) , ω ( t ) ] α V n + E τ n , n = 1 , 2 , , N ,
where, H D t ; [ ζ ( t ) , ω ( t ) ] α is the approximation of GCFD and | E τ n | = O ( τ 4 α ) .
For distinct values of n, the coefficients in (18) can be expressed as below.
For n = 1 ,
λ 0 = a 0 , λ 1 = a 0 .
For n = 2 ,
λ 0 = a 0 + b 0 , λ 1 = a 1 a 0 2 b 0 , λ 2 = a 1 + b 0 .
For n = 3 ,
λ 0 = A 1 , 0 , λ 1 = A 2 , 0 + a 1 + b 1 , λ 2 = A 3 , 0 + a 2 a 1 2 b 1 , λ 3 = A 4 , 0 a 2 + b 1 .
For n = 4 ,
λ 0 = A 1 , 0 , λ 1 = A 1 , 1 + A 2 , 0 , λ 2 = A 2 , 1 + A 3 , 0 + a 2 + b 2 , λ 3 = A 3 , 1 + A 4 , 0 + a 3 a 2 2 b 2 , λ 4 = A 4 , 1 a 3 + b 2 .
For n = 5 ,
λ 0 = A 1 , 0 , λ 1 = A 1 , 1 + A 2 , 0 , λ 2 = A 1 , 2 + A 2 , 1 + A 3 , 0 , λ 3 = A 2 , 2 + A 3 , 1 + A 4 , 0 + a 3 + b 3 , λ 4 = A 3 , 2 + A 4 , 1 + a 4 a 3 2 b 3 , λ 5 = A 4 , 2 a 4 + b 3 .
For n 6 ,
λ 0 = A 1 , 0 , λ 1 = A 1 , 1 + A 2 , 0 , λ 2 = A 1 , 2 + A 2 , 1 + A 3 , 0 , λ l = A 1 , l + A 2 , l 1 + A 3 , l 2 + A 4 , l 3 ( 3 l n 3 ) , λ n 2 = a n 2 + b n 2 + A 2 , n 3 + A 3 , n 4 + A 4 , n 5 , λ n 1 = a n 1 a n 2 2 b n 2 + A 3 , n 3 + A 4 , n 4 , λ n = a n 1 + b n 2 + A 4 , n 3 .
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, 1 l n ,
0 < a n 1 < < a n l < a n l 1 < < a 1 < a 0 .
(ii) 
The quadratic approximation coefficient satisfies, 2 l n ,
0 < b n 2 < < b n l < b n l 1 < < b 1 < b 0 .
Proof. 
(i) If scale function ζ is continuous on its respective domain then, by using mean-value theorem, there exist x ^ l ( t l 1 , t l ) such that
1 τ t l 1 t l [ ζ ( t n ) ζ ( s ) ] α d s = [ ζ ( t n ) ζ ( x ^ l ) ] α , 1 l n .
Since [ ζ ( t n ) ζ ( s ) ] α is a monotone increasing function, we easily get our required result.
(ii) Let η ( s ) = [ ζ ( t n ) ζ ( s ) ] 1 α ; suppose scale function ζ is sufficiently smooth on domain [ 0 , T ] then mean-value theorem yields
2 t l 1 t l η ( s ) d s η ( t l 1 ) + η ( t l ) = 2 η ( x ˜ l ) η ( t l 1 ) + η ( t l ) , x l ˜ ( t l 1 , t l ) , = θ η ( t l ) η ( t l 1 ) , 0 < θ < 1 = θ τ η ( ν l ) , ν l ( t l 1 , t l ) , = θ α ( 1 α ) τ [ ζ ( t l ) ζ ( t l 1 ) ] 2 [ ζ ( t n ) ζ ( ν l ) ] α 1 > 0 .
Using (21), we can easily get that b n l > 0 for positive strictly increasing weight function ω ( t ) . Since [ ζ ( t n ) ζ ( s ) ] α 1 is a monotone increasing function on temporal domain [ 0 , T ] , 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 α ( 0 , 1 ) , the following conditions hold for coefficients in (19)
(1) 
λ 0 > 0 , n 1 ,
(2) 
l = 0 n λ l = 0 .
Proof. 
(1) If n = 1 , then
λ 0 = a 0 = α 0 ( n ) [ ζ 1 ζ 0 ] α ,
as the scale function is strictly increasing, therefore ζ n 1 < ζ n , for n 1 . This implies that λ 0 = a 0 > 0 .
If n = 2 , then λ 0 = a 0 + b 0 ,
since
b 0 = α 0 ( n ) 1 2 α 1 2 [ ζ 2 ζ 1 ] α > 0 ,
and we have already shown that a 0 > 0 , therefore λ 0 = a 0 + b 0 > 0 .
If n 3 then, for 0 < α < 1 ,
λ 0 = A 1 , 0 = α 0 ( n ) 1 3 + 1 ( 2 α ) + 1 ( 2 α ) ( 3 α ) [ ζ n ζ n 1 ] α > 0 .
(2) If n = 1 , then λ 0 = λ 1 = a 0 > 0 , for 0 < α < 1 . This implies that λ 0 + λ 1 = 0 .
If n = 2 , then there exists an α ( 0 , 1 ) , by numerical analysis
λ 0 + λ 1 + λ 2 = ( a 0 + b 0 ) + ( a 1 a 0 + 2 b 0 ) + ( a 1 + b 0 ) = 0 .
If n 3 , then
l = 0 n λ l = A 1 , 0 + A 1 , 1 + A 2 , 0 + A 1 , 2 + A 2 , 1 + A 3 , 0 + l = 3 n 3 ( A 1 , l + A 2 , l 1 + A 3 , l 2 + A 4 , l 3 ) + a n 2 + b n 2 + A 2 , n 3 + A 3 , n 4 + A 4 , n 5 + a n 1 a n 2 2 b n 2 + A 3 , n 3 + A 4 , n 4 a n 1 + b n 2 + A 4 , n 3 = j = 0 n 3 ( A 1 , j + A 2 , j + A 3 , j + A 4 , j ) = 0 .
Lemma 4.
If scale function ζ is a Lipschitz function on interval [ t l 1 , t l ] with Lipschitz constant L, then
| ζ l ζ l 1 | L τ , 1 l n .

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 g ( t ) = ω ( t ) V ( t ) such that g ( t ) C 4 ( ( 0 , T ] ) and ζ ( t ) = v ; this implies t = ζ 1 ( v ) . Therefore, g ( v ) = ω ( ζ 1 ( v ) ) V ( ζ 1 ( v ) ) .
Theorem 1.
A triangle inequality gives the bound
| 0 C D t ; [ ζ ( t ) , ω ( t ) ] α g n H D t ; [ ζ ( t ) , ω ( t ) ] α g n | n = 1 N | E τ n | ,
with
E τ n = [ ω n ] 1 Γ ( 1 α ) ζ l 1 ζ l [ ζ n v ] α [ g ( v ) Π l g ( v ) ] d v , 1 l n .
( 1 ) | E τ 1 | α [ ω n ] 1 Γ ( 1 α ) 1 8 α + 1 2 ( 1 α ) ( 2 α ) max t 0 v s . t 1 | g ( 2 ) ( v ) | L 2 α ( τ ) 2 α , n = 1 , ( 2 ) | E τ 2 | α [ ω n ] 1 Γ ( 1 α ) { 1 12 max t 0 v s . t 1 | g ( 2 ) ( v ) | ( t 2 t 1 ) α 1 L 2 α ( τ ) 3 + [ 1 12 + 1 3 ( 1 α ) ( 2 α ) 1 2 + 1 ( 3 α ) ] max t 0 v s . t 2 | g ( 3 ) ( v ) | L 3 α ( τ ) 3 α } , n = 2 , ( 3 ) | E τ n | α [ ω n ] 1 Γ ( 1 α ) { 1 72 max t 0 x t 2 | g ( 3 ) ( x ) | ( t n t 2 ) α 1 L 3 α ( τ ) 4 + [ 3 128 α + 1 12 ( 1 α ) ( 2 α ) 1 + 3 ( 3 α ) + 3 ( 3 α ) ( 4 α ) ] max t 0 x 1 t n | g ( 4 ) ( x 1 ) | L 4 α ( τ ) 4 α } , n 3 .
Proof. 
( 1 ) Ref. [23] For n = 1 , the scheme (18) is a linear approximation of the GCFD and the convergence order for this approximation formula is O ( τ 2 α ) .
( 2 ) Ref. [23] For n = 2 , the scheme (18) is a quadratic approximation of GCFD and here the convergence order of the approximation formula is O ( τ 3 α ) .
( 3 ) For n 3 , by the Lagrange interpolation remainder theorem, we use quadratic interpolation function Π 2 g ( v ) to interpolate g ( v ) using node points ( ζ 0 , g 0 ) , ( ζ 1 , g 1 ) , ( ζ 2 , g 2 ) on interval [ ζ 0 , ζ 2 ] and cubic interpolation function Π l g ( v ) depends on ( ζ l 3 , g l 3 ) , ( ζ l 2 , g l 2 ) , ( ζ l 1 , g l 1 ) , ( ζ l , g l ) to interpolate g ( v ) on [ ζ l 3 , ζ l ] as follows
g ( v ) Π 2 g ( v ) = g ( 3 ) ( η 1 ) 3 ! ( v ζ 0 ) ( v ζ 1 ) ( v ζ 2 ) , v s . [ ζ 0 , ζ 2 ] , η 1 ( ζ 0 , ζ 2 ) ,
g ( v ) Π l g ( v ) = g ( 4 ) ( η l ) 4 ! ( v ζ l 3 ) ( v ζ l 2 ) ( v ζ l 1 ) ( v ζ l ) , v s . [ ζ l 3 , ζ l ] , η l ( ζ l 3 , ζ l ) ,
where 3 l n .
Now,
E τ n = [ ω n ] 1 Γ ( 1 α ) [ ζ 0 ζ 2 [ g ( v ) Π 2 g ( v ) ] [ ζ n v ] α d v + l = 3 n ζ l 1 ζ l [ g ( v ) Π l g ( v ) ] [ ζ n v ] α d v ] = [ ω n ] 1 Γ ( 1 α ) { [ g ( v ) Π 2 g ( v ) ] [ ζ n v ] α | ζ 0 ζ 2 α ζ 0 ζ 2 [ g ( v ) Π 2 g ( v ) ] [ ζ n v ] α 1 d v + l = 3 n [ g ( v ) Π l g ( v ) ] [ ζ n v ] α | ζ l 1 ζ l α ζ l 1 ζ l [ g ( v ) Π l g ( v ) ] [ ζ n v ] α 1 d v ] } = α [ ω n ] 1 Γ ( 1 α ) { ζ 0 ζ 2 [ g ( v ) Π 2 g ( v ) ] [ ζ n v ] α 1 d v + l = 3 n ζ l 1 ζ l [ g ( v ) Π l g ( v ) ] [ ζ n v ] α 1 d v }
Since, from (22) and (23),
[ g ( v ) Π 2 g ( v ) ] [ ζ n v ] α | ζ 0 ζ 2 = g ( 3 ) ( η 1 ) 3 ! ( v ζ 0 ) ( v ζ 1 ) ( v ζ 2 ) [ ζ n v ] α | ζ 0 ζ 2 = 0 , [ g ( v ) Π l g ( v ) ] [ ζ n v ] α | ζ l 1 ζ l = g ( 4 ) ( η l ) 4 ! ( v ζ l 3 ) ( v ζ l 2 ) ( v ζ l 1 ) ( v ζ l ) | ζ l 1 ζ l = 0 .
Consider the first integration of Equation (24)
| ζ 0 ζ 2 [ g ( v ) Π 2 ( v ) ] [ ζ n v ] α 1 d v | = | ζ 0 ζ 2 g ( 3 ) ( η 1 ) 3 ! ( v ζ 0 ) ( v ζ 1 ) ( v ζ 2 ) [ ζ n v ] α 1 d v | 1 6 max ζ 0 η 1 ζ 2 | g ( 3 ) ( η 1 ) | ( ζ n ζ 2 ) α 1 ( ζ 0 ζ 2 ) 3 ( ζ 0 2 ζ 1 + ζ 2 ) 12 1 6 max ζ 0 η 1 ζ 2 | g ( 3 ) ( η 1 ) | ( ζ n ζ 2 ) α 1 ( ζ 0 ζ 2 ) 3 ( ζ 0 ζ 1 ) 12 .
Using Lemma (4), we have the following inequality
| ζ 0 ζ 2 [ g ( v ) Π 2 g ( v ) ] [ ζ n v ] α 1 d v | 1 72 max t 0 x t 2 | g ( 3 ) ( x ) | ( t n t 2 ) α 1 L 3 α ( τ ) 4 .
Consider the second integration of Equation (24)
| l = 3 n ζ l 1 ζ l [ g ( v ) Π l g ( v ) ] [ ζ n v ] α 1 d v | | l = 3 n 1 ζ l 1 ζ l g ( 4 ) ( η l ) 4 ! ( v ζ l 3 ) ( v ζ l 2 ) ( v ζ l 1 ) ( v ζ l ) [ ζ n v ] α 1 d v | + | ζ n 1 ζ n g ( 4 ) ( η n ) 4 ! ( v ζ n 3 ) ( v ζ n 2 ) ( v ζ n 1 ) ( v ζ n ) [ ζ n v ] α 1 d v | .
Consider the first part of the RHS of Equation (26)
| l = 3 n 1 ζ l 1 ζ l g ( 4 ) ( η l ) 4 ! ( v ζ l 3 ) ( v ζ l 2 ) ( v ζ l 1 ) ( v ζ l ) [ ζ n v ] α 1 d v | 1 24 max ζ 2 η 2 ζ n 1 | g ( 4 ) ( η 2 ) | f ˜ ( ζ l 3 , ζ l 2 , ζ l 1 , ζ l ) ζ 2 ζ n 1 [ ζ n v ] α 1 d v 1 24 α max ζ 2 η 2 ζ n 1 | g ( 4 ) ( η 2 ) | f ˜ ( ζ l 3 , ζ l 2 , ζ l 1 , ζ l ) ( ζ n ζ n 1 ) α .
Here, f ˜ ( ζ l 3 , ζ l 2 , ζ l 1 , ζ l ) is the max ζ 2 v s . ζ n 1 [ ( v ζ l 3 ) ( v ζ l 2 ) ( v ζ l 1 ) ( v ζ l ) ] . Using Lemma 4, we get the following inequality
| l = 3 n 1 ζ l 1 ζ l g ( 4 ) ( η l ) 4 ! ( v ζ l 3 ) ( v ζ l 2 ) ( v ζ l 1 ) ( v ζ l ) [ ζ n v ] α 1 d v | 3 128 α max t 2 x 1 t n 1 | g ( 4 ) ( x 1 ) | L 4 α ( τ ) 4 α .
Remark 3.
If scale function ζ fulfills the condition of the Lipschitz function such that | ζ l ζ l 1 | L τ , then max ζ 2 v s . ζ n 1 | ( v ζ l 3 ) ( v ζ l 2 ) ( v ζ l 1 ) ( v ζ l ) | is obtained at v = ζ l 2 + L τ 2 , therefore
f ˜ ( ζ l 3 , ζ l 2 , ζ l 1 , ζ l ) = max ζ 2 v s . ζ n 1 [ ( v ζ l 3 ) ( v ζ l 2 ) ( v ζ l 1 ) ( v ζ l ) ] 9 16 L 4 ( τ ) 4 .
Consider the second part of the RHS of Equation (26)
| ζ n 1 ζ n g ( 4 ) ( η n ) 4 ! ( v ζ n 3 ) ( v ζ n 2 ) ( v ζ n 1 ) ( v ζ n ) [ ζ n v ] α 1 d v | 1 24 max ζ n 1 η 3 ζ n | g ( 4 ) ( η 3 ) | | ζ n 1 ζ n ( v ζ n 3 ) ( v ζ n 2 ) ( v ζ n 1 ) [ ζ n v ] α d v | 1 12 ( 1 α ) ( 2 α ) 1 + 3 ( 3 α ) + 3 ( 3 α ) ( 4 α ) max ζ n 1 η 3 ζ n | g ( 4 ) ( η 3 ) | ( ζ n ζ n 1 ) 4 α 1 12 ( 1 α ) ( 2 α ) 1 + 3 ( 3 α ) + 3 ( 3 α ) ( 4 α ) max t n 1 x 2 t n | g ( 4 ) ( x 2 ) | L 4 α ( τ ) 4 α .
Now combining (25), (27) and (28), then the error bound is
| E τ n | α [ ω n ] 1 Γ ( 1 α ) { 1 72 max t 0 x t 2 | g ( 3 ) ( x ) | ( t n t 2 ) α 1 L 3 α ( τ ) 4 + [ 3 128 α + 1 12 ( 1 α ) ( 2 α ) ( 1 + 3 ( 3 α ) + 3 ( 3 α ) ( 4 α ) ) ] max t 0 x 1 t n | g ( 4 ) ( x 1 ) | L 4 α ( τ ) 4 α } .

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 u ( x , t ) = ω ( t ) U ( x , t ) ω ( 0 ) U 0 ( x ) . We rewrite Equation (1) to a similar equation using a unity weight function and the same scale function ζ ( t ) . Then, Equation (1) converts into the following form:
0 C D t ; [ ζ ( t ) , 1 ] α u ( x , t ) = D 2 u ( x , t ) x 2 A u ( x , t ) x + f ( x , t ) , x Ω , t ( 0 , T ] , u ( x , 0 ) = 0 , x Ω ¯ = Ω Ω , u ( 0 , t ) = ϕ 1 ( t ) , u ( a , t ) = ϕ 2 ( t ) , t ( 0 , T ] ,
where ϕ 1 ( t ) = ω ( t ) ρ 1 ( t ) ω ( 0 ) ρ 1 ( 0 ) , ϕ 2 ( t ) = ω ( t ) ρ 2 ( t ) ω ( 0 ) ρ 2 ( 0 ) , and f ( x , t ) = D ω ( 0 ) U 0 ( x ) A ω ( 0 ) U 0 ( x ) + ω ( t ) g ( x , t ) .
For the uniform spatial mesh, let 0 = x 0 < x 1 < < x M = a of the interval [ 0 , a ] with step size h = a M where M denotes number of subintervals, the grid points x 0 + i h ( 0 i M ) , and τ = T N be the step size in the temporal direction with grids t n = n τ ( 0 = t 0 < t 1 < < t n = T ) .
Now, we discretize our problem (30) at ( x i , t n ) , then we get
0 C D t ; [ ζ ( t ) , 1 ] α u ( x i , t n ) = D u x x ( x i , t n ) A u x ( x i , t n ) + f ( x i , t n ) .
In Equation (31), for fixed t n and 1 i M 1 , the first- and second-order spatial derivatives are discretized by using the following central difference approximations:
u ( x i , t n ) x = u n i + 1 u n i 1 2 h + O ( h 2 ) ,
2 u ( x i , t n ) x 2 = u n i + 1 2 u n i + u n i 1 h 2 + O ( h 2 ) .
With the help of Equation (18), we get an approximation of the generalized Caputo-type fractional derivative term in (31) as follows:
H D t ; [ ζ ( t ) , 1 ] α u ( x i , t n ) = λ 0 u ( x i , t n ) + λ 1 u ( x i , t n 1 ) + λ 2 u ( x i , t n 2 ) + l = 3 n 3 λ n l u ( x i , t l ) + λ n 2 u ( x i , t 2 ) + λ n 1 u ( x i , t 1 ) + λ n u ( x i , t 0 ) + O ( τ 4 α ) .
where λ l is defined in Equation (19). Next, we use Equations (32)–(34) in discretized Equation (31), yielding
l = 0 n λ n l u ( x i , t l ) = D u x x ( x i , t n ) A u x ( x i , t n ) + f n i + e n i .
where | e n i | c ˜ ( τ 4 α + h 2 ) for some constant c ˜ .
Now, we use a numerical approximation u n i of u ( x i , t n ) to neglect the truncation error term in Equation (35), then we determine the following finite difference scheme:
λ 0 u n i + λ 1 u n 1 i + l = 3 n 2 λ n l u l i + λ n 2 u 2 i + λ n 1 u 1 i + λ n u 0 i = D u n i + 1 2 u n i + u n i 1 h 2 A u n i + 1 u n i 1 2 h + f n i ,
where 1 n N , 1 i M 1 . That is,
( D h 2 + A 2 h ) u 1 i 1 ( λ 0 + 2 D h 2 ) u 1 i + ( D h 2 A 2 h ) u 1 i + 1 = λ 1 u 0 i f 1 i , n = 1 , ( D h 2 + A 2 h ) u 2 i 1 ( λ 0 + 2 D h 2 ) u 2 i + ( D h 2 A 2 h ) u 2 i + 1 = λ 1 u 1 i + λ 2 u 0 i f 2 i , n = 2 , ( D h 2 + A 2 h ) u n i 1 ( λ 0 + 2 D h 2 ) u n i + ( D h 2 A 2 h ) u n i + 1 = λ 1 u n 1 i + l = 0 n 2 λ n l u l i f n i , n 3 . u 0 i = 0 , 0 i M , u n 0 = ϕ 1 ( t n ) , u n M = ϕ 2 ( t n ) , 0 n N .
We rewrite the matrix form of the above equation as follows:
K U 1 = λ 1 U 0 F 1 + H 1 , n = 1 , K U 2 = λ 1 U 1 + λ 2 U 0 F 2 + H 2 , n = 2 , K U n = λ 1 U n 1 + l = 0 n 2 λ n l U l F n + H n , n 3 ,
where the matrix as well as vectors of Equation (38) are defined as follows:
  • K = t r i D h 2 + A 2 h , λ 0 2 D h 2 , D h 2 A 2 h ( M 1 ) × ( M 1 ) ,
  • U n = ( u n 1 , u n 2 , , u n M 1 ) T ,
  • F n = ( f n 1 , f n 2 , , f n M 1 ) T ,
  • H n = ( D h 2 A 2 h ) u n 0 , 0 , , 0 , ( D h 2 + A 2 h ) u n M T , 1 n N .
Remark 4.
Since the coefficient matrix A is a tridiagonal and strictly diagonally dominant then d e t ( A ) 0 (Levy–Desplanques theorem). Therefore, at each time level t n , the proposed scheme (37) has a unique solution for 1 n N .
Theorem 2.
The local truncation error of difference scheme (36) at ( x i , t n ) , 1 i M 1 , 1 n N , is
| e n i | C ( h 2 + τ 4 α ) ,
where 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
e n i = λ 0 u ( x i , t n ) + λ 1 u ( x i , t n 1 ) + λ 2 u ( x i , t n 2 ) + l = 3 n 3 λ n l u ( x i , t l ) + λ n 2 u ( x i , t 2 ) + λ n 1 u ( x i , t 1 ) + λ n u ( x i , t 0 ) D u n i + 1 2 u n i + u n i 1 h 2 + A u n i + 1 u n i 1 2 h f n i = [ λ 0 u ( x i , t n ) + λ 1 u ( x i , t n 1 ) + λ 2 u ( x i , t n 2 ) + l = 2 n 3 λ n l u ( x i , t l ) + λ n 1 u ( x i , t 1 ) + λ n u ( x i , t 0 ) H D t ; [ ζ ( t ) , 1 ] α u ( x i , t n ) ] D u n i + 1 2 u n i + u n i 1 h 2 2 u ( x i , t n ) x 2 + A u n i + 1 u n i 1 2 h u ( x i , t n ) x = O ( τ 4 α ) D O ( h 2 ) + A O ( h 2 ) = O ( τ 4 α + h 2 ) .
In the next theorem, we use L 2 ( Ω ) -space with the norm . 2 and the inner product . , . .
Theorem 3.
(see [22]) If the tridiagonal matrix elements satisfy the inequality
λ 1 ( λ 0 ) M 3 λ 0 + D h 2 A 2 h λ 0 + D h 2 + A 2 h ,
then the finite difference scheme is stable.
Proof. 
Equation (38) can be rewritten as follows, for 1 n N
K U n = λ 1 U n 1 + l = 0 n 2 λ n l U l F n + H n ,
let,          V n = l = 0 n 2 λ n l U l F n + H n .
Since matrix K is invertible then the above equation can be rewritten as
U n = λ 1 K 1 U n 1 + K 1 V n .
Using the recurrence relation, we get the following equation
U n = ( λ 1 K 1 ) 2 U n 2 + ( λ 1 K 1 ) K 1 V n 1 + K 1 V n = ( λ 1 K 1 ) n U 0 + ( λ 1 K 1 ) n 1 K 1 V 1 + ( λ 1 K 1 ) n 2 K 1 V 2 + + K 1 V n .
Let U ˜ n be the approximate solution of Equation (38), then we define the error at grid point t n = n τ
e n = U n U ˜ n , 0 n N .
Then, we obtain
e n = ( λ 1 K 1 ) n e 0 , 1 n N .
From the definition of the compatible matrix norm
e n ( λ 1 K 1 ) n e 0 .
Since
( λ 1 K 1 ) n ( λ 1 K 1 ) ( λ 1 K 1 ) n 1 ( λ 1 K 1 ) n .
According to the Ostrowski therorem ([40], Theorem 3.1), let K = [ a i , j ] M 1 × M 1 , then the determinant of K satisfies
| d e t ( K ) | i = 1 M 1 | a i , i | j = 1 , j i M 1 | a i , j | = ( λ 0 ) M 3 λ 0 + D h 2 A 2 h λ 0 + D h 2 + A 2 h .
Using assumed inequality (40) into (42), we get | d e t ( K 1 ) | 1 λ 1 ;
this implies that
| λ 1 | K 1 1 .
Therefore,           | e n | | e 0 | .
Thus, the numerical scheme is stable. □
Theorem 4.
The solution U n i of the difference scheme (36) satisfies
max i , n | U ( x i , t n ) U ˜ ( x i , t n ) | C ( h 2 + τ 4 α )
for some constant C .
Proof. 
The truncation error for the difference scheme at ( x i , t n ) [ 0 , a ] × [ 0 , T ] is
| e n i | C ( h 2 + τ 4 α )
by Theorem 2. Now, using the theorem of stability, it implies
max i , n | U ( x i , t n ) U ˜ ( x i , t n ) | C | U ( x i , t 0 ) U ˜ ( x i , t 0 ) |
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 E and error E 2 corresponding to the L 2 -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.
E ( M , N ) = max 1 i M 1 | U N i u N i | ,
E 2 ( M , N ) = h i = 1 M 1 | U N i u N i | 2 1 / 2 ,
where { U n i } is the exact solution of advection diffusion equation and { u n i } is the approximate solution at the point ( x i , t n ) .
Moreover, the convergence order in the space and time directions for the described difference scheme corresponding to L -norm can be evaluated using the following formulas. R x is the order of convergence on the space side and R t on the temporal side.
R x = l o g ( E ( 2 M , N ) ) l o g ( E ( M , N ) ) l o g ( 2 ) ,
and
R t = l o g ( E ( M , 2 N ) ) l o g ( E ( M , N ) ) l o g ( 2 ) .
Example 1.
Ref. [39] Take function u ( t ) = t 4 + α , t [ 0 , 1 ] , for 0 < α < 1 . Determine the α-th order GCFD for u ( t ) at T = 1 numerically.
The maximum absolute error and rate of convergence for the scheme (18) to approximate the GCFD of function u ( t ) for α = 0.2 , 0.5 , 0.8 with uniform time steps 1 / 10 , 1 / 20 , 1 / 40 , 1 / 80 , 1 / 160 are calculated and shown in the following tables.
Table 1 shows the errors with respect to L -norm and convergence rates in time, and these data are found after calculating the classical Caputo derivative (i.e., taking ζ ( t ) = t , ω ( t ) = 1 ) of function u ( t ) 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 ( 4 α ) , while accuracy for the time derivative in [39] is ( 3 α ) . In Table 2, to compute the maximum errors E and order of convergence in the time direction, we take ω ( t ) = e t while the scale function is fixed with t. From Table 3, we validate the convergence for ζ ( t ) = t , ω ( t ) = t + 1 and the CPU time in seconds for α = 0.8 are discussed. Table 4 shows maximum errors and order of convergence for different choices of weight ω ( t ) = t 0.5 , t , t 4 , e 2 t ; scale is ζ ( t ) = t and α = 1 / 3 . In all cases, accuracy in time is obtained as ( 4 α ) 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 α = 0.5 and α = 0.8 , respectively.
Example 2.
Ref. [21] We take the following generalized fractional advection–diffusion equation:
0 C D t ; [ ζ ( t ) , ω ( t ) ] α u ( x , t ) = 2 u ( x , t ) x 2 u ( x , t ) x + f ( x , t ) , ( x , t ) ( 0 , 1 ) × ( 0 , 1 ) , u ( x , 0 ) = 0 , x ( 0 , 1 ) , u ( 0 , t ) = t 6 + α , u ( 1 , t ) = e t 6 + α , t ( 0 , 1 ] ,
where f ( x , t ) = e x t 6 Γ ( 7 + α ) 720 . When ζ ( t ) = t and ω ( t ) = 1 , then u ( x , t ) = e x t 6 + α is the exact solution.
To solve this example, we use the numerical scheme defined in (37). The maximum errors at time t = 1 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 h = 1 2000 and describe the numerical errors and convergence rates in time R t for different values of N. In Table 6, we fix τ = 1 500 and present the numerical errors and spatial convergence rates R x for different values of M. It is shown that our scheme (37) gives ( 4 α ) -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 ζ ( t ) 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 M = 2000 for particular choices of scale ζ ( t ) = t and weight function ω ( t ) = 1 .
Example 3.
Take the following fractional advection–diffusion equation:
0 C D t , [ ζ ( t ) , ω ( t ) ] α u ( x , t ) = 2 u ( x , t ) x 2 u ( x , t ) x + f ( x , t ) , ( x , t ) ( 0 , 1 ) × ( 0 , 1 ) , u ( x , 0 ) = 0 , x ( 0 , 1 ) , u ( 0 , t ) = 0 , u ( 1 , t ) = t 7 ( s i n 1 ) , t ( 0 , 1 ] ,
where f ( x , t ) = Γ 8 Γ ( 8 α ) t 7 α s i n ( x ) + t 7 ( s i n ( x ) + c o s ( x ) ) . When ζ ( t ) = t and ω ( t ) = 1 , the exact solution is u ( x , t ) = t 7 s i n ( x ) .
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 h = 1 / 12000 and display maximum absolute errors and errors E 2 with respect to norm L 2 , and also the rate of convergence for the temporal direction with different time steps N = 8 , 16 , 32 , 64 , 128 , 256 . In Table 9, we express maximum errors and convergence order in the space direction to set time steps fixed with τ = 1 / 500 and taking different M = 8 , 16 , 32 , 64 , 128 .
In Table 10, we validate the proposed scheme with [21] (Example 5.1). Firstly, we present the max error, the convergence order with fixed h = 1 6000 and varying N with 8, 16, 32, 64, 128 for α = 0.368 . 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 τ = 1 / 200 and changing M = 4 , 8 , 16 , 32 , 64 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 ζ ( t ) = t and ω ( t ) = 1 . Furthermore, we can compute numerical results for different suitable choices of scale and weight functions.

5. Conclusions

A high-order numerical scheme based on cubic interpolation formula is discussed for approximation of GCFD of α -th order. Properties of discretized coefficients are analyzed and local truncation error in approximation of GCFD is also discussed. Further, we establish a difference scheme for the generalized fractional advection–diffusion equation using the developed approximation formula for GCFD. The stability and convergence of this established scheme to approximate the time fractional generalized advection–diffusion equation are studied. Order of accuracy for the difference scheme is O ( τ 4 α + h 2 ) , with step sizes τ in time and h in space directions. The convergence order of the difference scheme is described by some numerical experiments.
Numerical calculations reveal that the proposed difference scheme has ( 4 α ) -th order convergence in time; and, in the space direction, it has second-order convergence. The temporal rate of convergence of the scheme for the generalized fractional advection–diffusion equation has the highest accuracy to date. In addition, the developed scheme can be directly used to get other Caputo-type time fractional advection–diffusion equations by selecting suitable scale ζ ( t ) and weight ω ( t ) functions in the generalized Caputo-type fractional derivative. The developed scheme is tested for the cases having smooth solutions of the considered fractional advection–diffusion equation. However, the case of the nonsmooth solutions will be presented in our future works.

Author Contributions

Conceptualization, S.K., R.K.P. and R.P.A.; methodology, S.K. and R.K.P.; software, S.K. and R.K.P.; validation, S.K. and R.K.P.; writing—original draft preparation, S.K., R.K.P. and R.P.A.; writing—review and editing, R.K.P. and R.P.A.; supervision, R.K.P.; funding acquisition, R.K.P. and R.P.A. All authors have equally contributed in preparation of the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Not applicable.

Acknowledgments

The authors sincerely thank the reviewers for their comments incorporated in the manuscript. The authors are thankful to K. Mustapha, Department of Mathematics and Statistics, King Fahd University of Petroleum and Minerals, Saudi Arabia for his valuable comments incorporated in the manuscript.

Conflicts of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References

  1. Agrawal, O.P. Some generalized fractional calculus operators and their applications in integral equations. Fract. Calc. Appl. Anal. 2012, 15, 700–711. [Google Scholar] [CrossRef]
  2. Benson, D.A.; Wheatcraft, S.W.; Meerschaert, M.M. Application of a fractional advection-dispersion equation. Water Resour. Res. 2000, 36, 1403–1412. [Google Scholar] [CrossRef] [Green Version]
  3. Zhou, L.; Selim, H. Application of the fractional advection-dispersion equation in porous media. Soil Sci. Soc. Am. J. 2003, 67, 1079–1084. [Google Scholar] [CrossRef]
  4. Baleanu, D.; Wu, G.C.; Zeng, S.D. Chaos analysis and asymptotic stability of generalized Caputo fractional differential equations. Chaos Solitons Fractals 2017, 102, 99–105. [Google Scholar] [CrossRef]
  5. Iqbal, M.S.; Yasin, M.W.; Ahmed, N.; Akgül, A.; Rafiq, M.; Raza, A. Numerical simulations of nonlinear stochastic Newell-Whitehead-Segel equation and its measurable properties. J. Comput. Appl. Math. 2023, 418, 114618. [Google Scholar] [CrossRef]
  6. Partohaghighi, M.; Mirtalebi, Z.; Akgül, A.; Riaz, M.B. Fractal–fractional Klein–Gordon equation: A numerical study. Results Phys. 2022, 42, 105970. [Google Scholar] [CrossRef]
  7. Dan, D.; Mueller, C.; Chen, K.; Glazier, J.A. Solving the advection-diffusion equations in biological contexts using the cellular Potts model. Phys. Rev. E 2005, 72, 041909. [Google Scholar] [CrossRef] [Green Version]
  8. Verwer, J.; Blom, J.; Hundsdorfer, W. An implicit-explicit approach for atmospheric transport-chemistry problems. Appl. Numer. Math. 1996, 20, 191–209. [Google Scholar] [CrossRef] [Green Version]
  9. Dehghan, M. Weighted finite difference techniques for the one-dimensional advection–diffusion equation. Appl. Math. Comput. 2004, 147, 307–319. [Google Scholar] [CrossRef]
  10. Mohebbi, A.; Dehghan, M. High-order compact solution of the one-dimensional heat and advection–diffusion equations. Appl. Math. Model. 2010, 34, 3071–3084. [Google Scholar] [CrossRef]
  11. Owolabi, K.M. High-dimensional spatial patterns in fractional reaction-diffusion system arising in biology. Chaos Solitons Fractals 2020, 134, 109723. [Google Scholar] [CrossRef]
  12. Podlubny, I. Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications; Elsevier: Amsterdam, The Netherlands, 1998. [Google Scholar]
  13. Anatoly, A.K. Hadamard-type fractional calculus. J. Korean Math. Soc. 2001, 38, 1191–1204. [Google Scholar]
  14. Gaboury, S.; Tremblay, R.; Fugère, B.J. Some relations involving a generalized fractional derivative operator. J. Inequalities Appl. 2013, 2013, 1–9. [Google Scholar] [CrossRef] [Green Version]
  15. Atangana, A.; Baleanu, D. New fractional derivatives with nonlocal and non-singular kernel: Theory and application to heat transfer model. arXiv 2016, arXiv:1602.03408. [Google Scholar] [CrossRef] [Green Version]
  16. Mustapha, K. An L1 approximation for a fractional reaction-diffusion equation, a second-order error analysis over time-graded meshes. SIAM J. Numer. Anal. 2020, 58, 1319–1338. [Google Scholar] [CrossRef]
  17. Alikhanov, A.A. A new difference scheme for the time fractional diffusion equation. J. Comput. Phys. 2015, 280, 424–438. [Google Scholar] [CrossRef] [Green Version]
  18. Abu Arqub, O. Numerical solutions for the Robin time-fractional partial differential equations of heat and fluid flows based on the reproducing kernel algorithm. Int. J. Numer. Methods Heat Fluid Flow 2018, 28, 828–856. [Google Scholar] [CrossRef]
  19. Li, Z.; Yan, Y. Error estimates of high-order numerical methods for solving time fractional partial differential equations. Fract. Calc. Appl. Anal. 2018, 21, 746–774. [Google Scholar] [CrossRef] [Green Version]
  20. Lv, C.; Xu, C. Error analysis of a high order method for time-fractional diffusion equations. SIAM J. Sci. Comput. 2016, 38, A2699–A2724. [Google Scholar] [CrossRef]
  21. Cao, J.; Li, C.; Chen, Y. High-order approximation to Caputo derivatives and Caputo-type advection-diffusion Equations (II). Fract. Calc. Appl. Anal. 2015, 18, 735–761. [Google Scholar] [CrossRef]
  22. Xu, Y.; Agrawal, O.P. Numerical solutions and analysis of diffusion for new generalized fractional Burgers equation. Fract. Calc. Appl. Anal. 2013, 16, 709–736. [Google Scholar] [CrossRef]
  23. Kumar, K.; Pandey, R.K.; Sultana, F. Numerical schemes with convergence for generalized fractional integro-differential equations. J. Comput. Appl. Math. 2021, 388, 113318. [Google Scholar] [CrossRef]
  24. Diethelm, K. An efficient parallel algorithm for the numerical solution of fractional differential equations. Fract. Calc. Appl. Anal. 2011, 14, 475–490. [Google Scholar] [CrossRef]
  25. Zheng, Y.; Li, C.; Zhao, Z. A note on the finite element method for the space-fractional advection diffusion equation. Comput. Math. Appl. 2010, 59, 1718–1726. [Google Scholar] [CrossRef] [Green Version]
  26. Mardani, A.; Hooshmandasl, M.R.; Heydari, M.H.; Cattani, C. A meshless method for solving the time fractional advection–diffusion equation with variable coefficients. Comput. Math. Appl. 2018, 75, 122–133. [Google Scholar] [CrossRef]
  27. Li, C.; Cai, M. High-order approximation to Caputo derivatives and Caputo-type advection–diffusion equations: Revisited. Numer. Funct. Anal. Optim. 2017, 38, 861–890. [Google Scholar] [CrossRef]
  28. Yadav, S.; Pandey, R.K.; Shukla, A.K.; Kumar, K. High-order approximation for generalized fractional derivative and its application. Int. J. Numer. Methods Heat Fluid Flow 2019, 29, 3515–3534. [Google Scholar] [CrossRef]
  29. Tian, W.; Deng, W.; Wu, Y. Polynomial spectral collocation method for space fractional advection–diffusion equation. Numer. Methods Partial Differ. Equations 2014, 30, 514–535. [Google Scholar] [CrossRef] [Green Version]
  30. Zhuang, P.; Liu, F.; Anh, V.; Turner, I. Numerical methods for the variable-order fractional advection-diffusion equation with a nonlinear source term. SIAM J. Numer. Anal. 2009, 47, 1760–1781. [Google Scholar] [CrossRef] [Green Version]
  31. Singh, D.; Sultana, F.; Pandey, R.K. Approximation of Caputo-Prabhakar derivative with application in solving time fractional Advection-Diffusion equation. Int. J. Numer. Methods Fluids 2022, 94, 896–919. [Google Scholar] [CrossRef]
  32. Kannan, R.; Wang, Z.J. LDG2: A variant of the LDG flux formulation for the spectral volume method. J. Sci. Comput. 2011, 46, 314–328. [Google Scholar] [CrossRef]
  33. Kannan, R.; Wang, Z.J. A study of viscous flux formulations for a p-multigrid spectral volume Navier Stokes solver. J. Sci. Comput. 2009, 41, 165–199. [Google Scholar] [CrossRef]
  34. Kumar, K.; Pandey, R.K.; Sharma, S.; Xu, Y. Numerical scheme with convergence for a generalized time-fractional Telegraph-type equation. Numer. Methods Partial Differ. Equ. 2019, 35, 1164–1183. [Google Scholar] [CrossRef]
  35. Cao, W.; Xu, Y.; Zheng, Z. Finite difference/collocation method for a generalized time-fractional KDV equation. Appl. Sci. 2018, 8, 42. [Google Scholar] [CrossRef] [Green Version]
  36. Xu, Y.; He, Z.; Agrawal, O.P. Numerical and analytical solutions of new generalized fractional diffusion equation. Comput. Math. Appl. 2013, 66, 2019–2029. [Google Scholar] [CrossRef]
  37. Ding, Q.; Wong, P.J. A higher order numerical scheme for generalized fractional diffusion equations. Int. J. Numer. Methods Fluids 2020, 92, 1866–1889. [Google Scholar] [CrossRef]
  38. Sultana, F.; Pandey, R.K.; Singh, D.; Agrawal, O.P. High order approximation on non-uniform meshes for generalized time-fractional telegraph equation. MethodsX 2022, 9, 101905. [Google Scholar] [CrossRef] [PubMed]
  39. Gao, G.; Sun, Z.; Zhang, H. A new fractional numerical differentiation formula to approximate the Caputo fractional derivative and its applications. J. Comput. Phys. 2014, 259, 33–50. [Google Scholar] [CrossRef]
  40. Xu, Y.; He, Z.; Xu, Q. Numerical solutions of fractional advection–diffusion equations with a kind of new generalized fractional derivative. Int. J. Comput. Math. 2014, 91, 588–600. [Google Scholar] [CrossRef]
Figure 1. Error plot of the numerical results for Example 1 at final time T = 1 for different values of α (left side (a) for α = 0.5 ; right side (b) for α = 0.8 ).
Figure 1. Error plot of the numerical results for Example 1 at final time T = 1 for different values of α (left side (a) for α = 0.5 ; right side (b) for α = 0.8 ).
Mathematics 11 01200 g001
Figure 2. Exact and approximate solutions of Example 2 for different α ’s with M = N = 200 , ω ( t ) = 1 , ζ ( t ) = t .
Figure 2. Exact and approximate solutions of Example 2 for different α ’s with M = N = 200 , ω ( t ) = 1 , ζ ( t ) = t .
Mathematics 11 01200 g002
Figure 3. Numerical solutions of Example 2 with different choices of weight functions ω ( t ) and scale functions ζ ( t ) , and M = N = 200 with α = 0.9 .
Figure 3. Numerical solutions of Example 2 with different choices of weight functions ω ( t ) and scale functions ζ ( t ) , and M = N = 200 with α = 0.9 .
Mathematics 11 01200 g003
Table 1. E errors and convergence rates R t for Example 1, when ζ ( t ) = t , ω ( t ) = 1 , with different α s.
Table 1. E errors and convergence rates R t for Example 1, when ζ ( t ) = t , ω ( t ) = 1 , with different α s.
α = 0.2 α = 0.5 α = 0.5  [39]
N E Error R t E Error R t E Error R t
101.6978  × 10 4 1.5401  × 10 3 1.3507  × 10 2
201.3130  × 10 5 3.69281.4383  × 10 4 3.42062.6121  × 10 3 2.3704
409.9792  × 10 7 3.71781.3116  × 10 5 3.45504.8618  × 10 4 2.4256
807.4966  × 10 8 3.73461.1811  × 10 6 3.47318.8645  × 10 5 2.4554
1605.5944  × 10 9 3.74421.0560  × 10 7 3.48351.5975  × 10 5 2.4722
Table 2. E errors and convergence rates R t for Example 1, when ζ ( t ) = t , ω ( t ) = e t , with different α s.
Table 2. E errors and convergence rates R t for Example 1, when ζ ( t ) = t , ω ( t ) = e t , with different α s.
   α = 0.2    α = 0.5    α = 0.8    α = 0.8
N E Error R t E Error R t E Error R t CPU Time (s)
107.5552  × 10 4 6.3075  × 10 3 3.2184  × 10 2 11.354861
207.1075  × 10 5 3.41016.7873  × 10 4 3.21624.2148  × 10 3 2.932833.306962
405.9370  × 10 6 3.58156.6677  × 10 5 3.34765.0375  × 10 4 3.064782.772072
804.6934  × 10 7 3.66106.2491  × 10 6 3.41555.7497  × 10 5 3.1312180.064565
1603.5650  × 10 8 3.71865.7027  × 10 7 3.45396.4090  × 10 6 3.1653419.978865
Table 3. E errors and convergence rates R t for Example 1, when ζ ( t ) = t , ω ( t ) = t + 1 , with different α s.
Table 3. E errors and convergence rates R t for Example 1, when ζ ( t ) = t , ω ( t ) = t + 1 , with different α s.
   α = 0.2    α = 0.5    α = 0.8    α = 0.8
N E Error R t E Error R t E Error R t CPU Time (s)
103.5019  × 10 4 3.1752  × 10 3 1.7251  × 10 2 10.194270
203.0621  × 10 5 3.51553.1426  × 10 4 3.33682.0904  × 10 3 3.044825.931157
402.4466  × 10 6 3.64572.9499  × 10 5 3.41322.3982  × 10 4 3.123865.689186
801.8844  × 10 7 3.69862.6976  × 10 6 3.45092.6799  × 10 5 3.1617155.623905
1601.4248  × 10 8 3.72522.4326  × 10 7 3.47112.9560  × 10 6 3.1805384.953845
Table 4. E errors and convergence rates R t for Example 1, when ζ ( t ) = t , α = 1 3 , with different weight functions.
Table 4. E errors and convergence rates R t for Example 1, when ζ ( t ) = t , α = 1 3 , with different weight functions.
   ω ( t ) = t 0.5     ω ( t ) = t     ω ( t ) = t 4     ω ( t ) = e 2 t
N E Error R t E Error R t E Error R t R t
109.5450  × 10 4 1.7465  × 10 3 3.2557  × 10 2
208.4374  × 10 5 3.49991.5027  × 10 4 3.53882.0348  × 10 3 4.00003.2332
407.0877  × 10 6 3.57341.2845  × 10 5 3.54831.2749  × 10 4 3.99643.4225
805.8166  × 10 7 3.60711.0650  × 10 6 3.59231.1169  × 10 5 3.51293.5220
1604.7123  × 10 8 3.62578.6813  × 10 8 3.61679.3941  × 10 7 3.57163.5764
Table 5. Errors E and E 2 with convergence rates in time for Example 2, when h = 1 / 2000 with different α s.
Table 5. Errors E and E 2 with convergence rates in time for Example 2, when h = 1 / 2000 with different α s.
α N E Error R t E 2 Error R t
0.8 8 1.6385  × 10 2 7.4195  × 10 3
16 2.2920  × 10 3 2.8377 1.3514  × 10 3 2.4569
32 2.8191  × 10 4 3.0233 1.8556  × 10 4 2.8645
64 3.2609  × 10 5 3.1119 2.2573  × 10 5 3.0392
128 3.6574  × 10 6 3.1564 2.5937  × 10 6 3.1215
256 4.0176  × 10 7 3.1864 1.2919  × 10 7 3.1721
0.5 8 3.7940  × 10 3 1.8168  × 10 3
16 4.2889  × 10 4 3.1451 2.5900  × 10 4 2.8104
32 4.3064  × 10 5 3.3160 2.8647  × 10 5 3.1765
64 4.0768  × 10 6 3.4010 2.8348  × 10 6 3.3371
128 3.7173  × 10 7 3.4551 2.6407  × 10 7 3.4243
256 3.0470  × 10 8 3.6088 9.8124  × 10 8 3.5948
0.2 8 5.4498  × 10 4 2.7293  × 10 4
16 5.1093  × 10 5 3.4150 3.1455  × 10 5 3.1172
32 4.2949  × 10 6 3.5724 2.8822  × 10 6 3.4480
64 3.3828  × 10 7 3.6663 2.3627  × 10 7 3.6087
128 2.2628  × 10 8 3.9020 1.6173  × 10 8 3.8687
256 1.8462  × 10 9 3.6155 5.7547  × 10 9 3.6574
Table 6. The maximum errors and convergence rates R x for Example 2, when τ = 1 / 500 with different α s.
Table 6. The maximum errors and convergence rates R x for Example 2, when τ = 1 / 500 with different α s.
α = 0.2 α = 0.5 α = 0.8
M E Error R x E Error R x E Error R x
82.3437  × 10 4 2.1275  × 10 4 1.8078  × 10 4
165.8814  × 10 5 1.99455.3334  × 10 5 1.99604.5242  × 10 5 1.9985
321.4733  × 10 5 1.99711.3373  × 10 5 1.99581.1319  × 10 5 1.9990
643.6832  × 10 6 2.00013.3408  × 10 6 2.00102.7940  × 10 6 2.0183
1289.2088  × 10 7 1.99998.3276  × 10 7 2.00426.6258  × 10 7 2.0762
Table 7. The maximum errors and convergence rates R t of Example 4.1 in [39] (page no 43) for different values of α with h = 1 / 2000 .
Table 7. The maximum errors and convergence rates R t of Example 4.1 in [39] (page no 43) for different values of α with h = 1 / 2000 .
Current SchemeL1−2 [39]
α N E R t E R t
0.9102.2536  × 10 3 1.8600  × 10 2
203.2727  × 10 4 2.78374.7228  × 10 3 1.9776
404.1906  × 10 5 2.96531.1488  × 10 3 2.0395
805.1072  × 10 6 3.03662.7367  × 10 4 2.0697
1606.1092  × 10 7 3.06356.4520  × 10 5 2.0846
0.5102.5470  × 10 4 2.4936  × 10 3
202.6005  × 10 5 3.29194.8366  × 10 4 2.3662
402.4615  × 10 6 3.40129.0163  × 10 5 2.4234
802.2827  × 10 7 3.43071.6457  × 10 5 2.4539
1602.3662  × 10 8 3.27012.9700  × 10 6 2.4701
Table 8. Errors E and E 2 with convergence rates in time for Example 3, when h = 1 / 12000 with different α s.
Table 8. Errors E and E 2 with convergence rates in time for Example 3, when h = 1 / 12000 with different α s.
α N E Error R t E 2 Error R t
0.8 8 5.2117  × 10 3 4.1672  × 10 1
16 7.4158  × 10 4 2.8159 5.9162  × 10 2 2.8163
32 9.1750  × 10 5 3.0148 7.2113  × 10 3 3.0364
64 1.0655  × 10 5 3.1062 8.3111  × 10 4 3.1171
128 1.1991  × 10 6 3.1515 9.3172  × 10 5 3.1571
256 1.3228  × 10 7 3.1803 1.0259  × 10 5 3.1830
0.5 8 1.5442  × 10 3 1.2679  × 10 1
16 1.7157  × 10 4 3.1700 1.3692  × 10 2 3.2111
32 1.7543  × 10 5 3.2898 1.3792  × 10 3 3.3113
64 1.6788  × 10 6 3.3854 1.3099  × 10 4 3.3964
128 1.5678  × 10 7 3.4207 1.2186  × 10 5 3.4262
256 1.4942  × 10 8 3.3913 1.1596  × 10 6 3.3935
0.2 8 2.8743  × 10 4 2.3613  × 10 2
16 2.5557  × 10 5 3.4914 2.0399  × 10 3 3.5330
32 2.2185  × 10 6 3.5261 1.7445  × 10 4 3.5476
64 1.7950  × 10 7 3.6275 1.4008  × 10 5 3.6385
128 1.5668  × 10 8 3.5181 1.2186  × 10 6 3.5230
256 2.2811  × 10 9 2.7800 1.7738  × 10 7 2.7803
Table 9. The maximum errors and convergence rates R x for Example 3, when τ = 1 / 500 with different α s.
Table 9. The maximum errors and convergence rates R x for Example 3, when τ = 1 / 500 with different α s.
α = 0.2 α = 0.5 α = 0.8
M  E Error R x   E Error R x   E Error R x
83.0308  × 10 4 2.7311  × 10 4 2.3102  × 10 4
167.6029  × 10 5 1.99516.8533  × 10 5 1.99465.8003  × 10 5 1.9938
321.9050  × 10 5 1.99681.7180  × 10 5 1.99611.4560  × 10 5 1.9941
644.7625  × 10 6 2.00004.2962  × 10 6 1.99963.6518  × 10 6 1.9953
1281.1909  × 10 6 1.99971.0752  × 10 6 1.99859.2444  × 10 7 1.9820
Table 10. The maximum errors and convergence rates R t for Example 5.1 in [21], when h = 1 / 6000 , and spatial convergence rates R x , when τ = 1 / 200 .
Table 10. The maximum errors and convergence rates R t for Example 5.1 in [21], when h = 1 / 6000 , and spatial convergence rates R x , when τ = 1 / 200 .
α = 0.368 Cao et al. [21] α = 0.368 Cao et al. [21]
N  E Error R t R t M  E Error R x R x
88.1539  × 10 4 48.8034  × 10 4
167.7720  × 10 5 3.39113.403182.2559  × 10 4 1.96431.9643
326.8845  × 10 6 3.49693.4972165.6583  × 10 5 1.99531.9953
645.8739  × 10 7 3.55103.5528321.4173  × 10 5 1.99721.9972
1284.9904  × 10 8 3.55713.5836643.5359  × 10 6 2.00302.0030
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Kumari, S.; Pandey, R.K.; Agarwal, R.P. High-Order Approximation to Generalized Caputo Derivatives and Generalized Fractional Advection–Diffusion Equations. Mathematics 2023, 11, 1200. https://doi.org/10.3390/math11051200

AMA Style

Kumari S, Pandey RK, Agarwal RP. High-Order Approximation to Generalized Caputo Derivatives and Generalized Fractional Advection–Diffusion Equations. Mathematics. 2023; 11(5):1200. https://doi.org/10.3390/math11051200

Chicago/Turabian Style

Kumari, Sarita, Rajesh K. Pandey, and Ravi P. Agarwal. 2023. "High-Order Approximation to Generalized Caputo Derivatives and Generalized Fractional Advection–Diffusion Equations" Mathematics 11, no. 5: 1200. https://doi.org/10.3390/math11051200

APA Style

Kumari, S., Pandey, R. K., & Agarwal, R. P. (2023). High-Order Approximation to Generalized Caputo Derivatives and Generalized Fractional Advection–Diffusion Equations. Mathematics, 11(5), 1200. https://doi.org/10.3390/math11051200

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop