1 Introduction

The calculus of derivatives and integrals of non-integer order (that is fractional order) is not a new concept. Its footprints of inspiration span three decades and has multidimensional patterns of its practical implementations in science and engineering. This connection to other sciences makes fractional calculus more vital and indispensable. The fractional calculus has a long track of exertion in various theories and their applications. For example, it has been successfully used in viscoplastic and viscoelastic flow [1], transport problems [2], control theory [3], tumor development [4], continuum mechanics [5], random walks [6] and turbulence [7, 8]. In contrary to the conventional integer order counterparts, the fractional models which are described by fractional differential or integral equations possess more reliable and certain characteristics and thus behave more appropriately. Consequently, there is a growing need to explore solution techniques to study these models. More often, the analytical solution of most of the fractional differential equations cannot be acquired. Therefore, the search for approximate and analytical approaches is a subject matter of many lately published research papers. In this paper, approximate solution of the following time fractional advection diffusion equation (TFADE) shall be studied:

$$\begin{aligned} \frac{\partial ^{\gamma }w(r,t)}{\partial t^{\gamma }}&=\Phi \frac{\partial ^{2}w(r,t)}{\partial r^2}\nonumber \\&\quad -\Psi \frac{\partial w(r,t)}{\partial r}+ q(r,t) , \ \ r\in [a,b],\ t\in [t_0,T], \end{aligned}$$
(1.1)

with initial condition:

$$\begin{aligned} w(r,t_0)=\phi (r) \end{aligned}$$
(1.2)

and the boundary conditions:

$$\begin{aligned} w(a,t)=\psi _1(t),\ \ w(b,t)=\psi _2(t), \end{aligned}$$
(1.3)

where \(\gamma \in (0,1)\) and q(r, t) is the source term. \(\Phi >0\) and \(\Psi \) are the diffusion coefficient and average velocity, respectively. \(\frac{\partial ^{\gamma }}{\partial t^{\gamma }}w(r,t)\) is the Atangana–Baleanu fractional derivative (ABFD) based on non-singular kernel. The ABFD contains the Mittag–Leffler function that is the cornerstone of the fractional calculus. Applications of the Atangana–Baleanu operator have been explored in fields as diverse as heat transfer [9], chaos theory [10], oscillators [11], the hepatitis E virus model [12], cancer chemotherapy [13], smoking model [14] and coronavirus [15].

As noted in the literature, significant work has been done on numerical and analytical solutions of the TFADE. Mardani et al. [16] presented a meshless scheme based on the moving least square approximation for solving the TFADE with variable coefficients involving Caputo derivative. Bu et al. [17] proposed a V-cycle multigrid method for multi-term TFADEs. Sarboland [18] employed a meshfree approach for approximate solution of the time fractional partial differential equations based on multiquadric quasi-interpolation operator. The polynomial spectral collocation scheme has been constructed by Tian et al. [19] for solving the space fractional advection diffusion equation (FADE). Zheng et al. [20] studied the space FADE by means of a finite element method. An explicit and an implicit difference approximations have been developed by Shen et al. [21] for the solution of space-time Riesz–Caputo FADE. Azin et al. [22] established a hybrid method based on the modified Legendre and the Chebyshev cardinal functions for solving TFADE. The analytical and numerical investigations of two-dimensional TFADE have been carried out by Ahmed et al. [23] involving memory and a concentrated source.

Mirza and Vieru [24] obtained the fundamental solutions of TFADE by employing the Laplace and Fourier transforms. Baleanu et al. [25] presented the variational homotopic perturbation and q-homotopic analysis methods for solving the fractional advection differential equation involving Caputo and Caputo–Fabrizio derivatives. Analytical solutions for FADE have been obtained by Rubbab et al. [26] utilizing the integral transforms technique based on time-dependent pulses on the boundary. Rubbab et al. [27] developed the finite difference and pseudo-spectral collocation schemes and analysed the numerical simulation of FADE with Caputo–Fabrizio operator in cylindrical domains. Korpinar et al. [28] presented the fractional homotopy perturbation transform scheme for the time fractional Gardner equation involving ABFD with Mittag–Leffler kernel. Owolabi [29, 30] proposed computational schemes for time fractional Kuramoto–Sivashinsky equation and multi-component systems with ABFD, respectively. A computational algorithm based on homotopic method has been employed by Kumar et al. [31] to investigate the fractional vibration equation. Hosseininia and Heydari [32] established a meshfree scheme by utilizing the moving least squares shape functions for the numerical solution of nonlinear 2D telegraph equation with variable-order ABFD. Inc et al. [33] analysed the logarithmic-KdV equation with Atangana–Baleanu operator. Bas and Ozarslan [34] obtained analytical solutions of fractional models (population growth, blood alcohol model, logistic equation, Newton’s law of cooling) involving ABFD using Laplace transform. Akgül [35] presented a novel technique for investigation of fractional differential equations involving the ABFD. Akgül and Modanli [36] proposed the Crank–Nicholson difference scheme method for the third order fractional partial differential equation with ABFD. Attia et al. [37] developed the computational approach to find an approximate solution of the TFADE.

In recent years, several authors worked on numerous numerical procedures using B-spline interpolation for solving fractional differential equations. In particular, due to rich geometrical features of the spline functions, they are found to be a powerful tool in curve approximation. Yaseen et al. [38] presented a computational method for the generalized nonlinear time-fractional Klein–Gordon equation involving Caputo operator via cubic trigonometric B-spline functions. Abbas et al. [39] proposed the new cubic B-spline approximation based method for numerical solution of non-linear third order Korteweg-de Vries equation. A redefined cubic B-spline functions based scheme has been developed by Khalid et al. [40] for numerical solution of time fractional Allen–Cahn equation. The authors in [41, 42] derived numerical algorithms from the extended cubic B-spline functions to investigate the computational solutions for the time-fractional fisher and telegraph equations. Iqbal et al. [43] investigated numerical simulation of Kuramoto–Sivashinsky equation using new quintic polynomial B-spline functions. Khalid et al. [44] presented a computational method using redefined extended B-spline approximation for solving TFADE.

This paper presents the numerical scheme for the TFADE based on CBS functions. Current proposed method used the \(\theta \)-weighted scheme and Atangana–Baleanu operator. The usage of non-singular kernel operator in B-spline methods is novel. The stability and convergence analysis of the scheme is also performed to avoid any false result. To the best of the author’s knowledge, the proposed algorithm for TFADE is novel and it has not been studied yet.

The current paper is structured in the following pattern: In Sect. 2, the ABFD, Parseval’s identity and basis functions are defined and currently developed scheme is presented in Sect. 3. Sections 4 and 5 contain the stability of algorithm and analysis of convergence, respectively. The validity and efficiency of the proposed method is analysed in Sect. 6 and finally, the conclusion is given in Sect. 7.

2 Preliminaries

Definition 1

The Atangana–Baleanu fractional derivative (ABFD) \(\frac{\partial ^{\gamma }}{\partial t^{\gamma }}w(r,t)\) without singular kernel of order \(\gamma \in (0,1)\) is given by

$$\begin{aligned} \frac{\partial ^{\gamma }}{\partial t^{\gamma }}w(r,t) =\frac{R(\gamma )}{1-\gamma } \int \limits _{0}^{t} \frac{\partial }{\partial \upsilon }w(r,\upsilon )E_{\gamma } \left[ -\frac{\gamma }{1-\gamma }(t-\upsilon )^{\gamma }\right] \mathrm{d}\upsilon , \end{aligned}$$
(2.1)

where \(R(\gamma )\) is a normalization function which satisfies \(R(0)=R(1)=1\). \(E_{\gamma ,\sigma }(r)\) is the Mittag–Leffler function satisfying \(E_{\gamma , 1}(r)=E_{\gamma }(r)\), which is defined as

$$\begin{aligned} E_{\gamma ,\sigma }(r)=\sum _{p=0}^{\infty }\frac{r^{p}}{\Gamma (\gamma p+\sigma )}. \end{aligned}$$

Definition 2

If \(\tilde{f}\in L^{2}[a,b]\), then Parseval’s identity is given as [45]:

$$\begin{aligned} \sum _{n=-\infty }^{\infty } |\hat{f}(n)|^2=\int _{a}^{b}|\tilde{f}(r)|^2 \mathrm{d}r, \end{aligned}$$
(2.2)

where \(\hat{f}(n)=\int _{a}^{b} \tilde{f}(r)e^{2\pi i nr} \mathrm{d}r\) for every integer n is its Fourier transform.

2.1 Basis functions

The spatial domain [a, b] such that \(a=r_0<r_1<\cdots <r_N=b\), where \(r_k=r_0+kh\), \(k=0(1)N\), be divided into N subintervals of equal length \(h=\frac{b-a}{N}\).

Now, we suppose that W(r, t) be the CBS approximation for w(r, t) s.t.

$$\begin{aligned} W(r,t)= \sum _{k={-1}}^{N+1}\mu _{k}^{m}(t)C_k(r), \end{aligned}$$
(2.3)

where \(\mu _{k}^{m}(t)\), the control points, are to be computed at each time level and \(C_k(r)\), the CBS functions, are defined as

$$\begin{aligned}&C_k(r)\nonumber \\&\quad =\frac{1}{6 h^3}\left\{ \begin{array}{ll} (r-r_{k-2})^3 , &{} \text {if }r\in [r_{k-2},r_{k-1}), \\ h^3+3h^2(r-r_{k-1})+3h(r-r_{k-1})^2-3(r-r_{k-1})^3 , &{} \text {if } r\in [r_{k-1},r_{k}), \\ h^3+3h^2(r_{k+1}-r)+3h(r_{k+1}-r)^2-3(r_{k+1}-r)^3 , &{} \text {if } r\in [r_{k},r_{k+1}), \\ (r_{k+2}-r)^3 , &{} \text {if }r\in [r_{k+1},r_{k+2}), \\ 0, &{} \text {otherwise}. \end{array}\right. \end{aligned}$$
(2.4)

The CBS functions retain abundant supply of geometrical properties, such as convex hull property, local support, geometric invariability, symmetry, non-negativity and the partition of unity [40]. Furthermore, \(C_{-1},C_0, \ldots , C_{N+1}\) have been organized in that manner that they can serve as basis for space of all third degree splines. The relations (2.3) and (2.4) yield the following approximations:

$$\begin{aligned} {\left\{ \begin{array}{ll} {(W)}_k^m= \left( \frac{1}{6}\right) \mu _{k-1}^m+\left( \frac{4}{6}\right) \mu _{k}^m+\left( \frac{1}{6}\right) \mu _{k+1}^m, \\ (W_r)_k^m=\left( \frac{1}{2h}\right) \mu _{k+1}^m+\left( -\frac{1}{2h}\right) \mu _{k-1}^m, \\ (W_{rr})_k^m= \left( \frac{1}{h^2}\right) \mu _{k-1}^m+\left( -\frac{2}{h^2}\right) \mu _{k}^m+\left( \frac{1}{h^2}\right) \mu _{k+1}^m. \end{array}\right. } \end{aligned}$$
(2.5)

3 Illustration of the scheme

Let the time domain [0, T] using the knots \(0=t_0<t_1<\cdots <t_M=T\), where \(t_m=m \Delta t\) and \(m=0,1,\ldots ,M\), be divided into M equal subintervals of the length \(\Delta t=\frac{T}{M}\). The ABFD involved in (1.1) is discretized at \(t=t_{m+1}\) as

$$\begin{aligned}&\frac{\partial ^{\gamma }}{\partial t^{\gamma }}w(r,t_{m+1}) \nonumber \\&\quad =\frac{R(\gamma )}{1-\gamma }\int \limits _{0}^{t_{m+1}} \frac{\partial }{\partial \upsilon }w(r,\upsilon )E_{\gamma }\nonumber \\&\quad \left[ -\frac{\gamma }{1-\gamma }(t_{m+1}-\upsilon )^{\gamma }\right] d\upsilon ,\ \ \ \ \ 0<\gamma <1, \nonumber \\&\quad =\frac{R(\gamma )}{1-\gamma } \sum _{s=0}^{m} \int \limits _{t_s}^{t_{s+1}} \frac{\partial }{\partial \upsilon }w (r,\upsilon )E_{\gamma }\left[ -\frac{\gamma }{1-\gamma }(t_{m+1} -\upsilon )^{\gamma }\right] \mathrm{d}\upsilon . \end{aligned}$$
(3.1)

Utilizing forward difference formulation, Eq. (3.1) becomes

$$\begin{aligned}&\frac{\partial ^{\gamma }}{\partial t^{\gamma }}w(r,t_{m+1}) \\&\quad =\frac{R(\gamma )}{1-\gamma } \sum _{s=0}^{m} \frac{w(r,t_{s+1})- w(r,t_s)}{\Delta t} \int \limits _{t_s}^{t_{s+1}} E_{\gamma }\\&\quad \left[ -\frac{\gamma }{1-\gamma }(t_{m+1}-\upsilon )^{\gamma }\right] \mathrm{d}\upsilon +\lambda _{\Delta t}^{m+1} \\&\quad =\frac{R(\gamma )}{1-\gamma } \sum _{s=0}^{m}\left[ w(r,t_{m-s+1}) -w(r,t_{m-s})\right] \left\{ (s+1)E_{\gamma ,2}\right. \\&\quad \left[ -\frac{\gamma }{1-\gamma } \left( (s+1)\Delta t\right) ^{\gamma }\right] \\&\quad \left. -sE_{\gamma ,2}\left[ -\frac{\gamma }{1-\gamma }(s\Delta t)^{\gamma }\right] \right\} +\lambda _{\Delta t}^{m+1}\\&\quad =\frac{R(\gamma )}{1-\gamma } \sum _{s=0}^{m}\left[ w(r,t_{m-s+1})- w(r,t_{m-s})\right] \\&\quad \left( (s+1)E_{s+1}-sE_{s}\right) +\lambda _{\Delta t}^{m+1}. \end{aligned}$$

Hence

$$\begin{aligned} \frac{\partial ^{\gamma }}{\partial t^{\gamma }}w(r,t_{m+1}) =\frac{R(\gamma )}{1-\gamma }\sum _{s=0}^{m} l_s \left[ w(r,t_{m-s+1}) -w(r,t_{m-s})\right] +\lambda _{\Delta t}^{m+1}, \end{aligned}$$
(3.2)

where \(E_{s}=E_{\gamma ,2}\left[ -\frac{\gamma }{1-\gamma }(s\Delta t)^{\gamma }\right] \) and \( l_s = (s+1)E_{s+1}-sE_{s}\). It is straightforward to observe that

  • \( l_s>0 \) and \( l_0=E_{1} \), \( s=1:1:m \),

  • \( l_0>l_1>l_2>\cdots > l_s, l_s \rightarrow 0 \) as \( s\rightarrow \infty \),

  • \(\sum _{s=0}^{m}(l_s-l_{s+1})+l_{m+1}=(E_{1}-l_1)+\sum _{s=1}^{m-1} (l_{s}-l_{s+1})+l_m=E_{1}\).

Moreover, the truncation error \( \lambda _{\Delta t}^{m+1}\) is given by [46]:

$$\begin{aligned} \lambda _{\Delta t}^{m+1}&=\frac{R(\gamma )}{1-\gamma }\sum _{s=0}^{m} \int \limits _{t_s}^{t_{s+1}} \frac{\Delta t}{2} \frac{\partial ^{2}w(r,t_s)}{\partial t^{2}} E_{\gamma } \left[ -\frac{\gamma }{1-\gamma }(t_{m+1}-\upsilon )^{\gamma }\right] \mathrm{d}\upsilon \\&=\frac{R(\gamma )}{1-\gamma }\frac{(\Delta t)^{2}}{2}\sum _{s=0}^{m} \frac{\partial ^{2}w(r,t_s)}{\partial t^{2}} \left\{ (m-s+1)E_{\gamma ,2}\right. \\&\quad \left[ -\frac{\gamma }{1-\gamma }\left( (m-s+1)\Delta t\right) ^{\gamma }\right] \\&\quad \left. -(m-s)E_{\gamma ,2}\left[ -\frac{\gamma }{1-\gamma }((m-s) \Delta t)^{\gamma }\right] \right\} \\&=\frac{R(\gamma )}{1-\gamma }\frac{(\Delta t)^{2}}{2}\sum _{s=0}^{m} \frac{\partial ^{2}w(r,t_s)}{\partial t^{2}} \\&\quad \left( (m-s+1)E_{m-s+1}-(m-s)E_{m-s}\right) \\&\le \frac{R(\gamma )}{1-\gamma }\frac{(\Delta t)^{2}}{2} \left[ \max _{0\le t\le t_{m}}\frac{\partial ^{2}w(r,t)}{\partial t^{2}}\right] c_{1}, \end{aligned}$$

where \(c_{1}\) is a constant:

$$\begin{aligned} | \lambda _{\Delta t}^{m+1} |\le \vartheta (\Delta t)^{2}, \end{aligned}$$
(3.3)

where \( \vartheta \) is a constant.

Using (3.2) and \(\theta \)-weighted scheme then Eq. (1.1) takes the following form:

$$\begin{aligned}&\frac{R(\gamma )}{1-\gamma }\sum _{s=0}^{m} l_s \left[ w(r,t_{m-s+1})- w(r,t_{m-s})\right] \nonumber \\&\quad = \theta \left( \Phi w_{rr}(r,t_{m+1})-\Psi w_{r}(r,t_{m+1})\right) \nonumber \\&\quad \quad + (1-\theta )\left( \Phi w_{rr}(r,t_{m})-\Psi w_{r}(r,t_{m})\right) +q(r,t_{m+1}). \end{aligned}$$
(3.4)

Discretizing (3.4) along spatial direction for \(\theta =1\), we acquire

$$\begin{aligned}&E_{1}w_{k}^{m+1}-\xi \Phi (w_{rr})_{k}^{m+1} +\xi \Psi (w_{r})_{k}^{m+1}\nonumber \\&\quad = E_{1}w_{k}^{m}-\sum _{s=1}^{m} l_s (w_{k}^{m-s+1}- w_{k}^{m-s})+\xi \, q_{k}^{m+1}, \end{aligned}$$
(3.5)

where \( \xi =\frac{1-\gamma }{R(\gamma )}\), \({w}_k^m=w(r_k,t_m)\) and \(q_{k}^{m+1}=q(r_k,t_{m+1})\).

Using (2.5) in (3.5), we get

$$\begin{aligned}&(E_{1}b_{1}-\xi \Phi b_{4}-\xi \Psi b_{3})\mu _{k-1}^{m+1}\nonumber \\&\qquad +(E_{1}b_{2}-\xi \Phi b_{5})\mu _{k}^{m+1} +(E_{1}b_{1}-\xi \Phi b_{4}+\xi \Psi b_{3})\mu _{k+1}^{m+1}\nonumber \\&\quad =E_{1}(b_{1}\mu _{k-1}^{m}+b_{2}\mu _{k}^{m}+b_{1}\mu _{k+1}^{m})\nonumber \\&\qquad -\sum _{s=1}^{m} l_s \left[ b_{1}(\mu _{k-1}^{m-s+1}-\mu _{k-1}^{m-s}) +b_{2}(\mu _{k}^{m-s+1}-\mu _{k}^{m-s})\nonumber \right. \\&\qquad \qquad \quad \left. +b_{1}(\mu _{k+1}^{m-s+1}-\mu _{k+1}^{m-s})\right] +\xi \,\, q_{k}^{m+1}, \end{aligned}$$
(3.6)

where \(b_1=\frac{1}{6}\), \(b_2=\frac{4}{6}\), \(b_3=\frac{1}{2h}\), \(b_4=\frac{1}{h^2}\), \(b_5=\frac{-2}{h^2}\) and \(\mu _{k}^m=\mu _{k}(t^{m})\).

The system (3.6) consists of \(N+1\) linear equations in \(N+3\) unknowns. For the unique solution, two additional equations can be achieved from the boundary conditions (1.3). Consequently, a matrix system of dimension \((N+3)\times (N+3)\) is obtained:

$$\begin{aligned} A\mu ^{m+1}=B\left( \sum _{s=0}^{m-1}(l_{s}-l_{s+1}) \mu ^{m-s}+l_m \mu ^0\right) +\xi J^{m+1}, \end{aligned}$$
(3.7)

where

$$\begin{aligned}&A=\begin{pmatrix} b_1 &{} b_2 &{} b_1 &{} &{} &{} &{} \\ y_{1} &{} y_{2} &{} y_{3} &{} &{} &{} &{} \\ &{} y_{1} &{} y_{2} &{}y_{3} &{} &{} &{} \\ &{} &{} \ddots &{} \ddots &{} \ddots &{} &{} \\ &{} &{} &{} y_{1} &{}y_{2} &{}y_{3} &{} \\ &{} &{} &{} &{} y_{1} &{} y_{2} &{}y_{3} \\ &{} &{} &{} &{} b_1 &{} b_2 &{} b_1 \end{pmatrix}, \\&B=\begin{pmatrix} 0 &{} 0 &{} 0 &{} &{} &{} &{} \\ b_1 &{} b_2 &{} b_1 &{} &{} &{} &{} \\ &{} b_1 &{} b_2 &{} b_1 &{} &{} &{} \\ &{} &{} \ddots &{} \ddots &{} \ddots &{} &{} \\ &{} &{} &{} b_1 &{} b_2 &{} b_1 &{} \\ &{} &{} &{} &{} b_1 &{} b_2 &{} b_1 \\ &{} &{} &{} &{} 0 &{} 0 &{} 0 \end{pmatrix},\\&\mu ^m=\begin{pmatrix} \mu ^{m}_{-1} \\ \mu ^{m}_0 \\ \mu ^{m}_1 \\ \vdots \\ \mu ^{m}_{N-1} \\ \mu ^{m}_{N} \\ \mu ^{m}_{N+1} \\ \end{pmatrix} \ \text {and}\ J^{m+1}=\begin{pmatrix} \psi _{1}^{m+1} \\ q^{m+1}_0 \\ q^{m+1}_1 \\ \vdots \\ q^{m+1}_{N-1} \\ q^{m+1}_{N} \\ \psi _{2}^{m+1} \\ \end{pmatrix}, \end{aligned}$$

where \(y_{1}= E_{1}b_{1}-\xi \Phi b_{4}-\xi \Psi b_{3}\), \(y_{2}=E_{1}b_{2}-\xi \Phi b_{5}\) and \(y_{3}=E_{1}b_{1}-\xi \Phi b_{4}+\xi \Psi b_{3}\). Before employing (3.7), the initial vector \(\mu ^{0} = [\mu ^{0}_{-1},\mu ^{0}_{0}, \ldots ,\mu ^{0}_{N+1}]^{T}\) is achieved by utilizing the initial conditions as

$$\begin{aligned} {\left\{ \begin{array}{ll} (W_r)_k^0=\phi '(r_k), &{} \ k=0, \\ (W)_k^0=\phi (r_k), &{} \ k=0:1:N, \\ (W_r)_k^0=\phi '(r_k), &{} \ k=N. \end{array}\right. } \end{aligned}$$
(3.8)

In matrix form, the above system of equations can be represented as

$$\begin{aligned} G \mu ^0=H, \end{aligned}$$
(3.9)

where

$$\begin{aligned}&G=\begin{pmatrix} -b_{3} &{} 0 &{} b_{3} &{} &{} &{} &{} \\ b_1 &{} b_2 &{} b_1 &{} &{} &{} &{} \\ &{} b_1 &{} b_2 &{} b_1 &{} &{} &{} \\ &{} &{} \ddots &{} \ddots &{} \ddots &{} &{} \\ &{} &{} &{} b_1 &{} b_2 &{} b_1 &{} \\ &{} &{} &{} &{} b_1 &{} b_2 &{} b_1 \\ &{} &{} &{} &{} -b_{3} &{} 0 &{} b_{3} \end{pmatrix}\ \text {and}\\&H=\begin{pmatrix} \phi '(r_0) \\ \phi (r_0) \\ \phi (r_1) \\ \vdots \\ \phi (r_{N-1}) \\ \phi (r_{N}) \\ \phi '(r_N) \\ \end{pmatrix}. \end{aligned}$$

Equation (3.9) is easily solvable for \(\mu ^{0}\) by means of a suitable numerical algorithm. All the numerical computations are performed in Mathematica 12.

4 The stability analysis

During the computational procedure, when the error does not amplify, the numerical scheme is presumed to be stable [47]. Here, the presented scheme is analysed for stability by implementing the Fourier method [41, 44]. For this, suppose \(\delta _k^m\) and \(\tilde{\delta }_k^m\) symbolize the growth factor and its estimation in Fourier mode. The error \(\chi _k^m\) can be presented as

$$\begin{aligned} \chi _k^m=\delta _k^m-\tilde{\delta }_k^m, \ \ m=0(1)M, \ \ k=1(1)N-1. \end{aligned}$$

For simplicity, we analyse the stability of the scheme given in (3.6) for force-free case (\(q=0\)) only, so that from Eq. (3.6), we obtain

$$\begin{aligned}&(E_{1}b_{1}-\xi \Phi b_{4}-\xi \Psi b_{3})\chi _{k-1}^{m+1} +(E_{1}b_{2}-\xi \Phi b_{5})\chi _{k}^{m+1}\nonumber \\&\qquad +(E_{1}b_{1}-\xi \Phi b_{4}+\xi \Psi b_{3})\chi _{k+1}^{m+1}\nonumber \\&\quad =E_{1}(b_{1}\chi _{k-1}^{m}+b_{2}\chi _{k}^{m}+b_{1}\chi _{k+1}^{m})\nonumber \\&\qquad -\sum _{s=1}^{m} l_s \left[ b_{1}(\chi _{k-1}^{m-s+1}-\chi _{k-1}^{m-s}) +b_{2}(\chi _{k}^{m-s+1}-\chi _{k}^{m-s})\right. \nonumber \\&\qquad \qquad \quad \left. +b_{1}(\chi _{k+1}^{m-s+1}-\chi _{k+1}^{m-s})\right] . \end{aligned}$$
(4.1)

From initial and boundary conditions, we can write

$$\begin{aligned} \chi _k^0 =\phi (r_k), \ \ k=1:1:N \end{aligned}$$
(4.2)

and

$$\begin{aligned} \chi _{0}^{m}=\psi _{1}(t_m), \ \ \chi _{N}^{m}=\psi _{2}(t_m), \ \ m=0:1:M. \end{aligned}$$
(4.3)

Define grid function as

$$\begin{aligned} \chi ^m={\left\{ \begin{array}{ll} \chi _{k}^{m}, &{} r_k-\frac{h}{2}<r\le r_k+\frac{h}{2}, \ k=1(1)N-1, \\ 0, &{} a\le r \le a+\frac{h}{2}\ or \ b-\frac{h}{2} \le r \le b. \end{array}\right. } \end{aligned}$$
(4.4)

The Fourier expansion for \(\chi ^m(r)\) can be represented as

$$\begin{aligned} \chi ^m(r)=\sum _{n=-\infty }^{\infty } \zeta ^{m}(n) \mathrm{e}^{\frac{2\pi z nr }{b-a}}, \end{aligned}$$
(4.5)

where

$$\begin{aligned}&\zeta ^{m}(n)=\frac{1}{b-a} \int \limits _{a}^{b} \chi ^m(r)\mathrm{e}^{\frac{-2\pi z nr }{b-a}} \mathrm{d}r, \nonumber \\&\quad m=0:1:M \ \ \ \text {and} \ \ \ \chi ^m=[\chi _1^m, \chi _2^m, \ldots , \chi _{N-1}^m]^T. \end{aligned}$$
(4.6)

Applying \( \Vert .\Vert _2\) norm, we acquire

$$\begin{aligned} \Vert \chi ^m\Vert _2&= \sqrt{ h \sum _{k=1}^{N-1}|\chi _k^m|^2 } \\&= \left( \int \limits _{a}^{a+\frac{h}{2}}|\chi ^m|^2 \mathrm{d}r\right. \\&\quad \left. +\sum _{k=1}^{N-1} \int \limits _{r_k-\frac{h}{2}}^{r_k+\frac{h}{2}} |\chi ^m|^2 \mathrm{d}r+ \int \limits _{b-\frac{h}{2}}^{b}|\chi ^m|^2 \mathrm{d}r \right) ^{\frac{1}{2}}\\&= \left( \int \limits _{a}^{b}|\chi ^m|^2 \mathrm{d}r \right) ^{\frac{1}{2}}. \end{aligned}$$

By Parseval’s identity (2.2), we have [45]

$$\begin{aligned} \int \limits _{a}^{b}|\chi ^m|^2 \mathrm{d}r = \sum _{n=-\infty }^{\infty } |\zeta ^{m}(n)|^2. \end{aligned}$$

Hence, we acquire

$$\begin{aligned} \Vert \chi ^m\Vert _2^2 = \sum _{n=-\infty }^{\infty } |\zeta ^{m}(n)|^2. \end{aligned}$$
(4.7)

Assume that (4.1)–(4.3) possess a solution in Fourier form as:

$$\begin{aligned} \chi _k^m=\zeta ^{m} \mathrm{e}^{i \alpha k h}, \end{aligned}$$
(4.8)

where \(\alpha \) is a real number and \(i=\sqrt{-1}\). Substituting (4.8) in (4.1) and simplifying, we acquire

$$\begin{aligned}&(E_{1}b_{1}-\xi \Phi b_{4}-\xi \Psi b_{3})\zeta ^{m+1}\mathrm{e}^{-i \alpha h} +(E_{1}b_{2}-\xi \Phi b_{5})\zeta ^{m+1}\nonumber \\&\quad +(E_{1}b_{1}-\xi \Phi b_{4}+\xi \Psi b_{3})\zeta ^{m+1}\mathrm{e}^{i \alpha h}\nonumber \\&\quad =(E_{1}b_{1}\zeta ^{m}\mathrm{e}^{-i \alpha h}+E_{1}b_{2}\zeta ^{m} +E_{1}b_{1}\zeta ^{m}\mathrm{e}^{i \alpha h})\nonumber \\&\quad -\sum _{s=1}^{m} l_s \left[ b_{1}(\zeta ^{m-s+1}\mathrm{e}^{-i \alpha h} -\zeta ^{m-s}\mathrm{e}^{-i \alpha h})\right. \nonumber \\&\quad + b_{2}(\zeta ^{m-s+1}-\zeta ^{m-s})\nonumber \\&\quad \left. + b_{1}(\zeta ^{m-s+1}\mathrm{e}^{i \alpha h}-\zeta ^{m-s}\mathrm{e}^{i \alpha h})\right] . \end{aligned}$$
(4.9)

Utilizing the relation \(\mathrm{e}^{i \alpha h }+\mathrm{e}^{-i \alpha h } = 2 \cos (\alpha h)\) and assembling the identical terms, we acquire the following equality:

$$\begin{aligned} \zeta ^{m+1}=\frac{1}{\epsilon }\left[ \zeta ^{m}-\frac{1}{E_{1}} \sum _{s=1}^{m} l_s\left( \zeta ^{m-s+1}- \zeta ^{m-s}\right) \right] , \end{aligned}$$
(4.10)

where \(\epsilon =1+\frac{12\xi \Phi \sin ^2(\alpha h/2)+3\xi \Psi kh\sin \alpha h}{h^2 E_{1}(3-2 \sin ^2 (\alpha h /2) )}\). Obviously \(\epsilon \ge 1\) and \(E_{1}>0\).

Lemma 4.1

If \(\zeta ^{m}\) is the solution for equation (4.10), then \(|\zeta ^m|\le |\zeta ^0|, \ \ m=0, 1, 2, \cdots , M \).

Proof

Mathematical induction will be availed. For \(m=0\), (4.10) gives

$$\begin{aligned} |\zeta ^{1}|=\frac{1}{\epsilon }|\zeta ^{0}|\le |\zeta ^{0}|, \ \ \ \ \epsilon \ge 1. \end{aligned}$$

Assume that \(|\zeta ^m|\le |\zeta ^0|\) for \(m=1,2,\ldots , M-1\), then

$$\begin{aligned} | \zeta ^{m+1}|&\le \frac{1}{\epsilon }|\zeta ^{m}|-\frac{1}{\epsilon E_{1}} \sum _{s=1}^{m} l_s\left( | \zeta ^{m-s+1}|- |\zeta ^{m-s}|\right) \\&\le \frac{1}{\epsilon }|\zeta ^{0}|-\frac{1}{\epsilon E_{1}} \sum _{s=1}^{m} l_s\left( | \zeta ^{0}|- |\zeta ^{0}|\right) \\&\le |\zeta ^0|. \end{aligned}$$

\(\square \)

Theorem 1

The present scheme (3.6) is unconditionally stable.

Proof

Utilizing Lemma 4.1 and the expression (4.7), we acquire

$$\begin{aligned} \Vert \chi ^m\Vert _2 \le |\chi ^0|_2, \ \ \ \forall m=0, 1, \ldots , M. \end{aligned}$$

Hence, the proposed numerical method is stable unconditionally. \(\square \)

5 The convergence analysis

We follow the procedure employed in [48] to study the convergence of the propounded method. First of all, the following theorem is brought forward [49, 50].

Theorem 2

Suppose that w(r, t), q belong to \(C^4[a,b]\) and \(C^2[a,b]\) , respectively, and \(\Upsilon =\{a=r_0,r_1, \ldots , r_N=b\}\) be the partition of [a, b] s.t. \(r_k=a+kh, k= 0,1, \ldots , N\). Let \(\tilde{W}(r,t)\) be the unique spline that interpolates the solution curve at knots \(r_k \in \Upsilon \), then there is a constant \(\varrho _k\) not depending on h, then for every \(t\ge 0\), we acquire

$$\begin{aligned} \Vert D^k\left( w(r,t)-\tilde{W}(r,t)\right) \Vert _\infty \le \varrho _k h^{4-k}, \ \ k=0,1,2. \end{aligned}$$
(5.1)

Lemma 5.1

The CBS set \(\{C_{-1}, C_{0}, \ldots , C_{N+1}\}\) presented in (2.4) satisfies the inequality

$$\begin{aligned} \sum _{k=-1}^{N+1}|C_k(r)| \le \frac{5}{3}, \ \ \ 0 \le r \le 1. \end{aligned}$$
(5.2)

Proof

Employing the triangular inequality, we are able to write

$$\begin{aligned} \left| \sum _{k=-1}^{N+1}C_k(r)\right| \le \sum _{k=-1}^{N+1}|C_k(r)|. \end{aligned}$$

At any knot \(r_k\), we achieve

$$\begin{aligned} \sum _{k=-1}^{N+1}|C_k(r)|&= |C_{k-1}(r_k)|+|C_{k}(r_{k})|+|C_{k+1}(r_k)| \\&= \frac{1}{6}+ \frac{4}{6}+\frac{1}{6} =1 < \frac{5}{3}. \end{aligned}$$

Moreover, for \(r\in [r_{k},r_{k+1}]\), we get

$$\begin{aligned} |C_{k}(r)|\le \frac{4}{6}\ \ \ \&\text {and} \ \ \ \ |C_{k-1}(r)|\le \frac{1}{6}, \\ |C_{k+1}(r)|\le \frac{4}{6}\ \ \ \&\text {and} \ \ \ \ |C_{k+2}(r)|\le \frac{1}{6}. \end{aligned}$$

Thus, for any point \( r_{k} \le r \le r_{k+1}\), we obtain

$$\begin{aligned} \sum _{k=-1}^{N+1}|C_k(r)| = |C_{k-1}(r)|+|C_{k}(r)|+|C_{k+1}(r)|+|C_{k+2}(r)| \le \frac{5}{3}. \end{aligned}$$

\(\square \)

Theorem 3

The computational approximation W(r, t) to the analytical solution w(r, t) for TFADE (1.1)–(1.3) exists. Furthermore, if q belongs to \(C^2[0,1]\), then

$$\begin{aligned} \Vert w(r,t)-W(r,t)\Vert _\infty \le \widetilde{\varrho } h^2, \ \ \forall \ \ t \ge 0, \end{aligned}$$
(5.3)

where h is appropriately small and \(\widetilde{\varrho }>0\) is free of h.

Proof

We assume that \(\tilde{W}(r,t)=\sum _{k=-1}^{N+1} u_k^{m}(t)C_k(r)\) be the approximated spline for W(r, t). By means of triangular inequality, we acquire

$$\begin{aligned} \Vert w(r,t)-W(r,t)\Vert _\infty \le \Vert w(r,t)-\tilde{W}(r,t)\Vert _\infty +\Vert \tilde{W}(r,t)-W(r,t)\Vert _\infty . \end{aligned}$$

With the aid of Theorem 2, we achieve

$$\begin{aligned} \Vert w(r,t)-W(r,t)\Vert _\infty \le \varrho _{0} h^{4} +\Vert \tilde{W}(r,t)-W(r,t)\Vert _\infty . \end{aligned}$$
(5.4)

The collocation conditions of proposed scheme are \(Lw(r_k,t)=LW(r_k,t)=q(r_k,t)\), \(k=0(1)N \). Suppose that

$$\begin{aligned} L\tilde{W}(r,t)=\tilde{q}(r_k,t), \ \ k=0(1)N. \end{aligned}$$

Thus, the difference \( L(\tilde{W}(r_k,t)-W(r_k,t))\) for any time stage m can be stated as

$$\begin{aligned}&(E_{1}b_{1}-\xi \Phi b_{4}-\xi \Psi b_{3})\Omega _{k-1}^{m+1}\nonumber \\&\qquad +(E_{1}b_{2}-\xi \Phi b_{5})\Omega _{k}^{m+1} +(E_{1}b_{1} -\xi \Phi b_{4}+\xi \Psi b_{3})\Omega _{k+1}^{m+1}\nonumber \\&\quad =(E_{1}b_{1}\Omega _{k-1}^{m}+E_{1}b_{2}\Omega _{k}^{m} +E_{1}b_{1}\Omega _{k+1}^{m})+\frac{\xi }{h^{2}}\,\, \rho _{k}^{m+1}\nonumber \\&\qquad -\sum _{s=1}^{m} l_s \left[ b_{1}(\Omega _{k-1}^{m-s+1} -\Omega _{k-1}^{m-s})+b_{2}(\Omega _{k}^{m-s+1}-\Omega _{k}^{m-s})\right. \nonumber \\&\qquad \qquad \quad \left. +b_{1}(\Omega _{k+1}^{m-s+1}-\Omega _{k+1}^{m-s})\right] . \end{aligned}$$
(5.5)

The boundary conditions can be described as

$$\begin{aligned} b_1\Omega _{k-1}^{m+1}+b_2\Omega _{k}^{m+1}+b_1\Omega _{k+1}^{m+1}=0, \ \ k=0,N, \end{aligned}$$

where

$$\begin{aligned} \Omega _k^m=\mu _k^m-u_k^m, \ \ k=-1:0:N+1 \end{aligned}$$

and

$$\begin{aligned} \rho _k^m=h^2[q_k^m-\tilde{q}_k^m], \ \ k=0:1:N. \end{aligned}$$

It is obvious from inequality (5.1) that

$$\begin{aligned} |\rho _k^m|=h^2|q_k^m-\tilde{q}_k^m|\le \varrho h^4. \end{aligned}$$

Define \(\rho ^m=max \{|\rho _k^m|;0 \le k \le N \}\), \(e_k^m=|\Omega _k^m|\) and \(e^m=max\{|e_k^m|;0 \le k \le N \}\).

When \(m=0\), Eq. (5.5) becomes

$$\begin{aligned}&(E_{1}b_{1}-\xi \Phi b_{4}-\xi \Psi b_{3})\Omega _{k-1}^{1}\\&\qquad +(E_{1}b_{2}-\xi \Phi b_{5})\Omega _{k}^{1} + (E_{1}b_{1} -\xi \Phi b_{4}+\xi \Psi b_{3})\Omega _{k+1}^{1}\\&\quad =(E_{1}b_{1}\Omega _{k-1}^{0}+E_{1}b_{2}\Omega _{k}^{0} +E_{1}b_{1}\Omega _{k+1}^{0})+\frac{\xi }{h^{2}}\,\, \rho _{k}^{1}, \end{aligned}$$

where \(k=0:1:N\). With the help of initial condition, \(e^0=0\):

$$\begin{aligned} (E_{1}b_{2}-\xi \Phi b_{5})\Omega _{k}^{1}= & {} -(E_{1}b_{1}-\xi \Phi b_{4}) (\Omega _{k-1}^{1}+\Omega _{k+1}^{1})\\&+\xi \Psi b_{3}(\Omega _{k-1}^{1} -\Omega _{k+1}^{1})+\frac{\xi }{h^{2}}\,\, \rho _{k}^{1}. \end{aligned}$$

Taking norms of \(\Omega _k^1\), \(\rho _k^1\) and adequately small mesh spacing h, we acquire

$$\begin{aligned} e_k^1\le \frac{3 \xi }{E_{1} h^2+12\xi \Phi -3\xi \Psi h}\varrho h^4, \ \ \ \ k=0:1:N. \end{aligned}$$

We obtain the values of \(e_{-1}^1\) and \( e_{N+1}^1\) from the boundary conditions:

$$\begin{aligned}&e_{-1}^1\le \frac{15 \xi }{E_{1} h^2+12\xi \Phi -3\xi \Psi h}\varrho h^4, \\&e_{N+1}^1\le \frac{15 \xi }{E_{1} h^2+12\xi \Phi -3\xi \Psi h}\varrho h^4, \end{aligned}$$

which implies

$$\begin{aligned} e^1\le \varrho _1 h^2, \end{aligned}$$
(5.6)

where \(\varrho _1\) is not depending on h.

To prove this theorem, mathematical induction on m is used. For this, we consider that \(e_k^f\le \varrho _f h^2\) is true for \(f=1,2,\ldots ,m\) and \(\varrho =max\{\varrho _f: f=0,1,\ldots ,m \}\), then from Eq. (5.5), we obtain

$$\begin{aligned}&(E_{1}b_{1}-\xi \Phi b_{4}-\xi \Psi b_{3})\Omega _{k-1}^{m+1} \\&\quad +(E_{1}b_{2}-\xi \Phi b_{5})\Omega _{k}^{m+1} + (E_{1}b_{1} -\xi \Phi b_{4}+\xi \Psi b_{3})\Omega _{k+1}^{m+1}\\&\quad =\left[ (l_0-l_1)(b_1 \Omega _{k-1}^{m}+b_2 \Omega _{k}^{m}+b_1 \Omega _{k+1}^{m})\right. \\&\quad +(l_1-l_2)(b_1 \Omega _{k-1}^{m-1}+b_2 \Omega _{k}^{m-1} +b_1 \Omega _{k+1}^{m-1})+ \cdots \\&\quad +(l_{m-1}-l_m)(b_1 \Omega _{k-1}^{1}+b_2 \Omega _{k}^{1}+b_1 \Omega _{k+1}^{1})\\&\quad \left. +l_m(b_1 \Omega _{k-1}^{0}+b_2 \Omega _{k}^{0} +b_1 \Omega _{k+1}^{0})\right] +\frac{\xi }{h^{2}}\,\, \rho _{k}^{m+1}. \end{aligned}$$

Again, employing the norms on \(\Omega _k^m\) and \(\rho _k^m\), we acquire

$$\begin{aligned} e_k^{m+1}\le \frac{3 \varrho h^4}{ E_{1}h^2+12\xi \Phi -3\xi \Psi h} \left( \xi + \sum _{s=0}^{m-1}( l_s -l_{s+1})\right) . \end{aligned}$$

Similarly, we get the values of \(e_{-1}^{m+1}\) and \( e_{N+1}^{m+1}\) from the boundary conditions:

$$\begin{aligned}&e_{-1}^{m+1}\le \frac{15 \varrho h^4}{E_{1}h^2+12\xi \Phi -3\xi \Psi h} \left( \xi + \sum _{s=0}^{m-1}( l_s -l_{s+1})\right) \\&e_{N+1}^{m+1}\le \frac{15 \varrho h^4}{E_{1}h^2+12\xi \Phi -3\xi \Psi h} \left( \xi + \sum _{s=0}^{m-1}( l_s -l_{s+1})\right) . \end{aligned}$$

Hence, for all m, we acquire

$$\begin{aligned} e^{m+1}\le \varrho h^2. \end{aligned}$$
(5.7)

In particular,

$$\begin{aligned} \tilde{W}(r,t)-W(r,t)=\sum _{k=-1}^{N+1}\left( u_k(t)-\mu _k(t)\right) C_k(r). \end{aligned}$$

Therefore, from Lemma 5.1 and inequality (5.7), we get

$$\begin{aligned} \Vert \tilde{W}(r,t)-W(r,t)\Vert _\infty \le \frac{5}{3} \varrho h^2. \end{aligned}$$
(5.8)

Using (5.8), the inequality (5.4) gives

$$\begin{aligned} \Vert w(r,t)-W(r,t)\Vert _\infty \le \varrho _{0} h^{4} +\frac{5}{3} \varrho h^2=\widetilde{\varrho } h^2, \end{aligned}$$

where \(\widetilde{\varrho }=\varrho _{0}h^2+\frac{5}{3} \varrho \). \(\square \)

Theorem 4

The TFADE with initial and boundary conditions is convergent.

Proof

Consider w(r, t) and W(r, t) be the solutions for TFADE analytically and numerically, respectively. Therefore, relation (3.3) and the above theorem validate that there are arbitrary constants \(\widetilde{\varrho }\) and \(\vartheta \), such that

$$\begin{aligned} \Vert w(r,t)-W(r,t)\Vert _\infty \le \widetilde{\varrho } h^2+\vartheta (\Delta t)^{2}. \end{aligned}$$

Hence, the presented scheme is second order convergent in spatial and time directions. \(\square \)

6 Illustration of numerical examples and discussion

Fig. 1
figure 1

Exact and numerical solutions for Example 6.1 at different time stages with \(\Delta t=0.001\)

Full size image
Fig. 2
figure 2

3D Exact and approximate solutions for Example 6.1, when \(N=1000\), \(t=1\), \(\gamma =0.4\), \(\Delta t=0.001\) and \(r\in [0,1]\)

Full size image
Fig. 3
figure 3

2D and 3D Error profiles for Example 6.1, when \(N=1000\), \(t=1\), \(\gamma =0.4\), \(\Delta t=0.001\) and \(r\in [0,1]\)

Full size image

In this section, numerical outcomes are reported to check the accuracy of the scheme through the error norms \(L_{2}(z)\), \(L_{\infty }(z)\) as

$$\begin{aligned} L_{2}(z)= & {} \Vert w(r_{k},t)-W(r_{k},t)\Vert _2\\= & {} \sqrt{h \sum _{r=0}^{N}| w(r_k,t)-W(r_k,t) |^2}, \\ L_{\infty }(z)= & {} \Vert w(r_{k},t)-W(r_{k},t)\Vert _\infty \\= & {} \max _{0\le k \le N}| w(r_k,t)-W(r_k,t) | \end{aligned}$$

and the convergence order [46] as

$$\begin{aligned} \log \left( \frac{L_{\infty }(z)}{L_{\infty }(z+1)}\right) \Bigg {/} \log \left( \frac{Q(z+1)}{Q(z)}\right) . \end{aligned}$$

All examples are solved by considering normalization function \(R(\gamma )=1\).

Example 6.1

Consider the TFADE [46]

$$\begin{aligned} \frac{\partial ^{\gamma }w(r,t)}{\partial t^{\gamma }}= & {} \frac{\partial ^{2}w(r,t)}{\partial r^2} -\frac{\partial w(r,t)}{\partial r}\\&+q(r,t) , \ \ \ r\in [0,1],\,\, t\in [0,1], \end{aligned}$$

with initial condition:

$$\begin{aligned} w(r,0)= 0 \end{aligned}$$

and boundary conditions:

$$\begin{aligned} w(0,t)=0, \ \ w(1,t)=0, \end{aligned}$$

where \(q(r,t)=2 \left( \frac{R(\gamma )}{1-\gamma }\right) r(r-1)t^{2}E_{\gamma ,3}\left[ \frac{-\gamma }{1-\gamma } t^{\gamma }\right] -2t^{2}+(2r-1)t^{2}\).

The exact solution is \(w(r,t)= r(r-1)t^{2}\). Tables 1 and 2 show the numerical outcomes and absolute errors of Example 6.1 for different choices of \(\gamma \) setting \(N=80, 10, 50\), \(\Delta t=0.001\) and \(t=1\). For different choices of \(\gamma \), error norms are reported in Table 3 at various time stages. A Comparison of error norm and convergence order with those obtained in [46] is given in Tables 4 and 5 along temporal and spatial directions. In Fig. 1, we can observe that close agreement between numerical results of the presented method and exact solutions at different time stages with \(\Delta t=0.001\). A 3-D plot of exact solution and computational results is illustrated in Fig. 2 by setting \(N=1000\), \(\Delta t=0.001\), \(\gamma =0.4\), \(t=1\) and \(r\in [0,1]\). In Fig. 3, the 2D and 3D error profiles are displayed at \(t=1\). Tables and figures represent that the proposed method is compatible with stated exact solution.

Example 6.2

Consider the TFADE [46]

$$\begin{aligned} \frac{\partial ^{\gamma }w(r,t)}{\partial t^{\gamma }}= & {} \frac{\partial ^{2}w(r,t)}{\partial r^2}-\frac{\partial w(r,t)}{\partial r}\\&+q(r,t) , \ \ \ r\in [0,1],\,\, t\in [0,1], \end{aligned}$$

with initial condition:

$$\begin{aligned} w(r,0)=0 \end{aligned}$$

and boundary conditions:

$$\begin{aligned} w(0,t)=0, \ \ w(1,t)=0, \end{aligned}$$

where \(q(r,t)=120 \left( \frac{R(\gamma )}{1-\gamma }\right) t^{5}\sin (\pi r)E_{\gamma ,6}\left[ \frac{-\gamma }{1-\gamma } t^{\gamma }\right] +\pi t^{5}\left( \pi \sin (\pi r)+\cos (\pi r)\right) \).

Table 1 Absolute errors for various choices of \(\gamma \), where \(\Delta t=0.001\) and \(N=80\) of Example 6.1 at \(t=1\)
Full size table
Table 2 Absolute errors for Example 6.1 at \(t=1\) and \(\Delta t=0.001\)
Full size table
Table 3 Error norms for various choices of \(\gamma \), where \(\Delta t=0.0005\), \(N=32\) and \(r \in [0,1]\) of Example 6.1
Full size table
Table 4 Comparison of error norm for Example 6.1 with different values of \(\Delta t=\frac{1}{m}\) when \(t=1\) and fixed \(h=\frac{1}{1000}\)
Full size table
Table 5 Comparison of error norm for Example 6.1 with different values of \(\Delta t=\frac{1}{m}\) and \(h=\frac{1}{n}\) when \(t=1\)
Full size table

The exact solution is \(w(r,t)=t^{5}\sin (\pi r)\). The absolute errors for Example 6.2 at numerous values of r setting \(\Delta t=0.0125, 0.05\), \(N=500, 70\), \(\gamma =0.2, 0.4\) and \(t=1\) are tabulated in Tables 6 and 7. Tables 8 and 9 display the error norms for different values of \(\gamma \) subject to \(N=85, 285\), \(\Delta t=0.01\) and \(r \in [0,1]\) at various time stages. Comparison of convergence order and errors are displayed in Tables 10 and 11 in temporal and spatial grids, respectively. It is found that the proposed scheme shows better efficiency and accuracy as compared to the scheme proposed in [46]. Figure 4 illustrates the behavior of computational results and exact values at different time stages. In Fig. 5, 3-D graphs of numerical outcomes and analytical solutions depict the accuracy of current scheme. Figure 6 shows the 2D and 3D error profiles, which exhibit exactness of the method.

Table 6 Absolute error for Example 6.2 at \(t=1\), \(\gamma = 0.2\), \(\Delta t=0.0125\) and \(N=500\)
Full size table
Table 7 Absolute error for Example 6.2 at \(t=1\), \(\gamma = 0.4\), \(\Delta t=0.05\) and \(N=70\)
Full size table
Table 8 Error norms for various choices of \(\gamma \), where \(\Delta t=0.01\), \(N=85\) and \(r \in [0,1]\) of Example 6.2
Full size table
Table 9 Error norms for Example 6.2 when \(\Delta t=0.01\), \(N=285\) and \(r \in [0,1]\)
Full size table
Table 10 Comparison of error norm for Example 6.2 with different values of \(\Delta t=\frac{1}{m}\) when \(t=1\) and fixed \(h=\frac{1}{1000}\)
Full size table
Table 11 Comparison of error norm for Example 6.2 with different values of \(h=\frac{1}{n}\) and \(\Delta t=\frac{1}{500}\) when \(t=1\)
Full size table
Fig. 4
figure 4

Exact and numerical solutions for Example 6.2 at different time stages

Full size image
Fig. 5
figure 5

3D exact and approximate solutions for Example 6.2, when \(N=1000\), \(t=1\), \(\gamma =0.4\), \(\Delta t=0.001\) and \(r\in [0,1]\)

Full size image
Fig. 6
figure 6

2D and 3D error profiles for Example 6.2, when \(N=1000\), \(t=1\), \(\gamma =0.4\), \(\Delta t=0.001\) and \(r\in [0,1]\)

Full size image

7 Conclusion

In current paper, an efficient solution of TFADE involving Atangana–Baleanu time derivative has attained through numerical scheme based on CBS functions. The standard finite difference formulation has been used to approximate the Atangana–Baleanu time fractional derivative, while interpolation of solution curve in spatial direction is obtained through CBS functions. The method proposed in current study is novel and provides reasonable accuracy when compared to those already developed in literature. The method is unconditionally stable along with second order spatial and temporal convergence. The implementation of the current algorithm on numerical examples reveal this method to be more efficient, simple and admissible.