Analysis of a Laplace Spectral Method for Time-Fractional Advection-Diffusion Equations Incorporating the Atangana-Baleanu Derivative

1  Introduction

Fractional calculus (FC) is a historical discipline in mathematics that traces its origins to the works of Leibniz and Euler, who explored integrals and derivatives of non-integer orders [1]. Despite significant advancements, this topic continues to captivate researchers due to its rich mathematical and numerical aspects. Today, FC extends beyond pure mathematics and has found applications in various scientific disciplines [2]. Replacing standard operators with fractional operators has significantly enhanced the precision and accuracy of numerous physical models and systems [3–5]. Researchers have explored numerous fractional operators and assessed various definitions of fractional derivatives, particularly those governed by power-law kernels, commonly known as singular or local fractional derivatives. Notable examples include the Grünwald-Letnikov, Riemann-Liouville, Caputo, Riesz, and Hadamard fractional derivatives [6]. However, these fractional derivatives can efficiently model the non-local and dissipative properties of physical processes. The non-locality of a fractional derivative refers to the fact that the value of a fractional derivative at a given point depends on the function’s values over an entire interval rather than just at the point itself.

In 2015, Caputo and Fabrizio [7] introduced a new fractional derivative, termed the Caputo-Fabrizio (CF) derivative, which employs an exponential kernel to overcome the limitations of fractional derivatives with singular kernels, such as those found in the Riemann-Liouville and Liouville-Caputo derivatives. The singularity in these earlier definitions often complicates their application to physical systems. The CF derivative’s exponential kernel, being nonsingular, provides a more natural transition and eliminates the issues associated with the Riemann-Liouville and Liouville-Caputo derivatives. However, despite its innovation, the CF derivative faced criticisms for a notable drawback: its kernel lacks nonlocality. To address these shortcomings, Atangana and Baleanu [8] introduced a novel fractional derivative in 2016, known as Atangana-Baleanu (AB) derivative, which is based on the Mittag-Leffler function. Unlike the CF derivative the AB derivative incorporates a nonlocal and nonsingular kernel, effectively combining the strengths of the Riemann-Liouville, Liouville-Caputo, and Caputo-Fabrizio derivatives. The AB derivative provides several advantages: (i) its nonlocal nature guarantees that it captures the full historical behavior of the function being differentiated; (ii) it has an adjustable parameter, that lets researchers change the fractional order to enhance data-fitting accuracy; (iii) the AB derivative exhibits greater flexibility than its predecessor’s derivatives, enabling accurate modeling of complex systems; and (iv) it offers a unifying structure that can refine and improve the existing models across various scientific domains by integrating the features of the Riemann-Liouville, Liouville-Caputo and extending their applicability. Due to these compelling attributes, the AB derivative has gained significant attention and has been successfully applied to a wide range of real-world problems [9].

Time-fractional advection-diffusion equations (TFADEs) are widely used models in applied mathematics to describe various physical systems. The advection term represents the movement of a fluid along a concentration gradient, while the diffusion term describes the process by which material spreads from regions of higher to lower concentration over time. TFADEs are applied to transport processes, including the long-range dispersion of air pollutants [10], turbulence [11], water transport in soil [12], dispersion in porous media [13], shallow water flow [14], ion transport in heterogeneous media [15], blood flow with chemical interactions [16], and contaminant transport in soil [17].

Analytical solutions for TFADEs have been derived by various researchers. For example, Sanskrityayn and Kumar [18] derived analytical solutions to TFADEs using Green’s function methods. Avci and Yetim [19] obtained analytical solutions for TFADEs incorporating the Atangana-Baleanu fractional derivative, while Mirza and Vieru [20] derived fundamental solutions for TFADEs using the Caputo-Fabrizio derivative. However, obtaining exact solutions for TFADEs is often challenging due to the involvement of complex functions, which can be difficult to handle analytically.

As a result, developing accurate and efficient numerical schemes has become essential. Many authors in the literature have proposed numerical solutions for TFADEs. Umer et al. [21] analyzed numerical solutions of advection-diffusion equations with the Atangana-Baleanu fractional derivative using an extended cubic B-spline technique. Fazio et al. [22] studied a finite difference method on non-uniform meshes for TFADEs with a source term. Ahmed et al. [23] applied a Haar wavelet-based numerical technique to solve TFADEs. Kamran et al. [24] investigated numerical inverse Laplace transform methods for approximating TFADEs.

Other contributions include the work by Pareek et al. [25], who developed the natural transform method for solving TFADEs, and Chawla et al. [26], who utilized extended one-step time integration schemes. Liu et al. [27] proposed a radial basis function (RBF)-based differential quadrature method for solving two-dimensional TFADEs. Nguyen and Reynen [28] devised a space-time least squares finite element scheme for approximating TFADEs, while Cunha et al. [29] employed the boundary element method to solve TFADEs. Sweilam et al. [30] simulated TFADEs using a spectral collocation method combined with a non-standard finite difference technique.

Many authors have developed and modified numerical methods for the approximation of fractional partial differential equations (FPDEs) from various perspectives, focusing on improving accuracy, stability, efficiency, consistency, and performance in terms of computational cost. In recent years, hybrid methods have gained significant attention due to their high accuracy, low computational cost, and ease of implementation for discretizing FPDEs.

Hybrid methods combine two or more approaches, enabling them to mitigate the limitations of individual methods. As a result, hybrid methods can approximate complex problems in a simple and effective manner. In the literature, several researchers have proposed hybrid methods. For example, Yin et al. [31] combined the Laplace transform and Legendre wavelet methods for the numerical simulation of Klein-Gordon equations. Soares and Mansur [32] coupled the boundary element method with the finite element method for solving acoustic elastodynamic problems. A hybrid method based on the Laplace transform and Legendre wavelet approaches was analyzed in [33] for the approximation of Lane-Emden equations. Khan et al. [34] combined the homotopy perturbation method with the Laplace transform method to solve fractional models. Joujehi et al. [35] developed a hybrid method based on Beta functions and fractional-order Bernoulli wavelets for approximating multi-term TFPDEs in fluid mechanics. A review of the Jacobi-Galerkin spectral method for linear partial differential equations is examined by Hafez and Youssri [36]. Lim and Li [37] coupled the boundary element method with the finite difference method to approximate fluid-structure interaction problems with dynamic analysis of outer hair cells. Kamran et al. [38] combined the Laplace transform with radial basis functions for the numerical approximation of the mobile-immobile advection-dispersion problem arising in solute transport. Sahu and Jena [39] developed an efficient technique for time fractional Klein-Gordon equation based on modified Laplace Adomian decomposition technique via hybridized Newton-Raphson Scheme arises in relativistic fractional quantum mechanics.

The main objective of this work is to develop and analyze a hybrid Laplace Transform-based Chebyshev Spectral Collocation Method (LT-CSCM) for the efficient numerical solution of time-fractional advection-diffusion equations (TFADEs) featuring a nonsingular kernel. By combining the LT for temporal discretization with the CSCM for spatial discretization, our method aims to achieve high accuracy and computational efficiency.

The CSCM is a subclass of spectral methods that has recently garnered significant attention due to its straightforward implementation for the spatial approximation of fractional partial differential equations (FPDEs). Within the numerical framework for FPDEs, CSCM belongs to the family of Weighted Residual Methods (WRMs) [40]. This family includes the Tau, Galerkin, and collocation methods, each employing distinct techniques to minimize residuals. In the Galerkin and Tau methods, residuals are projected onto a polynomial space and constrained to be zero, while the collocation method enforces zero residuals at specific grid points. For FPDEs, Chebyshev collocation points are highly effective as grid points due to their optimal distribution, which enhances numerical accuracy. CSCM utilizes basis functions, typically Lagrange interpolation polynomials, defined at these points [41]. This global approach achieves spectral convergence, delivering high accuracy for problems with simple geometries and smooth solutions [42,43]. Compared to finite difference or finite element methods, CSCM is easier to implement and significantly reduces computational cost.

The CSCMs have been widely adopted by numerous researchers for many applications. Such as Khader and Saad [44] utilized the CSCM for the solution of the fractional Fisher problems. In [45], the authors proposed approximating fractional-order diffusion problems using CSCM combined with a power-series method based on residuals. The authors of [46] obtained the solution of the time-fractional advection-diffusion equation using CSCM. Tohidi [47] developed a numerical scheme for finding approximate solutions of one-dimensional parabolic partial differential equations (PDEs) under non-classical boundary conditions. In [48], the authors proposed a spectral collocation method based on differentiated Chebyshev polynomials to obtain numerical solutions for various types of nonlinear partial differential equations. Li et al. [49] employed the CSCM to solve the transport equation with given initial and boundary conditions. Rongpei et al. [50] solved two-dimensional nonlinear reaction-diffusion equations with Neumann boundary conditions using a new highly accurate CSCM.

The LT is an efficient and error-free method for the temporal discretization of FPDEs, addressing stability issues often encountered with traditional time-marching methods. These methods are stable and accurate only if the error calculated in a one-time step does not amplify as computations progress. In other words, time-marching methods remain stable if the error diminishes or remains unchanged during computations. Moreover, achieving optimal accuracy with time-marching methods typically requires smaller time steps, which significantly increases computational cost [51], thus affecting the overall efficiency of the method. One of its key advantages is its ability to convert differential equations into algebraic equations, making complex problems more manageable. Additionally, it provides a systematic approach for handling initial conditions directly within the transformed domain, avoiding the need for numerical time-stepping methods that may suffer from stability issues [52]. The LT is particularly beneficial for solving linear time-invariant systems and fractional differential equations, as it allows for analytical solutions in many cases [53]. However, the method also has limitations. It is less effective for nonlinear problems, as transforming the nonlinear terms is not straightforward, often requiring approximations or numerical techniques [54]. By employing the LT for temporal discretization, the FPDEs are transformed into the Laplace domain. To obtain the solution in the time domain, the inverse Laplace transform (ILT) must be applied. However, the exact computation of Evaluating the Bromwich integral is computationally complex, prompting the adoption of numerical inverse Laplace transform methods (ILTMs). Several authors have developed numerical ILTMs. For instance: De Hoog et al. [55] employed the quotient-difference scheme to formulate an improved ILTM that accelerates Fourier series convergence. Stehfest [56] developed a linear acceleration method using Salzer’s approach for numerical inversion of the LT. Talbot [57] introduced an efficient ILTM. Weideman and Trefethen [58] utilized parabolic and hyperbolic contours to approximate the Bromwich integral. Weeks [59] employed Laguerre functions for numerical ILTM. Each numerical method for inverting the Laplace transform has specific applications and is best suited to a particular problem. Contour integration methods, such as those based on hyperbolic, parabolic, or Talbot contours, are particularly effective for partial differential equations (PDEs) due to their ease of implementation and high accuracy. These methods deform the integration path in the complex plane to optimize convergence, minimizing computational complexity while maintaining precision. In this work, we utilize the numerical ILTMs described in [58,60]. Consider the following 2-D TFADE:

{∂τσu(x¯,τ)=Lu(x¯,τ)+ξ(x¯,τ),(x¯,τ)∈Θ×[0,1],Bu(x¯,τ)=ζ(x¯,τ),x¯∈∂Θ,u(x¯,0)=u0,(1)

where, x¯=(x,y)∈Θ=[−1,1]2⊂R2, ∂Θ is the bound of Θ,  L is the linear differential operator  Lu(x¯,τ)=a1Δu(x¯,τ)−a2∇u(x¯,τ), with a1,a2 positive constants, B a linear boundary operator, ξ(x¯,τ), ζ(x¯,τ) continuous functions, and the Atangana-Baleanu-Caputo (ABC) time-fractional derivative of order σ∈[0,1] defined as [8]

∂τσu(x¯,τ)=B(σ)1−σ∫0τ∂su(x¯,s)Eσ(−σ(τ−s)σ1−σ)ds,(2)

with

B(σ)=σΓ(σ)+(1−σ),

and Eσ(.) the Mittag-Leffler function defined as

Eσ(τ)=∑r=0∞τrΓ(rσ+1),σ>0, τ∈(−∞,∞).(3)

2  Existence and Uniqueness of the Solution

In this section, we utilize the fixed-point theory to prove the existence and uniqueness of the solution to the fractional advection-diffusion model (1). Let us define a Banach space C(Λ,R) of continuous functions from Λ=Θ×[0,1] into R, equipped with the norm defined by ‖u‖∞:=sup{|u(x¯,τ)|;(x¯,τ)∈Λ}. By applying the AB-fractional integral operator to (1), we obtain

u(x¯,τ) = u0+1−σB(σ)(a1Δu(x¯,τ)−a2∇u(x¯,τ)+ξ(x¯,τ))+σΓ(σ)B(σ)∫0τ(τ−s)σ−1(a1Δu(x¯,τ)−a2∇u(x¯,τ)+ξ(x¯,τ))ds.(4)

Let define the operator X:C(Λ,R)→C(Λ,R), that reformulates problem (1) as a fixed-point problem, given by

Xu(x¯,τ) = u0+1−σB(σ)(a1Δu(x¯,τ)−a2∇u(x¯,τ)+ξ(x¯,τ))+σΓ(σ)B(σ)∫0τ(τ−s)σ−1(a1Δu(x¯,τ)−a2∇u(x¯,τ)+ξ(x¯,τ))ds.(5)

The fixed point of X, corresponds to the solution of problem (1). We assume the following hypotheses before proving our main results: for any (x¯,τ) ∈ Λ, there exist κk>0, k=1,2,3,4,5,6, such that

(A1)|Δu(x¯,τ)|≤κ1|u(x¯,τ)|,(A2)|∇u(x¯,τ)|≤κ2|u(x¯,τ)|,(A3)|u0|≤κ3,(A4)|ξ(x¯,τ)|≤κ4,(A5)|Δu1(x¯,τ)−Δu2(x¯,τ)|≤κ5|u1(x¯,τ)−u2(x¯,τ)|,(A6)|∇u1(x¯,τ)−∇u2(x¯,τ)|≤κ6|u1(x¯,τ)−u2(x¯,τ)|.

Theorem 1. If assumptions (A5,A6) hold then problem (1) has a unique solution.

Proof. The proof consists of several steps.

Step 1: Continuity of X. We will show that X is continuous. Suppose a sequence um→u, where u∈C(Λ,R). For (x¯,τ)∈Λ and using the bounds in (A5,A6) we can write

‖Xum(x¯,τ)−Xu(x¯,τ)‖∞=sup{|Xum(x¯,τ)−Xu(x¯,τ)|}=sup{|1−σB(σ)(a1Δum(x¯,τ)−a2∇um(x¯,τ))+σΓ(σ)B(σ)∫0τ(τ−s)σ−1(a1Δum(x¯,τ)−a2∇um(x¯,τ))ds−1−σB(σ)(a1Δu(x¯,τ)−a2∇u(x¯,τ))−σΓ(σ)B(σ)∫0τ(τ−s)σ−1(a1Δu(x¯,τ)−a2∇u(x¯,τ))ds|}≤ sup{1−σB(σ)(|a1||Δum(x¯,τ)−Δu(x¯,τ)|+|a2||∇um(x¯,τ)−∇u(x¯,τ)|)−σΓ(σ)B(σ)∫0τ(τ−s)σ−1(|a1||Δum(x¯,τ)−Δu(x¯,τ)|+|a2||∇um(x¯,τ)−∇u(x¯,τ)|)ds}≤ sup{1−σB(σ)+σΓ(σ)B(σ)∫0τ(τ−s)σ−1(κ5|a1||um(x¯,τ)−u(x¯,τ)|+κ6|a2||um(x¯,τ)−u(x¯,τ)|)ds}≤ sup{1−σB(σ)+τσΓ(σ)B(σ)(κ5|a1|+κ6|a2|)|um(x¯,τ)−u(x¯,τ)|}≤ τσ+(1−σ)Γ(σ)Γ(σ)B(σ)(κ5|a1|+κ6|a2|)‖um−u‖∞,

Since u is continuous. Hence, we obtain

‖um−u‖∞→0,as m→∞.

Therefore, X is continuous.

Step 2: Boundedness of X. Next, we show that X maps bounded sets to bounded sets. It is sufficient to prove that for γ>0 there exists ρ>0 such that if u∈Rγ={u∈C(Λ,R):‖u‖∞≤γ} then it is ‖Xu‖∞≤ϱ. For τ∈[0,T] and utilizing the assumptions (A1–A4) we have

|Xu(x¯,τ)|=|u0+1−σB(σ)(a1Δu(x¯,τ)−a2∇u(x¯,τ)+ξ(x¯,τ))+σΓ(σ)B(σ)∫0τ(τ−s)σ−1(a1Δu(x¯,τ)−a2∇u(x¯,τ)+ξ(x¯,τ))ds|≤ |u0|+1−σB(σ)(|a1||Δu(x¯,τ)|+|a2||∇u(x¯,τ)|+|ξ(x¯,τ)|)+σΓ(σ)B(σ)∫0τ(τ−s)σ−1(|a1||Δu(x¯,τ)|+|a2||∇u(x¯,τ)|+|ξ(x¯,τ)|)ds≤ κ3+1−σB(σ)(κ1|a1||u(x¯,τ)|+κ2|a2||u(x¯,τ)|+κ4)+σΓ(σ)B(σ)∫0τ(τ−s)σ−1(κ1|a1||u(x¯,τ)|+κ2|a2||u(x¯,τ)|+κ4)ds

≤κ3+1−σB(σ)(|a1|κ1‖u‖∞+|a2|κ2‖u‖∞+κ4)+σΓ(σ)B(σ)∫0τ(τ−s)σ−1(|a1|κ1‖u‖∞+|a2|κ2‖u‖∞+κ4)ds≤κ3+(1−σ)Γ(σ)B(σ)(|a1|κ1‖u‖∞+|a2|κ2‖u‖∞+κ4)+τσΓ(σ)B(σ)(κ1|a1|‖u‖∞+|a2|κ2‖u‖∞+κ4)≤κ3+(1−σ)Γ(σ)+τσΓ(σ)B(σ)(κ1|a1|‖u‖∞+κ2|a2|‖u‖∞+κ4)≤κ3+(1−σ)Γ(σ)+τσΓ(σ)B(σ)((κ1|a1|+κ2|a2|)‖u‖∞+κ4).

Defining

ρ=κ3+τσ+Γ(σ)(1−σ)Γ(σ)B(σ)((κ1|a1|+κ2|a2|)γ+κ4),

it follows that ‖Xu‖∞≤ρ proving that X maps bounded sets to bounded sets.

Step 3: Equicontinuity of X. Let consider u∈Rγ and (x¯1,τ1),(x¯2,τ2)∈Λ, such that x¯1<x¯2, τ1<τ2. We have

|Xu(x¯1,τ1)−Xu(x¯2,τ2)|=|1−σB(σ)(a1Δu(x¯1,τ1)−a2∇u(x¯1,τ1)+ξ(x¯1,τ1))+σΓ(σ)B(σ)∫0τ1(τ1−s1)σ−1(a1Δu(x¯1,τ1)−a2∇u(x¯1,τ1)+ξ(x¯1,τ1))ds1−1−σB(σ)(a1Δu(x¯2,τ2)−a2∇u(x¯2,τ2))−σΓ(σ)B(σ)∫0τ2(τ2−s2)σ−1(a1Δu(x¯2,τ2)−a2∇u(x¯2,τ2)+ξ(x¯2,τ2))ds2|≤1−σB(σ)(|a1||Δu(x¯1,τ1)|+|a2||∇u(x¯1,τ1)|+|ξ(x¯1,τ1)|)+σΓ(σ)B(σ)∫0τ1(τ1−s1)σ−1(|a1||Δu(x¯1,τ1)|+|a2||∇u(x¯1,τ1)|+|ξ(x¯1,τ1)|)ds1−1−σB(σ)(|a1||Δu(x¯2,τ2)|+|a2||∇u(x¯2,τ2)|+|ξ(x¯2,τ2)|)−σΓ(σ)B(σ)∫0τ2(τ2−s2)σ−1(|a1||Δu(x¯2,τ2)|+|a2||∇u(x¯2,τ2)|+|ξ(x¯2,τ2)|)ds2≤1−σB(σ)(κ1|a1||u(x¯1,τ1)|+κ2|a2||u(x¯1,τ1)|+κ4)+σΓ(σ)B(σ)∫0τ1(τ1−s1)σ−1(|a1|κ1|u(x¯1,τ1)|+|a2|κ2|u(x¯1,τ1)|+κ4)ds1−1−σB(σ){|a1|κ1|u(x¯2,τ2)|+|a2|κ2|u(x¯2,τ2)|+κ4}−σΓ(σ)B(σ)∫0τ1(τ2−s2)σ−1(|a1|κ1|u(x¯2,τ2)|+|a2|κ2|u(x¯2,τ2)|+κ4)ds2≤1−σB(σ)(|a1|κ1‖u‖∞+|a2|κ2‖u‖∞+κ4)+σΓ(σ)B(σ)∫0τ1(τ1−s1)σ−1(|a1|κ1‖u‖∞+|a2|κ2‖u‖∞+κ4)ds1−1−σB(σ)(|a1|κ1‖u‖∞+|a2|κ2‖u‖∞+κ4)−σΓ(σ)B(σ)∫0τ2(τ2−s2)σ−1(|a1|κ1‖u‖∞+|a2|κ2‖u‖∞+κ4)ds2

=σΓ(σ)B(σ)∫0τ1(τ1−s1)σ−1(|a1|κ1‖u‖∞+|a2|κ2‖u‖∞+κ4)ds1−σΓ(σ)B(σ)∫0τ2(τ2−s2)σ−1(|a1|κ1‖u‖∞+|a2|κ2‖u‖∞+κ4)ds2=(∫0τ1(τ1−s1)σ−1ds1−∫0τ2(τ2−s2)σ−1ds2)[σΓ(σ)B(σ)(|a1|κ1‖u‖∞+|a2|κ2‖u‖∞+κ4)]=τ1σ−τ2σΓ(σ)B(σ)(|a1|κ1‖u‖∞+|a2|κ2‖u‖∞+κ4).

It follows that

|Xu(x¯1,τ1)−Xu(x¯2,τ2)|→0as τ1 → τ2.

By the Arzelà-Ascoli Theorem [61], the operator X is completely continuous.

Step 4: A Priori Bound. Define χ={u∈C(Λ,R):u=εXu,ε∈(0,1)}. We prove that χ is bounded. If u∈χ, then u=εXu, with 0<ε<1. Then for τ∈[0,T], using Eq. (5) we have

|u|=|εXu|=|ε×(u0+1−σB(σ)(a1Δu(x¯,τ)−a2∇u(x¯,τ)+ξ(x¯,τ))+σΓ(σ)B(σ)∫0τ(τ−s)σ−1(a1Δu(x¯,τ)−a2∇u(x¯,τ)+ξ(x¯,τ))ds)|.

Now, using assumptions (A1−A4) we get

|u|=ε×(κ3+(1−σ)Γ(σ)+τσΓ(σ)B(σ)((κ1|a1|+κ2|a2|)‖u‖∞+κ4)),

which gives

‖u‖∞=ε×ρ,

where ρ is as defined in Step 2. This shows that χ is bounded. By the Schaefer Fixed Point Theorem [61], X has at least one fixed point. Consequently, the considered problem has at least one solution. □ (25A1)

Theorem 2. The problem define in Eq. (1) has a unique solution if the assumptions (A5,A6) hold and

(1−σ)Γ(σ)+τσΓ(σ)B(σ){κ5|a1|+κ6|a2|}<1.

Proof

‖Xu1(x¯,τ)−Xu2(x¯,τ)‖∞=sup{|Xu1(x¯,τ)−Xu2(x¯,τ)|}=sup{|1−σB(σ)(a1Δu1(x¯,τ)−a2∇u1(x¯,τ))+σΓ(σ)B(σ)∫0τ(τ−s)σ−1(a1Δu1(x¯,τ)−a2∇u1(x¯,τ))ds−1−σB(σ)(a1Δu2(x¯,τ)−a2∇u2(x¯,τ))−σΓ(σ)B(σ)∫0τ(τ−s)σ−1(a1Δu2(x¯,τ)−a2∇u2(x¯,τ))ds|}≤sup{1−σB(σ)(|a1||Δu1(x¯,τ)−Δu2(x¯,τ)|+|a2||∇u1(x¯,τ)−∇u2(x¯,τ)|)−σΓ(σ)B(σ)∫0τ(τ−s)σ−1(|a1||Δu1(x¯,τ)−Δu2(x¯,τ)|+|a2||∇u1(x¯,τ)−∇u2(x¯,τ)|)ds}≤sup{1−σB(σ)+σΓ(σ)B(σ)∫0τ(τ−s)σ−1(κ5|a1||u1(x¯,τ)−u2(x¯,τ)|+κ6|a2||u1(x¯,τ)−u2(x¯,τ)|)ds}≤sup{1−σB(σ)+τσΓ(σ)B(σ)(κ5|a1|+κ6|a2|)|u1(x¯,τ)−u2(x¯,τ)|}≤τσ+(1−σ)Γ(σ)Γ(σ)B(σ)(κ5|a1|+κ6|a2|)‖u1−u2‖∞.

Thus, we find that under the given assumptions X is a contraction. By the Banach Fixed-Point Theorem [61], X has a unique fixed point. Therefore, the problem defined in Eq. (1) has a unique solution. □ (25A1)

3  Methodology

In the LT-CSCM approach, the Laplace Transform (LT) transforms Eq. (1) into the Laplace domain, enabling efficient temporal discretization. The Chebyshev Spectral Collocation Method (CSCM) is then applied to discretize spatial variables with high accuracy. Finally, the time-domain solution is recovered using the Talbot method for numerical inverse Laplace transform, ensuring precision and computational efficiency.

3.1 Temporal Discretization

The Laplace Transform is used for the temporal discretization of the proposed problem defined in Eq. (1). The LT of u(x¯,τ) is defined as

ℒ{u(x¯,τ)}=u^(x¯,s)=∫0∞e−sτu(x¯,τ)dτ.

The LT of the ABC derivative, ∂τσu(x¯,τ), defined in Eq. (2), is given by [8]

ℒ{∂τσu(x¯,τ)}=B(σ)(sσu^(x¯,s)−sσ−1u0)sσ(1−σ)+σ.

Applying the LT to Eq. (1) yields

B(σ)(sσu^(x¯,s)−sσ−1u0)sσ(1−σ)+σ=Lu^(x¯,s)+ξ^(x¯,s),x¯∈Θ,(6)

and

Bu^(x¯,s)=ζ^(x¯,s),x¯∈∂Θ.(7)

Simplifying further, we obtain

{(B(σ)sσsσ(1−σ)+σ)I−L}u^(x¯,s)=S^(x¯,s),x¯∈Θ,(8)

and

Bu^(x¯,s)=ζ^(x¯,s),x¯∈∂Θ,(9)

where

S^(x¯,s)=B(σ)sσ−1u0sσ(1−σ)+σ+ξ^(x¯,s),

and I is the (N+1)×(N+1) identity operator. The time independent problem in Eqs. (8) and (9) is solved for each s in the Laplace space, using CSCM for spatial discretization. Finally, numerical ILTM is used to to recover the time-domain solution of (1). The next section explains the CSCM approach.

3.2 Spatial Discretization by Chebyshev Spectral Collocation Method

In CSCM, a global polynomial interpolant is utilized on specific nodes (Chebychev nodes) to approximate the unknown solution of a FPDEs. The spatial derivatives are computed using discrete derivative operators, also called differentiation matrices (DM) [41].

The solution is considered over [−1,1] and interpolates {(xm,u^(xm))}, by [42,46]

IN(x)=∑m=0Nlm(x)u^m,

where u^m=u^(xm), and the basic Lagrange polynomials are as follows:

lm(x)=(x−x0)…(x−xm−1)(x−xm+1)…(x−xN)(xm−x0)…(xm−xm−1)(xm−xm+1)…(xm−xN).(10)

For spatial discretization in [−1,1], the Chebyshev nodes are usually considered

xm=cos⁡{(lπN)}m=0N.(11)