1 Introduction

1.1 Problem statement

In this work, we consider the following one dimensional time fractional telegraph equation (TFTE) with reaction term [1]:

$$ \frac{\partial ^{2\alpha }u(x,t)}{\partial t^{2\alpha }}+2\lambda \frac{ \partial ^{\alpha }u(x,t)}{\partial t^{\alpha }}+\mu u(x,t)= \nu \frac{ \partial ^{2}u(x,t)}{\partial x^{2}}+f(x,t), \quad (x,t)\in [a,b]\times [0,T] $$
(1)

with initial and boundary conditions

$$ \textstyle\begin{cases} u(x,0)=f_{1}(x), \\ u_{t}(x,0)=f_{2}(x), \\ u(a,t)=g_{1}(t), \\ u(b,t)=g_{2}(t), \end{cases} $$
(2)

where \(0<\alpha <1\), and λ, μ, ν are arbitrary positive constants. \(f(x)\) is the forcing term and \(f_{1}(x)\), \(f_{2}(x)\), \(g_{1}(x)\), \(g_{2}(x)\) are sufficiently smooth prescribed functions. If \(\alpha =1\) Eq. (1) becomes the one dimensional hyperbolic telegraph equation. Time fractional derivatives \(\frac{\partial ^{2\alpha }u(x,t)}{\partial t^{2\alpha }}\) and \(\frac{\partial ^{\alpha }u(x,t)}{\partial t^{\alpha }}\) denote the Caputo fractional derivative of order 2α and α, respectively, which will be defined in Sect. 2.

This model has been developed to overcome the shortcomings of the classical telegraph equation which might not adequately model the abnormal diffusion phenomena during a long transmission process in a transmission line [2].

1.2 Application and literature review

Fractional differential equations have generated significant interest due to their appearance in various fields. Fractional differential equation models are more effective for the description of certain systems. For example, fractional order derivatives have been used successfully in diffusion processes, rheology, damping law visco-elasticity and fluid mechanics. They also appear in the modeling of many mathematical biology, chemical processes and a number of problems in engineering [3,4,5]. Our study will focus on the numerical solution of reaction diffusion models which contain fractional order derivatives. Telegraph equations are hyperbolic partial differential equations that are applicable in modeling the reaction diffusion processes. These models appear in the study of random walk theory, wave phenomena and wave propagation of electrical signal in the cable of a transmission line [6,7,8,9].

There has been much interest in TFTE as of lately. The well-posedness and asymptotical study about TFTE using Riemann–Liouville approach have been discussed by Cascaval et al. [10]. Analytical solution for the TFTE with three different nonhomogeneous boundary conditions using separation of variables has been derived by Chen et al. [11]. Approximate solutions of space and TFTE using Adomian decomposition method have been discussed by Momani [12]. Cauchy and signaling problems using Laplace and Fourier transforms and the boundary problem using spatial Sine transform have been solved by Huang [13] who derived analytical solution for three basic problems of TFTE. Dehghan and Shokri [14] presented a numerical method for solving hyperbolic telegraph equation using collocation points and approximated the solution via thin plate spline radial basic functions. Yousefi [15] solved the hyperbolic telegraph equation via the Legendre multiwavelet Galerkin method. Wang et al. [16] discussed and analyzed the Galerkin mixed finite element method for the numerical solution of TFTE. Li and Cao [17] presented a scheme based on a finite difference method for a kind of linear TFTE. Saadatmandi and Mohabbati [18] developed a computational technique for solving TFTE based on the Tau method and Legendre polynomials. Alkahtani et al. [19] studied the space-time fractional equation and obtained the solution via the Sumudu variational iteration method which is a combination of the Sumudu transform and the variational iteration method. Asgari et al. [1] discussed a method based on Bernstein polynomials operational matrices for solving TFTE. Hashemi and Baleanu [2] proposed a numerical method for the solution of TFTE using the Caputo fractional derivative in time direction, a combination of a group preserving scheme and the method of line for spatial direction. The reproducing kernel method has been presented for the solution of TFTE with initial boundary conditions by Wang et al. [20]. Wang and Mei [21] proposed a method for solving TFTE via the Legendre spectral Galerkin method and generalized finite difference scheme. Liu [22] discussed the Caputo fractional difference formula and Grünwald difference formula for the solution of TFTE.

Motivated by the success of B-splines in the numerical solution of differential equations, our focus is to use an appropriate B-splines for the numerical solution of TFTE. So far as we are aware there are no such studies on the use of splines for the fractional telegraph partial differential equation. In the literature, there are some studies, based on splines, for solving fractional partial differential equations. Tasbozan et al. [23] developed a numerical solution of fractional diffusion equation via the cubic B-spline collocation method. Akram and Tariq [24] presented a numerical scheme based on the quintic spline collocation method for the solution of fractional boundary value problems. The cubic B-spline collocation method has been used for the solution of fractional diffusion equation by Sayevand et al. [25]. Arshed [26] solved a time fractional super-diffusion fourth order differential equation using the quintic B-spline collocation method. Tasbozan and Esen [27] discussed numerical solutions of TFTE using the quadratic B-spline Galerkin method. Yaseen et al. [28] presented a scheme for the numerical solution of fractional diffusion equation using a finite difference method based on cubic trigonometric B-spline basis functions. Mohyud-Din et al. [29] constructed a fully implicit finite difference scheme for solving a time fractional diffusion equation by incorporating an extended cubic B-spline (ExCuBs) approach in its formulation. Because of the promising results obtained by this scheme, efforts are now being made in this work extend the formulation to solve the more complicated telegraph equation with fractional order derivatives. The purpose of this study is to describe a possible way to find the numerical solution of telegraph model which contains fractional order derivatives.

The structure of the paper is organized as follows: Sect. 2 provides the Caputo fractional derivative and basis functions. Solving the TFTE for the discretization of time by Caputo fractional derivative is presented in Sect. 3. The finite difference scheme by ExCuBs is discussed in Sect. 4. Initial state is presented in Sect. 5. The stability analysis and convergence are discussed in Sect. 6 and Sect. 7, respectively. In Sect. 8, some numerical experiments of TFTE are discussed to illustrate the reliability and capacity of proposed scheme. Finally a conclusion is discussed in Sect. 9.

2 Mathematical preliminaries

2.1 Fractional derivative

Definition 1

The Caputo’s time fractional derivative of order α is defined as [30]

$$ \frac{\partial ^{\alpha }u(x,t)}{\partial t^{\alpha }}= \textstyle\begin{cases} \frac{1}{\varGamma (n-\alpha )}\int ^{t}_{0}\frac{\partial u(x,\tau )}{ \partial \tau }\frac{d\xi }{(t-\tau )^{\alpha -n+1}}, & n-1< \alpha \leq n, n\in N, \\ \frac{\partial ^{n} u(x,t)}{\partial t^{n}}, & \alpha =n \in N. \end{cases} $$
(3)

where Γ is the Euler Gamma function.

By using Eq. (3), we can easily derive the fractional derivative of order 2α, which will be discussed in Sect. 3.

2.2 Basis functions

Suppose that the interval \([a,b]\) is divided into equal partitioning by the spatial knots \(x_{i}\) i.e. \(a=x_{0}< x_{1}<\cdots<x_{N-1}<x_{N}=b\) into N subintervals \([x_{i}, x_{i+1}]\) with equal length \(h= \frac{b-a}{N}\), \(i=0,1,\ldots,N\). The ExCuBs basis functions at the nodal point \(x_{i}\) over the interval \([a,b]\) can be presented as [31]

$$\begin{aligned}& E_{i}(x,\lambda ) \\ & \quad =\frac{1}{24 h^{4}} \textstyle\begin{cases} 4h(1-\lambda )(x-x_{i})^{3}+3\lambda (x-x_{i})^{4}, & x\in [x_{i}, x_{i+1}], \\ (4-\lambda )h^{4}+12h^{3}(x-x_{i+1})+6h^{2}(2+\lambda )(x-x_{i+1})^{2} \\ \quad {}-12h(x-x_{i+1})^{3}-3\lambda (x-x_{i+1})^{4}, & x\in [x_{i+1}, x_{i+2}], \\ (4-\lambda )h^{4}+12h^{3}(x_{i+3}-x)+6h^{2}(2+\lambda )(x_{i+3}-x)^{2} \\ \quad {}-12h(x_{i+3}-x)^{3}-3\lambda (x_{i+3}-x)^{4}, & x\in [x_{i+2}, x_{i+3}], \\ 4h(1-\lambda )(x_{i+4}-x)^{3}+3\lambda (x_{i+4}-x)^{4}, & x \in [x_{i+3}, x_{i+4}], \\ 0, & \text{otherwise,} \end{cases}\displaystyle \end{aligned}$$
(4)

where \(\lambda \in \mathbb{R}\) is a free parameter and \(x\in \mathbb{R}\) is a variable. For \(-8\leq \lambda \leq 1\), the ExCuBs basis and cubic B-spline possess the same properties. When \(\lambda =0\) it should be stated that the ExCuBs basis functions will be transformed into cubic B-spline basis. The spline \(\{E_{0},E_{0},\ldots,E_{N+1}\}\) forms a basis over the considered domain \([a,b]\). Coefficients of ExCuBs and its derivatives at different knots are given in Table 1.

Table 1 Coefficients of extended cubic basis \(E_{i}(x)\) and its derivatives at different knots
Full size table

3 Time discretization

The Caputo fractional derivative is employed to discretize the time of the given problem. Suppose \(t_{n}=n\tau \), \(n=0,1,\ldots,M\) in which \(\tau =\frac{T}{M}\) is the time step size. Forward finite difference technique is utilized for the discretization of Caputo fractional derivative. The Caputo derivatives \(\frac{\partial ^{\alpha }u(x,t)}{ \partial t^{\alpha }}\), \(\frac{\partial ^{2\alpha }u(x,t)}{\partial t ^{2\alpha }}\) of problem statement can be described as

$$\begin{aligned}& \frac{\partial ^{\alpha }u(x,t)}{\partial t^{\alpha }}=\frac{1}{\varGamma (1-\alpha )} \int _{0}^{t} \frac{\partial u(x,\tau )}{\partial \tau } \frac{d \tau }{(t-\tau )^{ \alpha }} , \end{aligned}$$
(5)
$$\begin{aligned}& \frac{\partial ^{2\alpha }u(x,t)}{\partial t^{2\alpha }}=\frac{1}{ \varGamma (2-2\alpha )} \int _{0}^{t} \frac{\partial ^{2} u(x,\tau )}{ \partial \tau ^{2}} \frac{d \tau }{(t-\tau )^{2\alpha -1}} , \end{aligned}$$
(6)

where \(0<\alpha <1\).

The discretized form of the Caputo derivative using a first order forward finite difference method [32] can be written as follows:

$$\begin{aligned} \frac{\partial ^{\alpha }(x,t_{n+1})}{\partial t^{\alpha }}=\frac{1}{ \varGamma (2-\alpha )}\sum ^{n}_{j=0}b^{\alpha }_{j} \frac{u(x,t_{n-j+1})-u(x,t _{n-j})}{\tau ^{\alpha }}+e^{n+1}_{\tau }, \end{aligned}$$
(7)

where \(b^{\alpha }_{j}=(j+1)^{1-\alpha }-j^{1-\alpha }\). The approximation of first order truncation error \(e^{n+1}_{\tau }\) bound is given in [33] as

$$ \bigl\vert e^{n+1}_{\tau } \bigr\vert \leq C\tau ^{2-\alpha }, $$
(8)

where C is a constant.

Discretization of the Caputo derivative using a second order forward finite difference method [16] can be written as

$$\begin{aligned} \frac{\partial ^{2\alpha }(x,t_{n+1})}{\partial t^{2\alpha }} =&\frac{1}{ \varGamma (3-2\alpha )}\sum^{n}_{j=0} \frac{u(x,t_{n-j+1})-2u(x,t_{n-j})+u(x,t _{n-j-1})}{\tau ^{2\alpha }} \\ &{}\times \bigl((j+1)^{2-2\alpha }-j^{2-2\alpha }\bigr)+r^{n+1}_{\tau }. \end{aligned}$$

The above equation can be rewritten as

$$\begin{aligned} \frac{\partial ^{2\alpha }(x,t_{n+1})}{\partial t^{2\alpha }} =&\frac{1}{ \varGamma (3-2\alpha )}\sum ^{n}_{j=0}b^{2\alpha }_{j} \frac{u(x,t_{n-j+1})-2u(x,t _{n-j})+u(x,t_{n-j-1})}{\tau ^{2\alpha }} \\ &{}+r^{n+1}_{\tau }, \end{aligned}$$
(9)

where \(b^{2\alpha }_{j}=(j+1)^{2-2\alpha }-j^{2-2\alpha }\). The approximation of second order truncation error \(r^{n+1}_{\tau }\) bound is given in [34] as

$$\begin{aligned} \bigl\vert r^{n+1}_{\tau } \bigr\vert \leq D\tau ^{2-\alpha }, \end{aligned}$$
(10)

where D is a constant. The following properties of coefficients \(b_{j}\) can easily be verified [29]:

  • \(b_{0}=1\)

  • \(b_{0}>b_{1}>b_{2}>\cdots>b_{j}\), \(b_{j}\rightarrow 0\) as \(j\rightarrow \infty \)

  • \(b_{j}>0\) for \(j=0,1,\ldots,n\)

  • \(\sum^{n}_{j=0}(b_{j}-b_{j+1})+b_{j+1}=(1-b_{1})+\sum^{n-1}_{s=1}(b _{j}-b_{j+1})+b_{n}=1\).

Substituting Eqs. (7) and (9) into Eq. (1), we obtain the following form of the time discretization:

$$\begin{aligned}& \beta _{1}\sum^{n}_{j=0}b^{2\alpha }_{j} \bigl[u(x,t_{n-j+1})-2u(x,t_{n-j})+u(x,t _{n-j-1}) \bigr] +\beta _{2}\sum^{n}_{j=0} \bigl[b^{\alpha }_{j}u(x,t_{n-j+1}) \\& \quad {}-u(x,t_{n-j})\bigr]+\mu u\bigl(x,t^{n+1}\bigr)=\nu \frac{\partial ^{2}u(x,t^{n+1})}{ \partial x^{2}}+f\bigl(x,t^{n+1}\bigr) \end{aligned}$$

where \(\beta _{1}=\frac{1}{\tau ^{2\alpha } \varGamma (3-2\alpha )}\) and \(\beta _{2}=\frac{2\lambda }{\tau ^{\alpha } \varGamma (2-\alpha )}\).

Suppose \(u^{n+1}=u(x,t_{n+1})\) and \(f^{n+1}=f(x,t_{n+1})\), the above equation can be rewritten as

$$\begin{aligned}& \beta _{1}\bigl(u^{n+1}-2u^{n}+u^{n-1} \bigr)+\beta _{2}\bigl(u^{n+1}-u^{n}\bigr)+ \beta _{1} \sum^{n}_{j=1}b^{2\alpha }_{j} \bigl(u^{n+1-j}-2u^{n-j} \\& \quad {}+u^{n-j-1}\bigr) +\beta _{2}\sum ^{n}_{j=1}b^{\alpha }_{j} \bigl(u^{n+1-j}-u^{n-j}\bigr)+ \mu u^{n+1} =\nu \frac{\partial ^{2}u^{n+1}}{\partial x^{2}}+f^{n+1} \end{aligned}$$
(11)

where \(n=0,1,\ldots,M\). It is noticed that \(u^{-1}\) will appear for \(j=0,n\). The initial velocity condition is used to calculate this term via a central difference formula. We obtain the following result:

$$ u^{-1}=u^{1}-2\tau f_{2}(x). $$
(12)

4 Description of technique

The approximated solution \(U(x,t)\) of given model using ExCuBs to the exact solution \(u(x,t)\) is described in the following form [35, 36]:

$$ U(x,t)=\sum^{N+1}_{i=-1}d_{i}(t)E_{i}(x, \lambda ), $$
(13)

where \(d_{i}(t)\) are the time dependent unknown coefficients which are to be required by some particular restrictions. Each subinterval \([x_{i},x_{i+1}]\) of basis function covers only three non-zero elements \(E_{i-1}\), \(E_{i}\), \(E_{i+1}\). The approximated solution \(u^{n}_{j}\) at the grid point \((x_{j},t_{n})\) at the nth time level to the exact solution is defined as

$$\begin{aligned} u^{n}_{j}=\sum^{i+1}_{j=i-1}d^{n}_{j}(t)E_{j}(x, \lambda ), \end{aligned}$$

where \(i=0,1,\ldots,N\). Using the above approximation and basis functions, the values \(u^{n}_{j}\) and their necessary derivatives up to second order as given below:

$$\begin{aligned}& u^{n}_{i} =c_{1}d^{n}_{i-1}+c_{2}d^{n}_{i}+c_{1}d^{n}_{i+1}, \\& (u_{x})^{n}_{i} =c_{3}d^{n}_{i+1}-c_{3}d^{n}_{i-1}, \\& (u_{xx})^{n}_{i} =c_{4}d^{n}_{i-1}-c_{5}d^{n}_{i}+c_{4}d^{n}_{i+1}, \end{aligned}$$

where \(c_{1}=\frac{4-\lambda }{24}\), \(c_{2}=\frac{8+\lambda }{12}\), \(c_{3}=\frac{1}{2h}\), \(c_{4}=\frac{2+\lambda }{2h^{2}}\), \(c_{5}=\frac{2+ \lambda }{h^{2}}\).

The Caputo derivatives and ExCuBs are used to discretize the model problem. Using the approximation and its derivatives in Eq. (11) and after simplification we obtain the recurrence relation in the following form:

$$\begin{aligned}& \bigl((\beta _{1}+\beta _{2}+\mu )c_{1}-\nu c_{4}\bigr)d^{n+1}_{j-1}+ \bigl((\beta _{1}+ \beta _{2}+\mu )c_{2}-\nu c_{5}\bigr)d^{n+1}_{j}+\bigl((\beta _{1}+\beta _{2}+ \mu )c_{1}-\nu c_{4}\bigr)d^{n+1}_{j+1} \\& \quad = \textstyle\begin{cases} (2\beta _{1}+\beta _{2})(c_{1}d^{n}_{j-1}+c_{2}d^{n}_{j}+c_{1}d^{n}_{j+1})- \beta _{1}(c_{1}d^{n-1}_{j-1}+c_{2}d^{n-1}_{j}+c_{1}d^{n-1}_{j+1}) \\ \quad {}-\beta _{1}\sum^{n}_{k=1}b^{2\alpha }_{k}[c_{1}(d^{n+1-k}_{j-1}-2d^{n-k} _{j-1}+d^{n-1-k}_{j-1})+c_{2}(d^{n+1-k}_{j}-2d^{n-k}_{j} \\ \quad {}+d^{n-1-k}_{j})+c_{1}(d^{n+1-k}_{j+1}-2d^{n-k}_{j+1}+d^{n-1-k}_{j+1})]- \beta _{2}\sum^{n}_{k=1}b^{\alpha }_{k}[c_{1}(d^{n+1-k}_{j-1} \\ \quad {}-d^{n-k}_{j-1}) +c_{2}(d^{n+1-k}_{j}-d^{n-k}_{j})+c_{1}(d^{n+1-k} _{j+1}-d^{n-k}_{j+1})]+f^{n+1}_{j}. \end{cases}\displaystyle \end{aligned}$$
(14)

The above system carries \((N+1)\times (N+3)\) dimensions. To solve the above system for unique solution we need two additional equations which will come from boundary conditions. Thus the system has \((N+3)\times (N+3)\) dimensions.

5 Initial case

In order to start the iterative process, it is necessary to find the initial vector \(d^{0}=[d^{0}_{0},d^{0}_{1},\ldots,d^{0}_{N}]\) which can be evaluated from initial conditions. We employ the initial condition with its derivatives explained below:

(i):

\(\frac{d }{dx}(u^{0}_{j})=\frac{d}{dx}(f_{1}(x_{j}))\), \(j=0,N\),

(ii):

\(u^{0}_{j}=f_{1}(x_{j})\), \(j=0,1,\ldots,N\).

Hence we obtain a system \(Au=b\), where A is the matrix of dimension \((N+3)\times (N+3)\) which can be written as

$$ A= \begin{bmatrix} \frac{-1}{2h} & 0 & \frac{1}{2h} & 0 & \ldots & \ldots & 0 \\ \frac{4-\lambda }{24} & \frac{8+\lambda }{12} & \frac{4-\lambda }{24} & 0 & \ldots & \ldots & 0 \\ 0 & \frac{4-\lambda }{24} & \frac{8+\lambda }{12} & \frac{4-\lambda }{24} & \ldots & \ldots & 0 \\ \vdots & \ldots & \ddots & \ddots & \ddots & \ldots & \vdots \\ \vdots & \ldots & \ldots & \ldots & \frac{4-\lambda }{24} & \frac{8+ \lambda }{12} & \frac{4-\lambda }{24} \\ 0& \ldots & \ldots & \ldots & \frac{-1}{2h} & 0 & \frac{1}{2h} \end{bmatrix} $$

and \(b=[f'_{1}(x_{0}),f_{1}(x_{0}),f_{1}(x_{1}),\ldots,f_{1}(x_{N}),f'_{1}(x _{N})]^{T}\).

6 Stability analysis

In this section, we use the Von Neumann stability analysis to investigate the stability of proposed scheme. Consider the growth factor in the form of a one Fourier mode as

$$ U^{n}_{j}=\xi ^{n}e^{i w h j}, $$
(15)

with no forcing term. Here \(i=\sqrt{-1}\), w and h are the mode number and the element size, respectively. We have

$$\begin{aligned}& \bigl[(\beta _{1}+\beta _{2}+\mu )c_{1}- \nu c_{4}\bigr]U^{n+1}_{j-1}+\bigl[(\beta _{1}+ \beta _{2}+\mu )c_{2}-\nu c_{5}\bigr]U^{n+1}_{j}\\& \qquad {} +\bigl[(\beta _{1}+\beta _{2}+ \mu )c_{1} -\nu c_{4}\bigr]U^{n+1}_{j+1} \\& \quad =(2\beta _{1}+\beta _{2})c_{1}U^{n}_{j-1}+(2 \beta _{1}+\beta _{2})c_{2}U^{n}_{j}+(2 \beta _{1}+\beta _{2})c_{1}U^{n} _{j+1} -\beta _{1}\bigl[c_{1}U^{n-1}_{j-1} \\& \qquad {}+c_{2}U^{n-1}_{j}+c_{1}U^{n-1}_{j} \bigr]-\beta _{1}\sum^{n}_{k=1}b^{2 \alpha }_{k} \bigl[c_{1}\bigl(U^{n-k+1}_{j-1}-2U^{n-k}_{j-1}+U^{n-k-1}_{j-1} \bigr) +c _{2}\bigl(U^{n-k+1}_{j} \\& \qquad {}-2U^{n-k}_{j}+U^{n-k-1}_{j} \bigr) +c_{1}\bigl(U^{n-k+1}_{j+1}-2U^{n-k}_{j+1}+U ^{n-k-1}_{j+1}\bigr)\bigr]-\beta _{2}\sum ^{n}_{k=1}b^{\alpha }_{k} \bigl[c_{1}\bigl(U^{n-k+1} _{j-1} \\& \qquad{} -U^{n-k}_{j-1}\bigr)+c_{2} \bigl(U^{n-k+1}_{j}-U^{n-k}_{j} \bigr) +c_{1}\bigl(U^{n-k+1} _{j+1}-U^{n-k}_{j+1} \bigr)\bigr]. \end{aligned}$$

The above equation shows a round off error equation. Consider Eq. (15) to be the solution, then the above equation becomes

$$\begin{aligned}& \bigl[(\beta _{1}+\beta _{2}+\mu )c_{1}- \nu c_{4}\bigr]\xi ^{n+1}e^{i w h (j-1)}+\bigl[( \beta _{1}+\beta _{2}+\mu )c_{2}-\nu c_{5}\bigr]\xi ^{n+1}e^{i w h j}\\& \qquad {} +\bigl[(\beta _{1}+\beta _{2} +\mu )c_{1}-\nu c_{4}\bigr]\xi ^{n+1}e^{i w h (j+1)} \\& \quad =(2\beta _{1}+\beta _{2})c_{1}\xi ^{n}e^{i w h (j-1)}+(2\beta _{1}+\beta _{2})c_{2}\xi ^{n}e ^{i w h j} +(2\beta _{1}+\beta _{2})c_{1}\xi ^{n}e^{i w h (j+1)}\\& \qquad {}-\beta _{1}\bigl[c_{1} \xi ^{n-1}e^{i w h (j-1)} +c_{2}\xi ^{n-1}e^{i w h j}+c_{1}\xi ^{n-1}e ^{i w h (j+1)}\bigr] \\& \qquad {}-\beta _{1}\sum^{n}_{k=1}b^{2\alpha }_{k} \bigl[c_{1}\bigl(\xi ^{n-k+1}-2\xi ^{n-k}+ \xi ^{n-k-1}\bigr) +c_{2}\bigl(\xi ^{n-k+1}-2\xi ^{n-k}+\xi ^{n-k-1}\bigr) \\& \qquad {}+c_{1}\bigl(\xi ^{n-k+1}-2\xi ^{n-k}+ \xi ^{n-k-1}\bigr)\bigr]\\& \qquad {}-\beta _{2}\sum ^{n}_{k=1}b ^{\alpha }_{k} \bigl[c_{1}\bigl(\xi ^{n-k+1}-\xi ^{n-k} \bigr)+c_{2}\bigl(\xi ^{n-k+1}-\xi ^{n-k}\bigr) \\& \qquad {}+c_{1}\bigl(\xi ^{n-k+1}-\xi ^{n-k} \bigr)\bigr]. \end{aligned}$$

Throughout dividing by \(e^{i w h j}\) and rearranging the terms, we obtain

$$\begin{aligned}& \bigl[(\beta _{1}+\beta _{2}+\mu ) \bigl(c_{2}+2c_{1} \cos (w h)\bigr)-\nu \bigl(c_{5}+2c _{4} \cos (w h)\bigr)\bigr]\xi ^{n+1} \\& \quad =(2\beta _{1}+\beta _{2}) \bigl(c_{2}+2c_{1} \cos (w h)\bigr)\xi ^{n}-\beta _{1}\bigl(c _{2}+2c_{1} \cos (w h)\bigr)\xi ^{n-1} \\& \qquad {}-\beta _{1}\bigl(c_{2}+2c_{1} \cos (w h)\bigr)\sum^{n}_{k=1}b^{2\alpha }_{k}[c _{1}\bigl(\xi ^{n-k+1}-2\xi ^{n-k}+\xi ^{n-k-1}\bigr) \\& \qquad {}-\beta _{2}\bigl(c_{2}+2c_{1} \cos (w h)\bigr)\sum^{n}_{k=1}b^{\alpha }_{k} \bigl[c _{1}\bigl(\xi ^{n-k+1}-\xi ^{n-k}\bigr) \bigr]. \end{aligned}$$

After some calculation, we get the following equality:

$$\begin{aligned} \xi ^{n+1} =&\frac{1}{\omega } \Biggl[(1+\alpha _{1})\xi ^{n}-\alpha _{1}\xi ^{n-1}-\alpha _{1}\sum^{n}_{k=1}b^{2\alpha }_{k} \bigl[\xi ^{n-k+1}-2\xi ^{n-k}+ \xi ^{n-k-1}\bigr] \\ &{}-\alpha _{2}\sum^{n}_{k=1}b^{\alpha }_{k} \bigl[\xi ^{n-k+1}-\xi ^{n-k}\bigr] \Biggr] \end{aligned}$$
(16)

where \(\alpha _{1}=\frac{\beta _{1}}{\beta _{1}+\beta _{2}}\), \(\alpha _{2}=\frac{ \beta _{2}}{\beta _{1}+\beta _{2}}\) and

$$ \omega =1+\frac{\mu }{\beta _{1}+\beta _{2}}+\frac{12\nu (2+\lambda ) \sin ^{2}w h/2}{h^{2}(\beta _{1}+\beta _{2}+\mu )(6+(\lambda -4)\sin ^{2}w h/2)}. $$

Obviously \(\omega \geq 1\), for all \(\lambda >-2\).

Proposition 6.1

Let \(\xi ^{n}, n=0,1,\ldots,T\times M\), be the solution of Eq. (16), we have

$$ \bigl\vert \xi ^{n} \bigr\vert \leq (1+\alpha _{1}) \bigl\vert \xi ^{0} \bigr\vert , \quad n=0,1,\ldots,T\times M $$
(17)

where \(\alpha _{1}\) is a positive constant.

Proof

We prove this proposition with the help of mathematical induction. For \(n=0\), in Eq. (16), we get the following relation:

$$ \xi ^{1}=\frac{1}{\omega }(1+\alpha _{1})\xi ^{0} \quad \Rightarrow\quad \bigl\vert \xi ^{1} \bigr\vert \leq (1+\alpha _{1}) \bigl\vert \xi ^{0} \bigr\vert ; \quad \omega \geq 1. $$

Suppose \(|\xi ^{n}| \leq (1+\alpha _{1})|\xi ^{0}|\) holds for \(n=0,1,\ldots,T \times (M-1)\), we have

$$\begin{aligned} \bigl\vert \xi ^{n+1} \bigr\vert =&\frac{1}{\omega }(1+ \alpha _{1}) \bigl\vert \xi ^{n} \bigr\vert - \frac{\alpha _{1}}{\omega } \bigl\vert \xi ^{n-1} \bigr\vert - \frac{\alpha _{1}}{\omega }\sum^{n}_{k=1}b ^{2\alpha }_{k} \bigl[ \bigl\vert \xi ^{n-k+1} \bigr\vert -2 \bigl\vert \xi ^{n-k} \bigr\vert \\ &{}+ \bigl\vert \xi ^{n-k-1} \bigr\vert \bigr]-\frac{\alpha _{2}}{\omega } \sum^{n}_{k=1}b^{ \alpha }_{k} \bigl[ \bigl\vert \xi ^{n-k+1} \bigr\vert - \bigl\vert \xi ^{n-k} \bigr\vert \bigr] \\ \leq & \frac{1}{\omega }(1+\alpha _{1})^{2} \bigl\vert \xi ^{0} \bigr\vert -\frac{\alpha _{1}(1+ \alpha _{1})}{\omega } \bigl\vert \xi ^{0} \bigr\vert -\frac{\alpha _{1}(1+\alpha _{1})}{ \omega }\sum ^{n}_{k=1}b^{2\alpha }_{k} \bigl[ \bigl\vert \xi ^{0} \bigr\vert \\ &{}-2 \bigl\vert \xi ^{0} \bigr\vert + \bigl\vert \xi ^{0} \bigr\vert \bigr]-\frac{\alpha _{2}(1+\alpha _{1})}{\omega }\sum ^{n}_{k=1}b^{\alpha }_{k} \bigl[ \bigl\vert \xi ^{0} \bigr\vert - \bigl\vert \xi ^{0} \bigr\vert \bigr] \\ \leq & (1+\alpha _{1})[1+\alpha _{1}-\alpha _{1}] \bigl\vert \xi ^{0} \bigr\vert , \\ \bigl\vert \xi ^{n+1} \bigr\vert \leq & (1+\alpha _{1}) \bigl\vert \xi ^{0} \bigr\vert . \end{aligned}$$

Thus \(|\xi ^{n+1}|=|U^{n+1}_{j}|\leq (1+\alpha _{1})|\xi ^{0}|=(1+\alpha _{1})|U^{0}_{j}|\), so that \(\|U^{n+1}_{j}\|_{2} \leq (1+\alpha _{1})\| \xi ^{0}\|_{2}\). Thus one concludes that the proposed numerical scheme is unconditionally stable. □

7 Convergence

In this part, we will investigate the convergence of proposed technique using the Lopez-Marcos [37] method, which plays a significant role in the theory of convergence analysis of fractional type equation. Here we take a few results and notations from [37]. Assume that \(\varOmega _{x}=\{x_{j}; 0\leq j\leq N\}\) and \(\varOmega _{t}=\{t_{n}; 0\leq n\leq M\}\) be the equidistant partitioning of intervals \([a,b]\) and \([0,T]\) with the step size h and τ, respectively. Consider \(u^{n}_{j}\) be the approximated solution at the grid point \((x_{j},t_{n})\) and \(V=\{v_{j}; 0\leq j\leq N\}\), \(W=\{w_{j}; 0\leq j\leq N\}\) be the two functions defined on \(\varOmega _{x}\). We define difference notation as follows:

$$\begin{aligned}& \delta ^{2}V=v_{j+1}-2v_{j}+v_{j-1}, \qquad \delta V=v_{j+1}-v_{j} \\& \Vert V \Vert ^{2}=(V,V), \qquad (V,W)=\sum ^{N}_{j=1}h v_{j} w_{j}, \\& (V_{xx},V)=-(V_{x},V_{x}), \qquad (V, W_{x})=-(V_{x},W). \end{aligned}$$

We also suppose that \(u_{tt}\), \(u_{xxxx}\) are continuous over the intervals \([a,b]\) and \([0,T]\), and that there is a positive constant F, such that

$$\begin{aligned} \vert u_{tt} \vert \leq F, \qquad \vert u_{xxxx} \vert \leq F. \end{aligned}$$
(18)

The above values are different at different locations and independent of j, n, h, τ and for \((x,t)\in \varOmega _{x}\times \varOmega _{t}\).

Proposition 7.1

([37])

Let \(\{z_{0}, z_{1},\ldots,z_{n},\ldots\}\) be a monotonically decreasing sequence with the properties \(z_{n}\geq 0\) and \(z_{n+1}+z_{n-1}\geq 2 z_{n}\). Then, for any positive integer S and for each vector \((v_{1},v_{2},\ldots,v_{S})\) with S real entries, we have

$$ \sum^{S-1}_{n=0} \Biggl(\sum^{n}_{m=0}z_{m}V_{n+1-m} \Biggr)V_{n+1} \geq 0. $$
(19)

So for the proposed scheme, we have

$$\begin{aligned}& \beta _{1}\sum^{n}_{j=0}b^{2\alpha }_{j} \bigl(u(x,t_{n-j+1})-2u(x,t_{n-j})+u(x,t _{n-j-1}) \bigr) +\beta _{2}\sum^{n}_{j=0}b^{\alpha }_{j} \bigl(u(x,t_{n-j+1}) \\& \quad {}-u(x,t_{n-j})\bigr) +\mu u\bigl(x,t^{n+1}\bigr) = \nu \frac{\partial ^{2}u(x,t^{n+1})}{ \partial x^{2}}+f\bigl(x,t^{n+1}\bigr)+O\bigl(\tau ^{2-\alpha }+h^{2}\bigr) \end{aligned}$$
(20)

and

$$\begin{aligned}& \beta _{1}\sum^{n}_{j=0}b^{2\alpha }_{j} \bigl(u^{n-j+1}-2u^{n-j}+u^{n-j-1}\bigr) +\beta _{2}\sum^{n}_{j=0}b^{\alpha }_{j} \bigl(u^{n-j+1}-u^{n-j}\bigr) +\mu u^{n+1} \\& \quad =\nu \frac{\partial ^{2}u^{n+1}}{\partial x^{2}}+f^{n+1} \end{aligned}$$
(21)

where \(u^{n}_{j}\) and \(u(x_{j},t_{n})\) are an approximated and the exact solution at point \((x_{j},t_{n})\), respectively.

Theorem 1

Suppose that \(u(x,t)\) and \(u^{n}_{j}\) be the solutions of given model and Eq. (20), respectively, and \(u(x,t)\) satisfies the smoothness condition (18), then we have

$$ \bigl\Vert E^{n+1} \bigr\Vert \leq O\bigl(\tau ^{2-\alpha }+h^{2}\bigr), $$
(22)

where \(E^{n+1}_{j}=u(x_{j},t_{n+1})-u^{n+1}_{j}\), for every \(t\geq 0\) and suitably small h and τ.

Proof

Subtract Eq. (20) from Eq. (21), we get the error equation as follows:

$$ \beta _{1}\sum^{n}_{k=0}b^{2\alpha }_{k} \delta ^{2}E^{n+1-k}_{j}+\beta _{2} \sum^{n}_{k=0}b^{\alpha }_{k} \delta E^{n+1-k}_{j}+\mu E^{n+1}_{j} =\nu \bigl(E^{n+1}_{j}\bigr)+O\bigl(\tau ^{2-\alpha }+h^{2}\bigr). $$

Consider \(p^{n+1}_{j}=O(\tau ^{2-\alpha }+h^{2})\) then multiply by \(h E^{n+1}_{j}\) on both sides of the above equation and sum up the terms for j, which varies from 1 to N, we obtain

$$\begin{aligned} \bigl\Vert E^{n+1} \bigr\Vert ^{2} =&- \frac{\beta _{1}}{\mu }\sum^{n}_{k=0}b^{2\alpha } _{k}\bigl(\delta ^{2}E^{n+1-k},E^{n+1} \bigr) -\frac{\beta _{2}}{\mu }\sum^{n}_{k=0}b ^{\alpha }_{k}\bigl(\delta E^{n+1-k},E^{n+1} \bigr) \\ &{}+\frac{1}{\mu }\bigl(E^{n+1},E^{n+1} \bigr)_{xx}+\frac{1}{\mu }\bigl(p^{n+1},E^{n+1} \bigr) \\ =&-\frac{\beta _{1}}{\mu }\sum^{n}_{k=0}b^{2\alpha }_{k} \bigl(\delta ^{2}E ^{n+1-k},E^{n+1}\bigr) - \frac{\beta _{2}}{\mu }\sum^{n}_{k=0}b^{\alpha } _{k}\bigl(\delta E^{n+1-k},E^{n+1}\bigr) \\ &{}-\frac{1}{\mu }\bigl(\bigl(E^{n+1}\bigr)_{x}, \bigl(E^{n+1}\bigr)_{x}\bigr)+\frac{1}{\mu } \bigl(p^{n+1},E ^{n+1}\bigr), \\ \bigl\Vert E^{n+1} \bigr\Vert ^{2} =&- \frac{\beta _{1}}{\mu }\sum^{n}_{k=0}b^{2\alpha } _{k}\bigl(\delta ^{2}E^{n+1-k},E^{n+1} \bigr) -\frac{\beta _{2}}{\mu }\sum^{n}_{k=0}b ^{\alpha }_{k}\bigl(\delta E^{n+1-k},E^{n+1} \bigr) \\ &{}-\frac{1}{\mu } \bigl\Vert \bigl(E^{n+1} \bigr)_{x} \bigr\Vert ^{2}+\frac{1}{\mu } \bigl(p^{n+1},E^{n+1}\bigr). \end{aligned}$$

After repositioning the terms, we obtain

$$\begin{aligned}& \bigl\Vert E^{n+1} \bigr\Vert ^{2} + \frac{\beta _{1}}{\mu }\sum^{n}_{k=0}b^{2\alpha } _{k}\bigl(\delta ^{2}E^{n+1-k},E^{n+1} \bigr) +\frac{\beta _{2}}{\mu }\sum^{n}_{k=0}b ^{\alpha }_{k}\bigl(\delta E^{n+1-k},E^{n+1} \bigr) \\& \quad {}+\frac{1}{\mu } \bigl\Vert \bigl(E^{n+1} \bigr)_{x} \bigr\Vert ^{2}=\frac{1}{\mu } \bigl(p^{n+1},E^{n+1}\bigr). \end{aligned}$$

Since \(\frac{1}{\mu }\|(E^{n+1})_{x}\|^{2}\geq 0\), we get

$$\begin{aligned} & \bigl\Vert E^{n+1} \bigr\Vert ^{2} + \frac{\beta _{1}}{\mu }\sum^{n}_{k=0}b^{2\alpha } _{k}\bigl(\delta ^{2}E^{n+1-k},E^{n+1} \bigr) +\frac{\beta _{2}}{\mu }\sum^{n}_{k=0}b ^{\alpha }_{k}\bigl(\delta E^{n+1-k},E^{n+1} \bigr) \\ &\quad \leq \frac{1}{\mu }\bigl(p^{n+1},E^{n+1} \bigr). \end{aligned}$$

The remaining terms are

$$\begin{aligned}& \bigl\Vert E^{n} \bigr\Vert ^{2}+ \frac{\beta _{1}}{\mu }\sum^{n-1}_{k=0}b^{2\alpha }_{k} \bigl( \delta ^{2}E^{n-k},E^{n}\bigr) + \frac{\beta _{2}}{\mu }\sum^{n-1}_{k=0}b^{ \alpha }_{k} \bigl(\delta E^{n-k},E^{n}\bigr)\leq \frac{1}{\mu } \bigl(p^{n},E^{n}\bigr), \\& \bigl\Vert E^{n-1} \bigr\Vert ^{2}+ \frac{\beta _{1}}{\mu }\sum^{n-2}_{k=0}b^{2\alpha } _{k}\bigl(\delta ^{2}E^{n-k-1},E^{n} \bigr) +\frac{\beta _{2}}{\mu }\sum^{n-2}_{k=0}b ^{\alpha }_{k}\bigl(\delta E^{n-k-1},E^{n-1} \bigr)\leq \frac{1}{\mu }\bigl(p^{n-1},E ^{n-1} \bigr), \\& \vdots \\& \bigl\Vert E^{2} \bigr\Vert ^{2}+ \frac{\beta _{1}}{\mu }\sum^{1}_{k=0}b^{2\alpha }_{k} \bigl( \delta ^{2}E^{2-k},E^{2}\bigr) + \frac{\beta _{2}}{\mu }\sum^{1}_{k=0}b^{ \alpha }_{k} \bigl(\delta E^{2-k},E^{2}\bigr)\leq \frac{1}{\mu } \bigl(p^{2},E^{2}\bigr), \\& \bigl\Vert E^{1} \bigr\Vert ^{2}+ \frac{\beta _{1}}{\mu }\sum^{0}_{k=0}b^{2\alpha }_{k} \bigl( \delta ^{2}E^{1-k},E^{1}\bigr) + \frac{\beta _{2}}{\mu }\sum^{0}_{k=0}b^{ \alpha }_{k} \bigl(\delta E^{1-k},E^{1}\bigr)\leq \frac{1}{\mu } \bigl(p^{1},E^{1}\bigr). \end{aligned}$$

Taking the sum of all the above inequalities, we have

$$\begin{aligned}& \sum^{n}_{k=0} \bigl\Vert E^{n+1} \bigr\Vert ^{2}+\frac{\beta _{1}}{\mu }\sum ^{n}_{m=0}\sum ^{m}_{k=0}b^{2\alpha }_{k} \bigl(\delta ^{2} E^{m+1-k},E^{m+1}\bigr) + \frac{\beta _{2}}{\mu }\sum^{n}_{m=0}\sum ^{m}_{k=0}b^{\alpha }_{k} \bigl(\delta E^{m+1-k},E ^{m+1}\bigr) \\& \quad \leq \frac{1}{\mu }\sum^{n}_{k=0} \bigl(p^{k+1},E^{k+1}\bigr). \end{aligned}$$
(23)

By Proposition 7.1, we can deduce that

$$ \frac{\beta _{1}}{\mu }\sum^{n}_{m=0}\sum ^{m}_{k=0}b^{2\alpha }_{k} \bigl( \delta ^{2} E^{m+1-k},E^{m+1}\bigr)\geq 0 $$

and

$$ \frac{\beta _{2}}{\mu }\sum^{n}_{m=0}\sum ^{m}_{k=0}b^{\alpha }_{k} \bigl( \delta E^{m+1-k},E^{m+1}\bigr)\geq 0. $$

Therefore Eq. (23) can be written as

$$ \sum^{n}_{k=0} \bigl\Vert E^{n+1} \bigr\Vert ^{2}\leq \frac{1}{\mu }\sum ^{n}_{k=0}\bigl(p^{k+1},E ^{k+1}\bigr).\vadjust{\goodbreak} $$

Hence

$$ \bigl\Vert E^{n+1} \bigr\Vert ^{2}\leq \frac{1}{\mu }\bigl(p^{k+1},E^{k+1}\bigr). $$

Using the Cauchy–Schwarz inequality, we get

$$ \bigl\Vert E^{n+1} \bigr\Vert ^{2}\leq \frac{1}{\mu }\bigl(p^{k+1},E^{k+1}\bigr) \leq \frac{1}{ \mu } \bigl\Vert p^{k+1} \bigr\Vert \bigl\Vert E^{k+1} \bigr\Vert . $$

From the above inequality, we can get the desired result,

$$ \bigl\Vert E^{n+1} \bigr\Vert \leq O\bigl(\tau ^{2-\alpha }+h^{2}\bigr). $$

 □

8 Numerical examples and discussions

In this section, some numerical experiments are discussed to demonstrate the feasibility of the proposed method. The calculated error norms are established by absolute \(L^{\infty }\) and Euclidean \(L^{2}\) norms, i.e.,

$$\begin{aligned}& L^{\infty } = \bigl\Vert U(x_{i},t)-u(x_{i},t) \bigr\Vert _{\infty }=\max_{0\leq i\leq N} \bigl\vert u(x _{i}, t)-u(x_{i},t) \bigr\vert , \\& L^{2} = \bigl\Vert U(x_{i},t)-u(x_{i},t) \bigr\Vert _{2}=\sqrt{h\sum^{N}_{i=0} \bigl\vert u(x_{i}, t)-u(x_{i},t) \bigr\vert ^{2}}. \end{aligned}$$

The following formula can be used to calculate the order of convergence [38] numerically:

$$ \mathrm{Order}=\frac{\log (L^{\infty }(N_{i}))-\log (L^{\infty }(N_{i+1}))}{ \log (N_{i+1})-\log (N_{i})}, $$

where \(L^{\infty }(N_{i})\) and \(L^{\infty }(N_{i+1})\) are the absolute errors at number of partitioning \(N_{i}\) and \(N_{i+1}\), respectively.

Problem 1

Consider the TFTE of the form

$$ \frac{\partial ^{2\alpha }u(x,t)}{\partial t^{2\alpha }}+\frac{ \partial ^{\alpha }u(x,t)}{\partial t^{\alpha }}= \frac{\partial ^{2}u(x,t)}{ \partial x^{2}}+f(x,t), $$
(24)

where \(f(x,t)\) is suitable with the exact solution \(u(x,t)=t^{2+ \alpha } \sin (2\pi x)\) [16], for all \((x,t)\in [0,1] \times [0,1]\), where the initial and boundary conditions are

$$ \textstyle\begin{cases} f_{1}(x)=0, \\ f_{2}(x)=0, \\ g_{1}(t)=0, \\ g_{2}(t)=0. \end{cases} $$
(25)

In Table 2, we calculate the \(L^{2}\)-norm for different spatial and temporal step size \(h=5\), \(\tau = \frac{1}{M}\), \((M=20, 40, 80)\). In Table 3, we determine the order of convergence [16, 29] from the computed data and present maximum absolute errors at different space-time step sizes. We give \(L^{2}\)-norm and maximum errors for \(\alpha =0.6, 0.7, 0.8,0.9\). Table 4 shows the maximum absolute error, the Euclidean norm, the order of convergence and the CPU time for \(\alpha =0.75\) and \(\tau =\frac{1}{100}\).

Table 2 A comparison of \(L^{2}\) norm with \(h=5\tau =\frac{1}{M}\) at \(T=1\) for Problem 1
Full size table
Table 3 A comparison of maximum error with \(h=5\tau =\frac{1}{M}\) at \(T=1\) for Problem 1
Full size table
Table 4 Maximum absolute errors and Euclidean norm \((L^{2})\) for Problem 1
Full size table

Problem 2

Consider the TFTE with \(\lambda =\frac{1}{2}\), \(\mu =0\) and \(\nu =\frac{1}{2}\):

$$ \frac{\partial ^{2\alpha }u(x,t)}{\partial t^{2\alpha }}+\frac{ \partial ^{\alpha }u(x,t)}{\partial t^{\alpha }}= \frac{1}{2} \frac{ \partial ^{2}u(x,t)}{\partial x^{2}}+f(x,t) $$
(26)

having the initial condition and boundary conditions

$$ \textstyle\begin{cases} f_{1}(x)=0, \\ f_{2}(x)=0, \\ g_{1}(t)=t^{2}, \\ g_{2}(t)=et^{2}, \end{cases} $$
(27)

where \(f(x,t)\) is appropriate for the exact solution \(u(x,t)=t^{2}e ^{x}\) [39] with \(0\leq t\leq 1\), \(0\leq x\leq 1\) and \(0<\alpha <1\).

The comparisons of maximum absolute errors are demonstrated in Table 5 for \(\alpha =0.64, 0.80, 0.96\). Here we choose step sizes \(h=\frac{1}{N}\) and \(\tau =\frac{1}{N^{2}}\) for \(N=4, 8, 12\). The \(L^{2}\)-norm and order of convergence are summarized in Table 6. In Table 7, we choose \(\alpha =0.7\), \(T=1\) and calculate the maximum absolute error, the Euclidean norm and the order of convergence. Figure 1 shows the two dimensional comparison graph of exact and approximated solutions for \(\alpha =0.7\), \(N=40\) and \(\tau =200\). Figure 2 depicts the error plot for Problem 1 and Problem 2 at different time levels. A high accuracy of the proposed method for these problems can be visualized in Fig. 2. Compatibility of exact and approximated solution for Problem 1 and Problem 2 can be viewed in Fig. 3 and Fig. 4, respectively.

Figure 1
figure 1

Comparison plots of exact solutions and approximated solutions for Problem 1 and Problem 2, respectively

Full size image
Figure 2
figure 2

Error plot at different time level for Problem 1 and Problem 2, respectively

Full size image
Figure 3
figure 3

3D plot for the exact and approximated solution of Problem 1

Full size image
Figure 4
figure 4

3D plot for the exact and approximated solution of Problem 2

Full size image
Table 5 A comparison of maximum absolute error with \(h=\frac{1}{N}\), \(\tau =\frac{1}{N^{2}}\) at \(T=1\) for Problem 2
Full size table
Table 6 A comparison of \(L^{2}\)-norm with \(h=\frac{1}{N}\), \(\tau =\frac{1}{N^{2}}\) at \(T=1\) for Problem 2
Full size table
Table 7 Maximum absolute errors and Euclidean norm \((L^{2})\) at \(T=1\) for Problem 2
Full size table

Problem 3

Consider TFTE of the form

$$ \frac{\partial ^{2\alpha }u(x,t)}{\partial t^{2\alpha }}+20\frac{ \partial ^{\alpha }u(x,t)}{\partial t^{\alpha }}+25u(x,t) =\frac{ \partial ^{2}u(x,t)}{\partial x^{2}}+f(x,t) $$
(28)

with the initial and boundary conditions

$$ \textstyle\begin{cases} f_{1}(x)=\sin (x), \\ f_{2}(x)=0, \\ g_{1}(t)=0, \\ g_{2}(t)=\cos (t)\sin (1),\vadjust{\goodbreak} \end{cases} $$
(29)

where \(f(x,t)=-20\sin (x)\sin (t)+25\sin (x)\cos (t)\) is relevant with the exact solution \(u(x,t)=\sin (x)\cos (t)\) [1], \(0\leq t\leq 1\), \(0\leq x\leq 1\) and \(0<\alpha <1\).

In Table 8, we present a comparison of maximum absolute errors at \(T=0.5\), for \(\alpha =0.925, 0.975\) with the results given by Asgari et al. [1]. In Table 9, we select \(h=10\tau =\frac{1}{M}\) at \(T=0.5\) for \(\alpha =0.9\) and present maximum norm, Euclidean norm and order of convergence.

Table 8 A comparison of maximum absolute errors at \(T=0.5\) for Problem 3
Full size table
Table 9 The \(L^{\infty }\)-norm and \(L^{2}\)-norm with \(h=10\tau =\frac{1}{M}\) at \(T=0.5\) for Problem 3
Full size table

Problem 4

Consider TFTE of the form

$$\begin{aligned}& \frac{\partial ^{2\alpha }u(x,t)}{\partial t^{2\alpha }}+40\frac{ \partial ^{\alpha }u(x,t)}{\partial t^{\alpha }}+100u(x,t) =\frac{ \partial ^{2}u(x,t)}{\partial x^{2}}+f(x,t) \end{aligned}$$
(30)
$$\begin{aligned}& \textstyle\begin{cases} f_{1}(x)=\sinh (x), \\ f_{2}(x)=-2\sinh (x), \\ g_{1}(t)=0, \\ g_{2}(t)=e^{-2t}\sinh (1), \end{cases}\displaystyle \end{aligned}$$
(31)

where \(f(x,t)=23e^{-2t}\sinh (x)\) is the appropriate forcing term with the exact solution \(u(x,t)=e^{-2t}\sinh (x)\) [1], \(0\leq t\leq 1\), \(0\leq x\leq 1\) and \(0<\alpha <1\).

Table 10 shows the maximum absolute errors of knots for time \(T=0.5\), \(\alpha =0.975\) and \(N=10\tau =\frac{1}{10}\) with the results given by Asgari et al. [1] for different degree polynomials. The \(L^{\infty }\)-norm, \(L^{2}\)-norm for different time step τ and corresponding rate of convergence are demonstrated in Table 11. Figure 5 depicts the error plots at different time levels for Problem 3 and Problem 4, respectively. The compatibility between the exact and the approximated solutions for the Problem 3 and Problem 4 can be seen in Fig. 6 and Fig. 7, respectively.

Figure 5
figure 5

Error plot at different time level for Problem 3 and Problem 4 respectively

Full size image
Figure 6
figure 6

3D plot for the exact and approximated solution of Problem 3

Full size image
Figure 7
figure 7

3D plot for the exact and approximated solution of Problem 4

Full size image
Table 10 A comparison of maximum absolute errors at \(T=0.5\) for Problem 4
Full size table
Table 11 The \(L^{\infty }\)-norm and \(L^{2}\)-norm with \(h=2\tau =\frac{1}{M}\) at \(T=0.5\) for Problem 4
Full size table

It is simple to notice that convergence rate obtained by the present method is compatible with the theoretical results. The proposed method needs a small storage and less CPU time, which shows the simplicity and strength of the proposed scheme. It is concluded that the present scheme has a great capacity to deal with the fractional order partial differential equations.

9 Conclusion

A finite difference scheme based on a combination of the ExCuBs method and Caputo’s fractional derivative for the numerical solution of TFTE has been presented. The proposed method is investigated and good compatibility was found with the exact solution. The proposed scheme is convergent of order \(O(\tau ^{2-\alpha }+h^{2})\) and unconditionally stable. Numerical experiments have been conducted using the proposed method and graphs show the feasibility and accuracy of the method. Studies on the order of convergence were also carried out in a numerical way.