Numerical Solution of Advection–Diffusion Equation of Fractional Order Using Chebyshev Collocation Method

Abstract

This work presents a highly accurate method for the numerical solution of the advection–diffusion equation of fractional order. In our proposed method, we apply the Laplace transform to handle the time-fractional derivative and utilize the Chebyshev spectral collocation method for spatial discretization. The primary motivation for using the Laplace transform is its ability to avoid the classical time-stepping scheme and overcome the adverse effects of time steps on numerical accuracy and stability. Our method comprises three primary steps: (i) reducing the time-dependent equation to a time-independent equation via the Laplace transform, (ii) employing the Chebyshev spectral collocation method to approximate the solution of the transformed equation, and (iii) numerically inverting the Laplace transform. We discuss the convergence and stability of the method and assess its accuracy and efficiency by solving various problems in two dimensions.

1. Introduction

Fractional calculus (FC) is the generalization of classical calculus, which includes integrals and derivatives in non-integer order [1,2]. In recent years, the concept of fractional derivatives has been applied with great success to model various real-life phenomena in many scientific fields [3]. Fractional order operators include the history of a physical phenomenon from the initial state to the current state. Therefore, fractional order operators are often applied to model systems that describe the influence of memory effects [4]. Fractional order operators have broadened the concept of integer order operators [5]. In the last ten years, PDEs of fractional order have attracted the remarkable attention of the research community due to their ability to model various problems, such as neural networks [6], optimal voltage controllers for electronic vehicle charging stations [7], robotics [8], optimal controls for impulsive systems [9], hydrology [10], medical sciences [11], and many others [12].

The time-fractional advection–diffusion equations (TFADEs) have been widely used in many fields of science because of their ability to model memory and non-local properties [13]. ADEs combine diffusion and advection terms to describe physical events in which particles, energy, or other physical quantities are moved inside a physical system by two processes, i.e., diffusion and advection. The concentration of substances for mass and heat transfer is represented by the advection and diffusion models. This equation has been used to model heat transfer in draining films [14], flow in porous media [15], and air pollution [16], etc. There have been many attempts to solve the TFADEs accurately. In [17], a second-order implicit difference method was proposed for solving a time–space-fractional ADE. Mardani et al. [18] studied a meshless method for the solution of ADEs with variable coefficients. A fully implicit finite difference scheme based on cubic B-splines was developed in [19] for solving TFADEs. Abbaszadeh and Amjadian [20] studied the spectral element method for the solution of TFADEs. In the current study, we investigate the numerical solution of two-dimensional TFADEs using the Chebyshev spectral collocation method (CSCM) coupled with the Laplace transform (LT).

Spectral collocation methods have attracted the attention of the research community due to their highly accurate solutions to partial differential equations. They converge exponentially fast as compared to the algebraic convergence rates for FEM and FDM. Spectral methods are global in nature. This means that the optimal accuracy can be obtained with few grid points. For the first time, spectral methods were introduced by Gottlieb and Orszag [21]. Their work paved the way for more advanced extensions of the method for a larger variety of problems [22,23]. The most important tool for spectral methods is the computation of differentiation matrices, discussed briefly by Welfert [24]. The spectral collocation methods have been used by many authors. For example, in [25], the authors discussed the Hermite spectral methods for unbounded domains. The solution of the fractional advection–dispersion equation by the Chebyshev spectral method was incorporated by Sweilam et al. [26]. In [27], the spectral Legendre–Chebyshev collocation method for variable coefficients was examined. Tian et al. [28] investigated the polynomial spectral collocation method for space-fractional ADEs. The spectral collocation method and the non-standard finite difference technique for the solution of TFADEs were introduced by [29].

In the aforementioned methods, the time discretization is carried out using the finite difference time-stepping method. The main shortcoming of the time-stepping technique is that it may not always result in a stable solution. A finite difference time-stepping scheme is stable if the errors calculated at a one-time step do not cause the errors to be increased during the calculations. In other words, the time-stepping scheme is stable whether the errors remain constant or decay during the computations. Furthermore, for optimal accuracy, we need to decrease the time step, which results in an increased computational time, and, thus, has low efficiency in simulating the time-fractional problems. To overcome this drawback of the finite difference time-stepping scheme, the LT has been used as one of the best alternative mathematical tools [30,31]. However, while using the LT, the main difficulty is its inversion. In the literature, many algorithms have been developed for the numerical inversion of the LT [32,33,34]. In the current study, we use Stehfest’s method [35] and Talbot’s method [32,36] for the numerical inversion of the LT.

2. Numerical Scheme

In our numerical scheme, a two-dimensional TFADE is taken into account. First, the LT is employed to reduce the TFADE to a time-independent problem in the LT domain. Then, the CSCM is employed to obtain the approximate solution to the transformed problem. Finally, the time domain solution is retrieved from the LT domain solution via the inverse LT.

2.1. Two-Dimensional TFADE

We consider a two-dimensional time-fractional advection–diffusion equation of the form

$$D_{\tau }^{\alpha }U(\overline{\mathbf{x}},\tau )={\lambda }_1\Delta U(\overline{\mathbf{x}},\tau )+(\overrightarrow{\mathbf{u}}·\nabla )U(\overline{\mathbf{x}},\tau )-{\lambda }_{2}U(\overline{\mathbf{x}},\tau )+\mathcal{Q}(\overline{\mathbf{x}},\tau ),\overline{\mathbf{x}}\in \Omega ,0\le \tau \le 1,$$

(1)

with boundary conditions

$$U(\overline{\mathbf{x}},\tau )=h(\overline{\mathbf{x}},\tau ),\overline{\mathbf{x}}\in \partial \Omega ,0\le \tau \le 1,$$

(2)

and the initial condition

$$U(\overline{\mathbf{x}},0)=U_0\left(\overline{\mathbf{x}}\right),\overline{\mathbf{x}}\in \Omega ,$$

(3)

where  $\overrightarrow{\mathbf{u}}=(u_1,u_2)$ is a known vector, ${\lambda }_1$ and ${\lambda }_2$ are arbitrary constants, $\mathcal{Q}(\overline{\mathbf{x}},\tau ),U_0\left(\overline{\mathbf{x}}\right),$$h(\overline{\mathbf{x}},\tau )$ are given continuous functions, $\Omega \subset {\mathbb{R}}^2$ is a bounded domain with a smooth boundary $\partial \Omega$, and $D_{\tau }^{\alpha }$ represents Caputo’s derivative, defined as [3]

$$D_{\tau }^{\alpha }U(\overline{\mathbf{x}},\tau )=\frac{1}{\Gamma (1-\alpha )}{\int }_0^{\tau }\frac{\partial U(\overline{\mathbf{x}},z)}{\partial z}{(\tau -z)}^{-\alpha }dz,0<\alpha \le 1.$$

2.2. Time Discretization

In this section, the LT is used to reduce the TFADE to an equivalent time-independent problem in the LT domain. The LT of the function $U(\overline{\mathbf{x}},\tau )$ is denoted by $\widehat{U}(\overline{\mathbf{x}},s)$ and is defined as

$$\widehat{U}(\overline{\mathbf{x}},s)={\int }_0^{\infty }{e}^{-s\tau }U(\overline{\mathbf{x}},\tau )d\tau ,$$

and the LT of Caputo’s derivative is defined as

$$\mathcal{L}\left\{D_{\tau }^{\alpha }U(\overline{\mathbf{x}},\tau )\right\}=s^{\alpha }\widehat{U}(\overline{\mathbf{x}},s)-s^{\alpha -1}{U}_0.$$

Applying the LT to (1)–(3) gives

$$s^{\alpha }\widehat{U}(\overline{\mathbf{x}},s)-s^{\alpha -1}{U}_0={\lambda }_1\Delta \widehat{U}(\overline{\mathbf{x}},s)+(\overrightarrow{\mathbf{u}}·\nabla )\widehat{U}(\overline{\mathbf{x}},s)-{\lambda }_2\widehat{U}(\overline{\mathbf{x}},s)+\widehat{\mathcal{Q}}(\overline{\mathbf{x}},s),$$

(4)

and

$$\widehat{U}(\overline{\mathbf{x}},s)=\widehat{h}(\overline{\mathbf{x}},s),$$

(5)

Equations (4) and (5) can be rearranged as

$$(s^{\alpha }I-\mathbf{L})\widehat{U}(\overline{\mathbf{x}},s)=\widehat{G}(\overline{\mathbf{x}},s),$$

(6)

and

$$\widehat{U}(\overline{\mathbf{x}},s)=\widehat{h}(\overline{\mathbf{x}},s).$$

(7)

where $\mathbf{L}={\lambda }_1\Delta +(\overrightarrow{\mathbf{u}}·\nabla )-{\lambda }_{2}I$ is a linear differential operator, $\widehat{G}(\overline{\mathbf{x}},s)=s^{\alpha -1}{U}_0+\widehat{\mathcal{Q}}(\overline{\mathbf{x}},s)$, and I denotes the identity operator. To obtain the approximate solution of the system defined in (6) and (7), first, the linear spatial operator $\mathbf{L}$ is discretized using the CSCM. Then, the fully discrete system is solved for each node s in the LT domain. Finally, the numerical inversion of the LT is performed to obtain the solution of the problem defined in Equations (1)–(3).

2.3. Chebyshev Spectral Collocation Method

In the CSCM, we interpolate the data $\left\{\left(x_j,\widehat{U}\left(x_j\right)\right)\right\}$ by using the Lagrange interpolation polynomial (LIP) $l_j\left(x\right)$ of a degree of at most n [23,37]

$${\mathcal{I}}_n\left(x\right)=\sum _{j=0}^n{l}_j\left(x\right){\widehat{U}}_j$$

where $l_j\left(x\right)$ is the LIP at the point $x_j(j=0,1,\dots ,n)$, given as

$$l_j\left(x\right)=\frac{(x-x_0)\dots (x-x_{j-1})(x-x_{j+1})\dots (x-x_n)}{(x_j-x_0)\dots (x_j-x_{j-1})(x_j-x_{j+1})\dots (x_j-x_n)}$$

(8)

where ${\widehat{U}}_j=\widehat{U}\left(x_j\right).$ In the CSCM, the Chebyshev nodes are selected as the interpolation nodes. The domain $[-1,1]$ is discretized using the Chebyshev points defined as

$$x_j=cos\left(\frac{j\pi }{n}\right),j=0,1,\dots ,n.$$

(9)

The approximation of the first derivative $\frac{\partial \widehat{U}\left(x\right)}{\partial x}$ on the Chebyshev points can be obtained as

$$\frac{\partial \widehat{U}\left(x\right)}{\partial x}\approx {\mathbf{D}}_n\widehat{U},$$

where the entries of the matrix ${\mathbf{D}}_n$ are

$${\left[{\mathbf{D}}_n\right]}_{i,j}=l_j^{\prime }\left(x_i\right),i,j=1,2,\dots ,n.$$

The off-diagonal entries of ${\left[{\mathbf{D}}_n\right]}_{i,j}$ are obtained as

$${\left[{\mathbf{D}}_n\right]}_{i,j}=\frac{{\alpha }_j}{{\alpha }_i(x_i-x_j)},i\ne j,$$

where ${\alpha }_j^{-1}={\prod }_{i\ne j}^n(x_i-x_j),$ and the diagonal elements are obtained as

$${\left[{\mathbf{D}}_n\right]}_{i,j}=-\sum _{j=0,j\ne i}^n{\left[{\mathbf{D}}_n\right]}_{i,j},i=0,1,2,..,n.$$

Next, the entries of the $\mu th$-order differentiation matrix ${\mathbf{D}}_n^{\mu }$ can be obtained analytically as

$${\left[{\mathbf{D}}_n^{\mu }\right]}_{i,j}=l_j^{\mu }\left(x_i\right),i,j=1,2,\dots ,n.$$

A more accurate and stable evaluation can be found in [24,38]. In [24], the author derived a useful recursion relation for the differentiation matrices, given as

$${\left[{\mathbf{D}}_n^{\mu }\right]}_{kj}=\frac{\mu }{x_k-x_j}\left(\frac{{\alpha }_j}{{\alpha }_k}{\left[{\mathbf{D}}_n^{(\mu -1)}\right]}_{kk}-{\left[{\mathbf{D}}_n^{(\mu -1)}\right]}_{kj}\right),k\ne j$$

For square domain ${[-1,1]}^2,$ assume the points ${\overline{\mathbf{x}}}_{ij}$ are defined as

$${\overline{\mathbf{x}}}_{ij}=\left(cos\left(\frac{i\pi }{n}\right),cos\left(\frac{j\pi }{n}\right)\right),i,j=0,1,2,\dots ,n$$

The LIP corresponding to the above Chebyshev points are given as

$$l_{ij}\left(\overline{\mathbf{x}}\right)=l_i\left(x\right)l_j\left(y\right),$$

(10)

where $l_{ij}\left({\overline{\mathbf{x}}}_{ij}\right)={\delta }_{ij},i,j=0,1,2,\dots ,n.$ In terms of x and y, the second derivatives of the LIPs (10) are given as

$$\frac{{\partial }^2{l}_{ij}\left({\overline{\mathbf{x}}}_{rs}\right)}{\partial x^2}=l_i^{″}\left(x_r\right)l_j\left(y_s\right)={\left[{\mathbf{D}}_n^2\right]}_{ri}{\delta }_{js},$$

$$\frac{{\partial }^2{l}_{ij}\left({\overline{\mathbf{x}}}_{rs}\right)}{\partial y^2}=l_i\left(x_r\right)l_j^{″}\left(y_s\right)={\delta }_{ri}{\left[{\mathbf{D}}_n^2\right]}_{sj},$$

where ${\mathbf{D}}_n^2$ denotes the second-order Chebyshev differentiation matrix. Applying the linear operator $\mathbf{L}$ on the LIP at the Chebyshev nodes, $\left\{{\overline{\mathbf{x}}}_{rs}\right\}$ gives

$$\mathbf{L}\left(l_{ij}\left({\overline{\mathbf{x}}}_{rs}\right)\right)={\lambda }_1\left({\left[{\mathbf{D}}_n^2\right]}_{ri}{\delta }_{js}+{\delta }_{ri}{\left[{\mathbf{D}}_n^2\right]}_{sj}\right)+\left(u_1\left({\left[{\mathbf{D}}_n\right]}_{ri}{\delta }_{js}\right)+u_2\left({\delta }_{ri}{\left[{\mathbf{D}}_n\right]}_{sj}\right)\right)-{\lambda }_2{\delta }_{ri}{\delta }_{sj},$$

(11)

Therefore, the discretized form of the linear differential operator $\mathbf{L}$ via the CSCM can be expressed as

$${\mathbf{L}}_{\mathbf{D}}={\lambda }_1\left(I_n\otimes {\mathbf{D}}_n^2+{\mathbf{D}}_n^2\otimes I_n\right)+\left(u_1\left(I_n\otimes {\mathbf{D}}_n\right)+u_2\left({\mathbf{D}}_n\otimes I_n\right)\right)-{\lambda }_2\left(I_n\otimes I_n\right).$$

(12)

By using Equation (12) in Equation (6), we have

$$(s^{\alpha }I-{\mathbf{L}}_{\mathbf{D}})\widehat{U}(\overline{\mathbf{x}},s)=\widehat{G}(\overline{\mathbf{x}},s),$$

(13)

where we incorporate the boundary conditions in Equation (7) by taking the matrix ${\mathbf{L}}_{\mathbf{D}}$, based on all collocation points $\overline{\mathbf{x}}$, and then replace the rows of ${\mathbf{L}}_{\mathbf{D}}$ corresponding to collocation at boundary points with unit vectors that have a one in the position corresponding to the diagonal of ${\mathbf{L}}_{\mathbf{D}}$. Thus, the boundary condition $U(\overline{\mathbf{x}},s)=\widehat{h}(\overline{\mathbf{x}},s)$ in (7) will be explicitly enforced at this point as soon as we set the right-hand side to the corresponding value of $\widehat{h}(\overline{\mathbf{x}},s)$ [23]. After incorporating the boundary conditions and solving Equation (13) for each node s, we obtain the approximate solution $\widehat{U}(\overline{\mathbf{x}},s)$ in the LT domain. In Section 2.4, we describe the numerical inverse LT method used to evaluate the approximate time domain solution $U(\overline{\mathbf{x}},\tau )$ for the original problem in Equations (1)–(3).

Stability and Error Analysis

Given the Chebyshev nodes in Equation (9) and the Lagrange interpolation polynomials defined in Equation (8), the interpolation operator is defined as follows [37]:

$${\mathcal{I}}_n:C(\Omega )\to {\mathcal{P}}_n,{\mathcal{I}}_n\left(U\right)=\sum _{j=0}^{n}U\left(x_j\right)l_j\left(x\right).$$

(14)

The steps proposed by Borm et al. [39] for constructing the error bound are utilized here. Let $M_n$ be a constant with the stability estimate given as

$$\parallel {\mathcal{I}}_n{\left(U\right)\parallel }_{\infty }\le M_n{\parallel U\parallel }_{\infty },forallU\in C[-1,1].$$

(15)

$${\mathcal{I}}_n\left(U\right)=U,forallU\in {\mathcal{P}}_n.$$

(16)

For the Chebyshev interpolation, we have

$$M_n=\frac{ln(n+1)}{\left(\frac{\pi }{2}\right)}+1\le (n+1).$$

(17)

Thus, the stability constant depends on n and develops extremely slowly. Additionally, the following approximation error bound [39] holds for a function $U\in C^{(n+1)}[-1,1],$

$$\parallel U-{\mathcal{I}}_n{\left(U\right)\parallel }_{\infty }\le \frac{2^{-n}}{(n+1)!}{\parallel U^{(n+1)}\parallel }_{\infty }.$$

(18)

Theorem 1

([39]). If (15) and (18) hold for $U\in C^{(n+1)}[-1,1],$ and considering $k\in \{0,1,2,\dots ,n\},$ then

$$\parallel U^{\left(k\right)}-{\mathcal{I}}_n{\left(U\right)}^{\left(k\right)}{\parallel }_{\infty }\le \frac{2(M_n^{\left(k\right)}+1)}{(n-k+1)!}{\left(\frac{1}{2}\right)}^{(n-k+1)}{\parallel U^{(n+1)}\parallel }_{\infty }$$

(19)

depending on the stability constant

$$M_n^{\left(k\right)}=\frac{M_n}{k!}\left(\frac{n!}{(n-k)!}\right).$$

Now, use (18) and (19) for the TFADE (1) in one dimension. For the $1-D$ case, the linear differential operator $\mathbf{L}$ is of the form $\mathbf{L}={\lambda }_1\frac{{\partial }^2}{\partial x^2}+u_1\frac{\partial }{\partial x}-{\lambda }_{2}I.$ Thus, we find the following error estimate

$$\begin{array}{cc}\hfill \mathcal{E}& =\parallel (D_t^{\alpha }U-\mathbf{L}U)-(D_t^{\alpha }{\mathcal{I}}_n\left(U\right)-\mathbf{L}{\mathcal{I}}_n\left(U\right)){\parallel }_{\infty }\hfill \\ \hfill & =\parallel (D_t^{\alpha }(U-{\mathcal{I}}_n\left(U\right))-\mathbf{L}(U-{\mathcal{I}}_n\left(U\right)){\parallel }_{\infty }\hfill \\ \hfill & \le \parallel \left(D_t^{\alpha }(U-{\mathcal{I}}_n\left(U\right))\right){\parallel }_{\infty }+{\parallel \mathbf{L}(U-{\mathcal{I}}_n\left(U\right))\parallel }_{\infty }\hfill \\ \hfill & \le \parallel D_t^{\alpha }(U-{\mathcal{I}}_n\left(U\right)){\parallel }_{\infty }+|{\lambda }_1|\parallel U_{xx}-{\mathcal{I}}_n{\left(U\right)}_{xx}{\parallel }_{\infty }\hfill \\ \hfill & +|u_1|\parallel U_x-{\mathcal{I}}_n{\left(U\right)}_x{\parallel }_{\infty }+|{\lambda }_2|\parallel U-{\mathcal{I}}_n{\left(U\right)\parallel }_{\infty }\hfill \\ \hfill \mathcal{E}& \le \parallel D_t^{\alpha }(U-{\mathcal{I}}_n\left(U\right)){\parallel }_{\infty }+|{\lambda }_1|\frac{2(M_n^{\left(2\right)}+1)}{(n-1)!}{\left(\frac{1}{2}\right)}^{n-1}{\parallel U^{(n+1)}\parallel }_{\infty }\hfill \\ \hfill & +|u_1|\frac{2(M_n^{\left(1\right)}+1)}{n!}{\left(\frac{1}{2}\right)}^n\parallel U^{(n+1)}{\parallel }_{\infty }+|{\lambda }_2|\frac{2^{-n}}{(n+1)!}{\parallel U^{(n+1)}\parallel }_{\infty }\hfill \end{array}$$

The time derivative is accurately evaluated. Therefore, the error bound of $\parallel \left(D_t^{\alpha }(U-{\mathcal{I}}_n\left(U\right))\right){\parallel }_{\infty }$ is the same order of $\parallel (U-{\mathcal{I}}_n\left(U\right)){\parallel }_{\infty }.$ Finally, we have [37]:

$$\mathcal{E}\le C\parallel U^{(n+1)}{\parallel }_{\infty }.$$

where C is constant, resulting from calculation of the coefficients of $\parallel U^{(n+1)}{\parallel }_{\infty }.$ We can use the tensor product interpolation operators for two-dimensional problems.

2.4. Numerical Inversion of Laplace Transform

In this section, we implement the inverse Laplace transform method to convert the CSCM solution $\widehat{U}(\overline{\mathbf{x}},s)$ from the Laplace domain to the time domain as follows:

$$U(\overline{\mathbf{x}},\tau )=\frac{1}{2\pi i}{\int }_{\rho -i\infty }^{\rho +i\infty }{e}^{s\tau }\widehat{U}(\overline{\mathbf{x}},s)ds,\rho >{\rho }_0.$$

(20)

Here the transform $\widehat{U}(\overline{\mathbf{x}},s)$ needs to be inverted, ${\rho }_0$ is the converging abscissa, and $\rho >{\rho }_0.$ This means that all of the singularities of $\widehat{U}(\overline{\mathbf{x}},s)$ lie in the open half-plane $Re\left(s\right)<\rho .$ We aim to approximate the integral defined in Equation (20). In most cases, it is quite difficult to evaluate the integral defined in Equation (20) analytically; thus, a numerical method must be used. There are several numerical algorithms in the literature that can be used to evaluate the integral defined in Equation (20). Among them, we can list the Fourier series method [40], the de Hoog method [41], Stehfest’s method [35,42], Talbot’s method [32,36], and the Weeks method [43,44], etc. Each approach has an identifiable use and is appropriate for a particular function. In this article, we use two popular inversion algorithms: the improved Talbot’s and Stehfest’s methods will be presented in the following sections.

2.4.1. TheImproved Talbot’s Method

Here, we use Talbot’s method to approximate the $U(\overline{\mathbf{x}},\tau )$

$$U(\overline{\mathbf{x}},\tau )=\frac{1}{2\pi i}{\int }_{\rho -i\infty }^{\rho +i\infty }{e}^{s\tau }\widehat{U}(\overline{\mathbf{x}},s)ds=\frac{1}{2\pi i}{\int }_C{e}^{s\tau }\widehat{U}(\overline{\mathbf{x}},s)ds,Re\left(s\right)>0.$$

(21)

where C is a suitably chosen line. It is very difficult to solve the above integral due to the highly oscillatory exponential function $e^{s\tau }$ and the slow-decaying transform $\widehat{U}(\overline{\mathbf{x}},s)$. Talbot [36] suggested that by utilizing contour integration, this problem can be resolved. He further said that the integration might be carried out so that the real parts of the contours begin and stop on the left side of the complex plane, which contains all of the singularities in the transformed function $\widehat{U}(\overline{\mathbf{x}},s)$. Due to the exponential element, the integrand decays quickly on such contours, making integral Equation (21) suitable for approximation with the midpoint rule. Consider a contour of the form [32]:

$$\Gamma :s=s\left(\eta \right),-\pi \le \eta \le \pi ,$$

(22)

where $Res(\pm \pi )=-\infty$, and $s\left(\eta \right)$ is given as

$$s\left(\eta \right)=\frac{M_T}{\tau }\zeta \left(\eta \right),\zeta \left(\eta \right)=-\delta +\varrho \eta cot\left(\mu \eta \right)+\gamma i\eta ,$$

(23)

where $\mu , \varrho , \mathrm{and} \gamma$ will be chosen precisely. Using Equation (23) in Equation (21), we obtain

$$U(\overline{\mathbf{x}},\tau )=\frac{1}{2\pi i}{\int }_{-\pi }^{\pi }{e}^{s\left(\eta \right)\tau }\widehat{U}(\overline{\mathbf{x}},s\left(\eta \right))s^{\prime }\left(\eta \right)d\eta .$$

(24)

The integral defined in (24) is approximated via the midpoint rule with spacing $h=\frac{2\pi }{M_T}$ as

$$U_{App}(\overline{\mathbf{x}},\tau )\approx \frac{1}{M_{T}i}\sum _{j=1}^{M_T}{e}^{s\left({\eta }_j\right)\tau }\widehat{U}(\overline{\mathbf{x}},s\left({\eta }_j\right))s^{\prime }\left({\eta }_j\right),{\eta }_j=-\pi +(j-\frac{1}{2})h.$$

(25)

Convergence

The approximate value of the defined integral in Equation (24) depends on the contour of integration $\Gamma ;$ different rates of convergence are attained for the suggested numerical scheme. Additionally, the quadrature step h determines that the suggested scheme will converge. In order to have optimal results, we need to search for the optimal contour of integration, which can be found by using the optimal values for the parameters involved in (23). In [32], the authors have proposed optimal values for the parameters as

$$\delta =0.6122,\gamma =0.2645,\mu =0.6407,\mathrm{and}\varrho =0.5017,$$

For the improved Talbot’s method, the error estimate is given as

$$Erro{r}_{est}=|U_{App}(\overline{\mathbf{x}},\tau )-U(\overline{\mathbf{x}},\tau )|=O\left(e^{-1.35800M_T}\right).$$

2.4.2. TheStehfest’s Method

The Gaver–Stehfest method is one of the most important techniques for the numerical inversion of the Laplace transform. It was designed in the second half of the 1960s. It has gained prominence in a number of disciplines, such as chemistry, finance, economics, and computational physics, due to its effectiveness and ease. The series of Gaver approximants, as determined by Gaver [42], is the foundation of the Gaver–Stehfest method. Acceleration was required because the Gaver approximants’ convergence was basically logarithmic. Stehfest [35] used the Salzer acceleration method to offer a linear acceleration method. The Gaver–Stehfest approach uses a series of functions to approximate $U(\overline{\mathbf{x}},\tau )$ as

$$U_{App}(\overline{\mathbf{x}},\tau )=\frac{ln2}{\tau }\sum _{j=1}^{M_S}{\alpha }_j\widehat{U}\left(\overline{\mathbf{x}},\frac{ln2}{\tau }j\right),$$

(26)

where the coefficients ${\alpha }_j$ are defined by

$${\alpha }_j={(-1)}^{\frac{M_S}{2}+j}\sum _{\ell =\lfloor \frac{j+1}{2}\rfloor }^{min(j,\frac{M_S}{2})}\frac{{\ell }^{\frac{M_S}{2}}(2\ell )!}{(\frac{M_S}{2}-\ell )!\ell !(\ell -1)!(j-\ell )!(2\ell -j)!}.$$

(27)

Solve Equation (6) for the corresponding Laplace parameters $s=\frac{ln2}{\tau }j,j=1,2,3,\dots ,$$M_S.$ The solution to the given problem in Equation (1) can be obtained by (26). There are a few appealing features of the Gaver–Stehfest algorithm: (i) $U_{App}(\overline{\mathbf{x}},\tau )$ are linear in terms of the values of $U(\overline{\mathbf{x}},s);$ (ii) the values of $U(\overline{\mathbf{x}},s)$ are needed only for the real value of s; (iii) the process of computing the coefficients is quite simple; and (iv) for constant functions, this approach yields highly accurate approximations. In the literature, this technique has been used by many authors in [45,46], where it has been demonstrated that this method converges very fast to $U_{App}(\overline{\mathbf{x}},\tau )$ (provided $U_{App}(\overline{\mathbf{x}},\tau )$ is non-oscillatory).

Convergence

In [45], the author has derived two conditions for $U(\overline{\mathbf{x}},\tau ),$ which guarantee the convergence of $U_{App}(\overline{\mathbf{x}},\tau ).$ The conditions are given in the following theorem.

Theorem 2.

Let $U:(0,\infty )\to \mathbb{R}$ be a locally integrable function and the approximate solution $U_{App}(\overline{\mathbf{x}},\tau )$ be defined by (26), and let its LT $\widehat{U}(\overline{\mathbf{x}},s)$ be defined for $s>0:$

1. 

$U_{App}(\overline{\mathbf{x}},\tau )$ converges for the values of $U_{App}(\overline{\mathbf{x}},\tau )$ in the neighborhood of τ.

2. 

Let for some real number m and some $0<\epsilon <1/4$

$${\int }_0^{\epsilon }|U(-\tau lo{g}_2(1/2+\xi ))+U(-\tau lo{g}_2(1/2-\xi ))-2m|{\xi }^{-1}d\xi <\infty .$$

Then, $U\to m$ as $M_S\to +\infty .$

3. 

Allow that $U_{App}(\overline{\mathbf{x}},\tau )$ has a bounded variation in the neighborhood of $\tau .$ Then,

$$U_{App}\left(\tau \right)\to \frac{U(\tau +0)+U(\tau -0)}{2}\mathrm{as}{M}_S\to +\infty .$$

Corollary 1.

Under the above assumptions, if

$$U(\tau +\xi )-U\left(\tau \right)=O\left({\left|\xi \right|}^{\upsilon }\right),$$

$\forall \xi$ and some υ in the neighborhood of $\tau ,$ then $U_{App}(\overline{\mathbf{x}},\tau )\to U(\overline{\mathbf{x}},\tau ),$ as $M_S\to +\infty .$

Additionally, the authors in [47] performed various experiments for the parameter effect on the accuracy of the numerical scheme, and they concluded in their findings that “If p significant digits are required, then let $M_S$ be the positive integer $\lceil 2.2p\rceil$. Set the system precision to $q=\lceil 1.1M_S\rceil$ and, for a given $M_S$, calculate ${\alpha }_i,1\le i\le M_S$, using (27). Then, for the given transform $\widehat{U}(\overline{\mathbf{x}},s)$ and the argument $\tau$, calculate the $U_{App}(\overline{\mathbf{x}},\tau )$ in (26)”. According to these conclusions, if the error of input data is ${10}^{-(q+1)}\le \frac{\parallel \widehat{U}-U\parallel }{\parallel U\parallel }\le {10}^{-q}$, with an even positive integer $M_S$ via $q=\lceil 1.1M_S\rceil ,$ then the final error is ${10}^{-(p+1)}\le \frac{\parallel \widehat{U}-U\parallel }{\parallel U\parallel }\le {10}^{-p},$ where $M_S=\lceil 2.2p\rceil$ [48]. Next, we present Algorithm 1 for the proposed numerical scheme.

Algorithm 1: Algorithm for the proposed numerical scheme

1: Input: The computational domain, the fractional order derivative, the final time, the optimal parameters for the inverse laplace transform methods, the inhomogeneous function, and other conditions.  2: Step i: Apply the Laplace transform to problems (1) and (2), and obtain the time-independent problems (6) and (7).   3: Step ii: Obtain ${\mathbf{L}}_{\mathbf{D}}$ by discretizing the linear differential operator $\mathbf{L}$ using the CSCM via (12).   4: Step iii: Solve Equation (13) with the boundary conditions given in Equation (7) for each point in the LT domain.   5: Step iv: Compute the approximate solution to the problem in Equations (1)–(3) using Equation (21) or Equation (26).   6: Step v: The approximate solution is $U_{App}(\overline{\mathbf{x}},\tau ).$

3. Numerical Results and Discussion

The numerical results and discussions are addressed in this section. Three examples are used to evaluate and validate the effectiveness and accuracy of the Laplace-transformed CSCM. We performed our experiments in MATLAB R2019a on a Windows 10 (64-bit) PC equipped with an Intel(R) Core TM i5-4310M CPU @ 2.70 GHz with 4 GB of RAM. Three error norms—the relative $L_2$ error, the maximum absolute error, and the RMS error—are used to evaluate the accuracy of the proposed method and are defined below:

$${er}_2=\sqrt{\frac{{\sum }_{j=1}^n{(U({\overline{\mathbf{x}}}_j,\tau )-U_{App}({\overline{\mathbf{x}}}_j,\tau ))}^2}{{\sum }_{j=1}^n{\left(U({\overline{\mathbf{x}}}_j,\tau )\right)}^2}},$$

$${er}_{\infty }=\underset{1\le j\le n}{max}|U({\overline{\mathbf{x}}}_j,\tau )-U_{App}({\overline{\mathbf{x}}}_j,\tau )|,$$

$${er}_{rms}=\sqrt{\frac{{\sum }_{j=1}^n{(U({\overline{\mathbf{x}}}_j,\tau )-U_{App}({\overline{\mathbf{x}}}_j,\tau ))}^2}{n}},$$

where $U({\overline{\mathbf{x}}}_j,\tau )$ and $U_{App}({\overline{\mathbf{x}}}_j,\tau )$ represent the analytic and numerical solutions, respectively.

Example 1.

Consider the 2D time-fractional advection–diffusion Equation (1) with an exact solution given as

$$U(x,y,\tau )=x^2+y^2+{\tau }^2,\Omega ={[-1,1]}^2,\tau \in [0,1],$$

where ${\lambda }_1=1,\overrightarrow{\mathbf{u}}={(1,1)}^T,$ and ${\lambda }_2=0.$ The forcing term $\mathcal{Q}(x,y,\tau )$, and the initial-boundary data are extracted from the exact solution. The numerical results obtained via Talbot’s and Stehfest’s methods are presented in

Table 1. The $e{r}_2,e{r}_{\infty },\mathrm{and}e{r}_{rms}$ obtained using Stehfest’s method for different values of $\alpha ,n,\mathrm{and}{M}_S$ at $\tau =1$ for Example 1.

$\mathit{\alpha }=0.30$	n	${\mathit{M}}_{\mathit{S}}$	${\mathit{er}}_2$	${\mathit{er}}_{\mathit{\infty }}$	${\mathit{er}}_{\mathit{rms}}$	CPU(s)
	21	12	9.7994 $\times {10}^{-4}$	6.0556 $\times {10}^{-5}$	4.6664 $\times {10}^{-5}$	0.233261
	23		1.1084 $\times {10}^{-3}$	8.2038 $\times {10}^{-5}$	4.8190 $\times {10}^{-5}$	0.200611
	25		1.1581 $\times {10}^{-3}$	7.1935 $\times {10}^{-5}$	4.6325 $\times {10}^{-5}$	0.266045
	25	08	1.0812$\times {10}^{-1}$	4.1584 $\times {10}^{-3}$	4.3246 $\times {10}^{-3}$	0.196162
		10	1.4553 $\times {10}^{-3}$	5.8142 $\times {10}^{-5}$	5.8213 $\times {10}^{-5}$	0.249590
		12	1.1581 $\times {10}^{-3}$	7.1935 $\times {10}^{-5}$	4.6325 $\times {10}^{-5}$	0.255773
$\alpha =0.50$
	21	12	9.8850 $\times {10}^{-4}$	5.6517 $\times {10}^{-5}$	4.7071 $\times {10}^{-5}$	0.167258
	23		1.0923 $\times {10}^{-3}$	7.1906 $\times {10}^{-5}$	4.7492 $\times {10}^{-5}$	0.204245
	25		1.2186 $\times {10}^{-3}$	1.3279 $\times {10}^{-4}$	4.8744 $\times {10}^{-5}$	0.302795
	25	08	1.0812$\times {10}^{-1}$	4.1584 $\times {10}^{-4}$	4.3246 $\times {10}^{-3}$	0.196037
		10	1.4539 $\times {10}^{-3}$	5.7437 $\times {10}^{-5}$	5.8156 $\times {10}^{-5}$	0.241286
		12	1.2186 $\times {10}^{-3}$	1.3279 $\times {10}^{-4}$	4.8744 $\times {10}^{-5}$	0.377014
$\alpha =0.80$
	21	12	1.0540 $\times {10}^{-3}$	9.2987 $\times {10}^{-5}$	5.0191 $\times {10}^{-5}$	0.142862
	23		9.8247 $\times {10}^{-4}$	6.2113 $\times {10}^{-5}$	4.2716 $\times {10}^{-5}$	0.235252
	25		1.2741 $\times {10}^{-3}$	1.5227 $\times {10}^{-4}$	5.0962 $\times {10}^{-5}$	0.287223
	25	08	1.0811$\times {10}^{-1}$	4.1585 $\times {10}^{-3}$	4.3246 $\times {10}^{-5}$	0.214762
		10	1.4514 $\times {10}^{-3}$	5.7719 $\times {10}^{-5}$	5.8056 $\times {10}^{-5}$	0.234821
		12	1.2741 $\times {10}^{-3}$	1.5227 $\times {10}^{-4}$	5.0962 $\times {10}^{-5}$	0.287223

and

Table 2. The $e{r}_2,e{r}_{\infty },\mathrm{and}e{r}_{rms}$ obtained via Talbot’s method for different values of $\alpha ,n$, and $M_T$ at $\tau =1$ for Example 1.

$\mathit{\alpha }=0.30$	n	${\mathit{M}}_{\mathit{T}}$	${\mathit{er}}_2$	${\mathit{er}}_{\mathit{\infty }}$	${\mathit{er}}_{\mathit{rms}}$	CPU(s)
	21	26	3.6788 $\times {10}^{-10}$	9.4858 $\times {10}^{-11}$	1.7518 $\times {10}^{-11}$	0.477112
	23		5.1649 $\times {10}^{-10}$	1.5900 $\times {10}^{-10}$	2.2456 $\times {10}^{-11}$	0.710826
	25		8.1265 $\times {10}^{-10}$	3.4762 $\times {10}^{-10}$	3.2506 $\times {10}^{-11}$	1.118694
	25	22	3.6227 $\times {10}^{-8}$	1.5662 $\times {10}^{-9}$	1.4491 $\times {10}^{-9}$	0.955487
		24	3.0969 $\times {10}^{-9}$	3.1978 $\times {10}^{-10}$	1.2388 $\times {10}^{-10}$	0.997716
		26	8.1265 $\times {10}^{-10}$	3.4762 $\times {10}^{-10}$	3.2506 $\times {10}^{-11}$	1.082287
$\alpha =0.50$
	21	26	4.5560 $\times {10}^{-10}$	1.7228 $\times {10}^{-10}$	2.1695 $\times {10}^{-11}$	0.468985
	23		5.0871 $\times {10}^{-10}$	1.8159 $\times {10}^{-10}$	2.2118 $\times {10}^{-11}$	0.741330
	25		8.4173 $\times {10}^{-10}$	3.3289 $\times {10}^{-10}$	3.3669 $\times {10}^{-11}$	1.042339
	25	22	3.6042 $\times {10}^{-8}$	1.4790 $\times {10}^{-9}$	1.4417 $\times {10}^{-9}$	0.952621
		24	2.8473 $\times {10}^{-9}$	3.0892 $\times {10}^{-10}$	1.1389 $\times {10}^{-10}$	0.987540
		26	8.4173 $\times {10}^{-10}$	3.3289 $\times {10}^{-10}$	3.3669 $\times {10}^{-11}$	1.042339
$\alpha =0.80$
	21	26	3.9392 $\times {10}^{-10}$	8.8192 $\times {10}^{-11}$	1.8758 $\times {10}^{-11}$	0.479675
	23		4.8220 $\times {10}^{-10}$	1.9319 $\times {10}^{-10}$	2.0965 $\times {10}^{-11}$	0.761089
	25		1.1098 $\times {10}^{-9}$	7.7091 $\times {10}^{-10}$	4.4393 $\times {10}^{-11}$	1.094016
	25	22	3.6116 $\times {10}^{-8}$	1.5182 $\times {10}^{-9}$	1.4446 $\times {10}^{-9}$	0.917780
		24	2.9420 $\times {10}^{-9}$	2.3727 $\times {10}^{-10}$	1.1768 $\times {10}^{-10}$	1.010576
		26	1.1098 $\times {10}^{-9}$	7.7091 $\times {10}^{-10}$	4.4393 $\times {10}^{-11}$	1.094016

by using different values of $\alpha ,M_S,M_T,$ and n at $\tau =1.$ The approximate solution obtained by the proposed method is shown in Figure 1a. The plots for the comparison of errors $e{r}_2,e{r}_{\infty }$, and $e{r}_{rms}$ using Stehfest’s method for several values of $\alpha ,M_S$, and τ at $n=25$ are shown in Figure 1b–d, respectively. Figure 2a–c show the plots of errors $e{r}_2,e{r}_{\infty }$, and $e{r}_{rms}$ obtained using Talbot’s method with different values of $\alpha ,M_T$, and τ at $n=25.$ Figure 2d presents the contour plot of the absolute error using Talbot’s method for different values of α with $M_T=26,\tau =1$, and $n=25.$ The accuracy of Stehfest’s method enhances for $M_S=8,10,12,14$ and then gradually diminishes for values exceeding $M_S=14$. From all of the obtained results presented in the figures and tables, we conclude that the proposed method is stable, accurate, and efficient. Clearly, the numerical results show that the accuracy of Talbot’s method is better than that of Stehfest’s method.

Example 2.

Consider the 2D time-fractional advection–diffusion Equation (1) with an exact solution given as

$$U(x,y,\tau )=sin(x+y)+{\tau }^2,\Omega ={[-1,1]}^2,\tau \in [0,1],$$

where ${\lambda }_1=1,\overrightarrow{\mathbf{u}}={(1,1)}^T,$ and ${\lambda }_2=0.$ The forcing term $\mathcal{Q}(x,y,\tau ),$ and the initial-boundary data are extracted from the exact solution. The numerical results obtained via Talbot’s and Stehfest’s methods are presented in

Table 3. The $e{r}_2,e{r}_{\infty },\mathrm{and}e{r}_{rms}$ obtained via Stehfest’s method for different values of $\alpha ,n$, and $M_S$ at $\tau =1$ of Example 2.

$\mathit{\alpha }=0.30$	n	${\mathit{M}}_{\mathit{S}}$	${\mathit{er}}_2$	${\mathit{er}}_{\mathit{\infty }}$	${\mathit{er}}_{\mathit{rms}}$	CPU(s)
	21	12	9.9430 $\times {10}^{-4}$	5.2055 $\times {10}^{-5}$	4.7348 $\times {10}^{-5}$	0.192349
	23		1.0940 $\times {10}^{-3}$	6.1845 $\times {10}^{-5}$	4.7567 $\times {10}^{-5}$	0.333385
	25		1.1758 $\times {10}^{-3}$	6.1712 $\times {10}^{-5}$	4.7031 $\times {10}^{-5}$	0.303723
	25	08	1.0812$\times {10}^{-1}$	4.1584 $\times {10}^{-3}$	4.3246 $\times {10}^{-3}$	0.271822
		10	1.4554 $\times {10}^{-3}$	5.8032 $\times {10}^{-5}$	5.8217 $\times {10}^{-5}$	0.257647
		12	1.1758 $\times {10}^{-3}$	6.1712 $\times {10}^{-5}$	4.7031 $\times {10}^{-5}$	0.303723
$\alpha =0.50$
	21	12	9.9254 $\times {10}^{-4}$	5.2330 $\times {10}^{-5}$	4.7264 $\times {10}^{-5}$	0.170534
	23		1.0888 $\times {10}^{-3}$	5.9124 $\times {10}^{-5}$	4.7338 $\times {10}^{-5}$	0.257113
	25		1.1700 $\times {10}^{-3}$	5.4659 $\times {10}^{-5}$	4.6800 $\times {10}^{-5}$	0.360627
	25	08	1.0812$\times {10}^{-1}$	4.1583 $\times {10}^{-3}$	4.3246 $\times {10}^{-3}$	0.197690
		10	1.4551 $\times {10}^{-3}$	5.6736 $\times {10}^{-5}$	5.8205 $\times {10}^{-5}$	0.383466
		12	1.1700 $\times {10}^{-3}$	5.4659 $\times {10}^{-5}$	4.6800 $\times {10}^{-5}$	0.360627
$\alpha =0.80$
	21	12	9.9170 $\times {10}^{-4}$	5.0410 $\times {10}^{-5}$	4.7224 $\times {10}^{-5}$	0.178846
	23		1.0836 $\times {10}^{-3}$	5.0672 $\times {10}^{-5}$	4.7115 $\times {10}^{-5}$	0.276298
	25		1.2056 $\times {10}^{-3}$	6.7927 $\times {10}^{-5}$	4.8224 $\times {10}^{-5}$	0.292158
	25	08	1.0812$\times {10}^{-1}$	4.1583 $\times {10}^{-3}$	4.3246 $\times {10}^{-5}$	0.274343
		10	1.4550 $\times {10}^{-3}$	5.6956 $\times {10}^{-5}$	5.8199 $\times {10}^{-5}$	0.269435
		12	1.2056 $\times {10}^{-3}$	6.7927 $\times {10}^{-5}$	4.8224 $\times {10}^{-5}$	0.292158

and

Table 4. The $e{r}_2,e{r}_{\infty },\mathrm{and}e{r}_{rms}$ obtained using Talbot’s method for different values of $\alpha ,n$, and $M_T$ at $\tau =1$ of Example 2.

$\mathit{\alpha }=0.30$	n	${\mathit{M}}_{\mathit{T}}$	${\mathit{er}}_2$	${\mathit{er}}_{\mathit{\infty }}$	${\mathit{er}}_{\mathit{rms}}$	CPU(s)
	21	26	2.3790 $\times {10}^{-10}$	6.3914 $\times {10}^{-11}$	1.1328 $\times {10}^{-11}$	0.508399
	23		2.6664 $\times {10}^{-10}$	5.2105 $\times {10}^{-11}$	1.1593 $\times {10}^{-11}$	0.829189
	25		5.6968 $\times {10}^{-10}$	3.8928 $\times {10}^{-10}$	2.2787 $\times {10}^{-11}$	1.151650
	25	22	3.6179 $\times {10}^{-8}$	1.4491 $\times {10}^{-9}$	1.4472 $\times {10}^{-9}$	0.977947
		24	2.9202 $\times {10}^{-9}$	2.3436 $\times {10}^{-10}$	1.1681 $\times {10}^{-10}$	1.044065
		26	5.6968 $\times {10}^{-10}$	3.8928 $\times {10}^{-10}$	2.2787 $\times {10}^{-11}$	1.151650
$\alpha =0.50$
	21	26	2.7669 $\times {10}^{-10}$	7.5912 $\times {10}^{-11}$	1.3176 $\times {10}^{-11}$	0.545475
	23		2.7187 $\times {10}^{-10}$	6.4024 $\times {10}^{-11}$	1.1820 $\times {10}^{-11}$	0.806244
	25		4.6716 $\times {10}^{-10}$	2.0896 $\times {10}^{-10}$	1.8687 $\times {10}^{-11}$	1.225075
	25	22	3.6229 $\times {10}^{-8}$	1.4880 $\times {10}^{-9}$	1.4492 $\times {10}^{-9}$	0.939474
		24	2.8311 $\times {10}^{-9}$	1.5034 $\times {10}^{-10}$	1.1324 $\times {10}^{-10}$	1.027609
		26	4.6716 $\times {10}^{-10}$	2.0896 $\times {10}^{-10}$	1.8687 $\times {10}^{-11}$	1.225075
$\mathit{\alpha }=0.80$
	21	26	2.4234 $\times {10}^{-10}$	4.0678 $\times {10}^{-11}$	1.1540 $\times {10}^{-11}$	0.512461
	23		2.4004 $\times {10}^{-10}$	7.5801 $\times {10}^{-11}$	1.0437 $\times {10}^{-11}$	0.820029
	25		4.3673 $\times {10}^{-10}$	2.3799 $\times {10}^{-10}$	1.7469 $\times {10}^{-11}$	1.162087
	25	22	3.6188 $\times {10}^{-8}$	1.4465 $\times {10}^{-9}$	1.4475 $\times {10}^{-9}$	0.988509
		24	2.8035 $\times {10}^{-9}$	1.4878 $\times {10}^{-10}$	1.1214 $\times {10}^{-10}$	1.047553
		26	4.3673 $\times {10}^{-9}$	2.3799 $\times {10}^{-10}$	1.7469 $\times {10}^{-11}$	1.162087

by using different values for $\alpha ,M_S,M_T,$ and n at $\tau =1.$ The approximate solution obtained by the proposed method is shown in Figure 3a. The plots of the comparison of errors $e{r}_2,e{r}_{\infty }$, and $e{r}_{rms}$ using Stehfest’s method for several values of $\alpha ,M_S$, and τ at $n=25$ are shown in Figure 3b–d, respectively. Figure 4a–c show the plots of errors $e{r}_2,e{r}_{\infty }$, and $e{r}_{rms}$ by using Talbot’s method with different values of $\alpha ,M_T$, and τ at $n=25.$ Figure 3e presents the contour plot of the absolute error obtained by using Stehfest’s method for different values of α with $M_S=12,\tau =1$, and $n=25.$ The accuracy of Stehfest’s method enhances for $M_S=8,10,12,14$ and then gradually diminishes for values exceeding $M_S=14$. From all of the obtained results presented in the figures and tables, we conclude that the proposed method is stable, accurate, and efficient. Clearly, the results show that the accuracy of Talbot’s method is better than that of Stehfest’s method.

Example 3.

Consider the 2D time-fractional advection–diffusion Equation (1) with an exact solution given as

$$U(x,y,\tau )={\tau }^2{e}^{(x+y)},\Omega ={[-1,1]}^2,\tau \in [0,1],$$

where ${\lambda }_1=1,\overrightarrow{\mathbf{u}}={(1,1)}^T,$ and ${\lambda }_2=0.$ The forcing term $\mathcal{Q}(x,y,\tau ),$ and the initial-boundary data are extracted from the exact solution. The numerical results obtained via Talbot’s and Stehfest’s methods are presented in

Table 5. The $e{r}_2,e{r}_{\infty },\mathrm{and}e{r}_{rms}$ obtained via Stehfest’s method for different values of $\alpha ,n$, and $M_S$ at $\tau =1$ of Example 3.

$\mathit{\alpha }=0.30$	n	${\mathit{M}}_{\mathit{S}}$	${\mathit{er}}_2$	${\mathit{er}}_{\mathit{\infty }}$	${\mathit{er}}_{\mathit{rms}}$	CPU(s)
	21	14	1.9041 $\times {10}^{-4}$	3.5446 $\times {10}^{-5}$	9.0670$\times {10}^{-6}$	0.142019
	23		1.7304 $\times {10}^{-4}$	2.7673 $\times {10}^{-5}$	7.5234$\times {10}^{-6}$	0.214226
	25		2.2300 $\times {10}^{-4}$	4.3858 $\times {10}^{-5}$	8.9200$\times {10}^{-6}$	0.278411
	25	10	3.3997 $\times {10}^{-3}$	4.1349 $\times {10}^{-4}$	1.3599 $\times {10}^{-4}$	0.234626
		12	2.7481 $\times {10}^{-3}$	3.3443 $\times {10}^{-4}$	1.0992 $\times {10}^{-4}$	0.275799
		14	2.2300 $\times {10}^{-4}$	4.3858 $\times {10}^{-5}$	8.9200$\times {10}^{-6}$	0.334381
$\alpha =0.50$
	21	14	1.7886 $\times {10}^{-4}$	3.4616 $\times {10}^{-5}$	8.5172$\times {10}^{-6}$	0.173237
	23		1.9875 $\times {10}^{-4}$	4.7601 $\times {10}^{-5}$	8.6413$\times {10}^{-6}$	0.238699
	25		1.8379 $\times {10}^{-4}$	2.8417 $\times {10}^{-5}$	7.3515$\times {10}^{-6}$	0.321624
	25	10	3.3995 $\times {10}^{-3}$	4.1349 $\times {10}^{-4}$	1.3598 $\times {10}^{-4}$	0.246361
		12	2.7522 $\times {10}^{-3}$	3.3443 $\times {10}^{-4}$	1.1009 $\times {10}^{-4}$	0.256780
		14	1.8379 $\times {10}^{-4}$	2.8417 $\times {10}^{-5}$	7.3515$\times {10}^{-6}$	0.320022
$\alpha =0.80$
	21	14	2.1292 $\times {10}^{-4}$	2.4907 $\times {10}^{-5}$	1.0139 $\times {10}^{-5}$	0.160703
	23		2.8307 $\times {10}^{-4}$	4.5577 $\times {10}^{-5}$	1.2307 $\times {10}^{-5}$	0.228405
	25		4.0222 $\times {10}^{-4}$	4.9547 $\times {10}^{-5}$	1.6089 $\times {10}^{-5}$	0.323644
	25	10	3.3995 $\times {10}^{-3}$	4.1349 $\times {10}^{-4}$	1.3598 $\times {10}^{-4}$	0.253768
		12	2.7525 $\times {10}^{-3}$	3.3443 $\times {10}^{-4}$	1.1010 $\times {10}^{-5}$	0.318281
		14	4.0222 $\times {10}^{-4}$	4.9547 $\times {10}^{-5}$	1.6089 $\times {10}^{-5}$	0.323644

and

Table 6. The $e{r}_2,e{r}_{\infty },\mathrm{and}e{r}_{rms}$ obtained using Talbot’s method for different values of $\alpha ,n$, and $M_T$ at $\tau =1$ for Example 3.

$\mathit{\alpha }=0.30$	n	${\mathit{M}}_{\mathit{S}}$	${\mathit{er}}_2$	${\mathit{er}}_{\mathit{\infty }}$	${\mathit{er}}_{\mathit{rms}}$	CPU(s)
	21	30	2.4416 $\times {10}^{-11}$	8.9271 $\times {10}^{-12}$	1.1627 $\times {10}^{-12}$	0.623102
	23		4.5926 $\times {10}^{-11}$	2.4259 $\times {10}^{-11}$	1.9968 $\times {10}^{-12}$	0.893166
	25		1.0575 $\times {10}^{-10}$	3.6040 $\times {10}^{-11}$	4.2300 $\times {10}^{-12}$	1.424708
	25	26	5.3511 $\times {10}^{-10}$	8.2624 $\times {10}^{-11}$	2.1404 $\times {10}^{-11}$	1.154428
		28	7.2756 $\times {10}^{-11}$	3.0986 $\times {10}^{-11}$	2.9102 $\times {10}^{-12}$	1.290241
		30	1.0575 $\times {10}^{-10}$	3.6040 $\times {10}^{-11}$	4.2300r $\times {10}^{-12}$	1.304660
$\alpha =0.50$
	21	30	2.4307 $\times {10}^{-11}$	8.8405 $\times {10}^{-12}$	1.1575 $\times {10}^{-12}$	0.590850
	23		4.5760 $\times {10}^{-11}$	1.3785 $\times {10}^{-11}$	1.9896 $\times {10}^{-12}$	0.895766
	25		6.7767 $\times {10}^{-11}$	3.3908 $\times {10}^{-11}$	2.7107 $\times {10}^{-12}$	1.342481
	25	26	5.0284 $\times {10}^{-10}$	7.5070 $\times {10}^{-11}$	2.0114 $\times {10}^{-11}$	1.118816
		28	8.2316 $\times {10}^{-11}$	3.8009 $\times {10}^{-11}$	3.2926 $\times {10}^{-12}$	1.221575
		30	6.7767 $\times {10}^{-11}$	3.3908 $\times {10}^{-11}$	2.7107 $\times {10}^{-12}$	1.282603
$\alpha =0.80$
	21	30	2.6976 $\times {10}^{-11}$	8.4350 $\times {10}^{-12}$	1.2846 $\times {10}^{-12}$	0.569593
	23		4.9832 $\times {10}^{-11}$	1.3129 $\times {10}^{-11}$	2.1666 $\times {10}^{-12}$	0.901900
	25		8.9002 $\times {10}^{-11}$	5.0524 $\times {10}^{-11}$	3.5601 $\times {10}^{-12}$	1.338635
	25	26	5.0864 $\times {10}^{-10}$	7.6800 $\times {10}^{-11}$	2.0346 $\times {10}^{-11}$	1.152948
		28	6.4788 $\times {10}^{-11}$	2.4306 $\times {10}^{-11}$	2.5915 $\times {10}^{-12}$	1.197204
		30	8.9002 $\times {10}^{-11}$	5.0524 $\times {10}^{-11}$	3.5601 $\times {10}^{-12}$	1.295479

by using different values of $\alpha ,M_S,and{M}_T,$ with n at $\tau =1.$ The approximate solution obtained by the proposed method is shown in Figure 5a. The plots of the comparison of the errors $e{r}_2,e{r}_{\infty }$, and $e{r}_{rms}$ using Stehfest’s method for several values of $\alpha ,M_S$, and τ at $n=25$ are shown in Figure 5b–d, respectively. Figure 6a–c show the plots of errors $e{r}_2,e{r}_{\infty }$, and $e{r}_{rms}$ via Talbot’s method with different values of $\alpha ,M_T$, and τ at $n=25.$ Figure 6d presents the contour plot of the absolute error by using Talbot’s method for different values of α with $M_T=26,\tau =1$, and $n=25.$ The accuracy of Stehfest’s method enhances for $M_S=8,10,12,14$ and then gradually diminishes for values exceeding than $M_S=14$. From all the obtained results presented in the figures and tables, we conclude that the proposed method is stable, accurate, and efficient. Clearly, the results show that the accuracy of Talbot’s method is better than that of Stehfest’s method.

4. Conclusions

In this article, a numerical method based on the LT coupled with the CSCM for the numerical solution of two-dimensional TFADEs in the Caputo sense has been developed. It first uses the LT to convert a TFADE into a time-independent inhomogeneous equation in the LT space. Then, it uses the CSCM to discretize the spatial derivatives of this Laplace-transformed inhomogeneous equation. Finally, it chooses Stehfest’s method and Talbot’s method for the numerical inversion of the LT to retrieve the numerical solutions of the TFADE from the corresponding CSCM solutions in the LT domain. Compared with the finite difference time-stepping scheme, the proposed method introduces the LT technique to avoid the calculation of a costly convolution integral in the approximation of time-fractional derivation and to avoid the effect of the time step on the stability and accuracy. Hence, it can successfully solve the long time history TFADEs. The proposed advection–diffusion equation of fractional order presents significant potential for extension into uncertainty propagation and sensitivity analysis, enabling a more comprehensive and nuanced understanding of system responses to parameter variations and initial conditions [49,50,51].