# Numerical Approaches of the Generalized Time-Fractional Burgers’ Equation with Time-Variable Coefficients

Authors: Dumitru Vieru; Constantin Fetecau; Nehad Ali Shah; Jae Dong Chung
Journal: Journal of Function Spaces
 (2021)

---

## Abstract

The generalized time-fractional, one-dimensional, nonlinear Burgers equation with time-variable coefficients is numerically investigated. The classical Burgers equation is generalized by considering the generalized Atangana-Baleanu time-fractional derivative. The studied model contains as particular cases the Burgers equation with Atangana-Baleanu, Caputo-Fabrizio, and Caputo time-fractional derivatives. A numerical scheme, based on the finite-difference approximations and some integral representations of the two-parameter Mittag-Leffler functions, has been developed. Numerical solutions of a particular problem with initial and boundary values are determined by employing the proposed method. The numerical results are plotted to compare solutions corresponding to the problems with time-fractional derivatives with different kernels.


---

## Body

## 1. Introduction

The nonlinear convective–diffusive partial differential equations describe various mathematical models in important fields such as heat and mass transfer, fluid mechanics, and engineering. During the last years, various solution methods of ordinary differential equations and partial differential equations have been elaborated.Burgers’ equation is one of the most important equations involving both nonlinear propagation effects and diffusive effects. A particular form of Burgers’ equation describes the nonlinear wave propagation (the inviscid Burgers’ equation).Burgers’ equation is a suitable tool for analysis in various fields such as turbulent flows, gas dynamics, shock wave theory, nonlinear wave propagation, longitudinal elastic waves in isotropic solids, sedimentation of polydispersive suspensions and colloids, growth of molecular interfaces, traffic flow, and cosmology [1].By studying Burgers’ equation with random initial conditions or random forcing, Bec and Khanin [2] explained Burgers’ turbulence. The study of random Lagrangian systems, stochastic partial differential equations, the applications of field theory to the understanding of dissipative anomalies, and of multiscaling in hydrodynamic turbulence are some fields that have significantly benefited from the progress in Burgers’ turbulence. Yu [3] analytically studied the stability and density waves for traffic flow using the perturbation method and shown that the triangular shock waves, soliton wave, and kink wave appear for the density waves.In the last years, researchers have proved that many phenomena in engineering, bioengineering, physics, and chemistry can be successfully described by mathematical models that use mathematical tools from fractional calculus, i.e., the theory of derivatives and integrals of noninteger order.Models of viscoelastic materials, Caputo and Mainardi [4]; the signal processing, Marks and Hall [5]; diffusion problems, Olmstead and Handelsman [6]; viscoplastic materials modeling, Diethelm and Freed [7]; mechanical systems subject to damping, Gaul et al. [8]; relaxation and reaction kinetics of polymers, Glockle and Nonnenmacher [9]; and heat conduction, Hristov [10, 11] are some of the important problems modeled with the help of fractional differential operators.A very useful collection of numerical algorithms for Caputo-type derivatives, Riemann-Liouville integral operator, and Mittag-Leffler functions is that of Diethelm et al. [12].In the literature, there are articles in which solutions of Burgers’ equation with different time-fractional derivatives have been determined. We recall some of them. The effects of fractional-order of derivatives on the wave solutions of the generalized Zakharov-Kuznetsov-Burgers equation have been investigated by Faraz et al. [13]. The analytical approximate wave solutions are obtained using the homotopy analysis method and the time-fractional Caputo’s derivative, while exact solutions are determined with the help of the first integral method and the fractional derivative in Jumarie’s modified Riemann–Liouville sense. Bira et al. [14] studied a nonlinear time-fractional system of Boussinesq-Burgers equations. Using Lie group analysis, the authors derived the infinitesimal groups of transformations, the system of optimal algebras for the symmetry group of transformations, and the similarity variables that reduce the system of fractional partial differential equations to a system of fractional ordinary differential equations. The exact solutions and the physical significance of the solutions are obtained under the invariance condition. Saad et al. [15] have extended the model of the Burgers equation to generalized models based on Liouville-Caputo, Caputo-Fabrizio, and Mittag-Leffler time-fractional derivatives. Using the homotopy analysis transform method, the authors obtained approximate solutions of the newly proposed models. Baleanu and Shiri [16] numerically solved a system of fractional differential equations involving nonsingular Mittag-Leffler kernel using the collocation methods on discontinuous piecewise polynomial space. The existence and regularity of solutions and convergence of the introduced methods are derived.Recently, Vieru et al. [17] have generalized the time-fractional Atangana-Baleanu derivative. The newly proposed definition contains as particular cases the time-fractional Caputo, Caputo-Fabrizio, and Atangana-Baleanu derivatives.It is important to note that the operators of fractional derivatives are nonlocal in time and therefore have the advantage of modeling phenomena with memory. Caputo fractional derivatives are nonlocal operators but their kernel is singular. This weakness could have a negative effect when modeling real-world problems.The fractional derivative operators with the Mittag-Leffler kernel have all the benefits of Caputo operators; in addition, the kernel is nonsingular. Also, their fractional integral operators are the fractional average of the Riemann–Liouville fractional integral of the given function and the function itself. Caputo derivative was conceived for a description of linear short time elastic responses of deformed solids. It was consequently applied to the field of linear viscoelasticity where the Riemann-Liouville derivative was already applied to describe viscoelastic effects. It is known that the asymptotic behaviors of derivative operators with Mittag-Leffler kernel match the power-law behavior. The new fractional operators based on Mittag-Leffler functions have stronger and complex memories allowing capturing behaviors combining simultaneously (crossover) classical diffusion and anomalous behavior. Therefore, to model more complex and nonlinear phenomena, the new operators could be useful tools [11].In this paper, a nonlinear, one-dimensional, generalized Burgers equation with time-variable coefficients is numerically studied. The generalization consists of considering the fractional differential Burgers’ equation with the generalized time-fractional Atangana-Baleanu fractional derivative with Mittag-Leffler kernel.A numerical scheme, based on the finite-difference approximations and some integral representations of the two-parameter Mittag-Leffler functions, has been developed along with the consistency, stability, and convergence of the proposed method.It is important to point out that the studied generalized model can be customized to generate solutions to the problems described by the time-fractional Atangana-Baleanu, Caputo-Fabrizio, and Caputo fractional derivatives.Numerical solutions of a particular problem with initial and boundary values are determined by employing the proposed method. The numerical results are plotted to compare solutions corresponding to the problems with time-fractional derivatives with different kernels.


## 2. Preliminary Mathematics

In this section, we present the basic mathematical elements regarding the two-parametric Mittag-Leffler functions and the generalized time-fractional Atangana-Baleanu derivatives. These mathematical notions are necessary for the next sections of this paper.


### 2.1. One-Parametric and Two-Parametric Mittag-Leffler Functions

The classical one-parametric Mittag-Leffler function is defined as [18, 19]
(1)Eαz=∑j=0∞zjΓαj+1,α>0,where Γς=∫0∞e−ττς−1dτ,Reς>0 is the Euler integral of the second kind.The two-parametric Mittag-Leffler function generalizes the functionEαz and is defined by
(2)Eα,βz=∑j=0∞zjΓαj+β,α>0,β∈ℂ.It is easy to notice that function (1) is a particular case of function (2), so we have
(3)Eαz=Eα,1z.Let us recall some properties of Mittag-Leffler functions.(4)Eα,βz=zEα,α+βz+1/Γβ,Eα,βz=βEα,1+βz+αzddzEα,1+βz,(5)tγE1,1+γat=tγ∑j=0∞atjΓj+1+γ=Etγ,a−theMiller−Rossfunction,(6)tγE1+γ,1+γat1+γ=tγ∑j=0∞ajt1+γjΓ1+γ1+j=Rtγ,a−theRobotnovfunction,(7)∫0ztβ−1Eα,βbtαdt=zβEα,β+1bzα,β>0,∫0ztβ−1z−tγ−1Eα,βbtαdt=zβ+γ−1Eα,β+γbzαΓα,β>0,γ>0,∫0zz−tβ−1ebtdt=zβE1,β+1bzΓβ,β>0.The following special form of the one-parametric Mittag-Leffler function [12]:
(8)Gt−τ=Eα−t−τγα,α∈0,1,γ>0,τ∈0,t,along with its derivative:
(9)Mt−τ=∂Gt−τ∂τ=−1t−τEα,0−t−τγα,α∈0,1,γ>0,τ∈0,t,has applications in the theory of fractional-order viscoelasticity and in some problems described by fractional differential equations with constant coefficients.Some numerical algorithms for determining numerical values of the Mittag-Leffler functions have been presented in the reference [12]. These algorithms are based on the integral representations of the Mittag-Leffler functions. We will use in this paper the following integral representations:Ifα∈0,1,β∈ℝ,0<ρ<z,argz>απ,z≠0, then
(10)Eα,βz=∫ρ∞Kα,β,x,zdx+∫−απαπPα,β,ρ,y,zdy,Kα,β,x,z=1παx1−βαexp−x1/αxsinπ1−β−zsinπ1−β+αx2−2xzcosπα+z2,Pα,β,ρ,y,z=ρ1+1−β/α2παexpρ1/αcosy/αexpiφρexpiy−z,φ=ρ1/αsiny/α+y1+1−β/α.The integral representation(11)∫0∞e−sttmαEαm±btαdt=m!sα−1sα∓bm+1,Res>0,Reα>0,m∈ℕ,along with the definition of the Laplace transform of a function ϑt, Lϑt=∫0∞ϑtexp−stdt give the following relationship:
(12)LtmαEαm±btα=m!sα−1sα∓bm+1,Res>0,Reα>0,m∈ℕ.In the particular casem=0, Equation (12) becomes
(13)LEα±btα=sα−1sα∓b.


### 2.2. Generalized Atangana-Baleanu Time-Fractional Derivative

The function(14)φt,α,β=11−αEβ−α1−αtβ,t≥0,α∈0,1,β>0,is called the generalized Atangana-Baleanu kernel.The Laplace transform of the kernel (14) is given by
(15)Lφt,α,β=sβ−11−αsβ+α.Using the Laplace transform, the following properties of the generalized Atangana-Baleanu kernel (14) are found:
(16)limα⟶0Lφt,α,β=Llimα⟶0φt,α,β=sβ−1sβ=L1,limα⟶1Lφt,α,β=Llimα⟶1φt,α,β=1s1−β=Lt−βΓ1−β=Lφ0t,β,limβ⟶0Lφt,α,β=Llimβ⟶0φt,α,β=1s=L1,limβ⟶1Lφt,α,β=Llimβ⟶1φt,α,β=11−αs+α=L11−αe−αt/1−α=Lφ1t,α,therefore,
(17)limα⟶0φt,α,β=limβ⟶0φt,α,β=1,limα⟶1φt,α,β=φ0t,β=t−αΓ1−α,limβ⟶1φt,α,β=φ1t,α=11−αexp−αt1−α,φt,α,α=φ2t,α=11−αEα−α1−αtα,limα⟶1β⟶1φt,α,β=δt.In the above relations, functionsφ0t,α,φ1t,α,φ2t,α, and δtare, respectively, Caputo kernel, Caputo-Fabrizio kernel, Atangana-Baleanu kernel, and the Dirac’s distribution.Definition 1.
The generalized Atangana-Baleanu fractional derivative in Caputo sense.
Iff∈H10,T,T>0,α∈0,1,β∈0,1, the generalized Atangana-Baleanu fractional derivative in Caputo sense, of order α of the function ft, is defined by the relation
(18)DGABtα,βft=φt,α,β∗f.t=∫0tφt−τ,α,βf.τdτ.Using Equations (17) and (18), we obtain the following properties of the generalized Atangana-Baleanu time-fractional derivative:
(19)DGABt0,βft=DGABtα,0ft=1∗f.t=∫0tf.τdτ=ft−f0,DGABt1,1ft=δt∗f.t=f.t=dftdt,DGABt1,βft=φ0t,β∗f.t=DCtβft,DGABtα,1ft=φ1t,β∗f.t=DCFtβft,DGABtα,αft=φ2t,β∗f.t=DABtαft,whereDCtβft denotes the time-fractional Caputo derivative, DCFtβft is time-fractional Caputo-Fabrizio derivative, and DABtαft denotes the time-fractional Atangana-Baleanu derivative.Associated with the generalized Atangana-Baleanu derivative, we define the following fractional integral operator:(20)Jtα,βft=1−αft+αψ0t,β∗ft,α∈0,1,β∈0,1,where the kernelψ0t,β is defined as
(21)ψ0t,β=tβ−1Γβ.It is observed thatLψ0t,β=1/sβ,limβ⟶0Lψ0t,β=1=Lδt; therefore,
(22)limβ⇒0ψ0t,β=δt.Using property (22), the fractional integral operator can be defined for β=0.The fractional integral operator (20) has the following properties:
(23)Jt1,0ft=δt∗ft=ft,Jt1,1ft=1∗ft=∫0tfτdτ.Regarding the generalized Atangana-Baleanu derivative and associated fractional integral operator, we remember the foloowingproposition.Proposition 2.
The following relationships are fulfilled:(24)DGABtα,βJtα,βft=ft−1−αf0φt,α,β,Jtα,βDGABtα,βft=ft−f0.The demonstration of the above proposition can be found in the reference [17].The generalized fractional integral operator (20) contains the following particular cases:
(25)10.α=1,β∈0,1,Jt1,βft=ψ0t,β∗ft=1Γβ∫0tt−τβ−1fτdτ,i.e., the well-known Riemann-Liouville fractional integral operator.(26)20.α∈0,1,β=1,Jtα,1ft=1−αft+α∫0tfτdτ,that is, the integral operator associated to the Caputo-Fabrizio derivative.
(27)30.α=β∈0,1,Jtα,αft=1−αft+αΓα∫0tt−τα−1fτdτ,that is, the fractional integral operator associated with the Atangana-Baleanu fractional derivative.


## 2.1. One-Parametric and Two-Parametric Mittag-Leffler Functions

The classical one-parametric Mittag-Leffler function is defined as [18, 19]
(1)Eαz=∑j=0∞zjΓαj+1,α>0,where Γς=∫0∞e−ττς−1dτ,Reς>0 is the Euler integral of the second kind.The two-parametric Mittag-Leffler function generalizes the functionEαz and is defined by
(2)Eα,βz=∑j=0∞zjΓαj+β,α>0,β∈ℂ.It is easy to notice that function (1) is a particular case of function (2), so we have
(3)Eαz=Eα,1z.Let us recall some properties of Mittag-Leffler functions.(4)Eα,βz=zEα,α+βz+1/Γβ,Eα,βz=βEα,1+βz+αzddzEα,1+βz,(5)tγE1,1+γat=tγ∑j=0∞atjΓj+1+γ=Etγ,a−theMiller−Rossfunction,(6)tγE1+γ,1+γat1+γ=tγ∑j=0∞ajt1+γjΓ1+γ1+j=Rtγ,a−theRobotnovfunction,(7)∫0ztβ−1Eα,βbtαdt=zβEα,β+1bzα,β>0,∫0ztβ−1z−tγ−1Eα,βbtαdt=zβ+γ−1Eα,β+γbzαΓα,β>0,γ>0,∫0zz−tβ−1ebtdt=zβE1,β+1bzΓβ,β>0.The following special form of the one-parametric Mittag-Leffler function [12]:
(8)Gt−τ=Eα−t−τγα,α∈0,1,γ>0,τ∈0,t,along with its derivative:
(9)Mt−τ=∂Gt−τ∂τ=−1t−τEα,0−t−τγα,α∈0,1,γ>0,τ∈0,t,has applications in the theory of fractional-order viscoelasticity and in some problems described by fractional differential equations with constant coefficients.Some numerical algorithms for determining numerical values of the Mittag-Leffler functions have been presented in the reference [12]. These algorithms are based on the integral representations of the Mittag-Leffler functions. We will use in this paper the following integral representations:Ifα∈0,1,β∈ℝ,0<ρ<z,argz>απ,z≠0, then
(10)Eα,βz=∫ρ∞Kα,β,x,zdx+∫−απαπPα,β,ρ,y,zdy,Kα,β,x,z=1παx1−βαexp−x1/αxsinπ1−β−zsinπ1−β+αx2−2xzcosπα+z2,Pα,β,ρ,y,z=ρ1+1−β/α2παexpρ1/αcosy/αexpiφρexpiy−z,φ=ρ1/αsiny/α+y1+1−β/α.The integral representation(11)∫0∞e−sttmαEαm±btαdt=m!sα−1sα∓bm+1,Res>0,Reα>0,m∈ℕ,along with the definition of the Laplace transform of a function ϑt, Lϑt=∫0∞ϑtexp−stdt give the following relationship:
(12)LtmαEαm±btα=m!sα−1sα∓bm+1,Res>0,Reα>0,m∈ℕ.In the particular casem=0, Equation (12) becomes
(13)LEα±btα=sα−1sα∓b.


## 2.2. Generalized Atangana-Baleanu Time-Fractional Derivative

The function(14)φt,α,β=11−αEβ−α1−αtβ,t≥0,α∈0,1,β>0,is called the generalized Atangana-Baleanu kernel.The Laplace transform of the kernel (14) is given by
(15)Lφt,α,β=sβ−11−αsβ+α.Using the Laplace transform, the following properties of the generalized Atangana-Baleanu kernel (14) are found:
(16)limα⟶0Lφt,α,β=Llimα⟶0φt,α,β=sβ−1sβ=L1,limα⟶1Lφt,α,β=Llimα⟶1φt,α,β=1s1−β=Lt−βΓ1−β=Lφ0t,β,limβ⟶0Lφt,α,β=Llimβ⟶0φt,α,β=1s=L1,limβ⟶1Lφt,α,β=Llimβ⟶1φt,α,β=11−αs+α=L11−αe−αt/1−α=Lφ1t,α,therefore,
(17)limα⟶0φt,α,β=limβ⟶0φt,α,β=1,limα⟶1φt,α,β=φ0t,β=t−αΓ1−α,limβ⟶1φt,α,β=φ1t,α=11−αexp−αt1−α,φt,α,α=φ2t,α=11−αEα−α1−αtα,limα⟶1β⟶1φt,α,β=δt.In the above relations, functionsφ0t,α,φ1t,α,φ2t,α, and δtare, respectively, Caputo kernel, Caputo-Fabrizio kernel, Atangana-Baleanu kernel, and the Dirac’s distribution.Definition 1.
The generalized Atangana-Baleanu fractional derivative in Caputo sense.
Iff∈H10,T,T>0,α∈0,1,β∈0,1, the generalized Atangana-Baleanu fractional derivative in Caputo sense, of order α of the function ft, is defined by the relation
(18)DGABtα,βft=φt,α,β∗f.t=∫0tφt−τ,α,βf.τdτ.Using Equations (17) and (18), we obtain the following properties of the generalized Atangana-Baleanu time-fractional derivative:
(19)DGABt0,βft=DGABtα,0ft=1∗f.t=∫0tf.τdτ=ft−f0,DGABt1,1ft=δt∗f.t=f.t=dftdt,DGABt1,βft=φ0t,β∗f.t=DCtβft,DGABtα,1ft=φ1t,β∗f.t=DCFtβft,DGABtα,αft=φ2t,β∗f.t=DABtαft,whereDCtβft denotes the time-fractional Caputo derivative, DCFtβft is time-fractional Caputo-Fabrizio derivative, and DABtαft denotes the time-fractional Atangana-Baleanu derivative.Associated with the generalized Atangana-Baleanu derivative, we define the following fractional integral operator:(20)Jtα,βft=1−αft+αψ0t,β∗ft,α∈0,1,β∈0,1,where the kernelψ0t,β is defined as
(21)ψ0t,β=tβ−1Γβ.It is observed thatLψ0t,β=1/sβ,limβ⟶0Lψ0t,β=1=Lδt; therefore,
(22)limβ⇒0ψ0t,β=δt.Using property (22), the fractional integral operator can be defined for β=0.The fractional integral operator (20) has the following properties:
(23)Jt1,0ft=δt∗ft=ft,Jt1,1ft=1∗ft=∫0tfτdτ.Regarding the generalized Atangana-Baleanu derivative and associated fractional integral operator, we remember the foloowingproposition.Proposition 2.
The following relationships are fulfilled:(24)DGABtα,βJtα,βft=ft−1−αf0φt,α,β,Jtα,βDGABtα,βft=ft−f0.The demonstration of the above proposition can be found in the reference [17].The generalized fractional integral operator (20) contains the following particular cases:
(25)10.α=1,β∈0,1,Jt1,βft=ψ0t,β∗ft=1Γβ∫0tt−τβ−1fτdτ,i.e., the well-known Riemann-Liouville fractional integral operator.(26)20.α∈0,1,β=1,Jtα,1ft=1−αft+α∫0tfτdτ,that is, the integral operator associated to the Caputo-Fabrizio derivative.
(27)30.α=β∈0,1,Jtα,αft=1−αft+αΓα∫0tt−τα−1fτdτ,that is, the fractional integral operator associated with the Atangana-Baleanu fractional derivative.


## 3. Problem Formulation

The classical one-dimensional Burgers equation with variable coefficients, defined fory,t∈0,L×0,T,L>0,T>0, is [20]
(28)∂uy,t∂t=at∂2uy,t∂y2+btuy,t∂uy,t∂y,whereat>0 and bt are differentiable and bounded functions of the variable t. For bt=const., Rizun and Engel’brekht [20] have determined the analytical solution of Equation (28).In the present paper, we consider a generalized form of Equation (28), namely,
(29)DGABtα,βuy,t=at∂2uy,t∂y2+btuy,t∂uy,t∂y,0<α,β<1,whereDGABtα,βuy,t is the generalized Atangana-Baleanu derivative defined by Equation (18). Along with Equation (29), we consider the initial and boundary conditions
(30)uy,0=ψy,y∈0,L,(31)u0,t=uL,t=0,t∈0,T.In the following, we will elaborate a numerical scheme for determining the solution of problem (29)–(31). The proposed scheme is based on the finite-difference method and on the properties of the Mittag-Leffler functions.Let us consider the discrete set of spatial nodesyk=kΔy,k=0,1,2,⋯,N1, respectively, the discrete set of the time nodes tj=jΔt,j=0,1,2,⋯,N2, where Δy=L/N1, Δt=T/N2 are the increment steps of y and t, respectively.


### 3.1. Numerical Evaluation of the Generalized Atangana-Baleanu Derivative

Using Equation (18), we have
(32)DGABtα,βuy,ty,t=yk,tn=11−α∫0tnEβ−α1−αtn−τβ∂uyk,τ∂τdτ=∑j=0n−111−α∫tjtj+1Eβ−α1−αtn−τβ∂uyk,τ∂τdτ,n=1,2,⋯,N2.The first-order time derivative∂uyk,τ/∂τ,τ∈tj,tj+1 is approximated by [21].
(33)∂uyk,τ∂τ≅uyk,tj+1−uyk,tjΔt−Δt2∂2uyk,τ∂τ2τ=tj+OΔt2.Using approximation (33), Equation (32) becomes
(34)DGABtα,βuy,ty,t=yk,tn=∑j=0n−1uyk,tj+1−uyk,tj1−αΔt∫tjtj+1Eβ−α1−αtn−τβdτ+Rkn=∑j=0n−1uyk,tj+1−uyk,tj1−αΔt∫tn−tj+1tn−tjEβ,1λxβdx+Rkn=∑j=0n−1uyk,tj+1−uyk,tj1−αΔt∫0tn−tjEβ,1λxβdx−∫0tn−tj+1Eβ,1λxβdx+Rkn,where λ=−α/1−α, and Rkn is the truncation error.Now, using the properties of Mittag-Leffler functions given in Equation (7), we obtain
(35)DGABtα,βuy,ty,t=yk,tn=∑j=0n−1uyk,tj+1−uyk,tj1−αΔttn−tjEβ,2λtn−tjβ−tn−tj+1Eβ,2λtn−tj+1β+Rkn=∑j=0n−1Anjukj+1−∑j=0n−1Bnjukj+Rkn,where
(36)Anj=tn−tjEβ,2λtn−tjβ1−αΔt,Bnj=tn−tj+1Eβ,2λtn−tj+1β1−αΔt,ukj=uyk,tj.Introducing notations(37)Dnj=n−1En1−nEn0,j=0,n−j+1Enj−1−2n−jEnjEnn−1,j=n,+n−j−1Enj+1,j∈1,2,⋯,n−1,where Enj=Eβ,2λtn−tjβ, Equation (35) can be written in the equivalent form
(38)DGABtα,βuy,ty,t=yk,tn=11−α∑j=0nDnjukj+Rkn.The truncation errorRkn is defined as
(39)Rkn=∑j=0n−111−α∫tjtj+1Δt2∂2uyk,τ∂τ2τ=tjEβ−α1−αtn−τβdτ=Δt221−α∑j=0n−1∂2uyk,τ∂τ2τ=tjn−jEnj−n−j−1Enj+1≤C1Δt221−αmax0≤t≤tn∂2uyk,τ∂τ2,where C1 is a constant coming from the bounded of the Mittag-Leffler function. The inequality (39) ensures the consistency of the proposed method because, assuming that the function uy,t is sufficiently smooth on its domain of definition, the truncation error Rkn tends to zero if the time step Δt tends to zero. So, at point yk,tn, the generalized Atangana-Baleanu time-fractional derivative is approximated by
(40)DGABtα,βuy,ty,t=yk,tn≅11−α∑j=0nDnjukj.In the above relations, the numerical values of the Mittag-Leffler functions are evaluated using the integral representation (10).An important property of the coefficientsDnj is given in the following.Lemma 3.
CoefficientsDnj given by (37) have negative values for j=0,1,⋯,n−1. For α∈0,1/2, Dnn>0 and Dnn⟶1 if Δt⟶0.
Proof. LetQz a function defined as
(41)Qz=zEβ,2λzΔtβ−z−1Eβ,2λz−1Δtβ.Using the definition of Mittag-Leffler function,Qz can be written as
(42)Qz=z∑k=0∞λkΓβk+2zΔtβk−z−1∑k=0∞λkΓβk+2z−1Δtβk=1Δt∑k=0∞λkΓβk+2zΔtβk+1−∑k=0∞λkΓβk+2z−1Δtβk+1.Using the formulaΓz+1=zΓz, the derivative of function Qzis given by
(43)dQzdz=∑k=0∞λkΓβk+1zΔtβk+1−∑k=0∞λkΓβk+1z−1Δtβk+1=Eβ,1λzΔtβ−Eβ,1λz−1Δtβ.It is known that Mittag-Leffler functionEβ,1z is an increasing function. Because, in this study, parameter λ<0, it results that function Qz is a decreasing function.The coefficientsDnj are written as
(44)Dnj=n−j+1Enj−1−n−jEnj−n−jEnj−n−j−1Enj+1=n−j+1Eβ,2λn−j+1Δtβ−n−jEβ,2λn−jΔtβ−n−jEβ,2λn−jΔtβ−n−j−1Eβ,2λn−j−1Δtβ=Qn−j+1−Qn−j.Since, the functionQz is decreasing we obtain that Dnj<0,j=1,2,⋯,n−1.Dnn=Enn−1=Eβ,2λΔtβ. Since α∈0,1/2, we obtain λ<1; therefore, λΔtβ⟶0 for Δtβ⟶0. Using the asymptotic expansion of Mittag-Leffler function [22], we have Eβ,2λΔtβ≃1+λΔtβ/Γβ+2 that proves the property in Lemma 3.In the following, two examples of the application of formula (40) are presented.Example 1.
The generalized Atangana-Baleanu derivative of functionft=t.
The time-fractional generalized Atangana-Baleanu derivative of functionft=t is given by
(45)DGABtα,βft=11−α∫0tEβ,1−α1−ατβdτ=t1−αEβ,2−α1−αtβ,t≥0.Table1 gives the values of this derivative determined with the analytical expression (45), respectively, with the numerical formula (40) for tk=0.02k,k=1,2,⋯,20,α=0.45,β=0.75.Table 1
Comparison between analytical and numerical results of (GAB) derivative of functionft=t.
tk=kΔtΔt=0.02Values given by Equation. (41)Values given by Equation (40)tk=kΔtΔt=0.02Values given by Equation (41)Values given by Equation (40)0.020.03540020.03540020.220.34214640.34214640.040.06953210.06953210.240.36952330.36952330.060.10267370.10267370.260.39642450.39642450.080.13495570.13495570.280.42287010.42287010.100.16646440.16646440.300.44887830.44887830.120.19526060.19526060.320.47446580.47446580.140.22740260.22740260.340.49964810.49964810.160.25693080.25693080.360.52443960.52443960.180.28587610.28587610.380.54885370.54885370.200.31427220.31427220.400.57290290.5729029Numerical results obtained by Equations (40) and (45) are graphically illustrated in Figure 1.Figure 1
Comparison between analytical and numerical (GAB)-fractional derivative of functionft=t.It can be seen in Table1 and Figure 1 that there is a very good accuracy of the numerical method given by equation (40).Example 2.
Find the solution of the fractional equation(46)DGABtα,βut=tγ,γ>0,α,β∈0,1.Using the Laplace transform, it is found that the analytical solution of Equation (46) is given by
(47)ut=1−αtλ+αΓ1+γΓ1+β+γtβ+λ.Table2 gives the values of the solution of Equation (46) determined with the analytical expression (47), respectively, with the numerical formula (40) for tk=0.01k,k=4,8,⋯,80,α=0.45,β=0.75,λ=2.Table 2
Comparison between analytical and numerical results of the solution of Equation (46).
tk=kΔtΔt=0.01Analytical solutionNumerical solutiontk=kΔtΔt=0.01Analytical solutionNumerical solution0.040.0912002×10−49.0985187×10−40.440.12776240.12776830.083.7158954×10−33.7171929×10−30.480.15375590.15376220.128.5174136×10−38.5192433×10−30.520.18241270.18241950.160.01539800.01540020.560.21378890.21379620.200.02443420.02443710.600.24793930.24794710.240.03569890.03570230.640.28491780.28492610.280.04926060.04926450.680.32477730.32478610.320.06518520.06518960.720.36757010.36757940.360.08353630.08354110.760.41334740.41335720.400.10437540.10438070.800.46216010.4621703Numerical results for the solution of Equation (46), obtained by Equations (40) and (47) are graphically illustrated in Figure 2.Figure 2
Comparison between analytical and numerical solutions of Equation (46).It can be seen in Table2 and Figure 2 that there is a very good accuracy of the numerical method given by Equation (40).


### 3.2. Particular Cases

Note that the expression (40) of the generalized Atangana-Baleanu time-fractional derivative can be easily customized in the following cases:
(a)
Ifα=β∈0,1, it is obtained the expression of the time-fractional Atangana-Baleanu derivative(b)
Ifα∈0,1 and β=1, it is obtained the expression of the time-fractional Caputo-Fabrizio derivative(c)
Ifα=1 and β∈0,1, it is obtained the expression of the time-fractional Caputo derivativeWe must note that in the first two cases, the formulas (34)–(40) which determine the numerical values of the fractional derivative remain valid, obviously with the corresponding particularizations of the parameters α and β.In the third case, forα=1, there is an indeterminacy because limα⟶1α/1−α=∞.To eliminate this indeterminacy, we use the following asymptotic expansion of the Mittag-Leffler function [22]:
(48)β∈0,1,z⟶∞,βπ≤argz≤π,Eβz=−∑r=1n01zrΓ1−βr+O1zn0+1,n0∈ℕ,n0>1.Using Equation (44) into (34), we obtain the following relation:
(49)limα⟶111−α∫tjtj+1Eβ−α1−αtn−τβdτ=∫tjtj+1limα⟶1−11−α∑r=1n01Γ1−βr1−αr−αtn−τβrdτ∫tjtj+1tn−τ−βΓ1−βdτ=−tn−τ1−β1−βΓ1−βtjtj+1=1Γ2−βtn−tj1−β−tn−tj+11−β.As expected, the integrand in (45) is the Caputo kernel φ0t,β. Replacing (45) into (34), we obtain the following approximate formula for the time-fractional Caputo derivative:
(50)DCtα,βuy,ty,t=yk,tn=∑j=0nCnjukj,where(51)Cnj=1ΔtβΓ2−βn−11−β−n1−β,j=0,1ΔtβΓ2−βn−j+11−β−2n−j1−β+n−j−11−β1ΔtβΓ2−β,j=n.,j=1,2,⋯,n−1,


### 3.3. Numerical Solution to the Burgers Equation with Generalized Atangana-Baleanu Time-Fractional Derivative

To determine the numerical solution of Equation (29), we approximate the first- and second-order derivative with respect to y by [23–25].
(52)∂uy,t∂yy,t=yk,tj≅uk+1j−uk−1j2Δy,k=1,2,⋯,N1−1,∂2uy,t∂y2y,t=yk,tj≅uk+1j−2ukj+uk−1jΔy2,k=1,2,⋯,N1−1.The nonlinear termuy,t∂uy,t/∂y is replaced by the equivalent term
(53)13uy,t∂uy,t∂y+∂u2y,t∂y.For simplicity of calculations and proving the stability of the scheme, the nonlinear term is approximated by [30, 31](54)uy,t∂uy,t∂yy,t=yk,tj=13ukj−1uk+1j−uk−1j2Δy+uk+1j−1uk+1j−uk−1j−1uk−1j2Δy.Using Equations (40), (52), and (54), we obtain the following numerical scheme for fractional Burgers equation (29):
(55)btn6Δyukn−1+uk−1n−1−atnΔy2uk−1n+Dnn1−α+2atnΔy2ukn−btn6Δyukn−1+uk+1n−1+atnΔy2uk+1n=−11−α∑j=0n−1Dnjukj.The initial and boundary conditions (30) and (31) are transformed in following discrete relationships:
(56)uk0=uyk,0=uyk,t0=ψyk=ψk,k=0,1,⋯,N1,(57)u0j=u0,tj=uy0,tj=0,j=1,2,⋯,N2,(58)uN1j=uL,tj=uyN1,tj=0,j=1,2,⋯,N2.Using (57) and (58), the numerical scheme (55) is written in the following metrical form:
(59)MnUn=Pn,n=1,2,⋯,N2,where the matrix Mn=mijni,j=1,2,⋯N1−1 is the following tridiagonal matrix:
(60)Mn=m11nm12n000⋯0m21nm22nm23n00⋯00m32nm33nm34n0⋯000m43nm44nm45n⋯0..⋯00⋯mN1−2N1−3nmN1−2N1−2nmN1−2N1−1n00⋯00mN1−1N1−2nmN1−1N1−1nN1−1,N1−1,with(61)mijn=−btn6Δyuin−1+ui+1n−1+atnΔy2δi+1,j+Dnn1−α+2atnΔy2δi,j+btn6Δyuin−1+ui−1n−1+atnΔy2δi,j+1,i,j=1,2,⋯,N1−1,Un=u1nu2nu3n⋯uN1−2nuN1−1nT,Pn=p1np2np3n⋯pN1−2npN1−1nT,where
(62)pkn=−11−α∑j=0n−1Dnjukj,k=1,2,⋯,N1−1.In the following, we assume thatα∈0,1/2,β0,1.Letux,t be the exact solution of Equation (29). The local truncation error of the numerical scheme (55) is
(63)rkn=Dnn1−α+2atnΔy2uyk,tn+btn6Δyuyk,tn−1+uyk−1,tn−1−atnΔy2uyk−1,tn−btn6Δyuyk,tn−1+uyk+1,tn−1+atnΔy2uyk+1,tn+11−α∑j=0n−1Dnjuyk,tj.According to the Taylor expansion, it is found that exists a constantk0 such that rkn≤k0Δt+Δy2; therefore, the discrete operator (with finite differences) converges towards the continuous operator (with derivatives) for Δt,Δy⟶0 (vanishing truncation error, so, the numerical scheme is consistent).Let us introduce the notations:(64)ukny=uk+1n−uknΔy,ukny1=ukn−uk−1nΔy,ukny2=uk+1n−uk−1n2Δy,uknyy1=uk+1n−2ukn+uk−1nΔy,Un=ukn,k=0,1,⋯,N1,Uny=ukny,k=0,1,⋯,N1,Un,Vn=Δy∑k=1N1−1uknvkn,Un2=Un,Un=∑k=1N1−1uknukn,n=1,2,⋯,N2,Un∞=lim1≤k≤N1−1ukn.The following boundedness theorem will be proved.Theorem 6.
Assume thatψ∈H020,L,α∈0,1/2,β∈0,1. There exists a constant K~ such that the numerical solution derived by the finite difference scheme (55) satisfies inequality Un≤K~,n=1,2,⋯,N2.Proof.
Equation (55) can be written in the equivalent form
(65)Dnn1−αukn+11−α∑j=1n−1Dnjukj+Dn01−αuk0−btn6Δyukn−1uk+1n−uk−1n+uk+1n−1uk+1n−uk−1n−1uk−1n=atnuk+1n−2ukn+uk−1nΔy2.
Multiplying Equation (65) by Δyukn, summing for k=1,2,⋯,N1−1, and using Equation (64), we obtain
(66)Dnn1−αUn2+11−α∑j=1n−1DnjUj,Un+Dn01−αU0,Un−Sn=atnΔy∑k=1N1−1uk+1n−2ukn+uk−1nΔy2ukn,
where(67)Sn=btn6∑k=1N1−1ukn−1uknuk+1n−uk−1n−1uk−1nukn+uknuk+1n−1uk+1n−uk−1nukn−1ukn.☐A straightforward calculus leads to(68)Sn=btn6uN1nuN1n−1uN1−1n+uN1nuN1−1n−1uN1−1n−u0nu1n−1u1n−u0nu0n−1u1n.Based on the boundary conditions (57) and (58), it results that Sn=0.Using the property(69)Δy∑k=1N1−1uk+1n−2ukn+uk−1nΔy2ukn=Unyy1,Un=−Uny,Uny=−Uny2from equality (66) and Lemma 3, we obtain
(70)Dnn1−αUn2+atnUny2=11−α∑j=1n−1−DnjUj,Un+−Dn01−αU0,Un≤121−α∑j=1n−1−DnjUj2+Un2+−Dn021−αU02+Un2=121−α∑j=1n−1−DnjUj2+−Dn021−αU02+121−α∑j=0n−1−DnjUnn=121−α∑j=1n−1−DnjUj2+−Dn021−αU02+Dnn21−αUnn.Using (70) and mathematical induction, we obtain that Uj2≤U02 for j=1,2,⋯,n−1. Therefore, we have
(71)Dnn21−αUn2+atnUny2≤121−α∑j=1n−1−DnjUj2+−Dn021−αU02≤121−α∑j=1n−1−DnjU02+−Dn021−αU02=121−α∑j=1n−1−DnjU02=Dnn21−αU02.Equation (71) implies that
(72)Dnn21−αUn2−U02≤−atnUny2≤0,therefore,(73)Un≤U0.But, there exist a constantK~ such that U0≤K~, respectively, Un≤K~. Theorem 6 shows that for a given function ψy=uy,0, the numerical solution of the problem is bounded as time increases.In order to prove a stability result, we consider a perturbed problem, i.e., the fractional equation (29) with a different initial condition ψ^y=uy,0.We say that the numerical method isglobally stable if there exists a constant K^ such that Un−U∧n≤K^ψy−ψ^y, where Un and U∧n are solutions corresponding to initial conditions ψyand ψ^y, respectively. Using the boundedness theorem one can prove that the numerical scheme (55) corresponding to the Equation (29) is globally stable.Theorem 7. (Convergence).
Assume thatukn,u¯kn are the solutions of the difference scheme (55) and (29) and ekn=ukn−u¯kn. There exists a positive constant C~ such that En∞≤C~Δt+Δy2.Proof.
Using the equalities (55) and (63), we obtain
(74)Dnn1−αekn−atneknyy1=rkn+11−α∑j=1n−1−Dnjekj+−Dn01−αek0+Tkn+Skn,where
(75)Tkn=btn3ukn−1ukny2−u¯kn−1u¯kny2,Skn=btn3ukn−1ukny2−u¯kn−1u¯kny2.
Multiplying Equation (74) by Δyekn and summing up from k=1 to N1−1, we get
(76)Dnn1−αEn2−atnEnyy1,En=Rn,En+11−α∑j=1n−1−DnjEj,En+−Dn01−αE0,En+Tn,En+Sn,En.☐Using Equation (69) and the Cauchy-Schwarz inequality ekjekn≤1/2ekj2+ekn2, we get
(77)Dnn1−αEn2+atnEny2≤121−α∑j=1n−1−DnjEj2+121−α∑j=0n−1−DnjEn2+−Dn021−αE02+Tn,En+Sn,En+Rn,En.We note thatE0=0, and ∑j=0n=1−Dnj=Dnn. Based on the Cauchy-Schwarz inequality, we obtain following inequalities:
(78)Rn,En≤γ1εnRn2+14γ1εnEn2,Tn,En=btnΔy3∑k=1N1−1ukn−1ukny2−u¯kn−1u¯kny2ekn=btnΔy3∑k=1N1−1ukn−1ekn+u¯kny2−u¯kn−1u¯kny2ekn=btnΔy3∑k=1N1−1ukn−1ekny2+ekn−1u¯kny2ekn≤btnΔyK~3∑k=1N1−1ekny2+ekn−1ekn≤btnK~312En−12+En2+εnEny2+14εnEn2,Sn,En≤btnK~3εnγ2+γ3Eny2+14γ2εnEn2+14γ3εnEn−12.Takingεn=3atn/4btnK~,γ1=1/3,γ2=γ3=1/2, we obtain the following inequality:
(79)δnEn2+atnEny2≤121−α∑j=1n−1−DnjEj2+πnEn−12+k1atnΔt+Δy22,where
(80)ρn=K~btn6atn6+2K~btn+atn,πn=K~btn18atn4K~btn+3atn,δn=Dnn21−α−ρn.Now, we will show that(81)δnEn2+atnEny2≤k1atnδnΔt+Δy22δn−πn−Dn0+Dnn/21−α,by mathematical induction. Forn=1, inequality (81) holds according to (79). Suppose that
(82)δnEj2+atnEjy2≤k1atnδnΔt+Δy22δn−πn−Dn0+Dnn/21−α,j=1,2,⋯n−1,then,
(83)Ej2≤k1atnΔt+Δy22δn−πn−Dn0+Dnn/21−α,j=1,2,⋯n−1.Introducing (83) into (79), we obtain
(84)δnEn2+atnEny2≤121−α∑j=1n−1−Dnjk1atnΔt+Δy22δn−πn−Dn0+Dnn/21−α+πnk1atnΔt+Δy22δn−πn−Dn0+Dnn/21−α+k1atnΔt+Δy22=Dnn+Dn0/21−α+πnδn−πn−Dn0+Dnn/21−α+1k1atnΔt+Δy22=k1atnδnΔt+Δy22δn−πn−Dn0+Dnn/21−αTherefore, it results in(85)En2≤k1atnΔt+Δy22δn−πn−Dn0+Dnn/21−α,Eny2≤k1δnΔt+Δy22δn−πn−Dn0+Dnn/21−α.Using the discrete Sobolev’s inequalityEn∞≤C1En+C2Eny, we obtain the result of Theorem 7.Some comparisons between our results and other results from the literature are presented at the end of this section.(a)
If we considerat=1,bt=−1,α=1,β∈0,1, then the problem given by (29)–(31) is identical to the problem studied by Li et al. [23], for p=1(b)
Ifat=V,bt=−1,α=β, our solutions are identical to that obtained by Yadav et al. [24] for gx,t=0(c)
Ifα=β, the numerical scheme (38) for the generalized Atangana-Baleanu derivative is identical with that given in [25], Equation (3.2).


## 3.1. Numerical Evaluation of the Generalized Atangana-Baleanu Derivative

Using Equation (18), we have
(32)DGABtα,βuy,ty,t=yk,tn=11−α∫0tnEβ−α1−αtn−τβ∂uyk,τ∂τdτ=∑j=0n−111−α∫tjtj+1Eβ−α1−αtn−τβ∂uyk,τ∂τdτ,n=1,2,⋯,N2.The first-order time derivative∂uyk,τ/∂τ,τ∈tj,tj+1 is approximated by [21].
(33)∂uyk,τ∂τ≅uyk,tj+1−uyk,tjΔt−Δt2∂2uyk,τ∂τ2τ=tj+OΔt2.Using approximation (33), Equation (32) becomes
(34)DGABtα,βuy,ty,t=yk,tn=∑j=0n−1uyk,tj+1−uyk,tj1−αΔt∫tjtj+1Eβ−α1−αtn−τβdτ+Rkn=∑j=0n−1uyk,tj+1−uyk,tj1−αΔt∫tn−tj+1tn−tjEβ,1λxβdx+Rkn=∑j=0n−1uyk,tj+1−uyk,tj1−αΔt∫0tn−tjEβ,1λxβdx−∫0tn−tj+1Eβ,1λxβdx+Rkn,where λ=−α/1−α, and Rkn is the truncation error.Now, using the properties of Mittag-Leffler functions given in Equation (7), we obtain
(35)DGABtα,βuy,ty,t=yk,tn=∑j=0n−1uyk,tj+1−uyk,tj1−αΔttn−tjEβ,2λtn−tjβ−tn−tj+1Eβ,2λtn−tj+1β+Rkn=∑j=0n−1Anjukj+1−∑j=0n−1Bnjukj+Rkn,where
(36)Anj=tn−tjEβ,2λtn−tjβ1−αΔt,Bnj=tn−tj+1Eβ,2λtn−tj+1β1−αΔt,ukj=uyk,tj.Introducing notations(37)Dnj=n−1En1−nEn0,j=0,n−j+1Enj−1−2n−jEnjEnn−1,j=n,+n−j−1Enj+1,j∈1,2,⋯,n−1,where Enj=Eβ,2λtn−tjβ, Equation (35) can be written in the equivalent form
(38)DGABtα,βuy,ty,t=yk,tn=11−α∑j=0nDnjukj+Rkn.The truncation errorRkn is defined as
(39)Rkn=∑j=0n−111−α∫tjtj+1Δt2∂2uyk,τ∂τ2τ=tjEβ−α1−αtn−τβdτ=Δt221−α∑j=0n−1∂2uyk,τ∂τ2τ=tjn−jEnj−n−j−1Enj+1≤C1Δt221−αmax0≤t≤tn∂2uyk,τ∂τ2,where C1 is a constant coming from the bounded of the Mittag-Leffler function. The inequality (39) ensures the consistency of the proposed method because, assuming that the function uy,t is sufficiently smooth on its domain of definition, the truncation error Rkn tends to zero if the time step Δt tends to zero. So, at point yk,tn, the generalized Atangana-Baleanu time-fractional derivative is approximated by
(40)DGABtα,βuy,ty,t=yk,tn≅11−α∑j=0nDnjukj.In the above relations, the numerical values of the Mittag-Leffler functions are evaluated using the integral representation (10).An important property of the coefficientsDnj is given in the following.Lemma 3.
CoefficientsDnj given by (37) have negative values for j=0,1,⋯,n−1. For α∈0,1/2, Dnn>0 and Dnn⟶1 if Δt⟶0.
Proof. LetQz a function defined as
(41)Qz=zEβ,2λzΔtβ−z−1Eβ,2λz−1Δtβ.Using the definition of Mittag-Leffler function,Qz can be written as
(42)Qz=z∑k=0∞λkΓβk+2zΔtβk−z−1∑k=0∞λkΓβk+2z−1Δtβk=1Δt∑k=0∞λkΓβk+2zΔtβk+1−∑k=0∞λkΓβk+2z−1Δtβk+1.Using the formulaΓz+1=zΓz, the derivative of function Qzis given by
(43)dQzdz=∑k=0∞λkΓβk+1zΔtβk+1−∑k=0∞λkΓβk+1z−1Δtβk+1=Eβ,1λzΔtβ−Eβ,1λz−1Δtβ.It is known that Mittag-Leffler functionEβ,1z is an increasing function. Because, in this study, parameter λ<0, it results that function Qz is a decreasing function.The coefficientsDnj are written as
(44)Dnj=n−j+1Enj−1−n−jEnj−n−jEnj−n−j−1Enj+1=n−j+1Eβ,2λn−j+1Δtβ−n−jEβ,2λn−jΔtβ−n−jEβ,2λn−jΔtβ−n−j−1Eβ,2λn−j−1Δtβ=Qn−j+1−Qn−j.Since, the functionQz is decreasing we obtain that Dnj<0,j=1,2,⋯,n−1.Dnn=Enn−1=Eβ,2λΔtβ. Since α∈0,1/2, we obtain λ<1; therefore, λΔtβ⟶0 for Δtβ⟶0. Using the asymptotic expansion of Mittag-Leffler function [22], we have Eβ,2λΔtβ≃1+λΔtβ/Γβ+2 that proves the property in Lemma 3.In the following, two examples of the application of formula (40) are presented.Example 1.
The generalized Atangana-Baleanu derivative of functionft=t.
The time-fractional generalized Atangana-Baleanu derivative of functionft=t is given by
(45)DGABtα,βft=11−α∫0tEβ,1−α1−ατβdτ=t1−αEβ,2−α1−αtβ,t≥0.Table1 gives the values of this derivative determined with the analytical expression (45), respectively, with the numerical formula (40) for tk=0.02k,k=1,2,⋯,20,α=0.45,β=0.75.Table 1
Comparison between analytical and numerical results of (GAB) derivative of functionft=t.
tk=kΔtΔt=0.02Values given by Equation. (41)Values given by Equation (40)tk=kΔtΔt=0.02Values given by Equation (41)Values given by Equation (40)0.020.03540020.03540020.220.34214640.34214640.040.06953210.06953210.240.36952330.36952330.060.10267370.10267370.260.39642450.39642450.080.13495570.13495570.280.42287010.42287010.100.16646440.16646440.300.44887830.44887830.120.19526060.19526060.320.47446580.47446580.140.22740260.22740260.340.49964810.49964810.160.25693080.25693080.360.52443960.52443960.180.28587610.28587610.380.54885370.54885370.200.31427220.31427220.400.57290290.5729029Numerical results obtained by Equations (40) and (45) are graphically illustrated in Figure 1.Figure 1
Comparison between analytical and numerical (GAB)-fractional derivative of functionft=t.It can be seen in Table1 and Figure 1 that there is a very good accuracy of the numerical method given by equation (40).Example 2.
Find the solution of the fractional equation(46)DGABtα,βut=tγ,γ>0,α,β∈0,1.Using the Laplace transform, it is found that the analytical solution of Equation (46) is given by
(47)ut=1−αtλ+αΓ1+γΓ1+β+γtβ+λ.Table2 gives the values of the solution of Equation (46) determined with the analytical expression (47), respectively, with the numerical formula (40) for tk=0.01k,k=4,8,⋯,80,α=0.45,β=0.75,λ=2.Table 2
Comparison between analytical and numerical results of the solution of Equation (46).
tk=kΔtΔt=0.01Analytical solutionNumerical solutiontk=kΔtΔt=0.01Analytical solutionNumerical solution0.040.0912002×10−49.0985187×10−40.440.12776240.12776830.083.7158954×10−33.7171929×10−30.480.15375590.15376220.128.5174136×10−38.5192433×10−30.520.18241270.18241950.160.01539800.01540020.560.21378890.21379620.200.02443420.02443710.600.24793930.24794710.240.03569890.03570230.640.28491780.28492610.280.04926060.04926450.680.32477730.32478610.320.06518520.06518960.720.36757010.36757940.360.08353630.08354110.760.41334740.41335720.400.10437540.10438070.800.46216010.4621703Numerical results for the solution of Equation (46), obtained by Equations (40) and (47) are graphically illustrated in Figure 2.Figure 2
Comparison between analytical and numerical solutions of Equation (46).It can be seen in Table2 and Figure 2 that there is a very good accuracy of the numerical method given by Equation (40).


## 3.2. Particular Cases

Note that the expression (40) of the generalized Atangana-Baleanu time-fractional derivative can be easily customized in the following cases:
(a)
Ifα=β∈0,1, it is obtained the expression of the time-fractional Atangana-Baleanu derivative(b)
Ifα∈0,1 and β=1, it is obtained the expression of the time-fractional Caputo-Fabrizio derivative(c)
Ifα=1 and β∈0,1, it is obtained the expression of the time-fractional Caputo derivativeWe must note that in the first two cases, the formulas (34)–(40) which determine the numerical values of the fractional derivative remain valid, obviously with the corresponding particularizations of the parameters α and β.In the third case, forα=1, there is an indeterminacy because limα⟶1α/1−α=∞.To eliminate this indeterminacy, we use the following asymptotic expansion of the Mittag-Leffler function [22]:
(48)β∈0,1,z⟶∞,βπ≤argz≤π,Eβz=−∑r=1n01zrΓ1−βr+O1zn0+1,n0∈ℕ,n0>1.Using Equation (44) into (34), we obtain the following relation:
(49)limα⟶111−α∫tjtj+1Eβ−α1−αtn−τβdτ=∫tjtj+1limα⟶1−11−α∑r=1n01Γ1−βr1−αr−αtn−τβrdτ∫tjtj+1tn−τ−βΓ1−βdτ=−tn−τ1−β1−βΓ1−βtjtj+1=1Γ2−βtn−tj1−β−tn−tj+11−β.As expected, the integrand in (45) is the Caputo kernel φ0t,β. Replacing (45) into (34), we obtain the following approximate formula for the time-fractional Caputo derivative:
(50)DCtα,βuy,ty,t=yk,tn=∑j=0nCnjukj,where(51)Cnj=1ΔtβΓ2−βn−11−β−n1−β,j=0,1ΔtβΓ2−βn−j+11−β−2n−j1−β+n−j−11−β1ΔtβΓ2−β,j=n.,j=1,2,⋯,n−1,


## 3.3. Numerical Solution to the Burgers Equation with Generalized Atangana-Baleanu Time-Fractional Derivative

To determine the numerical solution of Equation (29), we approximate the first- and second-order derivative with respect to y by [23–25].
(52)∂uy,t∂yy,t=yk,tj≅uk+1j−uk−1j2Δy,k=1,2,⋯,N1−1,∂2uy,t∂y2y,t=yk,tj≅uk+1j−2ukj+uk−1jΔy2,k=1,2,⋯,N1−1.The nonlinear termuy,t∂uy,t/∂y is replaced by the equivalent term
(53)13uy,t∂uy,t∂y+∂u2y,t∂y.For simplicity of calculations and proving the stability of the scheme, the nonlinear term is approximated by [30, 31](54)uy,t∂uy,t∂yy,t=yk,tj=13ukj−1uk+1j−uk−1j2Δy+uk+1j−1uk+1j−uk−1j−1uk−1j2Δy.Using Equations (40), (52), and (54), we obtain the following numerical scheme for fractional Burgers equation (29):
(55)btn6Δyukn−1+uk−1n−1−atnΔy2uk−1n+Dnn1−α+2atnΔy2ukn−btn6Δyukn−1+uk+1n−1+atnΔy2uk+1n=−11−α∑j=0n−1Dnjukj.The initial and boundary conditions (30) and (31) are transformed in following discrete relationships:
(56)uk0=uyk,0=uyk,t0=ψyk=ψk,k=0,1,⋯,N1,(57)u0j=u0,tj=uy0,tj=0,j=1,2,⋯,N2,(58)uN1j=uL,tj=uyN1,tj=0,j=1,2,⋯,N2.Using (57) and (58), the numerical scheme (55) is written in the following metrical form:
(59)MnUn=Pn,n=1,2,⋯,N2,where the matrix Mn=mijni,j=1,2,⋯N1−1 is the following tridiagonal matrix:
(60)Mn=m11nm12n000⋯0m21nm22nm23n00⋯00m32nm33nm34n0⋯000m43nm44nm45n⋯0..⋯00⋯mN1−2N1−3nmN1−2N1−2nmN1−2N1−1n00⋯00mN1−1N1−2nmN1−1N1−1nN1−1,N1−1,with(61)mijn=−btn6Δyuin−1+ui+1n−1+atnΔy2δi+1,j+Dnn1−α+2atnΔy2δi,j+btn6Δyuin−1+ui−1n−1+atnΔy2δi,j+1,i,j=1,2,⋯,N1−1,Un=u1nu2nu3n⋯uN1−2nuN1−1nT,Pn=p1np2np3n⋯pN1−2npN1−1nT,where
(62)pkn=−11−α∑j=0n−1Dnjukj,k=1,2,⋯,N1−1.In the following, we assume thatα∈0,1/2,β0,1.Letux,t be the exact solution of Equation (29). The local truncation error of the numerical scheme (55) is
(63)rkn=Dnn1−α+2atnΔy2uyk,tn+btn6Δyuyk,tn−1+uyk−1,tn−1−atnΔy2uyk−1,tn−btn6Δyuyk,tn−1+uyk+1,tn−1+atnΔy2uyk+1,tn+11−α∑j=0n−1Dnjuyk,tj.According to the Taylor expansion, it is found that exists a constantk0 such that rkn≤k0Δt+Δy2; therefore, the discrete operator (with finite differences) converges towards the continuous operator (with derivatives) for Δt,Δy⟶0 (vanishing truncation error, so, the numerical scheme is consistent).Let us introduce the notations:(64)ukny=uk+1n−uknΔy,ukny1=ukn−uk−1nΔy,ukny2=uk+1n−uk−1n2Δy,uknyy1=uk+1n−2ukn+uk−1nΔy,Un=ukn,k=0,1,⋯,N1,Uny=ukny,k=0,1,⋯,N1,Un,Vn=Δy∑k=1N1−1uknvkn,Un2=Un,Un=∑k=1N1−1uknukn,n=1,2,⋯,N2,Un∞=lim1≤k≤N1−1ukn.The following boundedness theorem will be proved.Theorem 6.
Assume thatψ∈H020,L,α∈0,1/2,β∈0,1. There exists a constant K~ such that the numerical solution derived by the finite difference scheme (55) satisfies inequality Un≤K~,n=1,2,⋯,N2.Proof.
Equation (55) can be written in the equivalent form
(65)Dnn1−αukn+11−α∑j=1n−1Dnjukj+Dn01−αuk0−btn6Δyukn−1uk+1n−uk−1n+uk+1n−1uk+1n−uk−1n−1uk−1n=atnuk+1n−2ukn+uk−1nΔy2.
Multiplying Equation (65) by Δyukn, summing for k=1,2,⋯,N1−1, and using Equation (64), we obtain
(66)Dnn1−αUn2+11−α∑j=1n−1DnjUj,Un+Dn01−αU0,Un−Sn=atnΔy∑k=1N1−1uk+1n−2ukn+uk−1nΔy2ukn,
where(67)Sn=btn6∑k=1N1−1ukn−1uknuk+1n−uk−1n−1uk−1nukn+uknuk+1n−1uk+1n−uk−1nukn−1ukn.☐A straightforward calculus leads to(68)Sn=btn6uN1nuN1n−1uN1−1n+uN1nuN1−1n−1uN1−1n−u0nu1n−1u1n−u0nu0n−1u1n.Based on the boundary conditions (57) and (58), it results that Sn=0.Using the property(69)Δy∑k=1N1−1uk+1n−2ukn+uk−1nΔy2ukn=Unyy1,Un=−Uny,Uny=−Uny2from equality (66) and Lemma 3, we obtain
(70)Dnn1−αUn2+atnUny2=11−α∑j=1n−1−DnjUj,Un+−Dn01−αU0,Un≤121−α∑j=1n−1−DnjUj2+Un2+−Dn021−αU02+Un2=121−α∑j=1n−1−DnjUj2+−Dn021−αU02+121−α∑j=0n−1−DnjUnn=121−α∑j=1n−1−DnjUj2+−Dn021−αU02+Dnn21−αUnn.Using (70) and mathematical induction, we obtain that Uj2≤U02 for j=1,2,⋯,n−1. Therefore, we have
(71)Dnn21−αUn2+atnUny2≤121−α∑j=1n−1−DnjUj2+−Dn021−αU02≤121−α∑j=1n−1−DnjU02+−Dn021−αU02=121−α∑j=1n−1−DnjU02=Dnn21−αU02.Equation (71) implies that
(72)Dnn21−αUn2−U02≤−atnUny2≤0,therefore,(73)Un≤U0.But, there exist a constantK~ such that U0≤K~, respectively, Un≤K~. Theorem 6 shows that for a given function ψy=uy,0, the numerical solution of the problem is bounded as time increases.In order to prove a stability result, we consider a perturbed problem, i.e., the fractional equation (29) with a different initial condition ψ^y=uy,0.We say that the numerical method isglobally stable if there exists a constant K^ such that Un−U∧n≤K^ψy−ψ^y, where Un and U∧n are solutions corresponding to initial conditions ψyand ψ^y, respectively. Using the boundedness theorem one can prove that the numerical scheme (55) corresponding to the Equation (29) is globally stable.Theorem 7. (Convergence).
Assume thatukn,u¯kn are the solutions of the difference scheme (55) and (29) and ekn=ukn−u¯kn. There exists a positive constant C~ such that En∞≤C~Δt+Δy2.Proof.
Using the equalities (55) and (63), we obtain
(74)Dnn1−αekn−atneknyy1=rkn+11−α∑j=1n−1−Dnjekj+−Dn01−αek0+Tkn+Skn,where
(75)Tkn=btn3ukn−1ukny2−u¯kn−1u¯kny2,Skn=btn3ukn−1ukny2−u¯kn−1u¯kny2.
Multiplying Equation (74) by Δyekn and summing up from k=1 to N1−1, we get
(76)Dnn1−αEn2−atnEnyy1,En=Rn,En+11−α∑j=1n−1−DnjEj,En+−Dn01−αE0,En+Tn,En+Sn,En.☐Using Equation (69) and the Cauchy-Schwarz inequality ekjekn≤1/2ekj2+ekn2, we get
(77)Dnn1−αEn2+atnEny2≤121−α∑j=1n−1−DnjEj2+121−α∑j=0n−1−DnjEn2+−Dn021−αE02+Tn,En+Sn,En+Rn,En.We note thatE0=0, and ∑j=0n=1−Dnj=Dnn. Based on the Cauchy-Schwarz inequality, we obtain following inequalities:
(78)Rn,En≤γ1εnRn2+14γ1εnEn2,Tn,En=btnΔy3∑k=1N1−1ukn−1ukny2−u¯kn−1u¯kny2ekn=btnΔy3∑k=1N1−1ukn−1ekn+u¯kny2−u¯kn−1u¯kny2ekn=btnΔy3∑k=1N1−1ukn−1ekny2+ekn−1u¯kny2ekn≤btnΔyK~3∑k=1N1−1ekny2+ekn−1ekn≤btnK~312En−12+En2+εnEny2+14εnEn2,Sn,En≤btnK~3εnγ2+γ3Eny2+14γ2εnEn2+14γ3εnEn−12.Takingεn=3atn/4btnK~,γ1=1/3,γ2=γ3=1/2, we obtain the following inequality:
(79)δnEn2+atnEny2≤121−α∑j=1n−1−DnjEj2+πnEn−12+k1atnΔt+Δy22,where
(80)ρn=K~btn6atn6+2K~btn+atn,πn=K~btn18atn4K~btn+3atn,δn=Dnn21−α−ρn.Now, we will show that(81)δnEn2+atnEny2≤k1atnδnΔt+Δy22δn−πn−Dn0+Dnn/21−α,by mathematical induction. Forn=1, inequality (81) holds according to (79). Suppose that
(82)δnEj2+atnEjy2≤k1atnδnΔt+Δy22δn−πn−Dn0+Dnn/21−α,j=1,2,⋯n−1,then,
(83)Ej2≤k1atnΔt+Δy22δn−πn−Dn0+Dnn/21−α,j=1,2,⋯n−1.Introducing (83) into (79), we obtain
(84)δnEn2+atnEny2≤121−α∑j=1n−1−Dnjk1atnΔt+Δy22δn−πn−Dn0+Dnn/21−α+πnk1atnΔt+Δy22δn−πn−Dn0+Dnn/21−α+k1atnΔt+Δy22=Dnn+Dn0/21−α+πnδn−πn−Dn0+Dnn/21−α+1k1atnΔt+Δy22=k1atnδnΔt+Δy22δn−πn−Dn0+Dnn/21−αTherefore, it results in(85)En2≤k1atnΔt+Δy22δn−πn−Dn0+Dnn/21−α,Eny2≤k1δnΔt+Δy22δn−πn−Dn0+Dnn/21−α.Using the discrete Sobolev’s inequalityEn∞≤C1En+C2Eny, we obtain the result of Theorem 7.Some comparisons between our results and other results from the literature are presented at the end of this section.(a)
If we considerat=1,bt=−1,α=1,β∈0,1, then the problem given by (29)–(31) is identical to the problem studied by Li et al. [23], for p=1(b)
Ifat=V,bt=−1,α=β, our solutions are identical to that obtained by Yadav et al. [24] for gx,t=0(c)
Ifα=β, the numerical scheme (38) for the generalized Atangana-Baleanu derivative is identical with that given in [25], Equation (3.2).


## 4. Example

In this section, we investigate the problem described by the fractional differential equation (29), with at=0.01+0.1sin2πt/4,bt=−10sinπt/4, along with the following initial and boundary conditions:
(86)uy,0=ψy=21−exp−y,y∈0,1,u0,t=u0.5,t=0.The discussions in this section are mainly focused on the comparison between the solutions corresponding to the different types of fractional derivatives with the non-singular kernels. Because the kernels of fractional derivatives play an essential role in describing the effects of memory, Figure3 shows the profiles of the kernels corresponding to the generalized Atangana-Baleanu, Atangana-Baleanu, and Caputo-Fabrizio derivatives, versus fractional parameter α, for t=0.3. It is observed that for α<0.8, the values of the generalized kernel (GAB) are higher than the (AB) ones but smaller than the (CF) ones; therefore, the different damping effects will be transferred on the values of the function uy,t.Figure 3
The profiles of the (GAB), (AB), and (CF) kernels of fractional derivatives.The numerical solutions of Burgers equation with generalized Atangana-Baleanu time-fractional derivatives (GAB) are plotted in Figure4 for different fractional parameters α and β, for y∈0,0.5, and for different values of the time t. The computational results are taken at the time step size Δt=10−2 and spatial step size Δy=10−2.Figure 4
Profiles of the solutionuy,t for Burgers equation with generalized Atangana-Baleanu time-fractional derivative.Figure4 shows that, for the considered cases, the function uy,t is decreasing in relation to time tand it is increasing in relation to the fractional parameters α. The curves in Figure 4 indicate that the numerical solutions are bounded, and the numerical method is stable.A comparison between the solutions of the Burgers equation with the generalized Atangana-Baleanu derivative, the Atangana-Baleanu derivative and the Caputo-Fabrizio derivative is presented in Tables3–5, for x=0.4, and for different values of time t.Table 3
Comparison between numerical solutions corresponding to different fractional derivatives with the nonsingular kernel.
Values of timet(GAB)-fractional derivative forα=0.1,β=0.5(AB)-fractional derivative forα=0.1,β=0.1(CF)-fractional derivative forα=0.1,β=1.00.100.0220000.0210000.0220000.200.0190000.0180000.0190000.300.0150000.0140000.0150000.400.0120000.0110000.0120000.500.0120000.0110000.0120000.600.0074550.0070550.0076120.700.0061070.0057800.0062350.800.0051060.0048330.0052130.900.0043550.0041220.0059291.000.0037840.0035820.004913Table 4
Comparison between numerical solutions corresponding to different fractional derivatives with the nonsingular kernel.
Values of timet(GAB)-fractional derivative forα=0.3,β=0.5(AB)-fractional derivative forα=0.3,β=0.3(CF)-fractional derivative forα=0.3,β=1.00.100.0260000.0250000.0290000.200.0220000.0210000.0240000.300.0180000.0170000.0410000.400.0140000.0130000.0150000.500.0140000.0130000.0150000.600.0089210.0082990.0095370.700.0073060.0067980.0078260.800.0061080.0056840.0065530.900.0052100.0048480.0055951.000.0045270.0042120.004865Table 5
Comparison between numerical solutions corresponding to different fractional derivatives with the nonsingular kernel.
Values of timet(GAB)-fractional derivative forα=0.4,β=0.5(AB)-fractional derivative forα=0.4,β=0.4(CF)-fractional derivative forα=0.4,β=1.00.100.0290000.0280000.0330000.200.0250000.0240000.0550000.300.0200000.0190000.0220000.400.0160000.0150000.0170000.500.0160000.0150000.0170000.600.0098860.0094330.0110000.700.0080960.0077260.0090240.800.0067680.0064600.0075550.900.0057730.0055090.0064511.000.0050160.0047870.005610As expected, due to the properties of the three nonsingular kernels, the values of the solution corresponding to the Atangana-Baleanu derivative are the lowest, respectively, and the values of the solution corresponding to the Caputo-Fabrizio derivative are the highest.


## 5. Conclusion

A numerical scheme to solve the fractional Burgers equation with variable coefficients has been developed.The proposed algorithm is a linear implicit finite difference scheme.Unconditional stability and convergence of the method are proved for the fractional parametersα∈0.1/2,β∈0,1.The proposed scheme is suitable to solve fractional Burgers’ equation with four fractional derivatives, namely, generalized Atangana-Baleanu, Atangana-Baleanu, Caputo-Fbrizio, and Caputo derivatives.



---

*Source: 8803182-2021-12-08.xml*