A Fitted Numerical Scheme with Algorithm for Singularly Perturbed Parabolic Partial Differential Equations with Large Negative Shift

1. Introduction

A class of differential equations, which involve a perturbation parameter $\epsilon$ as a coefficient of the term with highest order derivative and and involves one or more shift arguments are known as singularly perturbed delay differential equations. Such equations have found applications in various areas of sciences and engineering where time scales and delays play a crucial role. For instance, singularly perturbed differential equations with delay are used to study the stability and transient behavior of fluid flows [1]. In biofluid dynamics, singularly perturbed differential equations with delay are employed to model physiological or pathological processes that involve delays in fluid dynamics [2]. In bistable devices, these problems are used to model the dynamics of the circuits, accounting for the time delays in feedback or input signal propagation [3]. They are also employed to model chemical reaction-diffusion systems in fluid mechanics which involve delayed reactions where the reaction rates depend on the concentrations of reactants at previous time points [4]. Singularly perturbed differential equations with delay play an important role in the analysis and design of control systems for fluid dynamics [5]. These problems are utilized in modeling neuronal networks, control systems and signal processing applications. The models are account for the delays in signal transmission, feedback loops or information processing in biological or artificial systems aiding in understanding system behaviors and designing control strategies [6].

The solution involves boundary layers because of the influence of $\epsilon$ and inner layer due the negative shift term. The behavior of layers causes difficulties to find analytical solution. On the other hand, the standard numerical techniques do not yield adequate solution as they do not consider the behavior of layers. As a result, such problems require parameter independent numerical techniques.

To overcome such computational difficulties, various numerical schemes are developed for different forms of SP-PDEs. For instance, Kumar and Kumari [7] treated singular perturbation problems with time-lag by constructing a numerical scheme using the extended cubic B-spline basis functions involving a free parameter $\lambda$ on a fitted mesh. In [8], an exponentially fitted scheme is formulated for a singularly perturbed partial differential equations involving small delay on the space variable. In the articles [9,10,11], various numerical schemes are developed for singularly perturbed partial differential equations consisting of delay in the spatial direction, and obtained that the methods are parameter-uniformly convergent. Sharma and Kaushik [12] solved singular perturbation problems with a large shift and integral boundary condition by developing a higher-order numerical scheme on a specially designed non-uniform grids. In [13], a nonstandard difference scheme is formulated for singularly perturbed partial differential equations involving shift, and found that the scheme is convergent regardless of $\epsilon$ and it is of order two both in time and space variables.

In this article, we presented an exponentially fitted tension spline method with computational algorithms for a SP-PDE involving large shift. The scheme is constructed applying implicit Euler finite difference technique in the temporal variable and an exponentially fitted tension spline function is fitted in the space variable on equally spaced grids. A computational algorithm for the developed scheme is also included in this work. The convergence analysis of the method is examined and proved that the scheme is uniformly convergent of first order in time and second order in space. The developed scheme and technique are confirmed with various forms of numerical experiments.

The remaining sections of the article are outlined below: Section 2 presents the preliminaries and a priori estimates of the continuous problem. Numerical formulation and computational algorithm of the study are presented in Section 3. To test the developed scheme, numerical experiments and discussions are provided in Section 4. The finding of this study is concluded in Section 5.

Notation: In this work, we represented C as a generic constant independent of $\epsilon$ and the mesh numbers. And we used the norm $\parallel .\parallel$ to be $\parallel .\parallel ={max}_{\overline{\Xi }}|.|$ for any function in $\overline{\Xi }$.

2. Preliminaries and a priori estimates

Consider a SP-PDE on $\Xi =[0,2]\times [0,T]$ given by

$$\begin{array}{cc}\hfill L_{\epsilon }\Im (n,t)& ={\Im }_t(n,t)-\epsilon {\Im }_{nn}(n,t)+\mu \left(n\right)\Im (n,t)+\nu \left(n\right)\Im (n-1,t)=\gamma (n,t),\hfill \end{array}$$

(2.1)

$$\begin{array}{cc}\hfill & \left\{\begin{array}{c}\Im (n,0)=\Im \left(n\right),n\in [0,2],\hfill \\ \Im (n,t)=\omega (n,t),(n,t)\in {\Xi }^{-},\hfill \\ \Im (2,t)=\lambda \left(t\right),(2,t)\in {\Xi }^{+},\hfill \end{array}\right.\hfill \end{array}$$

(2.2)

where $L_{\epsilon }$ is a differential operator, ${\Xi }^{-}=\{(n,t):(n,t)\in [-1,0]\times [0,T]\}$, ${\Xi }^{+}=\{(2,t),t\in [0,T]\}$, $0<\epsilon \ll 1$ and T is a fixed finite time. The functions $\mu \left(n\right)$, $\nu \left(n\right)$, $\gamma (n,t)$, $\omega (n,t)$ and $\lambda \left(t\right)$ are considered to be continuous. Moreover, to obtain non-oscillatory solution, we assume that the coefficient functions satisfy

$$\mu \left(n\right)+\nu \left(n\right)\ge \beta >0\mathrm{and}\nu \left(n\right)<0\mathrm{for}n\in [0,2].$$

(2.3)

To obtain a unique solution of (2.1)-(2.2), we assume the Holder continuity of the given data and impose compatibility conditions given as

$${\Im }_0(0,0)=\omega (0,0),{\Im }_0(2,0)=\lambda (2,0)$$

(2.4)

and

$$\left\{\begin{array}{c}{\omega }_t(0,0)-\epsilon {\left({\omega }_0\right)}_{nn}\left(0\right)+\mu \left(0\right){\Im }_0\left(0\right)+\nu \left(0\right)\varphi (-1,0)=\gamma (0,0),\hfill \\ {\lambda }_t(2,0)-\epsilon {\left({\Im }_0\right)}_{nn}\left(2\right)+\mu \left(2\right){\Im }_0\left(2\right)+\nu \left(2\right){\Im }_0(1,0)=\gamma (2,0).\hfill \end{array}\right.$$

(2.5)

Considering the above basic assumptions, the following estimates hold true for the solutions of the continuous problem.

Lemma 2.1. 

Let $\Im \in C^{2,1}\left(\overline{\Xi }\right)$ satisfying $\Im (n,t)\ge 0,(n,t)\in \partial \Xi$ and $L_{\epsilon }\Im (n,t)\ge 0,(n,t)\in \Xi$. Then $\Im (n,t)\ge 0,(n,t)\in \overline{\Xi }$.

Proof. 

Suppose that $(n^{*},t^{*})\in \overline{\Xi }$ such that $\Im (n^{*},t^{*})={min}_{\overline{\Xi }}\Im (n,t)<0$. Then,

Case 1: If $n\in (0,1]$, then we get

$$\begin{array}{cc}\hfill L_{\epsilon ,1}\Im (n^{*},t^{*})=& {\Im }_t(n^{*},t^{*})-\epsilon {\Im }_{nn}(n^{*},t^{*})+\mu \left(n^{*}\right)\Im (n^{*},t^{*})\hfill \\ \hfill & -\epsilon {\Im }_{nn}(n^{*},t^{*})+\mu \left(n^{*}\right)\Im (n^{*},t^{*})\le 0.\hfill \end{array}$$

Case 2: If $n\in (1,2)$, then we obtain

$$\begin{array}{cc}\hfill L_{\epsilon ,2}\Im (n^{*},t^{*})=& {\Im }_t(n^{*},t^{*})-\epsilon {\Im }_{nn}(n^{*},t^{*})+\mu \left(n^{*}\right)\Im (n^{*},t^{*})+\nu \left(n^{*}\right)\Im (n^{*}-1,t^{*})\hfill \\ \hfill & \le 0\hfill \end{array}$$

A contradiction is occurred from the two cases. Thus, $\Im (n^{*},t^{*})\ge 0$ and consequently $\Im (n,t)\ge 0,(n,t)\in \overline{\Xi }$.    □

Lemma 2.2. 

The stability estimate of the solution $\Im (n,t)$ of (2.1) is given as

$$|\Im (n,t)|\le {\beta }^{-1}\parallel L_{\epsilon }\Im \parallel +max\left\{|{\Im }_0\left(n\right)|,|\omega (0,t)|,|\Im (2,t)\left|\right\}\right\}.$$

Proof. 

Introduce ${\varpi }^{\pm }(n,t)={\beta }^{-1}\parallel L_{\epsilon }\Im \parallel +max\left\{\right|{\Im }_0\left(n\right)|,|\omega (0,t)|,|\Im (2,t)\left|\right\}\pm \Im (n,t)$. Then, obviously we find that ${\varpi }^{\pm }(0,t)\ge 0$ and ${\varpi }^{\pm }(2,t)\ge 0$.

$$\begin{array}{cc}\hfill \mathrm{For}0<n\le 1,& L_{\epsilon }{\varpi }^{\pm }(n,t)={\varpi }_t^{\pm }-\epsilon {\varpi }_{nn}^{\pm }+\mu \left(n\right){\varpi }^{\pm }(n,t)+\nu \left(n\right){\omega }^{\pm }(n-1,t)\hfill \\ \hfill & \ge \beta max\left\{\right|{\Im }_0\left(n\right)|,|\omega (0,t)|,|\Im (2,t)\}\pm \gamma (n,t)\ge 0.\hfill \\ \hfill \mathrm{For}1<n<2,& L_{\epsilon }{\varpi }^{\pm }(n,t)={\varpi }_t^{\pm }-\epsilon {\varpi }_{nn}^{\pm }+\mu \left(n\right){\varpi }^{\pm }(n,t)+\nu \left(n\right){\varpi }^{\pm }(n-1,t)\hfill \\ \hfill & \ge 2\beta max\left\{\right|{\Im }_0\left(n\right)|,|\omega (0,t)|,|\Im (2,t)\}\pm \gamma (n,t)\ge 0.\hfill \end{array}$$

Applying Lemma 2.1 yields ${\varpi }^{\pm }(n,t)\ge 0$, $(n,t)\in \overline{\Xi }$, which implies the required estimate.    □

Lemma 2.3. 

Let k and l are non-negative integers satisfying $0\le k+2l\le 4$. Then, the derivatives of $\Im (n,t)$ is bounded as

$$\left|\frac{{\partial }^{k+l}\Im (n,t)}{\partial n^k\partial t^l}\right|\le C{\epsilon }^{-\frac{k}{2}}.$$

Proof. 

Refer the proof of Lemma 4 of [14].    □

3. Numerical formulation and computational procedures

In this section, the fitted difference scheme of the continuous problem is derived and its convergence analysis is shown on a uniform meshes.

Step 1: Consider a uniform grids ${\overline{\Xi }}_t^M=\{t_j=j\tau ,\tau =\frac{T}{M},j=0\left(1\right)M\}$ on $[0,T]$. Applying implicit Euler technique, (2.1) becomes

$$\frac{{\Im }^{j+1}\left(n\right)-{\Im }^j\left(n\right)}{\tau }-\epsilon {\Im }_{nn}^{j+1}\left(n\right)+\mu \left(n\right){\Im }^{j+1}\left(n\right)+\nu \left(n\right){\Im }^{j+1}(n-1)={\gamma }^{j+1}\left(n\right).$$

On rearrangement and simplification, we obtain

$$\begin{array}{cc}\hfill L_{\epsilon }^M{\Im }^{j+1}\left(n\right)=& -{\zeta }_{\epsilon }{\Im }_{nn}^{j+1}\left(n\right)+p\left(n\right){\Im }^{j+1}\left(n\right)\hfill \\ \hfill =& \left\{\begin{array}{c}-q\left(n\right){\omega }^{j+1}(n-1)+{\varphi }^{j+1}\left(n\right)+{\Im }^j\left(n\right),n\in (0,1],\hfill \\ -q\left(n\right){\Im }^{j+1}(n-1)+{\varphi }^{j+1}\left(n\right)+{\Im }^j\left(n\right),n\in (1,2),\hfill \end{array}\right.\hfill \end{array}$$

(3.1)

subject to the considered conditions, where ${\zeta }_{\epsilon }=\epsilon \tau$, $p\left(n\right)=1+\tau \mu \left(n\right)$, $q\left(n\right)=\tau \nu \left(n\right)$ and ${\varphi }^{j+1}\left(n\right)=\tau {\gamma }^{j+1}\left(n\right)$. The stability estimate of the semi-discrete scheme (3.1) can be easily established following the procedures in Lemma 2.3 on a uniform time mesh, and its error estimate is determined as follow.

Lemma 3.1. 

Assume that $\left|\frac{{\partial }^l\Im (n,t)}{\partial t^l}\right|\le C$, $(n,t)\in [0,2]\times [0,T]$, $l=0,1,2$. The local truncation is given by $\parallel E_L^{j+1}\parallel \le C{\tau }^2$, and the global truncation error is given by $\parallel E_G^{j+1}\parallel \le C\tau$.

Proof. 

The estimate of the local truncation error is referred from Lemma 7 of [14]. Using the estimates of the point-wise error up to ${(j+1)}^{th}$ level, we have

$$\begin{array}{c}\hfill \left|\right|E_G^{j+1}\left|\right|=\parallel \sum _{\varsigma =1}^j{E}_L^{\varsigma }\parallel \le \parallel e^1\parallel +...+\parallel e^j\parallel \le C\tau ,j\tau \le T.\end{array}$$

Hence, the scheme (3.1) is first order convergent in temporal direction.    □

Step 2: Set $\epsilon =0$ to obtain the reduced solution with the appropriate interval conditions. Let ${\Im }_0^{j+1}\left(n\right)$ be the solution of the reduced problem given as

$$p\left(n\right){\Im }_0^{j+1}\left(n\right)=\left\{\begin{array}{c}-q\left(n\right){\omega }_0^{j+1}(n-1)+{\varphi }^{j+1}\left(n\right)+{\Im }_0^j\left(n\right),n\in (0,1],\hfill \\ -q\left(n\right){\Im }_0^{j+1}(n-1)+{\varphi }^{j+1}\left(n\right)+{\Im }_0^j\left(n\right),n\in (1,2)\hfill \end{array}\right.$$

(3.2)

with ${\Im }_0^{j+1}\left(n\right)=\omega \left(n\right)$, $n\in [-1,0]$. Using (3.2), we approximate ${\Im }^{j+1}\left(1\right)$.

Step 3: Setup a uniform spatial mesh as ${\overline{\Xi }}_n^N=\{n_i=ih$, $n_0=0,n_N=1,n_{2N}=2,i=0,1,...,n_{2N}\}$, where $h=\frac{1}{N}$. Then, define an exponential tension function as

$$E_i\left(n\right)=a_i{e}^{\ell (n-n_i)}+b_i{e}^{-\ell (n-n_i)}+c_i(n-n_i)+d_i,$$

(3.3)

where $a_i$, $b_i$, $c_i$ and $d_i$ are constants and ℓ is a tension parameter, which handles the continuity of the function and enhances the accuracy of the method. After certain manipulations, we obtain

$$\begin{array}{cc}\hfill & a_i=\frac{M_{i+1}-M_i{e}^{-\ell h}}{2{\ell }^{2}sinh(\ell h)},b_i=\frac{e^{M_i\ell h}-M_i+1}{2{\ell }^{2}sinh(\ell h)},\hfill \\ \hfill & c_i={\Im }_i^{j+1}-\frac{M_i}{{\ell }^2}\mathrm{and}{d}_i=\frac{{\Im }_{i+1}^{j+1}-{\Im }_i^{j+1}}{h}-\frac{M_{i+1}-M_i}{{\ell }^{2}h}.\hfill \end{array}$$

From the continuity conditions of the first derivative of $E_i\left(n\right)$ at $n=n_i$, we have

$$\frac{1}{h^2}[{\Im }_{i+1}^{j+1}-2{\Im }_i^{j+1}+{\Im }_{i-1}^{j+1}]=\alpha M_{i+1}+2\beta M_i+\alpha M_{i-1},i=0,1,...,2N,$$

(3.4)

where $\alpha =\frac{1}{{\left(\tau h\right)}^2}-\frac{1}{\tau hsinh\left(\tau h\right)}$ and $\beta =\frac{coth\left(\tau h\right)}{\tau h}-\frac{1}{{\left(\tau h\right)}^2}$.

From (3.1), we have

$$\begin{array}{cc}\hfill & {\zeta }_{\epsilon }{M}_i=\left\{\begin{array}{c}{p}_i{\Im }_i^{j+1}+q_i{\omega }_{i-N}^{j+1}+{\varphi }_i^{j+1}+{\Im }_i^j,i=0,1,...,N,\hfill \\ p_i{\Im }_i^{j+1}+q_i{\Im }_{i-N}^{j+1}+{\varphi }_i^{j+1}+{\Im }_i^j,i=N+1,...,2N,\hfill \end{array}\right.\hfill \end{array}$$

(3.5)

$$\begin{array}{cc}\hfill & {\zeta }_{\epsilon }{M}_{i\pm 1}=\left\{\begin{array}{c}{p}_{i\pm 1}{\Im }_{i\pm 1}^{j+1}+q_{i\pm 1}{\omega }_{i\pm 1-N}^{j+1}+{\varphi }_{i\pm 1}^{j+1}+{\Im }_{i\pm 1}^j,i=0,1,...,N,\hfill \\ p_{i\pm 1}{\Im }_{i\pm 1}^{j+1}+q_{i\pm 1}{\Im }_{i\pm 1-N}^{j+1}+{\varphi }_{i\pm 1}^{j+1}+{\Im }_{i\pm 1}^j,i=N+1,...,2N,\hfill \end{array}\right.\hfill \end{array}$$

(3.6)

Substituting (3.5) and (3.6) into (3.4), introducing a fitting factor $\sigma$ and taking the limit on both sides as $h\to 0$ yields

$$\begin{array}{cc}\hfill & {\sigma }_0=\frac{{\rho }^{2}p\left(0\right)(\alpha +\beta )}{4\Delta t{sinh}^2\left(\rho \sqrt{\frac{p\left(0\right)}{4\Delta t}}\right)}on[0,1]\hfill \\ \hfill & and\hfill \end{array}$$

(3.7a)

$$\begin{array}{cc}& {\sigma }_2=\frac{{\rho }^{2}p\left(2\right)(\alpha +\beta )}{4\Delta t{sinh}^2\left(\rho \sqrt{\frac{p\left(2\right)}{4\Delta t}}\right)}on(1,2].\hfill \end{array}$$

(3.7b)

From (3.4)-(3.7), we obtain

$$\begin{array}{cc}\hfill L_{\epsilon }^{N,M}{\Im }_i^{j+1}=& \frac{-\epsilon \Delta t\sigma }{h^2}[{\Im }_i^{j+1}-2{\Im }_i^{j+1}+{\Im }_{i-1}^{j+1}]+p_{i+1}{\Im }_{i+1}^{j+1}+p_i{\Im }_i^{j+1}+p_{i-1}{\Im }_{i-1}^{j+1}\hfill \\ \hfill =& \left\{\begin{array}{c}-q_{i+1}{\omega }_{i+1-N}-q_i{\omega }_{i-N}-q_{i-1}{\omega }_{i-1-N}+{\varphi }_{i+1}^{j+1}+{\varphi }_i^{j+1}+{\varphi }_{i-1}^{j+1}\hfill \\ +{\Im }_{i+1}^j+{\Im }_i^j+{\Im }_{i-1}^j,i=1,2,...,N\hfill \\ -q_{i+1}{\Im }_{i+1-N}^{j+1}-q_i{\Im }_{i-N}^{j+1}-q_{i-1}{\Im }_{i-1-N}^{j+1}+{\varphi }_{i+1}^{j+1}+{\varphi }_i^{j+1}+{\varphi }_{i-1}^{j+1}\hfill \\ +{\Im }_{i+1}^j+{\Im }_i^j+{\Im }_{i-1}^j,i=N+1,N+2,...,2N.\hfill \end{array}\right.\hfill \end{array}$$

(3.8)

From (3.8), we obtain system of equation as

$${\wp }_1^{+}{\Im }_{i+1}^{j+1}+{\wp }_1^0{\Im }_i^{j+1}+{\wp }_1^{-}{\Im }_{i-1}^{j+1}=D_1{\Im }_{i+1}^j+E_1{\Im }_i^j+D_1{\Im }_{i-1}^j+F_i,i=1,2,...,N-1,$$

(3.9)

with the boundary conditions ${\Im }_0^{j+1}=\omega \left(0\right)$ and ${\Im }_N^{j+1}={\Im }^{j+1}\left(1\right)$, where

$$\begin{array}{cc}\hfill & {\wp }_1^{+}=\frac{\epsilon \Delta t{\sigma }_0}{h^2}-\alpha p_{i+1},{\wp }_1^0=\frac{-2\epsilon \Delta t{\sigma }_0}{h^2}-2\beta p_i,{\wp }_1^{-}=\frac{\epsilon \Delta t{\sigma }_0}{h^2}-\alpha p_{i-1},\hfill \\ \hfill & D_1=-\alpha ,E_1=-2\beta ,\hfill \\ \hfill & F_i=\alpha q_{i+1}{\omega }_{i+1-N}+2\beta q_i{\omega }_{i-N}+\alpha q_{i-1}{\omega }_{i-1-N}-\alpha {\varphi }_{i+1}^{j+1}-2\beta {\varphi }_i^{j+1}-\alpha {\varphi }_{i-1}^{j+1}.\hfill \end{array}$$

Similarly, we have

$${\wp }_2^{+}{\Im }_{i+1}^{j+1}+{\wp }_2^0{\Im }_i^{j+1}+{\wp }_2^{-}{\Im }_{i-1}^{j+1}=D_2{\Im }_{i+1}^j+E_2{\Im }_i^j+D_2{\Im }_{i-1}^j+G_i,i=N+1,,...,2N,$$

(3.10)

with the boundary conditions ${\Im }_N^{j+1}={\Im }^{j+1}\left(1\right)$ and ${\Im }_{2N}^{j+1}=\lambda (2,t)$, where

$$\begin{array}{cc}\hfill & {\wp }_2^{+}=\frac{\epsilon \Delta t{\sigma }_2}{h^2}-\alpha p_{i+1},{\wp }_1^0=\frac{-2\epsilon \Delta t{\sigma }_2}{h^2}-2\beta p_i,{\wp }_1^{-}=\frac{\epsilon \Delta t{\sigma }_2}{h^2}-\alpha p_{i-1},\hfill \\ \hfill & D_1=-\alpha ,E_1=-2\beta ,\hfill \\ \hfill & G_i=\alpha q_{i+1}{\Im }_{i+1-N}^{j+1}+2\beta q_i\Im *{j+1}_{i-N}+\alpha q_{i-1}{\Im }_{i-1-N}^{j+1}-\alpha {\varphi }_{i+1}^{j+1}-2\beta {\varphi }_i^{j+1}-\alpha {\varphi }_{i-1}^{j+1}.\hfill \end{array}$$

From (3.9) and (3.10), we have system of equations which can be solved using suitable solvers following the computational algorithm given below.

Algorithm 1:

Discrete stability and convergence analyses

Lemma 3.2. 

Let ${\phi }_i^{j+1}$ is a mesh function which satisfies ${\phi }_0^{j+1}\ge 0$, ${\phi }_{2N}^{j+1}\ge 0$ and $L_{\epsilon }^{N.M}{\phi }_i^{j+1}\ge 0$, $i=1\left(1\right)2N-1$. Then, ${\phi }_i^{j+1}\ge 0$, $i=0\left(1\right)2N$.

Proof. 

We refer Lemma 3.6 of [15]. □

Lemma 3.3. 

The stability of the solution ${\Im }_i^{j+1}$, $i=0,1,...,2N$ of (3.8) is estimated as

$$|{\Im }_i^{j+1}|\le {\beta }^{-1}\parallel L_{\epsilon }^{N,M}\Im \parallel +max\left\{\right|{\omega }^{j+1}\left(0\right)|,|{\lambda }^{j+1}\left(2\right)\left|\right\}.$$

Proof. 

Consider barrier functions ${\varpi }_{i,\pm }^{j+1}$ given by

$${\varpi }_{i,\pm }^{j+1}=\frac{1}{\beta }\parallel L_{\epsilon }^{N,M}\Im \parallel +sup\left\{\right|{\omega }^{j+1}\left(0\right)|,|{\lambda }^{j+1}\left(2\right)\left|\right\}\pm {\Im }_i^{j+1}.$$

For $i=0$, we have ${\varpi }_{0,\pm }^{j+1}\ge 0$ and for $i=2N$, we have ${\varpi }_{2N,\pm }^{j+1}\ge 0$.

For $i=1\left(1\right)N$, by (3.8) we obtain

$$\begin{array}{cc}\hfill L_{\epsilon ,1}^{j+1}{\varpi }_{i,\pm }^{j+1}=& \frac{-\epsilon \Delta t{\sigma }_0}{h^2}[{\varpi }_{i,\pm }^{j+1}-2{\varpi }_{i,\pm }^{j+1}+{\varpi }_{i,\pm }^{j+1}]+p_{i+1}{\varpi }_{i+1,\pm }^{j+1}+p_i{\varpi }_{i,\pm }^{j+1}+p_i{\varpi }_{i-1,\pm }^{j+1}\hfill \\ \hfill & +q_{i+1}{\omega }_{i+1-N}+q_i{\omega }_{i-N}+q_{i-1}{\omega }_{i-1-N}\hfill \\ \hfill =& \pm L_{\epsilon ,1}^{j+1}{\Im }_i^{j+1}+p_{i+1}[{\beta }^{-1}\parallel L_{\epsilon ,1}^{N,M}{\Im }_{i+1}^{j+1}\parallel +sup\left\{\right|{\omega }^{j+1}\left(0\right)|,|{\lambda }^{j+1}\left(2\right)\left|\right\}]\hfill \\ \hfill & +p_i[{\beta }^{-1}\parallel L_{\epsilon ,1}^{N,M}{\Im }_i^{j+1}\parallel +sup\left\{\right|{\omega }^{j+1}\left(0\right)|,|{\lambda }^{j+1}\left(2\right)\left|\right\}]\hfill \\ \hfill & +p_{i-1}[{\beta }^{-1}\parallel L_{\epsilon ,1}^{N,M}{\Im }_{i-1}^{j+1}\parallel +sup\left\{\right|{\omega }^{j+1}\left(0\right)|,|{\lambda }^{j+1}\left(2\right)\left|\right\}]\hfill \\ \hfill \ge & (P_{i+1}+p_i+p_{i-1})sup\left\{\right|{\omega }^{j+1}\left(0\right)|,|{\lambda }^{j+1}\left(2\right)\left|\right\}\ge 0.\hfill \end{array}$$

Similarly, for $i=1\left(1\right)N$, one can obtain

$$\begin{array}{cc}\hfill L_{\epsilon ,2}^{j+1}{\varpi }_{i,\pm }^{j+1}=& \frac{-\epsilon \Delta t{\sigma }_2}{h^2}[{\varpi }_{i,\pm }^{j+1}-2{\varpi }_{i,\pm }^{j+1}+{\varpi }_{i,\pm }^{j+1}]+p_{i+1}{\varpi }_{i+1,\pm }^{j+1}+p_i{\varpi }_{i,\pm }^{j+1}+p_i{\varpi }_{i-1,\pm }^{j+1}\hfill \\ \hfill & +q_{i+1}{\Im }_{i+1-N}^{j+1}+q_i{\Im }_{i-N}^{j+1}+q_{i-1}{\Im }_{i-1-N}^{j+1}\hfill \\ \hfill =& \pm L_{\epsilon ,1}^{j+1}{\Im }_i^{j+1}+p_{i+1}[{\beta }^{-1}\parallel L_{\epsilon ,1}^{N,M}{\Im }_{i+1}^{j+1}\parallel +sup\left\{\right|{\omega }^{j+1}\left(0\right)|,|{\lambda }^{j+1}\left(2\right)\left|\right\}]\hfill \\ \hfill & +p_i[{\beta }^{-1}\parallel L_{\epsilon ,1}^{N,M}{\Im }_i^{j+1}\parallel +sup\left\{\right|{\omega }^{j+1}\left(0\right)|,|{\lambda }^{j+1}\left(2\right)\left|\right\}]\hfill \\ \hfill & +p_{i-1}[{\beta }^{-1}\parallel L_{\epsilon ,1}^{N,M}{\Im }_{i-1}^{j+1}\parallel +sup\left\{\right|{\omega }^{j+1}\left(0\right)|,|{\lambda }^{j+1}\left(2\right)\left|\right\}]\hfill \\ \hfill & +q_{i+1}[{\beta }^{-1}\parallel L_{\epsilon ,1}^{N,M}{\Im }_{i+1-N}^{j+1}\parallel +sup\left\{\right|{\omega }^{j+1}\left(0\right)|,|{\lambda }^{j+1}\left(2\right)\left|\right\}]\hfill \\ \hfill & +q_i[{\beta }^{-1}\parallel L_{\epsilon ,1}^{N,M}{\Im }_{i-N}^{j+1}\parallel +sup\left\{\right|{\omega }^{j+1}\left(0\right)|,|{\lambda }^{j+1}\left(2\right)\left|\right\}]\hfill \\ \hfill & +q_{i-1}[{\beta }^{-1}\parallel L_{\epsilon ,1}^{N,M}{\Im }_{i-1-N}^{j+1}\parallel +sup\left\{\right|{\omega }^{j+1}\left(0\right)|,|{\lambda }^{j+1}\left(2\right)\left|\right\}]\hfill \\ \hfill \ge & \beta sup\left\{\right|{\omega }^{j+1}\left(0\right)|,|{\lambda }^{j+1}\left(2\right)\left|\right\}\ge 0.\hfill \end{array}$$

Applying Lemma 3.1, the stability estimate holds true. □

Theorem 3.1. 

Let ${\Im }^{j+1}\left(n_i\right)$ and ${\Im }_i^{j+1}$ are the solutions of (3.1) and (3.8), respectively. Then, the point-wise error in the spatial direction is bounded as

$$|{\Im }^{j+1}\left(n_i\right)-{\Im }_i^{j+1}|\le C{N}^{-2},i=1,2,...,2N.$$

for constant C independent of ε and N.

Proof. 

By local truncation error of finite difference methods [16], we have

$$\begin{array}{cc}\hfill |L_{\epsilon }^{N,M}({\Im }^{j+1}\left(n_i\right)-{\Im }_i^{j+1})|=& |L_{\epsilon }^M{\Im }^{j+1}\left(n_i\right)-L_{\epsilon }^{N,M}{\Im }_i^{j+1}|\hfill \\ \hfill =& |-c_{\epsilon }{\Im }_{nn}^{j+1}\left(n_i\right)+\frac{c_{\epsilon }\sigma }{h^2}[{\Im }_{i+1}^{j+1}-2{\Im }_i^{j+1}+{\Im }_{i-1}^{j+1}]|.\hfill \end{array}$$

(3.11)

Taylor’s series expansion on ${\Im }_{i-1}^{j+1}$ and ${\Im }_{i+1}^{j+1}$ gives

$$\begin{array}{cc}\hfill & {\Im }_{i+1}^{j+1}={\Im }_i^{j+1}+h{\Im }_n^{j+1}+\frac{h^2}{2!}{\Im }_{nn}^{j+1}+\frac{h^3}{3!}{\Im }_{nnn}^{j+1}+O\left(h^4\right),\hfill \end{array}$$

(3.12)

$$\begin{array}{cc}\hfill and& {\Im }_{i-1}^{j+1}={\Im }_i^{j+1}-h{\Im }_n^{j+1}+\frac{h^2}{2!}{\Im }_{nn}^{j+1}-\frac{h^3}{3!}{\Im }_{nnn}^{j+1}+O\left(h^4\right).\hfill \end{array}$$

(3.13)

Plugging (3.12) and (3.13) into (3.11), we obtain

$$\begin{array}{cc}\hfill |L_{\epsilon }^{N,M}({\Im }^{j+1}\left(n_i\right)-{\Im }_i^{j+1})|=& |-\epsilon \Delta t{\Im }_{nn}^{j+1}+\epsilon \Delta t\sigma [{\Im }_{nn}^{j=1}+\frac{2h^2}{4!}{\Im }_{nnnn}^{j+1}+O\left(h^4\right)]\hfill \\ \hfill \le & C|\epsilon \sigma -\epsilon |\parallel {\Im }_{nn}^{j+1}\left({\xi }_i\right)\parallel +C\epsilon \sigma h^2\parallel {\Im }_{nnnn}^{j+1}\left({\xi }_i\right)\parallel ,{\xi }_i\in [n_{i-1},n_{i+1}].\hfill \end{array}$$

As estimated in [17], we have $|\epsilon \sigma -\epsilon |\le C{h}^2$ and hence for sufficiently small h, we get

$$|L_{\epsilon }^{N,M}({\Im }^{j+1}\left(n_i\right)-{\Im }_i^{j+1})|\le C{N}^{-2}.$$

Hence, applying Lemma 3.1 gives the required estimate. □

Theorem 3.2. 

Let $\Im (n,t)$ and ${\Im }_i^{J+1}$ are the solutions of (2.1) and 3.8, respectively. Then, the uniform error estimation is given by

$$\parallel \Im (n_i,t_{j+1})-{\Im }_i^{j+1}\parallel \le C(\tau +N^{-2}).$$

Proof. 

By the triangular inequality, we have

$$\parallel \Im (n_i,t_{j+1})-{\Im }_i^{j+1}\parallel \le \parallel \Im (n_i,t_{j+1})-{\Im }_i^{j+1}\parallel +\parallel {\Im }^{j+1}\left(n_i\right)-{\Im }_i^{j+1}\parallel$$

(3.14)

Using Lemma 3.1 and Theorem 3.1 into (3.14), yields required the uniform error estimate. □

4. Numerical experiments and discussions

In this section, to validate the developed numerical numerical techniques, we solved model examples. As it is tedious to compute the exact solutions of the examples, we apply a variant of the double mesh principle [16,18] as given in Algorithm 1.

Example 4.1. 

[9,10]. Consider a SP-PDE given by ${\Im }_t(n,t)-\epsilon {\Im }_{nn}(n,t)+3\Im (n,t)-\Im (n-1,t)=1$, $\Im (n,0)=0$, $n\in [0,2]$, $\Im (n,t)=0$, $(n,t)\in \left\{\right(n,t):-1\le n\le 0,0\le t\le 2\}$ and $\Im (2,t)=0$, $(2,t)\in \left\{\right(2,t):t\in [0,2\left]\right\}$.

Example 4.2. 

[10]. Consider a SP-PDE given by ${\Im }_t(n,t)-\epsilon {\Im }_{nn}(n,t)+(n+6)\Im (n,t)-(n^2+1)\Im (n-1,t)=3$, $\Im (n,0)=0$, $n\in [0,2]$, $\Im (n,t)=0$, $(n,t)\in \left\{\right(n,t):-1\le n\le 0,0\le t\le 2\}$ and $\Im (2,t)=0$, $(2,t)\in \left\{\right(2,t):t\in [0,2\left]\right\}$ .

Example 4.3. 

[9]. Consider a SP-PDE given by ${\Im }_t(n,t)-\epsilon {\Im }_{nn}(n,t)+5\Im (n,t)-2\Im (n-1,t)=-2$, $\Im (n,0)=sin\left(\pi n\right)$, $n\in [0,2]$, $\Im (n,t)=0$, $(n,t)\in \left\{\right(n,t):-1\le n\le 0,0\le t\le 2\}$ and $\Im (2,t)=0$, $(2,t)\in \left\{\right(2,t):t\in [0,2\left]\right\}$.

The examples are treated applying the developed numerical algorithm on a Matlab software package. We computed the maximum absolute errors (MAEs) and convergences rates (CRs) to validate the method and the numerical solutions are simulated graphically for each example. In Table 1, Table 2 and Table 3, the MAEs and CRs are shown for each examples respectively. From the tables, one can see that reducing the number of meshes reduces the error and reducing $\epsilon$ yields stabled Maximum point-wise errors, which indicate that the technique is uniformly convergent. Table 4 shows comparisons of the present result with other published results.

Simulations of the numerical solutions and error profiles are displayed graphically for each examples. In the line plots given in Figure 1, Figure 4 and Figure 7, the effect of $\epsilon$ on the solution profile are shown for Examples 4.1, 4.2 and 4.3, respectively. The influence of $\epsilon$ is also simulated in the surface plots given in Figure 2, Figure 5 and Figure 8 for each examples. From such figures, one can visualize that reducing the value of $\epsilon$ reduces the width of the layer regions yielding more steep gradients. From the maximum error profiles given in Figure 3, Figure 6 and Figure 9, one can observe that maximum error is committed in the interior layers raised due to the large negative shift. Figure 10 is log-log plots of the MAEs versus N for the examples, confirming the convergence of the scheme uniformly.

5. Conclusion

In this article, a numerical method is formulated and its computational algorithm is provided to solve singularly perturbed parabolic partial differential equations involving a negative shift. The influence of $\epsilon$ causes boundary layers and the negative shift term causes interior layer in the solutions. We treated the continuous problem by transformed it into a finite difference method applying implicit Euler technique in time variable and using an exponentially fitted tension spline function on a uniform grids in space variable. The stability and convergence analyses of the formulated method are investigated and proved. The step-wise computational algorithm is also provided, and to test the formulated technique, numerical experiments are carried out. For the considered examples, the numerical solutions and the effect of $\epsilon$ and the large negative shift on the solution profile are simulated graphically. MAEs and CRs of the numerical solutions are computed for each example. From the uniform convergence analysis and numerical experiments, we conclude that the formulated technique convergent uniformly with order one in time and second order in the spatial direction.