A second-order numerical scheme for the time-fractional partial differential equations with a time delay

Abstract

This work proposes a numerical scheme for a class of time-fractional convection–reaction–diffusion problems with a time lag. Time-fractional derivative is considered in the Caputo sense. The numerical scheme comprises the discretization technique given by Crank and Nicolson in the temporal direction and the spline functions with a tension factor are used in the spatial direction. Through the von Neumann stability analysis, the scheme is shown conditionally stable. Moreover, a rigorous convergence analysis is presented through the Fourier series. Two test problems are solved numerically to verify the effectiveness of the proposed numerical scheme.

1 Introduction

Though fractional calculus has a long history and is as old as classical calculus, people did not apply it to the problems arising in science and engineering for a long time. However, having wide applications, fractional calculus began to attract the increasing attention of researchers and scientists in the last few decades. To the best of the authors’ knowledge, Leibnitz was the first who has introduced the non-integer derivatives in 1695. However, Riemann and Liouville were the first who build the foundation of fractional calculus. The first book published in 1974 by Oldham and Spanier (1974) is solely devoted to fractional calculus. They have observed that the derivatives/integrals of fractional orders are more useful than the integer derivatives/integrals. It motivates them towards the fundamental mathematical properties of the fractional differential/integral operators and the applications of arbitrary order operators. Later, several authors such as Kilbas et al. (2006), Li and Zeng (2015), Miller and Ross (1993), and Podlubny (1999) have published several books on the theory and applications of fractional integrals/derivatives. For a good survey paper on the application of fractional derivatives in modern mechanics, the readers are referred to Xu and Tan (2006).

Fractional calculus has many applications in various branches of mathematics, science, and engineering, such as in the oil industry, neural networks, mass transport in fluids, control systems, finance, viscoelasticity, rheology, anomalous transport in biology, fluid flow in porous mediums, underground environmental problems, and acoustics (Bagley and Torvik (1983); Benson et al. (2000, 2001); Fitt et al. (2009); Galue et al. (2007); Hatano and Hatano (1998); Höfling and Franosch (2013); Machado (2001); Oeser and Freitag (2009); Raberto et al. (2002); Scalas et al. (2000); Arafa et al. (2013)). One cannot model these problems appropriately without using fractional calculus. As a result, the theories and applications of the fractional differential and integral operators are well established. Their applicabilities to the problems mentioned above are being considered extensively. It has been shown that the fractional-order models are better than the integer models for the complex systems with memory and genetic effects. The problems with fractional derivatives are mainly classified into two categories: space-fractional and time-fractional. This paper considers the problems with the fractional derivative falling in the latter category.

In contrast to partial differential equations (PDEs), the physical problems for which the evolution not only depends on the present state of the system but also on the history are modeled by the delay partial differential equations (DPDEs) and time-fractional delay partial differential equations (TF-DPDEs). Often, the future state is determined by their history in the systems with subdiffusive processes. One cannot ignore the delay terms despite the whole history being considered by the fractional differential equations (FDEs). These equations frequently arise in many areas of science and engineering, such as in time to maturity and incubation time, delayed feedback, time to transport, and the time lag for getting information.

The three well-known derivatives of fractional order are widely used in the literature, namely, Caputo derivative, Riemann–Liouville derivative, and the Grünwald–Letnikov derivative. However, recently two new fractional derivatives, namely, Atangana-Baleanu–Caputo derivative (Aliyu et al. (2018)), and Caputo–Fabrizio derivative (Owolabi and Atangana (2017)) have been introduced. For reference, we define the Caputo and Riemann-Liouville derivatives as follows:

Definition 1.1

(Caputo fractional derivative) The Caputo fractional derivative of order $\gamma$ (denoted by ${}_0^CD_t^{\gamma }$) of w(t) is defined as

$$\begin{aligned}{}_0^CD_t^{\gamma }w(t) = {\left\{ \begin{array}{ll} \displaystyle \frac{1}{\Gamma (n-\gamma )}\int _0^t\frac{w^{(n)}(s)}{(t-s)^{\gamma -n+1}}\,ds, &{} n-1\le \gamma < n,\,n\in {\mathbb {N}},\\ w^{(n)}(t), &{} \gamma =n\in {\mathbb {N}}. \end{array}\right. }\end{aligned}$$

Definition 1.2

(Riemann-Liouville derivative) The Riemann-Liouville derivative (denoted by ${}_0^{RL}D_t^{\gamma }$) of order $\gamma$ of w(t) is defined as

$$\begin{aligned} {}_0^{RL}D_t^{\gamma }w(t) = {\left\{ \begin{array}{ll} \displaystyle \frac{1}{\Gamma (n-\gamma )}\frac{d^n}{dt^n}\int _0^t\frac{w(s)}{(t-s)^{\gamma -n+1}}\,ds, &{} n-1 < \gamma \le n,\,n\in {\mathbb {N}},\\ w^{(n)}(t), &{} \gamma =n\in {\mathbb {N}}, \end{array}\right. } \end{aligned}$$

where $\Gamma$ represents the gamma function and defined for a complex number z with non-negative real part as

$$\begin{aligned} \Gamma (z)=\int _{0}^{\infty }e^{-t}t^{z-1}\,dt. \end{aligned}$$

Fractional partial differential equations (FPDEs) are where a fractional derivative replaces the classical derivative. Because of the precise and powerful descriptions of a large variety of natural phenomena, mathematicians and scientists are analyzing FPDEs with delay analytically and numerically (Yan and Kou (2012); Ouyang (2011); Zhou et al. (2009); Zhang et al. (2017); Sakara et al. (2016); Rodríguez et al. (2012); Mohebbi (2019); Nicaise and Pignotti (2006)). The models containing the time delay can be seen in the automatic control systems (Lazarević and Spasić (2009); Si-Ammour et al. (2009)), random walk (Ohira and Milton (2009)), and the modeling human immunodeficiency virus (HIV) infection of $CD4^+\,T$-cells (Culshaw et al. (2003); Hao et al. (2016); Yan and Kou (2012)). Furthermore, for the fractional models on HIV and corona virus disease, the readers are referred to the recent articles Chen et al. (2021); Chu et al. (2021); Pandey et al. (2021). An extensive study of the delay differential equations (DDEs) in the context of ordinary differential equations (ODEs) can be seen in Bellen and Zennaro (2003); Guglielmi (2006); Guglielmi and Hairer (2004, 2001); Rihan (2000) and the references therein. Furthermore, the analytical solutions, including the spectral methods and the integral transform methods such as Laplace transforms and Mellin transforms for FPDEs, can be seen in Langlands et al. (2011); Jiang et al. (2012); Kilbas et al. (2006); Ding and Jiang (2013); Ding and Nieto (2015, 2017). On the other hand, the DPDEs and TF-DPDEs are less explored than DDEs and FPDEs.

Because the evolution of a dependent variable of TF-DPDEs at any time t depends not only on its value at $t-\tau$ (for some time delay $\tau$) but also on all previous solutions, it is a difficult task to solve TF-PDEs with delay effectively and accurately. However, some analytical methods have been investigated to find the solutions to these problems. For instance, to obtain the actual solution to the TF-DPDEs, Prakash et al. (2020) proposed the invariant subspace approach. For the existence of solutions to DDEs, the readers are referred to Lakshmikantham (2008). Moreover, the asymptotic properties of the solution can be found in Hale and Verduyn Lunel (1993); Küchler and Mensch (1992). Furthermore, the stability bounds on the solution for the class of fractional delay difference equations (FDDEs) can be seen in Chen and Moore (2002); Deng et al. (2007); Hwang and Cheng (2006). For the finite-time stability of robotic systems where a time delay appears in $PD^{\gamma }$ fractional control system, the readers are referred to Lazarević (2006); Lazarević and Spasić (2009). As a consequence, we look for numerical methods for TF-DPDEs. Some numerical methods have also been developed for these problems; for instance, for a numerical scheme for space-fractional diffusion equations with a delay in time, the readers are referred to Hao et al. (2016). They have used the Taylor series expansion to linearize the nonlinear term and weighted shifted Grünwald–Letnikov formula to approximate the space-fractional derivative. Rihan (2009) extended the $\theta$ method to solve the TF-DPDEs of parabolic type in the Caputo sense. To test the BIBO stability of the system of TF-DDEs, Hwang and Cheng (2006) presented an effective numerical algorithm based on Cauchy’s integral theorem. Mohebbi (2019) constructed an efficient numerical method to solve the TF-DPDEs of convection-reaction-diffusion type with a nonlinear source term. The Chebyshev spectral collocation method and second-order finite difference approximations are used in spatial and temporal directions. Sakara et al. (2016) presented a homotopy perturbation method for nonlinear TF-DPDEs. Zhang et al. (2017) proposed a linearized compact finite difference method (FDM) for the semilinear TF-DPDEs. We use the compact finite difference approximation and $L_1$ approximation to discretize the spatial and temporal derivatives. Through rigorous analysis, the method is shown fourth-order convergent in space and of order $2-\gamma$ in time. Chen et al. (2020) suggested an approach comprising the use of global radial basis functions and a collocation method to solve a class of time-delay fractional optimal control problems. Chen et al. (2021) proposed a direct scheme comprising the use of a set of Dickson polynomials as basis functions and a collocation method to solve a class of time-delay fractional optimal control problems.

Motivated by the initial-boundary value problems for singularly perturbed delay parabolic PDEs (Kumar and Kumari (2020)), and the model considered by Kumar et al. (2021) for the singularly perturbed fractional delay PDEs, we analyze and develop an efficient numerical algorithm for the regular perturbation case. Furthermore, the following time-fractional diffusion model with a time-lag considered by Rihan (2009)

$$\begin{aligned}\frac{\partial ^\gamma y(x,t)}{\partial t^{\gamma }}=a\frac{\partial ^2 y(x,t)}{\partial x^2}+by(x,t-\tau ),\end{aligned}$$

motivates us to extend this model for the convection-reaction-diffusion fractional order PDEs.

In this work, we consider the following TF-DPDEs:

$$\begin{aligned}&\frac{\partial ^\gamma y(x,t)}{\partial t^{\gamma }}-\frac{\partial ^2 y(x,t)}{\partial x^2}+b(x,t)\frac{\partial y(x,t)}{\partial x}+c(x,t)y(x,t)\nonumber \\&\quad = -d(x,t)y(x,t-\tau )+f(x,t), \end{aligned}$$

(1.1a)

for all $(x,t)\in {\mathfrak {D}}\equiv \varOmega \times \varLambda =(0,1)\times (0,{\mathfrak {T}}]$, with the following interval and boundary conditions

$$\begin{aligned} y(x,t) = \psi _b (x,t),\quad (x,t)\in [0,1]\times [-\tau ,0], \end{aligned}$$

(1.1b)

$$\begin{aligned} y(0,t) = \psi _l (t), \quad y(1,t) = \psi _r(t), \quad t\in [0,{\mathfrak {T}}], \end{aligned}$$

(1.1c)

where the fractional derivative of order $\gamma \in (0,1)$ is estimated in the Caputo sense given by

$$\begin{aligned} \frac{\partial ^\gamma y(x,t)}{\partial t^\gamma } =\frac{1}{\Gamma (1-\gamma )}\int _0^t(t-s)^{-\gamma }\frac{\partial y(x,s)}{\partial t}\,ds. \end{aligned}$$

A hybrid scheme and a compact FDM have been investigated in Marzban and Tabrizidooz (2012); Li et al. (2015) for the classical integer order DPDEs that is the simplified form of the Eq. (1.1) for $\gamma =1$. Several numerical schemes have been developed when (1.1) is equivalent to the wave equation (for $\gamma = 2$) with delay, see e.g. Nicaise and Pignotti (2006); Garrido-Atienza and Real (2003); Rodríguez et al. (2012). Du et al. (2010) developed a high order difference method for the fractional diffusion-wave equation associated with $\gamma \in (1,2)$ in (1.1) with constant coefficients and without the convective term. Li et al. (2019) investigates a linearized FDM for the problem similar to (1.1) without the reaction term and non-linear source term. First, the original problem is transformed into an equivalent semilinear fractional delay reaction-diffusion equation. Then they used the central finite difference formula for the space derivative and $L_1$ approximation for the Caputo derivative. Finally, the inverse exponential recovery method is used to obtain the numerical solution.

This remainder article is outlined as follows. In Sect. 2, we semi-discretize the problem (1.1) in the temporal direction followed by the full-discretization in the spatial direction. The local truncation error is computed in Sect. 3, whereas the proof of the stability of the proposed scheme is given in Sect. 4. The rigorous convergence analysis is provided in Sect. 5. In Sect. 6, the numerical results of two test problems are reported to confirm the accuracy and efficiency of the proposed scheme. Finally, the paper ends with some observations, concluding remarks, and the proposed method’s future scope in Sect. 7.

Remark 1

Throughout the manuscript, we take C as a positive generic constant that takes different values at different places.

2 Discretization of the problem

2.1 The time semi-discretization

To discretize the time domain, we use an equidistant mesh. Let $\varLambda ^{k}_\tau =\{-\tau =t_{-k}<t_{-k+1}<\ldots<t_{-1}<t_0=0\}$ be the partition of $[-\tau ,0]$ divided into k sub-intervals using step size $\Delta t=\tau /k$. Then, using the mesh size $\Delta t$, we discretize $[0,{\mathfrak {T}}]$ into $M_t={\mathfrak {T}}/\Delta t$ sub-intervals and denote by $\varLambda ^{M_{t}}$, the collection of all nodal points of the domain $[0,{\mathfrak {T}}]$. Then

$$\begin{aligned} \varLambda ^{M_{t}}=\{0=t_0<t_1<\ldots<t_k=\tau<\ldots<t_{M_{t}-1}<t_{M_{t}}={\mathfrak {T}}\}.\end{aligned}$$

Now we semi-discretize the problem (1.1) on $\varLambda ^{M_{t}}$ as

$$\begin{aligned}&-\frac{\partial ^2 {\widetilde{y}}(x,t_{n+\frac{1}{2}})}{\partial x^2}+b(x,t_{n+\frac{1}{2}})\frac{\partial {\widetilde{y}}(x,t_{n+\frac{1}{2}})}{\partial x}+c(x,t_{n+\frac{1}{2}}) {\widetilde{y}}(x,t_{n+\frac{1}{2}})\nonumber \\&\quad =-d(x,t_{n+\frac{1}{2}}) {\widetilde{y}}(x,t_{n+\frac{1}{2}-k})+f(x,t_{n+\frac{1}{2}})-\frac{\partial ^\gamma \widetilde{y}(x,t_{n+\frac{1}{2}})}{\partial t^\gamma }, \end{aligned}$$

(2.1a)

for $n=1,2,\ldots ,M_{t}-1$ with

$$\begin{aligned} {\widetilde{y}}(x,t_{n+1})= & {} \psi _b(x,t_{n+1}),\ \text {for}\, (x,t_{n+1})\in [0,1]\times \varLambda ^{k}_\tau , \end{aligned}$$

(2.1b)

$$\begin{aligned} {\widetilde{y}}(0,t_{n+1})= & {} \psi _l(t_{n+1}),\quad \widetilde{y}(1,t_{n+1})=\psi _r(t_{n+1}), \ t_{n+1}\in \varLambda ^{M_{t}}, \end{aligned}$$

(2.1c)

where ${\widetilde{y}}$ represents the approximate solution to the problem (1.1). We approximate $\frac{\partial ^{\gamma }{\widetilde{y}}(x,t_{n+1/2})}{\partial t^{\gamma }}$ using the following approach:

$$\begin{aligned}&\frac{\partial ^{\gamma }{\widetilde{y}}(x,t_{n+1/2})}{\partial t^{\gamma }}\\&\quad =\ \frac{1}{\Gamma (1-\gamma )}\int _{0}^{t_{n+1/2}}{\widetilde{y}}_t(x,t_{n+1/2})(t_{n+1/2}-s)^{-\gamma }ds\\&\quad =\ \frac{1}{\Gamma (1-\gamma )}\int _{0}^{t_{n}}{\widetilde{y}}_t(x,t_{n+1/2})\left( \left( n+\frac{1}{2}\right) \Delta t-s\right) ^{-\gamma }ds\\&\qquad +\frac{1}{\Gamma (1-\gamma )}\int _{t_{n}}^{t_{n+1/2}}\left[ \frac{{\widetilde{y}}^{n+1}-{\widetilde{y}}^{n}}{\Delta t}+{\mathcal {O}}((\Delta t)^{2})\right] \left( \left( n+\frac{1}{2}\right) \Delta t-s\right) ^{-\gamma }ds\\&\quad = \frac{1}{\Gamma (1-\gamma )}\sum _{i=1}^{n}\int _{(i-1)\Delta t}^{i \Delta t}\left[ \frac{{\widetilde{y}}^{i}-{\widetilde{y}}^{i-1}}{\Delta t}+{\mathcal {O}}((\Delta t)^2)\right] \left( \left( n+\frac{1}{2}\right) \Delta t-s\right) ^{-\gamma }ds\\&\qquad +\frac{1}{\Gamma (1-\gamma )}\int _{n \Delta t}^{(n+1/2)\Delta t}\left[ \frac{{\widetilde{y}}^{n+1}-{\widetilde{y}}^{n}}{\Delta t}+{\mathcal {O}}((\Delta t)^{2})\right] \left( \left( n+\frac{1}{2}\right) \Delta t-s\right) ^{-\gamma }ds.\\&\quad =\ \frac{1}{\Gamma (1-\gamma )}\sum _{i=1}^{n}\left( \frac{{\widetilde{y}}^{i}-{\widetilde{y}}^{i-1}}{\Delta t}\right) \int _{(i-1) \Delta t}^{i \Delta t}\left( \left( n+\frac{1}{2}\right) \Delta t-s\right) ^{-\gamma }ds\\&\qquad +\frac{{\mathcal {O}}((\Delta t)^2)}{\Gamma (1-\gamma )}\sum _{i=1}^{n}\int _{(i-1)\Delta t}^{i \Delta t}\left( \left( n+\frac{1}{2}\right) \Delta t-s\right) ^{-\gamma }ds\\&\qquad +\frac{1}{\Gamma (1-\gamma )}\int _{n \Delta t}^{(n+1/2) \Delta t}\left[ \frac{{{\widetilde{y}}}^{n+1}-{{\widetilde{y}}}^{n}}{\Delta t}+{\mathcal {O}}((\Delta t)^{2})\right] \left( \left( n+\frac{1}{2}\right) \Delta t-s\right) ^{-\gamma }ds\\&\quad =\ \frac{1}{\Gamma (1-\gamma )}\frac{1}{(\Delta t)^{\gamma }(1-\gamma )}\sum _{i=1}^{n}\left[ {\widetilde{y}}^{i}-{\widetilde{y}}^{i-1}\right] \left[ \left( n-i+\frac{3}{2}\right) ^{1-\gamma }-\left( n-i+\frac{1}{2}\right) ^{1-\gamma }\right] \\&\qquad +\frac{1}{\Gamma (1-\gamma )}\frac{1}{(1-\gamma )}\sum _{i=1}^{n}\left[ \left( n-i+\frac{3}{2}\right) ^{1-\gamma }-\left( n-i+\frac{1}{2}\right) ^{1-\gamma }\right] {\mathcal {O}}((\Delta t)^{3-\gamma })\\&\qquad +\frac{1}{\Gamma (1-\gamma )}\frac{1}{(\Delta t)^{\gamma }(1-\gamma )}({\widetilde{y}}^{n+1}-\widetilde{y}^{n})\frac{1}{2^{1-\gamma }}+\frac{1}{\Gamma (1-\gamma )}\frac{1}{(1-\gamma )}\frac{1}{2^{1-\gamma }}{\mathcal {O}}((\Delta t)^{3-\gamma }). \end{aligned}$$

Thus, we get

$$\begin{aligned}&\frac{\partial ^{\gamma }{\widetilde{y}}(x,t_{n+1/2})}{\partial t^{\gamma }}\nonumber \\&\quad =\ \frac{1}{\Gamma (2-\gamma )}\frac{1}{(\Delta t)^{\gamma }}({\widetilde{y}}^{n+1}-{\widetilde{y}}^{n})\frac{1}{2^{1-\gamma }}\nonumber \\&\qquad +\nonumber \frac{1}{\Gamma (2-\gamma )}\frac{1}{(\Delta t)^{\gamma }}\sum _{i=1}^{n}\left[ {\widetilde{y}}^{i}-{\widetilde{y}}^{i-1}\right] \left[ \left( n-i+\frac{3}{2}\right) ^{1-\gamma }-\left( n-i+\frac{1}{2}\right) ^{1-\gamma }\right] \nonumber \\&\qquad +\nonumber \frac{1}{\Gamma (2-\gamma )}\sum _{i=1}^{n}\left[ \left( n-i+\frac{3}{2}\right) ^{1-\gamma }-\left( n-i+\frac{1}{2}\right) ^{1-\gamma }\right] {\mathcal {O}}((\Delta t)^{3-\gamma })\nonumber \\&\qquad +\frac{1}{\Gamma (2-\gamma )}\frac{1}{2^{1-\gamma }}{\mathcal {O}}((\Delta t)^{3-\gamma }). \end{aligned}$$

(2.2)

Let $\sigma = \frac{1}{\Gamma (2-\gamma )}\frac{1}{(\Delta t)^{\gamma }}$, and $w_{j} = \left( j+\frac{1}{2}\right) ^{1-\gamma }-\left( j-\frac{1}{2}\right) ^{1-\gamma }$, so that $\sum _{j=1}^{n_1}w_{j}= \left( n_1+\frac{1}{2}\right) ^{1-\gamma }-\left( \frac{1}{2}\right) ^{1-\gamma }$. Substituting $\sigma$ and $w_{j}$ into (2.2), we get

$$\begin{aligned} \frac{\partial ^{\gamma }{\widetilde{y}}(x,t_{n+1/2})}{\partial t^{\gamma }}&=\sigma \left[ w_{1}\widetilde{y}^{n}+\sum _{i=1}^{n-1}(w_{n-i+1}-w_{n-i})\widetilde{y}^{i}-w_{n}{\widetilde{y}}^{0}+\frac{\widetilde{y}^{n+1}-{\widetilde{y}}^{n}}{2^{1-\gamma }}\right] +R_{1}+R_{2}, \end{aligned}$$

where

$$\begin{aligned} R_{1}+R_{2}=\frac{1}{\Gamma (2-\gamma )}\sum _{i=1}^{n}w_{n-i+1}{\mathcal {O}}((\Delta t)^{3-\gamma })+\frac{1}{\Gamma (2-\gamma )}\frac{1}{2^{1-\gamma }}{\mathcal {O}}((\Delta t)^{3-\gamma }). \end{aligned}$$

Let $j = n-i+1$, then

$$\begin{aligned} R_{1}+R_{2}&=\frac{1}{\Gamma (2-\gamma )}\left[ \sum _{j=1}^{n}w_{j}{\mathcal {O}}((\Delta t)^{3-\gamma })+\frac{1}{2^{1-\gamma }}{\mathcal {O}}((\Delta t)^{3-\gamma })\right] \\&=\frac{1}{\Gamma (2-\gamma )}\left[ \biggl (n+\frac{1}{2}\biggr )^{1-\gamma }-\biggl (\frac{1}{2}\biggr )^{1-\gamma }+\frac{1}{2^{1-\gamma }}\right] {\mathcal {O}}((\Delta t)^{3-\gamma })\\&=\frac{1}{\Gamma (2-\gamma )}\biggl (n+\frac{1}{2}\biggr )^{1-\gamma }{\mathcal {O}}((\Delta t)^{3-\gamma })\\&=\frac{1}{\Gamma (2-\gamma )}\biggl (\frac{t_{n}}{\Delta t}+\frac{1}{2}\biggr )^{1-\gamma }{\mathcal {O}}((\Delta t)^{3-\gamma })\\&\le C(\Delta t)^2. \end{aligned}$$

Thus, we have the following approximation at $(x,t_{n+\frac{1}{2}})$

$$\begin{aligned} \frac{\partial ^{\gamma }{\widetilde{y}}(x,t_{n+1/2})}{\partial t^{\gamma }}= & {} \sigma \left[ w_{1}\widetilde{y}^{n}+\sum _{i=1}^{n-1}(w_{n-i+1}-w_{n-i})\widetilde{y}^{i}-w_{n}{\widetilde{y}}^{0}\right. \nonumber \\&\left. +\frac{\widetilde{y}^{n+1}-\widetilde{y}^{n}}{2^{1-\gamma }}\right] +{\mathcal {O}}((\Delta t)^{2}). \end{aligned}$$

(2.3)

Thus, we obtained a second-order approximation to $\frac{\partial ^{\gamma }{\widetilde{y}}(x,t_{n+1/2})}{\partial t^{\gamma }}$.

2.2 The spatial discretization

Let $\varOmega ^{N_{x}}=\{0=x_0<x_1<\ldots<x_{N_{x}-1}<x_{N_{x}}=1\}$ be the partition of [0, 1] divided into $N_x$ sub-intervals each of width $\Delta x=\frac{1}{N_x}$. Let the function ${\mathcal {S}}(x,t,\delta )$ belongs to the class $C^2[0,1]$, which gives interpolation of y(x, t) at the mesh point $(x_m,t_n)$, where $\delta$ is termed as tension factor. As $\delta \rightarrow 0$ the function ${\mathcal {S}}(x,t,\delta )$ is turned to be parametric cubic spline in [0, 1]. Let ${\mathcal {S}}^n_m$ be an approximation of ${\widetilde{y}}_m^n$ obtained by the segment of the functions passing through the points $(x_m,{\mathcal {S}}_m^n)$ and $(x_{m+1},{\mathcal {S}}_{m+1}^n)$. For $\delta >0$, ${\mathcal {S}}_{\pi }(x,t_n,\delta )$ satisfies the following differential equation in the sub-interval $[x_m,x_{m+1}]$

$$\begin{aligned} {\mathcal {S}}''_\pi (x,t_n,\delta )-\delta {\mathcal {S}}_\pi (x,t_n,\delta )=&\ [{\mathcal {S}}''_\pi (x_m,t_n,\delta )-\delta {\mathcal {S}}_\pi (x_m,t_n,\delta )](x_{m+1}-x)/\Delta x\nonumber \\&+[{\mathcal {S}}''_\pi (x_{m+1},t_n,\delta )-\delta {\mathcal {S}}_\pi (x_{m+1},t_n,\delta )](x-x_m)/\Delta x, \end{aligned}$$

(2.4)

with the following interpolation conditions

$$\begin{aligned} {\mathcal {S}}_\pi (x_m,t_n,\delta )={\widetilde{y}}_m^n,\, {\mathcal {S}}_\pi (x_{m+1},t_n,\delta )={\widetilde{y}}_{m+1}^n. \end{aligned}$$

(2.5)

The function derivative $y''(x_m,t_n)$ has the spline derivative approximation given by

$$\begin{aligned} {\mathcal {S}}''_\pi (x_m,t_n,\delta )={\mathcal {M}}(x_m,t_n)\, ,\, {\mathcal {S}}''_\pi (x_{m+1},t_n,\delta )={\mathcal {M}}(x_{m+1},t_n). \end{aligned}$$

(2.6)

The solution of (2.4) with the interpolatory conditions (2.5) can be written as

$$\begin{aligned}&{\mathcal {S}}_\pi (x,t_n,\delta )\nonumber \\&\quad =\ \frac{(\Delta x)^2}{\nu ^2\sinh (\nu )}\left[ {\mathcal {M}}_{m+1}^{n}\sinh \frac{\nu (x-x_m)}{\Delta x}+{\mathcal {M}}_m^n\sinh \frac{\nu (x_{m+1}-x)}{\Delta x}\right] \nonumber \\&\qquad -\frac{(\Delta x)^2}{\nu ^2}\left[ \frac{(x-x_m)}{\Delta x}\left( {\mathcal {M}}_{m+1}^n-\frac{\nu ^2}{(\Delta x)^2}\widetilde{y}_{m+1}^n\right) +\frac{(x_{m+1}-x)}{\Delta x}\left( {\mathcal {M}}_{m}^n-\frac{\nu ^2}{(\Delta x)^2}\widetilde{y}_{m}^n\right) \right] . \end{aligned}$$

(2.7)

Rewrite the Eq. (2.7) as follows

$$\begin{aligned} {\mathcal {S}}_\pi (x,t_n,\delta )=\mu {{\widetilde{y}}}_{m+1}^n+{\overline{\mu }} {{\widetilde{y}}}_m^n +(\Delta x)^2[g(\mu ){\mathcal {M}}_{m+1}^n+g(\overline{\mu }){\mathcal {M}}_{m}^n]/\nu ^2, \end{aligned}$$

(2.8)

where $\mu =(x-x_m)/\Delta x$, ${\overline{\mu }}=1-\mu$, $\nu =\Delta x\sqrt{\delta }$, $g(\mu )=\frac{\sinh (\nu \mu )}{\sinh (\nu )}-\mu$. Differentiating (2.8) and let $x\rightarrow x_m$ to get

$$\begin{aligned} {\mathcal {S}}'_\pi (x_m^+,t_n,\delta )=\frac{1}{\Delta x}(\widetilde{y}_{m+1}^n-{{\widetilde{y}}}_m^n)+\frac{\Delta x}{\nu ^2}\left[ \left( \frac{\nu }{\sinh (\nu )}-1\right) {\mathcal {M}}_{m+1}^n+\left( 1-\frac{\nu \cosh (\nu )}{\sinh (\nu )}\right) {\mathcal {M}}_m^n\right] . \nonumber \\ \end{aligned}$$

(2.9)

Similarly, we proceed for the interval $[x_{m-1},x_m]$, to obtain

$$\begin{aligned} {\mathcal {S}}'_\pi (x_m^-,t_n,\delta )=\frac{1}{\Delta x}(\widetilde{y}_{m}^n-{{\widetilde{y}}}_{m-1}^n)+\frac{\Delta x}{\nu ^2}\left[ \left( 1-\frac{\nu }{\sinh (\nu )}\right) {\mathcal {M}}_{m-1}^n+\left( \frac{\nu \cosh (\nu )}{\sinh (\nu )}-1\right) {\mathcal {M}}_m^n\right] . \nonumber \\ \end{aligned}$$

(2.10)

The first derivative of ${\mathcal {S}}_\pi (x,t_n,\delta )$ is continuous at $x=x_m$, so ${\mathcal {S}}'_\pi (x_m^-,t_n,\delta )={\mathcal {S}}'_\pi (x_m^+,t_n,\delta )$ for $m=1,2,\ldots , N_x-1,\, n=1,2,\ldots ,M_t$, which gives

$$\begin{aligned} \alpha {\mathcal {M}}_{m-1}^n+\beta {\mathcal {M}}_m^n+\alpha {\mathcal {M}}_{m+1}^n=\frac{1}{(\Delta x)^2}\left[ {{\widetilde{y}}}_{m-1}^n-2{{\widetilde{y}}}_m^n+\widetilde{y}_{m+1}^n\right] , \end{aligned}$$

(2.11)

where $\alpha =\frac{1}{\nu ^2}\left( 1-\frac{\nu }{\sinh (\nu )}\right) ,\ \beta =\frac{2}{\nu ^2}\left( \frac{\nu \cosh (\nu )}{\sinh (\nu )}-1\right)$. Using Eq. (2.3), the Eq. (2.1a) at $(x_m,t_{n+1/2})$ can be written as

$$\begin{aligned}&-\left( \frac{({{\widetilde{y}}}_{xx})_m^n+({{\widetilde{y}}}_{xx})_m^{n+1}}{2}\right) +b_m^{n+1/2}\left( \frac{({{\widetilde{y}}}_{x})_m^n+({{\widetilde{y}}}_{x})_m^{n+1}}{2}\right) +c_m^{n+1/2}\left( \frac{{{\widetilde{y}}}_m^n+{{\widetilde{y}}}_m^{n+1}}{2}\right) \nonumber \\&\quad =-d_m^{n+1/2}\left( \frac{{{\widetilde{y}}}_m^{n-k}+\widetilde{y}_m^{n+1-k}}{2}\right) +f_m^{n+1/2}-\sigma \left( \frac{{\widetilde{y}}_m^{n+1}-{\widetilde{y}}_m^{n}}{2^{1-\gamma }}\right) \nonumber \\&\qquad -\sigma \left[ w_{1}\widetilde{y}_m^{n}+\sum _{i=1}^{n-1}(w_{n-i+1}-w_{n-i})\widetilde{y}_m^{i}-w_{n}{\widetilde{y}}_m^{0}\right] . \end{aligned}$$

(2.12)

Now on replacing the second-order derivative in (2.12) by the tension spline function defined as

$$\begin{aligned} (y_{xx})_m^n={\mathcal {S}}''(x_m,t_n)={\mathcal {M}}_m^n+{\mathcal {O}}((\Delta x)^2), \end{aligned}$$

(2.13)

we get

$$\begin{aligned} \frac{{\mathcal {M}}_m^n+{\mathcal {M}}_m^{n+1}}{2}=&\ \sigma \left[ w_{1}{\widetilde{y}}_m^{n}+\sum _{i=1}^{n-1}(w_{n-i+1}-w_{n-i}){\widetilde{y}}_m^{i}-w_{n}{\widetilde{y}}_m^{0}+\frac{{\widetilde{y}}_m^{n+1}-{\widetilde{y}}_m^{n}}{2^{1-\gamma }}\right] \nonumber \\&+b_m^{n+1/2}\left( \frac{({{\widetilde{y}}}_{x})_m^n+({{\widetilde{y}}}_{x})_m^{n+1}}{2}\right) +c_m^{n+1/2}\left( \frac{{{\widetilde{y}}}_m^n+{{\widetilde{y}}}_m^{n+1}}{2}\right) \nonumber \\&+d_m^{n+1/2}\left( \frac{{{\widetilde{y}}}_m^{n-k}+\widetilde{y}_m^{n+1-k}}{2}\right) -f_m^{n+1/2}+{\mathfrak {R}}_m^{n+1/2}. \end{aligned}$$

(2.14)

Rewrite (2.11) at $(n+1)$-th time-level as

$$\begin{aligned} \alpha {\mathcal {M}}_{m-1}^{n+1}+\beta {\mathcal {M}}_m^{n+1}+\alpha {\mathcal {M}}_{m+1}^{n+1}= & {} \frac{1}{(\Delta x)^2}\left[ {{\widetilde{y}}}_{m-1}^{n+1}-2{{\widetilde{y}}}_m^{n+1}+\widetilde{y}_{m+1}^{n+1}\right] ,\nonumber \\&m=1,2,\ldots , N_x-1,\, n=0,1,\ldots ,M_t-1. \end{aligned}$$

(2.15)

On adding (2.11) and (2.15), we get

$$\begin{aligned}&\alpha \left( \frac{{\mathcal {M}}_{m-1}^n+{\mathcal {M}}_{m-1}^{n+1}}{2}\right) +\beta \left( \frac{{\mathcal {M}}_m^n+{\mathcal {M}}_m^{n+1}}{2}\right) +\alpha \left( \frac{{\mathcal {M}}_{m+1}^n+{\mathcal {M}}_{m+1}^{n+1}}{2}\right) \nonumber \\&\quad =\frac{1}{(\Delta x)^2}\left[ \frac{{{\widetilde{y}}}_{m-1}^n+\widetilde{y}_{m-1}^{n+1}}{2}-{{\widetilde{y}}}_m^n-\widetilde{y}_m^{n+1}+\frac{{{\widetilde{y}}}_{m+1}^n+\widetilde{y}_{m+1}^{n+1}}{2}\right] . \end{aligned}$$

(2.16)

A use of (2.14) in (2.16) yields

$$\begin{aligned}&\left[ \alpha \left( \frac{\sigma }{2^{1-\gamma }}-\frac{b_{m-1}^{n+1/2}}{2\Delta x}+\frac{c_{m-1}^{n+1/2}}{2}\right) -\frac{\beta b_m^{n+1/2}}{4\Delta x}-\frac{1}{2(\Delta x)^2}\right] {{\widetilde{y}}}_{m-1}^{n+1}\nonumber \\&\qquad +\bigg [\frac{\alpha b_{m-1}^{n+1/2}}{2(\Delta x)}+\beta \left( \frac{\sigma }{2^{1-\gamma }}+\frac{c_{m}^{n+1/2}}{2}\right) -\frac{\alpha b_{m+1}^{n+1/2}}{2\Delta x}+\frac{1}{(\Delta x)^2}\bigg ]{{\widetilde{y}}}_m^{n+1}\nonumber \\&\qquad +\left[ \alpha \left( \frac{\sigma }{2^{1-\gamma }}+\frac{b_{m+1}^{n+1/2}}{2\Delta x}+\frac{c_{m+1}^{n+1/2}}{2}\right) +\frac{\beta b_m^{n+1/2}}{4\Delta x}-\frac{1}{2(\Delta x)^2}\right] {{\widetilde{y}}}_{m+1}^{n+1}\nonumber \\&\quad =\ \left[ \alpha \left( \frac{\sigma }{2^{1-\gamma }}+\frac{b_{m-1}^{n+1/2}}{2\Delta x}-\frac{c_{m-1}^{n+1/2}}{2}\right) +\frac{\beta b_m^{n+1/2}}{4\Delta x}+\frac{1}{2(\Delta x)^2}\right] {{\widetilde{y}}}_{m-1}^{n}\nonumber \\&\qquad +\bigg [-\frac{\alpha b_{m-1}^{n+1/2}}{2\Delta x}+\beta \left( \frac{\sigma }{2^{1-\gamma }}-\frac{c_{m}^{n+1/2}}{2}\right) +\frac{\alpha b_{m+1}^{n+1/2}}{2\Delta x}-\frac{1}{(\Delta x)^2}\bigg ]{{\widetilde{y}}}_m^{n}\nonumber \\&\qquad +\left[ \alpha \left( \frac{\sigma }{2^{1-\gamma }}-\frac{b_{m+1}^{n+1/2}}{2\Delta x}-\frac{c_{m+1}^{n+1/2}}{2}\right) -\frac{\beta b_m^{n+1/2}}{4\Delta x}+\frac{1}{2(\Delta x)^2}\right] {{\widetilde{y}}}_{m+1}^{n}\nonumber \\&\quad -\sigma w_1\left( \alpha {{\widetilde{y}}}_{m-1}^n+\beta {{\widetilde{y}}}_{m}^n+\alpha {{\widetilde{y}}}_{m+1}^n\right) +\sigma w_n\left( \alpha {{\widetilde{y}}}_{m-1}^0+\beta {{\widetilde{y}}}_{m}^0+\alpha {{\widetilde{y}}}_{m+1}^0\right) \nonumber \\&\quad -\sigma \sum _{q=1}^{n-1}(w_{n-q+1}-w_{n-q})\left( \alpha {{\widetilde{y}}}_{m-1}^q+\beta {{\widetilde{y}}}_{m}^q+\alpha {{\widetilde{y}}}_{m+1}^q\right) -\frac{1}{2}\bigg [\alpha d_{m-1}^{n+1/2}({{\widetilde{y}}}_{m-1}^{n+1-k}+{{\widetilde{y}}}_{m-1}^{n-k})\nonumber \\&\quad +\beta d_{m}^{n+1/2}({{\widetilde{y}}}_{m}^{n+1-k}+{{\widetilde{y}}}_{m}^{n-k})+\alpha d_{m+1}^{n+1/2}({{\widetilde{y}}}_{m+1}^{n+1-k}+{{\widetilde{y}}}_{m+1}^{n-k})\bigg ]\nonumber \\&\quad +\left( \alpha f_{m-1}^{n+1/2}+\beta f_{m}^{n+1/2}+\alpha f_{m+1}^{n+1/2}\right) +{\mathfrak {T}}_m^{n+1/2}, \end{aligned}$$

(2.17)

where ${\mathfrak {T}}_m^{n+1/2}=\alpha {\mathfrak {R}}_{m-1}^{n+1/2}+\beta {\mathfrak {R}}_{m}^{n+1/2}+\alpha {\mathfrak {R}}_{m+1}^{n+1/2}$. We assume that $Y_m^n$ is the approximate solution of ${{\widetilde{y}}}_m^n$ and omitting the truncation error from (2.17), we find the following numerical scheme for (1.1)

$$\begin{aligned}&\left[ \alpha \left( \frac{\sigma }{2^{1-\gamma }}-\frac{b_{m-1}^{n+1/2}}{2\Delta x}+\frac{c_{m-1}^{n+1/2}}{2}\right) -\frac{\beta b_m^{n+1/2}}{4\Delta x}-\frac{1}{2(\Delta x)^2}\right] Y_{m-1}^{n+1}\nonumber \\&\qquad +\bigg [\frac{\alpha b_{m-1}^{n+1/2}}{2\Delta x}+\alpha \left( \frac{\sigma }{2^{1-\gamma }}+\frac{c_{m}^{n+1/2}}{2}\right) -\frac{\alpha b_{m+1}^{n+1/2}}{2\Delta x}+\frac{1}{(\Delta x)^2}\bigg ]Y_m^{n+1}\nonumber \\&\qquad +\left[ \alpha \left( \frac{\sigma }{2^{1-\gamma }}+\frac{b_{m+1}^{n+1/2}}{2\Delta x}+\frac{c_{m+1}^{n+1/2}}{2}\right) +\frac{\beta b_m^{n+1/2}}{4\Delta x}-\frac{1}{2(\Delta x)^2}\right] Y_{m+1}^{n+1}\nonumber \\&\quad =\ \left[ \alpha \left( \frac{\sigma }{2^{1-\gamma }}+\frac{b_{m-1}^{n+1/2}}{2\Delta x}-\frac{c_{m-1}^{n+1/2}}{2}\right) +\frac{\beta b_m^{n+1/2}}{4\Delta x}+\frac{1}{2(\Delta x)^2}\right] Y_{m-1}^{n} \nonumber \\&\qquad +\bigg [-\frac{\alpha b_{m-1}^{n+1/2}}{2\Delta x}+\beta \left( \frac{\sigma }{2^{1-\gamma }}-\frac{c_{m}^{n+1/2}}{2}\right) +\frac{\alpha b_{m+1}^{n+1/2}}{2\Delta x}-\frac{1}{(\Delta x)^2}\bigg ]Y_m^{n}\nonumber \\&\quad +\left[ \alpha \left( \frac{\sigma }{2^{1-\gamma }}-\frac{b_{m+1}^{n+1/2}}{2\Delta x}-\frac{c_{m+1}^{n+1/2}}{2}\right) -\frac{\beta b_m^{n+1/2}}{4\Delta x}+\frac{1}{2(\Delta x)^2}\right] Y_{m+1}^{n}\nonumber \\&\quad -\sigma w_1\left( \alpha Y_{m-1}^n+\beta Y_{m}^n+\alpha Y_{m+1}^n\right) +\sigma w_n\left( \alpha Y_{m-1}^0+\beta Y_{m}^0+\alpha Y_{m+1}^0\right) \nonumber \\&\quad -\sigma \sum _{q=1}^{n-1}(w_{n-q+1}-w_{n-q})\left( \alpha Y_{m-1}^q+\beta Y_{m}^q+\alpha Y_{m+1}^q\right) \nonumber \\&\quad -\frac{1}{2}\Bigg [\alpha d_{m-1}^{n+1/2}(Y_{m-1}^{n+1-k}+Y_{m-1}^{n-k})+\beta d_{m}^{n+1/2}(Y_{m}^{n+1-k}+Y_{m}^{n-k})+\alpha d_{m+1}^{n+1/2}(Y_{m+1}^{n+1-k}\nonumber \\&\quad +Y_{m+1}^{n-k})\Bigg ]+\left( \alpha f_{m-1}^{n+1/2}+\beta f_{m}^{n+1/2}+\alpha f_{m+1}^{n+1/2}\right) ,\nonumber \\&\qquad m=1,2,\ldots , N_x-1, \quad n=0,1,\ldots ,M_t-1, \end{aligned}$$

(2.18a)

with the interval and boundary conditions

$$\begin{aligned} Y_m^n= & {} \psi _b(x_m,t_n),\quad 0\le m\le N_x, \quad -k\le n\le 0, \end{aligned}$$

(2.18b)

$$\begin{aligned} Y_0^{n+1}= & {} \psi _l(t_{n+1}), \quad Y_{N_x}^{n+1}=\psi _r(t_{n+1}), \quad 0\le n\le M_{t}-1. \end{aligned}$$

(2.18c)

Let

$$\begin{aligned} E(m,n)&=\alpha \left( \frac{\sigma }{2^{1-\gamma }}-\frac{b_{m-1}^{n+1/2}}{2\Delta x}+\frac{c_{m-1}^{n+1/2}}{2}\right) -\frac{\beta b_m^{n+1/2}}{4\Delta x}-\frac{1}{2(\Delta x)^2},\nonumber \\ F(m,n)&=\frac{\alpha b_{m-1}^{n+1/2}}{2\Delta x}+\beta \left( \frac{\sigma }{2^{1-\gamma }}+\frac{c_{m}^{n+1/2}}{2}\right) -\frac{\alpha b_{m+1}^{n+1/2}}{2\Delta x}+\frac{1}{(\Delta x)^2},\nonumber \\ G(m,n)&=\alpha \left( \frac{\sigma }{2^{1-\gamma }}+\frac{b_{m+1}^{n+1/2}}{2\Delta x}+\frac{c_{m+1}^{n+1/2}}{2}\right) +\frac{\beta b_m^{n+1/2}}{4\Delta x}-\frac{1}{2(\Delta x)^2}, \end{aligned}$$

(2.19)

then we can write the system (2.18) in matrix form as

$$\begin{aligned}AC^{n+1}=B^n,\end{aligned}$$

with the coefficient matrix

$$\begin{aligned}A=\begin{bmatrix} F(1,n) &{} G(1,n) &{} 0 \\ E(2,n)&{} F(2,n) &{} G(2,n) \\ \ddots &{} \ddots &{} \ddots \\ &{}&{}\ddots &{}\ddots &{}\ddots \\ &{}&{} E(N_x-2,n) &{} F(N_x-2,n) &{} G(N_x-2,n)\\ &{}&{}0 &{} E(N_x-1,n)&{} F(N_x-1,n)\\ \end{bmatrix}_{(N_x-1)\times (N_x-1)}.\end{aligned}$$

3 Local truncation error

To find the local truncation error ${\mathfrak {e}}_m^n$ of the proposed scheme, we replace ${\mathcal {M}}_{m-1}^n$, ${\mathcal {M}}_{m}^n$, and ${\mathcal {M}}_{m+1}^n$ by $\frac{\partial ^2{{\widetilde{y}}}_{m-1}^n}{\partial x^2}$, $\frac{\partial ^2{{\widetilde{y}}}_{m}^n}{\partial x^2}$, and $\frac{\partial ^2{{\widetilde{y}}}_{m+1}^n}{\partial x^2}$, respectively in (2.11), to get

$$\begin{aligned} {{\mathfrak {e}}_m^n}=\alpha \frac{\partial ^2{{\widetilde{y}}}_{m-1}^n}{\partial x^2}+\beta \frac{\partial ^2{{\widetilde{y}}}_{m}^n}{\partial x^2}+\alpha \frac{\partial ^2{{\widetilde{y}}}_{m+1}^n}{\partial x^2}-\frac{1}{(\Delta x)^2}\left[ {{\widetilde{y}}}_{m-1}^n-2{{\widetilde{y}}}_m^n+{{\widetilde{y}}}_{m+1}^n\right] +{\mathcal {O}}((\Delta x)^p),\ p\ge 0. \end{aligned}$$

Using Taylor series expansion for ${{\widetilde{y}}}_{m-1}^n$ and ${{\widetilde{y}}}_{m+1}^n$ in term of ${{\widetilde{y}}}_{m}^n$ and its derivatives, we get

$$\begin{aligned} {{\mathfrak {e}}_m^n}= & {} ((2\alpha +\beta )-1)\frac{\partial ^2{{\widetilde{y}}}_{m}^n}{\partial x^2}+ (\Delta x)^2\left( \alpha -\frac{1}{12}\right) \frac{\partial ^4\widetilde{y}_{m}^n}{\partial x^4}+(\Delta x)^4\left( \frac{1}{360}-\frac{\alpha }{12}\right) \frac{\partial ^6\widetilde{y}_{m}^n}{\partial x^6}\nonumber \\&\quad +{\mathcal {O}}((\Delta x)^p). \end{aligned}$$

(3.1)

It is clear from (3.1) that for suitable arbitrary values of $\alpha$ ($\alpha \ne 1/12$) and $\beta$ such that $2\alpha +\beta =1$, the value of p is 2. Now, from (2.12) the local truncation error ${\mathfrak {R}}_m^{n+1/2}$ is

$$\begin{aligned} {\mathfrak {R}}_m^{n+1/2}=&-\sigma \left[ w_{1}{\widetilde{y}}_m^{n}+\sum _{i=1}^{n-1}(w_{n-i+1}-w_{n-i}){\widetilde{y}}_m^{i}-w_{n}{\widetilde{y}}_m^{0}+\frac{{\widetilde{y}}_m^{n+1}-{\widetilde{y}}_m^{n}}{2^{1-\gamma }}\right] \\&+b_m^{n+1/2}\left( \frac{({{\widetilde{y}}}_{x})_m^n+({{\widetilde{y}}}_{x})_m^{n+1}}{2}\right) +c_m^{n+1/2}\left( \frac{{{\widetilde{y}}}_m^n+{{\widetilde{y}}}_m^{n+1}}{2}\right) \\&-d_m^{n+1/2}\left( \frac{{{\widetilde{y}}}_m^{n-k}+\widetilde{y}_m^{n+1-k}}{2}\right) + {\mathcal {M}}_m^{n+1/2}+f_m^{n+1/2}\\ =&\ {\mathcal {O}}((\Delta x)^2+(\Delta t)^2). \end{aligned}$$

Thus, the local truncation error ${\mathfrak {e}}_m^{n+1/2}$ in (2.17) is given by

$$\begin{aligned} {{\mathfrak {e}}_m^{n+1/2}}&=\ \alpha {\mathfrak {R}}_{m-1}^{n+1/2}+\beta {\mathfrak {R}}_{m}^{n+1/2}+\alpha {\mathfrak {R}}_{m+1}^{n+1/2}\nonumber \\&= -\alpha \sigma \left[ w_{1}{\widetilde{y}}_{m-1}^{n}+\sum _{i=1}^{n-1}(w_{n-i+1}-w_{n-i}){\widetilde{y}}_{m-1}^{i}-w_{n}{\widetilde{y}}_{m-1}^{0}+\frac{{\widetilde{y}}_{m-1}^{n+1}-{\widetilde{y}}_{m-1}^{n}}{2^{1-\gamma }}\right] \nonumber \\&\quad -\beta \sigma \left[ w_{1}{\widetilde{y}}_{m}^{n}+\sum _{i=1}^{n-1}(w_{n-i+1}-w_{n-i}){\widetilde{y}}_{m}^{i}-w_{n}{\widetilde{y}}_{m}^{0}+\frac{{\widetilde{y}}_{m}^{n+1}-{\widetilde{y}}_{m}^{n}}{2^{1-\gamma }}\right] \nonumber \\&\quad -\alpha \sigma \left[ w_{1}{\widetilde{y}}_{m+1}^{n}+\sum _{i=1}^{n-1}(w_{n-i+1}-w_{n-i}){\widetilde{y}}_{m+1}^{i}-w_{n}{\widetilde{y}}_{m+1}^{0}+\frac{{\widetilde{y}}_{m+1}^{n+1}-{\widetilde{y}}_{m+1}^{n}}{2^{1-\gamma }}\right] \nonumber \\&\quad +\alpha \Bigg [b_{m-1}^{n+1/2}\left( \frac{({{\widetilde{y}}}_{x})_{m-1}^n+({{\widetilde{y}}}_{x})_{m-1}^{n+1}}{2}\right) +c_{m-1}^{n+1/2}\left( \frac{{{\widetilde{y}}}_{m-1}^n+{{\widetilde{y}}}_{m-1}^{n+1}}{2}\right) \nonumber \\&\quad -d_{m-1}^{n+1/2}\left( \frac{{{\widetilde{y}}}_{m-1}^{n-k}+{{\widetilde{y}}}_{m-1}^{n+1-k}}{2}\right) \Bigg ]\nonumber \\&\quad +\beta \Bigg [b_m^{n+1/2}\left( \frac{({{\widetilde{y}}}_{x})_m^n+({{\widetilde{y}}}_{x})_m^{n+1}}{2}\right) +c_m^{n+1/2}\left( \frac{{{\widetilde{y}}}_m^n+{{\widetilde{y}}}_m^{n+1}}{2}\right) \nonumber \\&\quad -d_m^{n+1/2}\left( \frac{{{\widetilde{y}}}_m^{n-k}+{{\widetilde{y}}}_m^{n+1-k}}{2}\right) \Bigg ]\nonumber \\&\quad +\alpha \Bigg [b_{m+1}^{n+1/2}\left( \frac{({{\widetilde{y}}}_{x})_{m+1}^n+({{\widetilde{y}}}_{x})_{m+1}^{n+1}}{2}\right) +c_{m+1}^{n+1/2}\left( \frac{{{\widetilde{y}}}_{m+1}^n+{{\widetilde{y}}}_{m+1}^{n+1}}{2}\right) \nonumber \\&\quad -d_{m+1}^{n+1/2}\left( \frac{{{\widetilde{y}}}_{m+1}^{n-k}+{{\widetilde{y}}}_{m+1}^{n+1-k}}{2}\right) \Bigg ]\nonumber \\&\quad +\left( (\alpha {\mathcal {M}}_{m-1}^{n+1/2}+\beta {\mathcal {M}}_m^{n+1/2}+\alpha {\mathcal {M}}_{m+1}^{n+1/2}\right) +\left( \alpha f_{m-1}^{n+1/2}+\beta f_{m}^{n+1/2}+\alpha f_{m+1}^{n+1/2}\right) \nonumber \\&=\ (2\alpha +\beta )({\mathcal {O}}(\Delta x)^2+{\mathcal {O}}(\Delta t)^2)\nonumber \\&\le \ C((\Delta x)^2+(\Delta t)^2). \end{aligned}$$

(3.2)

4 Stability analysis

In this section, we will discuss the stability of the numerical scheme (2.18) using the von Neumann stability method. We write the numerical scheme in the following form

$$\begin{aligned}&E(m,n)Y_{m-1}^{n+1}+F(m,n)Y_m^{n+1}+G(m,n)Y_{m+1}^{n+1}=E_1(m,n)Y_{m-1}^{n}+F_1(m,n)Y_m^{n}\nonumber \\&\quad +G_1(m,n)Y_{m+1}^{n}-\sigma w_1\left( \alpha Y_{m-1}^n+\beta Y_{m}^n+\alpha Y_{m+1}^n\right) +\sigma w_n\left( \alpha Y_{m-1}^0+\beta Y_{m}^0+\alpha Y_{m+1}^0\right) \nonumber \\&\quad -\sigma \sum _{q=1}^{n-1}(w_{n-q+1}-w_{n-q})\left( \alpha Y_{m-1}^q+\beta Y_{m}^q+\alpha Y_{m+1}^q\right) -\frac{1}{2}\Bigg [\alpha d_{m-1}^{n+1/2}(Y_{m-1}^{n+1-k}+Y_{m-1}^{n-k})\nonumber \\&\quad +\beta d_{m}^{n+1/2}(Y_{m}^{n+1-k}+Y_{m}^{n-k})+\alpha d_{m+1}^{n+1/2}(Y_{m+1}^{n+1-k}+Y_{m+1}^{n-k})\Bigg ]\nonumber \\&\quad +\left( \alpha f_{m-1}^{n+1/2}+\beta f_{m}^{n+1/2}+\alpha f_{m+1}^{n+1/2}\right) , \quad m=1,2,\ldots , N_x-1, \quad n=0,1,\ldots ,M_t-1, \end{aligned}$$

(4.1)

where

$$\begin{aligned}&E_1(m,n)=\alpha \left( \frac{\sigma }{2^{1-\gamma }}+\frac{b_{m-1}^{n+1/2}}{2\Delta x}-\frac{c_{m-1}^{n+1/2}}{2}\right) +\frac{\beta b_m^{n+1/2}}{4\Delta x}+\frac{1}{2(\Delta x)^2},\nonumber \\&F_1(m,n)=-\frac{\alpha b_{m-1}^{n+1/2}}{2\Delta x}+\beta \left( \frac{\sigma }{2^{1-\gamma }}-\frac{c_{m}^{n+1/2}}{2}\right) +\frac{\alpha b_{m+1}^{n+1/2}}{2\Delta x}-\frac{1}{(\Delta x)^2},\nonumber \\&G_1(m,n)=\alpha \left( \frac{\sigma }{2^{1-\gamma }}-\frac{b_{m+1}^{n+1/2}}{2\Delta x}-\frac{c_{m+1}^{n+1/2}}{2}\right) -\frac{\beta b_m^{n+1/2}}{4\Delta x}+\frac{1}{2(\Delta x)^2}. \end{aligned}$$

Let ${\hat{Y}}_m^n$ be the approximate solution of the system (2.18) and $\xi _m^n=Y_m^n-{\hat{Y}}_m^n$. The error equation of (4.1) is

$$\begin{aligned}&E(m,n)\xi _{m-1}^{n+1}+F(m,n)\xi _m^{n+1}+G(m,n)\xi _{m+1}^{n+1}=E_1(m,n)\xi _{m-1}^{n}+F_1(m,n)\xi _m^{n}\nonumber \\&\quad +G_1(m,n)\xi _{m+1}^{n}-\sigma w_1(\alpha \xi _{m-1}^n+\beta \xi _{m}^n+\alpha \xi _{m+1}^n)+\sigma w_n(\alpha \xi _{m-1}^0+\beta \xi _{m}^0+\alpha \xi _{m+1}^0) \nonumber \\&\quad -\sigma \sum _{q=1}^{n-1}(w_{n-q+1}-w_{n-q})\left( \alpha \xi _{m-1}^q+\beta \xi _{m}^q+\alpha \xi _{m+1}^q\right) \nonumber \\&\quad -\frac{1}{2}\Bigg [\alpha d_{m-1}^{n+1/2}(\xi _{m-1}^{n+1-k}+\xi _{m-1}^{n-k})+\beta d_{m}^{n+1/2}(\xi _{m}^{n+1-k}+\xi _{m}^{n-k})\nonumber \\&\quad +\alpha d_{m+1}^{n+1/2}(\xi _{m+1}^{n+1-k}+\xi _{m+1}^{n-k})\Bigg ],\nonumber \\&\quad m=1,2,\ldots ,N_x-1,\quad n=0,1,\ldots ,M_t-1, \end{aligned}$$

(4.2)

along with the conditions

$$\begin{aligned} \xi _0^n=\xi _{N_x}^n=0,\,n=0,1,\ldots ,M_t. \end{aligned}$$

The grid functions

$$\begin{aligned} \xi ^n(x)= {\left\{ \begin{array}{ll} \xi _m^n, &{}x_m-\frac{\Delta x}{2}<x<x_m+\frac{\Delta x}{2},\, m=1,2,\ldots , N_x-1,\\ 0, &{} 0\le x\le \frac{\Delta x}{2} \, \text {or} \, 1-\frac{\Delta x}{2}\le x\le 1, \end{array}\right. } \end{aligned}$$

have the following Fourier series expansion

$$\begin{aligned} \xi ^n(x)=\sum _{j=-\infty }^{\infty }\zeta ^n(j)e^{i(2j\pi x)},\ n=0,1,\ldots ,M_t, \end{aligned}$$

where

$$\begin{aligned} \zeta ^n(j)=\int _0^1 \xi ^n(x)e^{-i(2j\pi x)}dx. \end{aligned}$$

Using Parseval’s identity in $L_2-$norm, we find

$$\begin{aligned} ||\xi ^n||_2^2=\sum _{m=1}^{N_x-1}\Delta x|\xi _m^n|^2=\int _0^1|\xi ^n(x)|^2dx=\sum _{j=-\infty }^{\infty }|\zeta ^n(j)|^2. \end{aligned}$$

According to above analysis, we can assume that

$$\begin{aligned} \xi _m^n=\zeta ^ne^{i\theta m\Delta x}, \end{aligned}$$

where $i=\sqrt{-1}$ and $\theta$ is wave number, then from (4.2)

$$\begin{aligned}&E(m,n)\zeta ^{n+1} e^{i \theta (m-1)\Delta x}+F(m,n)\zeta ^{n+1} e^{i \theta m\Delta x}+G(m,n)\zeta ^{n+1} e^{i \theta (m+1)\Delta x}\nonumber \\&\quad =\ E_1(m,n)\zeta ^{n} e^{i \theta (m-1)\Delta x}+F_1(m,n)\zeta ^{n} e^{i \theta m\Delta x}+G_1(m,n)\zeta ^{n} e^{i \theta (m+1)\Delta x}\nonumber \\&\qquad -\sigma w_1\left( \alpha \zeta ^{n} e^{i \theta (m-1)\Delta x}+\beta \zeta ^{n} e^{i \theta m\Delta x}+\alpha \zeta ^{n} e^{i\theta (m+1)\Delta x}\right) \nonumber \\&\qquad +\sigma w_n\left( \alpha \zeta ^{0} e^{i \theta (m-1)\Delta x}+\beta \zeta ^{0} e^{i\theta m\Delta x}+\alpha \zeta ^{0} e^{i \theta (m+1)\Delta x}\right) \nonumber \\&\qquad -\sigma \sum _{q=1}^{n-1}(w_{n-q+1}-w_{n-q})\left( \alpha \zeta ^{q} e^{i\theta (m-1)\Delta x}+\beta \zeta ^{q} e^{i \theta m\Delta x}+\alpha \zeta ^{q} e^{i \theta (m+1)\Delta x}\right) \nonumber \\&\qquad -\frac{1}{2}\Bigg [\alpha d_{m-1}^{n+1/2}(\zeta ^{n+1-k} e^{i \theta (m-1)\Delta x}+\zeta ^{n-k} e^{i\theta (m-1)\Delta x})+\beta d_{m}^{n+1/2}(\zeta ^{n+1-k} e^{i \theta m\Delta x}\nonumber \\&\qquad +\zeta ^{n-k} e^{i \theta m\Delta x})\Bigg ]-\frac{1}{2}\alpha d_{m+1}^{n+1/2}(\zeta ^{n+1-k} e^{i\theta (m+1)\Delta x}+\zeta ^{n-k} e^{i \theta (m+1)\Delta x}). \end{aligned}$$

(4.3)

Dividing both sides by $e^{i \theta m\Delta x}$, we obtain

$$\begin{aligned}&\zeta ^{n+1}\Bigg [E(m,n) e^{-i \theta \Delta x}+F(m,n)+G(m,n) e^{i \theta \Delta x}\Bigg ]=\zeta ^{n}\Bigg [E_1(m,n) e^{-i \theta \Delta x}+F_1(m,n)\nonumber \\&\quad +G_1(m,n) e^{i \theta \Delta x}\Bigg ]-\sigma w_1\zeta ^{n}\left( \alpha e^{-i \theta \Delta x}+\beta +\alpha e^{i \theta \Delta x}\right) +\sigma w_n\zeta ^{0}\nonumber \\&\quad \left( \alpha e^{-i \theta \Delta x}+\beta +\alpha e^{i \theta \Delta x}\right) -\sigma \sum _{q=1}^{n-1}(w_{n-q+1}-w_{n-q})\zeta ^{q}\left( \alpha e^{-i \theta \Delta x}+\beta +\alpha e^{i \theta \Delta x}\right) \nonumber \\&\quad -\frac{1}{2}(\zeta ^{n+1-k} +\zeta ^{n-k} )\left( \alpha d_{m-1}^{n+1/2}e^{-i \theta \Delta x}+\beta d_{m}^{n+1/2}-\alpha d_{m+1}^{n+1/2} e^{i \theta \Delta x}\right) . \end{aligned}$$

(4.4)

After simplification for $n=0,1,\ldots ,k-1$, we get the following form

$$\begin{aligned} \zeta ^{n+1}=\left( \frac{U_1-w_1U}{U_2}\right) \zeta ^n-\frac{U}{U_2}\left( \sum _{q=1}^{n-1}(w_{n-q+1}-w_{n-q})\zeta ^{q}+w_n\zeta ^{0}\right) , \end{aligned}$$

(4.5a)

and for $n=k,k+1,\ldots ,M_t-1$, we have

$$\begin{aligned} \zeta ^{n+1}&=\left( \frac{U_1-w_1U}{U_2}\right) \zeta ^n-\frac{U}{U_2}\left( \sum _{q=1}^{n-1}(w_{n-q+1}-w_{n-q})\zeta ^{q}+w_n\zeta ^{0}\right) \nonumber \\&\quad -\frac{1}{2U_2}(\zeta ^{n+1-k} +\zeta ^{n-k})\left( \alpha d_{m-1}^{n+1/2}e^{-i \theta \Delta x}+\beta d_{m}^{n+1/2}-\alpha d_{m+1}^{n+1/2} e^{i \theta \Delta x}\right) , \end{aligned}$$

(4.5b)

where

$$\begin{aligned} U&=\sigma (2\alpha \cos (\theta \Delta x)+\beta ),\\ U_1&=E_1(m,n) e^{-i \theta \Delta x}+F_1(m,n)+G_1(m,n) e^{i \theta \Delta x},\\ U_2&=E(m,n) e^{-i \theta \Delta x}+F(m,n)+G(m,n) e^{i \theta \Delta x}. \end{aligned}$$

For stability of numerical scheme (2.18), we prove that $|\zeta ^{n+1}|\le |\zeta ^0|$ by use of mathematical induction. For $n=0$, (4.5a) yields

$$\begin{aligned} |\zeta ^{1}|&= \left| \left( \frac{U_1-w_1U}{U_2}\right) \zeta ^0\right| \\&= \left| \frac{(2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}-\sigma w_1\right) +(1-\cos (\theta \Delta x))(-X_1)-X_2-i X_3\sin (\theta \Delta x)}{(2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}\right) +(1-\cos (\theta \Delta x))(X_1)+X_2+iX_3\sin (\theta \Delta x)}\right| |\zeta ^0|\\&\le \frac{\left| (2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}-\sigma w_1\right) +(1-\cos (\theta \Delta x))(-X_1)-X_2\right| +X_3\sin (\theta \Delta x)}{\left| (2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}\right) +(1-\cos (\theta \Delta x))(X_1)+X_2+i X_3\sin (\theta \Delta x)\right| }|\zeta ^0|\\&=R|\zeta ^0|, \end{aligned}$$

where

$$\begin{aligned} X_1&= \frac{1}{(\Delta x)^2}+\frac{\alpha }{2\Delta x}b_{m-1}^{1/2}-\frac{\alpha }{2\Delta x}b_{m+1}^{1/2},\\ X_2&= \frac{\beta }{2}c_m^{1/2}+\frac{\alpha \cos (\theta \Delta x)}{2}(c_{m-1}^{1/2}+c_{m+1}^{1/2}),\\ X_3&= \frac{\beta }{2\Delta x}b_{m}^{1/2}+\frac{\alpha }{2\Delta x}b_{m-1}^{1/2}+\frac{\alpha }{2\Delta x}b_{m+1}^{1/2}-\frac{\alpha }{2}c_{m-1}^{1/2}+\frac{\alpha }{2}c_{m+1}^{1/2},\\ R&= \frac{\left| (2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}-\sigma w_1\right) -(1-\cos (\theta \Delta x))X_1-X_2\right| +X_3\sin (\theta \Delta x)}{\left| (2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}\right) +(1-\cos (\theta \Delta x))X_1+X_2+i X_3\sin (\theta \Delta x)\right| }. \end{aligned}$$

To show that $R\le 1$, we consider two cases.

Case 1. If $(2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}-\sigma w_1\right) -(1-\cos (\theta \Delta x))X_1-X_2\ge 0$. Then $R\le 1$ gives,

$$\begin{aligned} \frac{\left| (2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}-\sigma w_1\right) -(1-\cos (\theta \Delta x))X_1-X_2\right| +X_3\sin (\theta \Delta x)}{\left| (2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}\right) +(1-\cos (\theta \Delta x))X_1+X_2+ iX_3\sin (\theta \Delta x)\right| }\le 1, \end{aligned}$$

which gives

$$\begin{aligned}&(2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}-\sigma w_1\right) -(1-\cos (\theta \Delta x))X_1-X_2\\&\qquad \le (2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}\right) +(1-\cos (\theta \Delta x))X_1+X_2, \end{aligned}$$

if and only if

$$\begin{aligned} (2\alpha \cos (\theta \Delta x)+\beta )\left( \sigma w_1\right) +2(1-\cos (\theta \Delta x))X_1+2X_2\ge 0, \end{aligned}$$

which is true for all $\theta$ and $\Delta x$.

Case 2. If $(2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}-\sigma w_1\right) -(1-\cos (\theta \Delta x))X_1-X_2\le 0$. Then $R\le 1$ gives

$$\begin{aligned} \frac{\left| (2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}-\sigma w_1\right) -(1-\cos (\theta \Delta x))X_1-X_2\right| +X_3\sin (\theta \Delta x)}{\left| (2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}\right) +(1-\cos (\theta \Delta x))X_1+X_2+i X_3\sin (\theta \Delta x)\right| }\le 1, \end{aligned}$$

which gives

$$\begin{aligned}&-(2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}-\sigma w_1\right) +(1-\cos (\theta \Delta x))X_1+X_2\\&\qquad \le (2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}\right) +(1-\cos (\theta \Delta x))X_1+X_2, \end{aligned}$$

which conclude that $w_1\le \frac{2}{2^{1-\gamma }}$, true for all values of $\gamma$. Now we assume that

$$\begin{aligned} |\zeta ^j|\le |\zeta ^0|,\ \text {for } j=2,3,\ldots , n. \end{aligned}$$

(4.6)

We will prove that (4.6) holds for $j=n+1$. The equation (4.5a) yields

$$\begin{aligned} |\zeta ^{n+1}|\le&\ \left( \frac{|U_1-w_1U|+w_1|U|}{|U_2|}\right) |\zeta ^0|\\ =&\ \left| \frac{(2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}-\sigma w_1\right) -(1-\cos (\theta \Delta x))X_1-X_2-iX_3\sin (\theta \Delta x)}{(2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}\right) +(1-\cos (\theta \Delta x))X_1+X_2+iX_3\sin (\theta \Delta x)}\right| |\zeta ^0|\\&+\frac{\sigma w_1(2\alpha \cos (\theta \Delta x)+\beta )}{\left| (2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}\right) +(1-\cos (\theta \Delta x))X_1+X_2+i X_3\sin (\theta \Delta x)\right| }|\zeta ^0|\\ \le&\ \frac{\left| (2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}-\sigma w_1\right) -(1-\cos (\theta \Delta x))X_1-X_2\right| +X_3\sin (\theta \Delta x)}{\left| (2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}\right) +(1-\cos (\theta \Delta x))X_1+X_2+i X_3\sin (\theta \Delta x)\right| }|\zeta ^0|\\&+\frac{\sigma w_1(2\alpha \cos (\theta \Delta x)+\beta )}{\left| (2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}\right) +(1-\cos (\theta \Delta x))X_1+X_2+i X_3\sin (\theta \Delta x)\right| }|\zeta ^0|. \end{aligned}$$

Again, there arise two cases.

Case 1. If $(2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}-\sigma w_1\right) -(1-\cos (\theta \Delta x))X_1-X_2\ge 0$, then

$$\begin{aligned} |\zeta ^{n+1}|&\le \frac{(2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}\right) -(1-\cos (\theta \Delta x))X_1-X_2+X_3\sin (\theta \Delta x)}{\left| (2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}\right) +(1-\cos (\theta \Delta x))X_1+X_2+i X_3\sin (\theta \Delta x)\right| }|\zeta ^0|\\&\le |\zeta ^0|. \end{aligned}$$

Case 2. If $(2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}-\sigma w_1\right) -(1-\cos (\theta \Delta x))X_1-X_2\le 0$, then

$$\begin{aligned} |\zeta ^{n+1}|\le \frac{-(2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}-2\sigma w_1\right) +(1-\cos (\theta \Delta x))X_1+X_2+X_3\sin (\theta \Delta x)}{\left| (2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}\right) +(1-\cos (\theta \Delta x))X_1+X_2+i X_3\sin (\theta \Delta x)\right| }|\zeta ^0|. \end{aligned}$$

Now for $|\zeta ^{n+1}|\le |\zeta ^0|$, we have

$$\begin{aligned} \frac{-(2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}-2\sigma w_1\right) +(1-\cos (\theta \Delta x))X_1+X_2+X_3\sin (\theta \Delta x)}{\left| (2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}\right) +(1-\cos (\theta \Delta x))X_1+X_2+i X_3\sin (\theta \Delta x)\right| }\le 1, \end{aligned}$$

that gives

$$\begin{aligned} -(2\alpha \cos (\theta \Delta x)&+\beta )\left( \frac{\sigma }{2^{1-\gamma }}-2\sigma w_1\right) +(1-\cos (\theta \Delta x))X_1+X_2\\&\le (2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}\right) +(1-\cos (\theta \Delta x))X_1+X_2. \end{aligned}$$

Last inequality holds true when $\frac{1}{2^{1-\gamma }}\ge w_1$ or $3^\gamma \ge \frac{3}{2}$. Now for $n=k,k+1,\ldots ,M_t-1$, the Eq. (4.5b) becomes

$$\begin{aligned} |\zeta ^{n+1}|\le&\ \left( \frac{|U_1-w_1U|+w_1|U|+|\alpha d_{m-1}^{n+1/2}e^{-i \theta \Delta x}+\beta d_{m}^{n+1/2}-\alpha d_{m+1}^{n+1/2} e^{i \theta \Delta x}|}{|U_2|}\right) |\zeta ^0|\\ =&\ \left| \frac{(2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}-\sigma w_1\right) -(1-\cos (\theta \Delta x))X_1-X_2-i X_3\sin (\theta \Delta x)}{(2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}\right) +(1-\cos (\theta \Delta x))X_1+X_2+iX_3\sin (\theta \Delta x)}\right| |\zeta ^0|\\&+\frac{\sigma w_1(2\alpha \cos (\theta \Delta x)+\beta )+|\alpha d_{m-1}^{n+1/2}e^{-i \theta \Delta x}+\beta d_{m}^{n+1/2}-\alpha d_{m+1}^{n+1/2} e^{i \theta \Delta x}|}{\left| (2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}\right) +(1-\cos (\theta \Delta x))X_1+X_2+i X_3\sin (\theta \Delta x)\right| }|\zeta ^0|. \end{aligned}$$

It is clear from previous analysis and mathematical induction that as $\Delta x,\Delta t\rightarrow 0$, $|\zeta ^{n+1}|\le |\zeta ^0|$ for each $n=0,1,\ldots , M_t-1$. Hence, the numerical scheme (2.18) is conditionally stable with the condition $3^\gamma \ge \frac{3}{2}$.

5 Convergence analysis

To prove the convergence of the numerical scheme (2.18), we assume

$$\begin{aligned} E_m^n = {\widetilde{y}}_m^n-Y_m^n,\, m=0,1,\ldots ,N_x,\, n=0,1,\ldots ,M_t. \end{aligned}$$

(5.1)

From (2.17) and (2.18), we obtain

$$\begin{aligned} E(m,n)&E_{m-1}^{n+1}+F(m,n)E_m^{n+1}+G(m,n)E_{m+1}^{n+1}=E_1(m,n)E_{m-1}^{n}+F_1(m,n)E_m^{n}\nonumber \\&+G_1(m,n)E_{m+1}^{n}-\sigma w_1\left( \alpha E_{m-1}^n+\beta E_{m}^n+\alpha E_{m+1}^n\right) +\sigma w_n\nonumber \\&\times \left( \alpha E_{m-1}^0+\beta E_{m}^0+\alpha E_{m+1}^0\right) -\sigma \sum _{q=1}^{n-1}(w_{n-q+1}-w_{n-q})\left( \alpha E_{m-1}^q+\beta E_{m}^q+\alpha E_{m+1}^q\right) \nonumber \\&-\frac{1}{2}\Bigg [\alpha d_{m-1}^{n+1/2}(E_{m-1}^{n+1-k}+E_{m-1}^{n-k})+\beta d_{m}^{n+1/2}(E_{m}^{n+1-k}+E_{m}^{n-k})+\alpha d_{m+1}^{n+1/2}\nonumber \\&\times (E_{m+1}^{n+1-k}+E_{m+1}^{n-k})\Bigg ]+\left( \alpha {\mathfrak {R}}_{m-1}^{n+1/2}+\beta {\mathfrak {R}}_{m}^{n+1/2}+\alpha {\mathfrak {R}}_{m+1}^{n+1/2}\right) , \end{aligned}$$

(5.2)

along with the conditions

$$\begin{aligned} E_0^n=E_{N_x}^n=0,\, \text {and} \, E_m^0=0, \end{aligned}$$

where

$$\begin{aligned} |{\mathfrak {R}}_m^{n+1/2}|\le C{\mathfrak {e}}_m^{n+1/2}\le C ((\Delta x)^2+(\Delta t)^2). \end{aligned}$$

(5.3)

The grid functions

$$\begin{aligned} E^n(x)= {\left\{ \begin{array}{ll} E_m^n, &{}x_m-\frac{\Delta x}{2}<x<x_m+\frac{\Delta x}{2},\, m=1,2,\ldots , N_x-1,\\ 0, &{} 0\le x\le \frac{\Delta x}{2} \, \text {or} \, 1-\frac{\Delta x}{2}\le x\le 1, \end{array}\right. } \end{aligned}$$

and

$$\begin{aligned} {\mathfrak {R}}^{n+1/2}(x)= {\left\{ \begin{array}{ll} {\mathfrak {R}}_m^{n+1/2}, &{}x_m-\frac{\Delta x}{2}<x<x_m+\frac{\Delta x}{2},\, m=1,2,\ldots , N_x-1,\\ 0, &{} 0\le x\le \frac{\Delta x}{2} \, \text {or} \, 1-\frac{\Delta x}{2}\le x\le 1, \end{array}\right. } \end{aligned}$$

have the following Fourier series expansions

$$\begin{aligned} E^n(x)=\sum _{j=-\infty }^{\infty }\chi ^n(j)e^{i(2j\pi x)},\, n=0,1,\ldots ,M_t, \end{aligned}$$

and

$$\begin{aligned} {\mathfrak {R}}^{n+1/2}(x)=\sum _{j=-\infty }^{\infty }r^{n+1/2}(j)e^{i(2j\pi x)},\,n=0,1,\ldots ,M_t, \end{aligned}$$

respectively, where

$$\begin{aligned} \chi ^n(j)=\int _0^1 E^n(x)e^{-i(2j\pi x)}dx, \end{aligned}$$

and

$$\begin{aligned}r^{n+1/2}(j)=\int _0^1 {\mathfrak {R}}^{n+1/2}(x)e^{-i(2j\pi x)}dx.\end{aligned}$$

Using Parseval’s identity in $L_2-$norm, we find

$$\begin{aligned} ||E^n||_2^2=\sum _{m=1}^{N_x-1}\Delta x|E_m^n|^2=\int _0^1|E^n(x)|^2dx=\sum _{j=-\infty }^{\infty }|\chi ^n(j)|^2, \end{aligned}$$

(5.4)

and

$$\begin{aligned} ||{\mathfrak {R}}^{n+1/2}||_2^2=\sum _{m=1}^{N_x-1}\Delta x|{\mathfrak {R}}_m^{n+1/2}|^2=\int _0^1|{\mathfrak {R}}^{n+1/2}(x)|^2dx=\sum _{j=-\infty }^{\infty }|r^{n+1/2}(j)|^2. \end{aligned}$$

(5.5)

Now, we assume that

$$\begin{aligned} E_m^n=\chi ^ne^{i\theta m\Delta x}, \end{aligned}$$

(5.6)

and

$$\begin{aligned} {\mathfrak {R}}_m^{n+1/2}=r^{n+1/2}e^{i\theta m \Delta x}, \end{aligned}$$

(5.7)

where $\theta =2\pi j$. Using Eqs. (5.2) and (5.6) in (5.7) for $n=0,1,\ldots ,k-1$, we have

$$\begin{aligned} \chi ^{n+1}=\left( \frac{U_1-w_1U}{U_2}\right) \chi ^n-\frac{U}{U_2}\left( \sum _{q=1}^{n-1}(w_{n-q+1}-w_{n-q})\chi ^{q}+w_n\chi ^{0}\right) +\frac{Ur^{n+1/2}}{\sigma U_2}, \end{aligned}$$

(5.8a)

and for $n=k,k+1,\ldots ,M_t-1$, we have

$$\begin{aligned} \chi ^{n+1}=&\ \left( \frac{U_1-w_1U}{U_2}\right) \chi ^n-\frac{U}{U_2}\left( \sum _{q=1}^{n-1}(w_{n-q+1}-w_{n-q})\chi ^{q}+w_n\chi ^{0}\right) \nonumber \\&-\frac{1}{2U_2}(\chi ^{n+1-k} +\chi ^{n-k} )\left( \alpha d_{m-1}^{n+1/2}e^{-i \theta \Delta x}+\beta d_{m}^{n+1/2}-\alpha d_{m+1}^{n+1/2} e^{i \theta \Delta x}\right) +\frac{Ur^{n+1/2}}{\sigma U_2}. \end{aligned}$$

(5.8b)

Proposition 5.1

There exist a positive constant C such that $|\chi ^{n+1}|\le (1+C\Delta t)^{n+1}|r^{1/2}|,\, n=0,1,\ldots ,M_t-1$.

Proof

Clearly $\chi ^0=\chi ^0(j)=0$ (as $E^0=0$). From the convergence of the series in right hand side of the Eq. (5.5), there exist a positive constant C such that $|r^{n+1/2}|\le C\Delta t|r^{1/2}|,\, n=0,1,\ldots ,M_t-1$. From (5.3) and (5.5), we get

$$\begin{aligned} ||{\mathfrak {R}}^{n+1/2}||_2=C((\Delta x)^2+(\Delta t)^2), \, n=0,1,\ldots ,M_t-1,\end{aligned}$$

by the use of mathematical induction, for $n=0$, we have

$$\begin{aligned} |\chi ^1|&=\left| \frac{(2\alpha \cos (\theta (\Delta x))+\beta )r^{1/2}}{(2\alpha \cos (\theta (\Delta x))+\beta )\left( \frac{\sigma }{2^{(1-\gamma )}}\right) +(1-\cos (\theta (\Delta x)))X_1+X_2+iX_3\sin (\theta (\Delta x))}\right| \\&\le C\Delta t|r^{1/2}|\le (1+C\Delta t)|r^{1/2}|. \end{aligned}$$

Now assuming that

$$\begin{aligned} |\chi ^\wp |\le (1+C\Delta t)^\wp |r^{1/2}|,\,\wp =2,3,\ldots ,n, \end{aligned}$$

(5.9)

we will prove that (5.9) is true for $\wp =n+1$. Equation (5.8a) yields

$$\begin{aligned} |\chi ^{n+1}|\le&\left( \frac{|U_1-w_1U|+w_1|U|}{|U_2|}\right) (1+C\Delta t)^n|r^{1/2}|\\&+\left| \frac{(2\alpha \cos (\theta \Delta x)+\beta )r^{1/2}}{(2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}\right) +(1-\cos (\theta \Delta x))X_1+X_2+i X_3\sin (\theta \Delta x)}\right| , \end{aligned}$$

on the basis of proof of stability it can be easily prove that

$$\begin{aligned}|\chi ^{n+1}|\le (1+C\Delta t)^{n+1}|r^{1/2}|\ \text {for}\, n=0,1,\ldots , k-1.\end{aligned}$$

For $n=k,k+1,\dots ,M_t-1$, Eq. (5.8b) yields

$$\begin{aligned} |\chi ^{n+1}|\le&\left( \frac{|U_1-w_1U|+w_1|U|+|\alpha d_{m-1}^{n+1/2}e^{-i \theta \Delta x}+\beta d_{m}^{n+1/2}-\alpha d_{m+1}^{n+1/2} e^{i \theta \Delta x}|}{|U_2|}\right) \\&\qquad (1+C\Delta t)^n|r^{1/2}|\\&+\left| \frac{(2\alpha \cos (\theta \Delta x)+\beta )r^{1/2}}{(2\alpha \cos (\theta \Delta x)+\beta )\left( \frac{\sigma }{2^{1-\gamma }}\right) +(1-\cos (\theta \Delta x))X_1+X_2+i X_3\sin (\theta \Delta x)}\right| . \end{aligned}$$

Similarly, as in the stability, we can prove that

$$\begin{aligned} |\chi ^{n+1}|\le (1+C\Delta t)^{n+1}|r^{1/2}|\ \text {for}\, {n=k,k+1,\ldots ,M_t-1}. \end{aligned}$$

Hence, by the mathematical induction we can say that $|\chi ^{n}|\le (1+C\Delta t)^{n}|r^{1/2}|\ \text {for}\, n=1,2,\ldots ,M_t$. $\square$

Theorem 5.1

The numerical scheme (2.18) is convergent and the solution satisfies

$$\begin{aligned} ||{\widetilde{y}}_m^n-Y_m^n||_2\le C((\Delta x)^2+(\Delta t)^2),\, m=1,2,\ldots ,N_x,\, n=1,2,\ldots ,M_t, \end{aligned}$$

where C is a positive constant independent of $\Delta x$ and $\Delta t$.

Proof

From Eqs. (5.4) and (5.5) and proposition (5.1), we obtain

$$\begin{aligned} ||E^n||_2\le (1+C\Delta t)^n||R^{1/2}||_2\le Ce^{(Cn\Delta t)}((\Delta x)^2+(\Delta t)^2), \,n=1,2,\ldots ,M_t, \end{aligned}$$

then

$$\begin{aligned} ||E^n||_2\le C((\Delta x)^2+(\Delta t)^2),\,n=1,2,\ldots ,M_t. \end{aligned}$$

Thus, the numerical scheme (2.18) is second-order convergent. $\square$

6 Numerical illustrations

In this section, the efficiency and effectiveness of the proposed technique are demonstrated by implementing our numerical method on two test examples for the TF-DPDEs. The exact solutions to these problems are not available, so we employ the double mesh principle to determine the errors in $L_\infty$ and $L_2-$norms as follows:

$$\begin{aligned} {\mathcal {E}}^{N_{x},M_{t}}_{\infty }&=\max _{1\le n\le M_t}\left( \max _{1 \le m\le N_{x}-1}|{\widetilde{y}}_{i}^{j}-{\widetilde{y}}_{2i}^{2j}|\right) ,\\ {\mathcal {E}}^{N_{x},M_{t}}_{2}&=\max _{1\le n\le M_t}\sqrt{\Delta x\sum _{m=1}^{N_{x}-1}({\widetilde{y}}_{i}^{j}-{\widetilde{y}}_{2i}^{2j})^{2}}, \end{aligned}$$

where ${\widetilde{y}}_{i}^{j}$ and ${\widetilde{y}}_{2i}^{2j}$ are the numerical solutions obtained using ($N_{x}+1$,$M_{t}+1$) and ($2N_{x}+1$,$2M_{t}+1$) points, respectively. We also compute the corresponding orders of convergence using the formula

$$\begin{aligned} {\mathcal {Q}}_{q}^{N_x,M_t}=\log _2\left( \frac{{\mathcal {E}}_{q}^{N_x,M_t}}{{\mathcal {E}}_{q}^{2N_x,2M_t}}\right) , \quad q = 2,\,\infty . \end{aligned}$$

We consider the following TF-DPDEs with $\tau =1$ and suitable interval and boundary conditions.

Example 6.1

$$\begin{aligned}&{}_0^CD^{\gamma }_t y-y_{xx}+(2-x^2)y_x+xy = y(x,t-\tau )+10t^2e^{-t}(x-x^2), \, (x,t)\in (0,1)\times (0,2],\\&\psi _b(x,t)=0, \ \psi _l(t)=0, \ \psi _r(t)=0. \end{aligned}$$

Table 1 $L_2$ and $L_{\infty }$ errors, and respective orders of convergence for Example 6.1

Fig. 1

Numerical solution for Example 6.1 for different values of $\gamma$

Fig. 2

Graph of numerical solution for Example 6.1 at different time levels

Fig. 3

Absolute error for Example 6.1

Fig. 4

Numerical solution for Example 6.2 for different values of $\gamma$

Example 6.2

$$\begin{aligned}&{}_0^CD^{\gamma }_t y-y_{xx}+e^xy_x+(x+1)(t+1)y = y(x,t-\tau )+10t^2e^{-t}(2x-x^2),\\&\, (x,t)\in (0,1)\times (0,2],\\&\psi _b(x,t)=0, \ \psi _l(t)=0, \ \psi _r(t)=0. \end{aligned}$$

Replacing the $\gamma$ order fractional derivative by first-order integer derivative, we define the following problem

$$\begin{aligned} y_{t}-y_{xx}+b(x,t)y_{x}+c(x,t)y(x,t) = -d(x,t)y(x,t-\tau )+f(x,t), \end{aligned}$$

(6.1)

where b(x, t), c(x, t), d(x, t), f(x, t), and the interval and boundary conditions are the same as defined in Examples 6.1 and 6.2 .

We have shown the errors in $L_2$ and $L_\infty$ norms for different values of $\gamma$, in Tables 1 and 3, for Examples 6.1 and 6.2 , respectively. From these tables, it can be observed that as we increase $N_{x}$ and $M_{t}$, the errors decrease and are in good agreement with the theoretical results. It confirms our theoretical estimates that the suggested scheme is second-order convergent. From Tables 2 and 4 , one can observe that the errors in $L_2$ and $L_{\infty }$ norms for Examples 6.1 and 6.2 approach to the errors for the problem (6.1) as $\gamma$ approaches 1. In Figures 1 and 4, we portray the numerical solution profiles of Examples 6.1 and 6.2 , respectively for different values of $\gamma$. Figures 2 and 5 exhibit the graphs of the numerical solutions of Examples 6.1 and 6.2 at different time-levels. Keeping $\gamma$ fixed, Figures 3 ($\gamma = 0.9$) and 6 ($\gamma = 0.8$) show that the error decrease when we increase the number of points in both directions.

Remark 2

To draw the Figs. 1, 2, 4, and 5 , we have used $N_x=M_t=64$.

Table 2 Comparison of $L_2$ and $L_{\infty }$ errors between Example 6.1 and problem (6.1)

Table 3 $L_2$ and $L_{\infty }$ errors, and respective orders of convergence for Example 6.2

Table 4 Comparison of $L_2$ and $L_{\infty }$ errors between Example  6.2 and problem (6.1)

Fig. 5

Graph of numerical solution for Example 6.2 at different time levels

Fig. 6

Absolute error for Example 6.2

7 Concluding remarks

This paper suggests a numerical method comprising Crank-Nicolson scheme with tension spline for TF-DPDEs. The Caputo fractional derivative is used to discretize the fractional-order time derivative. Using the Fourier series analysis, the method is proved conditionally stable. Moreover, through rigorous analysis, the method is shown second-order convergent for arbitrary suitable choices of $\alpha \, (\alpha \ne 1/12)$ and $\beta$ with $2\alpha + \beta = 1$. The method is easy to apply to the problem (1.1) and can easily be implemented through MATLAB software. Numerical examples demonstrate the effectiveness and adaptability of the suggested method. The suggested technique can be extended to the nonlinear TF-DPDEs and system of FPDEs.