Skip to main content

Fitted computational method for solving singularly perturbed small time lag problem



An accurate exponentially fitted numerical method is developed to solve the singularly perturbed time lag problem. The solution to the problem exhibits a boundary layer as the perturbation parameter approaches zero. A priori bounds and properties of the continuous solution are discussed.


The backward-Euler method is applied in the time direction and the higher order finite difference method is employed for the spatial derivative approximation. An exponential fitting factor is induced on the difference scheme for stabilizing the computed solution. Using the comparison principle, the stability of the method is examined and analyzed. It is proved that the method converges uniformly with linear order of convergence. To validate the theoretical findings and analysis, two test examples are given. Comparison is made with the results available in the literature. The proposed method has better accuracy than the schemes in the literature.


A differential equation in which its highest order derivative term is multiplied by a small parameter is known as a singularly perturbed differential equation (SPDE). It commonly occurs in the modeling of chemical processes, fluid flows, water quality problems in river networks, mechanical systems and simulation of oil extraction from underground reservoirs [1]. The solution of such type of equation possesses a multi-scale character that varies quickly in the boundary layer region and it slows in the outer layer regions.

Due to the multi-scale character of the solution, the classical numerical methods fail to give an accurate result. Currently, it becomes interesting to develop a numerical method which gives accurate results; and its convergence does not depend on the perturbation parameter. For solving the considered problem in Kumar et al. [2] proposed an adaptive mesh method using the concept of entropy function. Hybrid scheme of the midpoint upwind in the outer region and the central difference method in the boundary layer region are used in [3, 4]. Gowrisankar and Natesan [5] used the upwind finite difference method on a piecewise uniform mesh. The upwind finite difference method on Shishkin mesh is used in [6]. Podila and Kumar [7] used a stable finite difference scheme, which works on a uniform and an adaptive mesh. An exponentially fitted scheme is discussed in [8, 9]. The non-standard finite difference method is used in [10]. The numerical schemes developed in [11,12,13,14,15,16] works for both the large and small delay cases.

In this paper, we proposed a numerical scheme using higher order finite difference scheme fitted by the exponential fitting factor. Moreover, the main aim of this study is to develop a more accurate, stable and uniformly convergent numerical scheme for solving singularly perturbed convection-diffusion problem having small time lag.

In this paper, C is a generic positive constant, which does not depend on the mesh parameters and the perturbation parameter. The norm \(\Vert .\Vert\), defined by \(\Vert g \Vert =max_{(s,t) \in \Omega }\vert g(s,t) |\) is the maximum norm.

Continuous problem

Consider a class of singularly perturbed convection-diffusion problem of the form

$$\begin{aligned} \left\{ \begin{array}{lll} z_{t}(s,t) + {\mathcal {L}}z(s,t) = -c(s,t)z(s, t- \tau ) + g(s,t) , (s,t) \in \Omega =(0,1)\times (0,T], \\ z(s,t) = \psi _{b}(s,t) \ on \ \eta _{b}:= [0,1] \times [-\tau , 0], \\ z(0,t)= \psi _{l}(t) \ on \ \eta _{l} := \{ (0,t):0\le t\le T \}, \\ z(1,t) = \psi _{r}(t) \ on \ \eta _{r} := \{(1,t):0\le t \le T\} , \end{array}\right. \end{aligned}$$

where \({\mathcal {L}}z(s,t)=- \varepsilon z_{ss}(s,t) + b(s,t)z_{s}(s,t) + a(s)z(s,t)\), \(0 < \varepsilon \ll 1\) is the perturbation parameter and \(\tau > 0\) is the delay parameter. The coefficients \(b(s,t),\ a(s), \ c(s,t), \ g(s,t)\) on \(\bar{\Omega }\) and \(\psi _{b}(s,t),\ \psi _{l}(t),\ \psi _{r}(t)\) on \(\eta = \eta _{l}\cup \eta _{r} \cup \eta _{b}\) are assumed to be smooth and bounded functions that satisfy the conditions \(a(s) \ge \alpha>0,\ c(s,t) \ge \gamma >0 , \ b(s,t) \le \beta <0 ,\ a(s)+ c(s,t)\ge 0\) on \(\bar{\Omega }\). Under these conditions, the problem exhibits a boundary layer along \(s=0\) [2].

Properties of the analytical solution

In this part, we present the analytical aspects of the solution and its derivatives. The existence and uniqueness of a solution of (1).

Lemma 2.1

[17, 18] The solution z(st) of (1) satisfies \(\vert z(s,t) - \psi _{b}(s,0)|\le Ct,\; (s,t)\in \bar{\Omega }=[0,1]\times [0,T],\) where \(C >0\) is a constant that does not depends on \(\varepsilon\).

The operator \(L(s,t)= z_{t}(s,t) + {\mathcal {L}}z(s,t)\) of (1) satisfies the following minimum principle.

Lemma 2.2

[19] Let \(\nu (s,t) \in C^{2}(\Omega ) \cap C^{0}(\bar{\Omega })\), satisfies \(\nu (s,t) \ge 0\)  \((s,t) \in \partial \Omega=\bar{\Omega}-\Omega\). If \(L\nu (s,t) \le 0, ~(s,t) \in \Omega\), then \(\nu (s,t) \ge 0, (s,t)\in \bar{\Omega }\).

Lemma 2.3

[18, 19] [Stability result] The solution of (1) satisfies the bound

$$\begin{aligned} \vert z(s,t)|\le \alpha ^{-1}\Vert g\Vert + \max \{\vert \psi _{l}(t)|, \vert \psi _{b}(s,t)\vert , \vert \phi _{r}(t)|\}, \end{aligned}$$

where \(\alpha\) is the lower bound of a(s) and \(\Vert g \Vert = \max _{s,t\in \Omega }\vert g(s,t)|\).

Lemma 2.4

[20, 21] The following bounds are satisfied for the derivative of the solution z(st) of (1) with respect to s and t

$$\begin{aligned} \bigg \vert \frac{\partial ^{n}z(s,t)}{\partial s^{n}}\bigg |\le C\left(1 + \varepsilon ^{-n}\exp \left(-\frac{\beta }{\varepsilon }s\right)\right), \ (s,t)\in \bar{\Omega }, \ n= 0,1,2,3,4, \end{aligned}$$


$$\begin{aligned} \bigg \vert \frac{\partial ^{l}z(s,t)}{\partial t^{l}}\bigg |\le C, \ (s,t) \in \bar{\Omega } , \ l = 0,1,2, \end{aligned}$$

where \(\beta\) is lower bound of b(st).

Main text: numerical scheme

Temporal semi-discretization

Using Taylor’s series expansion for the delay term; using a uniform mesh in t-direction with step size \(\Delta t = T/M\) given by \(\Omega _{t}^{M} = \{ t_{n} = n\Delta t, n = 0,1,2,...,M, t_{M} = T\}\), where M is the number of mesh points in [0, T]. Note that \(T = r\tau\) for some positive integer r. Using the backward-Euler formula, we get

$$\begin{aligned} \begin{aligned} (1 - \tau c(s,t_{n}))\frac{z^{n}(s) - z^{n-1}(s)}{\Delta t}&-\varepsilon \frac{d^{2}z^{n}(s)}{ds^{2}} + b(s,t_{n})\frac{d z^{n}(s)}{ds} \\ {}&+(a(s)+c(s,t_{n}))z^{n}(s) = g(s,t_{n}). \end{aligned} \end{aligned}$$

Equivalently, we write

$$\begin{aligned} -\varepsilon \frac{\ d^{2}Z^{n}(s)}{\ ds^{2}} + b(s,t_{n})\frac{\ d Z^{n}(s)}{\ ds} + P(s,t_{n})Z^{n}(s) = R(s,t_{n}), \end{aligned}$$

where \(P(s,t_{n}) = a(s)+c(s,t_{n})+\frac{1 - \tau c(s,t_{n})}{\Delta t},\) and \(R(s,t_{n}) = g(s,t_{n})+ \frac{1 - \tau c(s,t_{n})}{\Delta t}Z^{n-1}(s)\) with the boundary

$$\begin{aligned} Z^{n}(0) = \psi _{l}(t_{n}),~~~~~~ Z^{n}(1) = \psi _{r}(t_{n}),~~~0\le n\le M. \end{aligned}$$

Now, we rewrite (6) as

$$\begin{aligned} -\varepsilon \frac{\ d^{2}Z(s)}{\ ds^{2}} + b(s,t_{n})\frac{\ dZ(s)}{\ d s} + P(s,t_{n})Z(s) = R(s,t_{n}), \end{aligned}$$

where \(Z(s) = Z^{n}(s) = Z(s,t_{n})\). At each time step the local error is defined as \(e_{n}(s): = z(s,t_{n}) - Z^{n}(s), \ 0\le n \le M.\)

Lemma 3.1

The local error estimate in the temporal direction satisfies the bound

$$\begin{aligned} \Vert e_{n}\Vert \le C_{1}(\Delta t)^{2}, \end{aligned}$$

and the global error at \(n^{th}\) time level satisfies the bound

$$\begin{aligned} \Vert E_{n}\Vert \le C(\Delta t). \end{aligned}$$


Refer from the Appendix section. \(\square\)

Lemma 3.2

For \(0 \le n \le M-1\), the solution of (6)–(7), satisfies the bound

$$\begin{aligned} \bigg |\frac{d^{r}Z^{n}(s)}{ds^{r}} \bigg |\le C\left(1 + \varepsilon ^{-r}\exp \left(-\frac{\beta }{\varepsilon }(s)\right)\right),~~~s\in \bar{\Omega }_{s},~~0\le r\le 4. \end{aligned}$$


For the proof see [22]. \(\square\)

Spatial discretization

The spatial domain [0, 1] into N equal number of sub-intervals with the length of h is given by \(0= s_{0},s_{1},...,s_{N}=1\) and \(s_{i} = ih, i = 0,1,2,...,N\). Assume smooth function \(Z(s) = Z^{n}(s) = Z(s,t_{n})\) in the interval [0, 1]. From Taylor’s series expansion, we have

$$\begin{aligned} \begin{aligned} Z_{i+1} \approx\,&Z_{i}+hZ^{'}_{i}+\frac{h^{2}}{2!}Z^{''}_{i}+\frac{h^{3}}{3!}Z^{'''}_{i}+\frac{h^{4}}{4!}Z^{(4)}_{i} +\frac{h^{5}}{5!}Z^{(5)}_{i}+\frac{h^{6}}{6!}Z^{(6)}_{i}+\frac{h^{7}}{7!}Z^{(7)}_{i} +\frac{h^{8}}{8!}Z^{(8)}_{i},\\ Z_{i-1} \approx\,&Z_{i}-hZ^{'}_{i}+\frac{h^{2}}{2!}Z^{''}_{i}-\frac{h^{3}}{3!}Z^{'''}_{i}+\frac{h^{4}}{4!}Z^{(4)}_{i} -\frac{h^{5}}{5!}Z^{(5)}_{i}+\frac{h^{6}}{6!}Z^{(6)}_{i}-\frac{h^{7}}{7!}Z^{(7)}_{i} +\frac{h^{8}}{8!}Z^{(8)}_{i}. \end{aligned} \end{aligned}$$

Combining the result in (11), we arrive at

$$\begin{aligned} \begin{aligned} Z_{i-1} - 2Z_{i} + Z_{i+1} =&\frac{2h^{2}}{2!}Z^{''}_{i}+\frac{2h^{4}}{4!}Z^{(4)}_{i} +\frac{2h^{6}}{6!}Z^{(6)}_{i} +\frac{2h^{8}}{8!}Z^{(8)}_{i}+O(h^{10}) \\ Z^{''}_{i-1} - 2Z^{''}_{i} + Z^{''}_{i+1} =&\frac{2h^{2}}{2!}Z^{(4)}_{i}+\frac{2h^{4}}{4!}Z^{(6)}_{i} +\frac{2h^{6}}{6!}Z^{(8)}_{i} +\frac{2h^{8}}{8!}Z^{(10)}_{i}+O(h^{12}). \end{aligned} \end{aligned}$$

Using (12), we get \(\frac{h^{4}}{12}Z^{(6)}_{i}\) and, we obtain

$$\begin{aligned} \begin{aligned} Z_{i-1} - 2Z_{i} + Z_{i+1} =&\frac{h^{2}}{30}(Z^{''}_{i-1}+28Z^{''}_{i} +Z^{''}_{i+1}) + \zeta , \end{aligned} \end{aligned}$$

where \(\zeta = \frac{h^{4}}{20}Z^{(4)}_{i}- \frac{13h^{8}}{302,400}Z^{(8)}_{i} + O(h^{10})\). Using (8), we get

$$\begin{aligned} \begin{aligned} -\varepsilon Z^{''}_{i+1} =&-b(s_{i+1},t_{n})Z^{'}_{i+1} - P(s_{i+1},t_{n})Z_{i+1}+R(s_{i+1},t_{n}) ,\\ -\varepsilon Z^{''}_{i} =&-b(s_{i},t_{n})Z^{'}_{i} -P(s_{i},t_{n})Z_{i}+R(s_{i},t_{n}), \\ -\varepsilon Z^{''}_{i-1} =&-b(s_{i-1},t_{n})Z^{'}_{i-1} - P(s_{i-1},t_{n})Z_{i-1}+R(s_{i-1},t_{n}). \end{aligned} \end{aligned}$$

Next, using the non-symmetric finite difference schemes we approximate \(Z^{'}_{i+1}, Z^{'}_{i}\) and \(Z^{'}_{i+1}\) as

$$\begin{aligned} \begin{aligned} Z^{'}_{i} =&\frac{Z_{i+1} - Z_{i-1}}{2h}+O(h^2), \;\; Z^{'}_{i+1} = \frac{3Z_{i+1} -4Z_{i}+ Z_{i-1}}{2h} - hZ^{''}_{i}+O(h^{2}), \\ Z^{'}_{i-1} =&\frac{-Z_{i+1} +4Z_{i}- 3Z_{i-1}}{2h} + hZ^{''}_{i}+O(h^{2}). \end{aligned} \end{aligned}$$

Substituting (15) into (14) and substituting the resulting equation into (13) after simplifying, we obtain

$$\begin{aligned}&-{\left(\varepsilon -{\frac{hb(s_{i-1},t_{n})}{30}}+\frac{hb(s_{i+1},t_{n})}{30}\right)}\left({\frac{Z_{i-1} - 2Z_{i} + Z_{i+1}}{h^2}}\right) + {\frac{b(s_{i-1},t_{n})}{60h}}(-3Z_{i-1}+4Z_{i}-Z_{i+1})\\&+{\frac{28b(s_{i},t_{n})}{60h}(Z_{i+1}-Z_{i-1})} +{\frac{b(s_{i+1},t_{n})}{60h}}(Z_{i-1}-4Z_{i}+3Z_{i+1}) + {\frac{P(s_{i-1},t_{n})}{30}}Z_{i-1}\\ {}&+{\frac{28P(s_{i},t_{n})}{30}Z_{i}} +{\frac{P(s_{i+1},t_{n})}{30}Z_{i+1}} = {\frac{1}{30}}(R(s_{i-1},t_{n}) +28R(s_{i},t_{n})+R(s_{i+1},t_{n})). \end{aligned}$$

Computing the exponential fitting factor

In this part, we introduce the fitting factor \(\sigma\) and for the obtained scheme of (6)–(7) at (in)th level. As the theory of singular perturbation given in [9, 23], the zero order asymptotic solution of the problem of the form

$${\left\{ \begin{array}{ll} -\varepsilon Z^{''}(s) +b(s)Z^{'}+p(s)Z(s) =q(s),~~~ s\in \Omega _{s} = (0,1),\\ Z(0) = \alpha _{l},~~~Z(1) = \alpha _{r}, \end{array}\right. }$$

is given by

$$\begin{aligned} Z(s) \approx Z_{0}(s) + \frac{b(0)}{b(s)}(\alpha _{l} - Z_{0}(0))\exp\bigg(-{\int _{0}^{s}\left(\frac{b(s)}{\varepsilon }-\frac{p(s)}{b(s)}\right)ds\bigg) + O(\varepsilon )}. \end{aligned}$$

From Taylor’s series expansion for b(s) and p(s) restricting to their first terms about \(s =0\) and the simplified form gives

$$\begin{aligned} Z(s) = Z_{0}(s)+(\alpha _{l}-Z_{0}(0))\exp \left(-\frac{b(0)}{\varepsilon }s\right), \end{aligned}$$

where \(Z_{0}\) is the reducible problem solution. Considering h is fairly small and solving the result in (19) at \(s_{i}\) which gives

$$\begin{aligned} Z(s_{i}) =Z(ih) = Z_{0}(0)+(\alpha _{l}-Z_{0}(0))\exp (-b(0)i\rho ), \end{aligned}$$

where \(\rho = \frac{h}{\varepsilon }.\) We present an exponentially fitting factor \(\sigma\) of the scheme (16)

$$\begin{aligned}&-{\left(\varepsilon \sigma -{\frac{hb(s_{i-1},t_{n})}{30}}+{\frac{hb(s_{i+1},t_{n})}{30}}\right)}\left({\frac{Z_{i-1} - 2Z_{i} + Z_{i+1}}{h^2}}\right) + {\frac{b(s_{i-1},t_{n})}{60h}}(-3Z_{i-1}+4Z_{i}-Z_{i+1})\\&+{\frac{28b(s_{i},t_{n})}{60h}}(Z_{i+1}-Z_{i-1}) +{\frac{b(s_{i+1},t_{n})}{60h}}(Z_{i-1}-4Z_{i}+3Z_{i+1}) + {\frac{P(s_{i-1},t_{n})}{30}}Z_{i-1}\\ {}&+{\frac{28P(s_{i},t_{n})}{30}}Z_{i}+{\frac{P(s_{i+1},t_{n})}{30}}Z_{i+1} = {\frac{1}{30}}(R(s_{i-1},t_{n}) +28R(s_{i},t_{n})+R(s_{i+1},t_{n})). \end{aligned}$$

Putting \(\rho = \frac{h}{\varepsilon }\) and after multiplying both side by h and letting \(h\rightarrow 0\) which gives

$$\begin{aligned}&-{\frac{\sigma }{\rho }}\lim _{h \rightarrow 0} (Z_{i-1} -2Z_{i}+Z_{i+1}) + {\frac{b_{0}(t_{n})}{60}}\lim _{h \rightarrow 0}(-3Z_{i-1}+4Z_{i}-Z_{i+1}) \\&+{\frac{28b_{0}(t_{n})}{60}}\lim _{h \rightarrow 0}(Z_{i+1}-Z_{i-1}) +{\frac{b_{0}(t_{n})}{60}}\lim _{h \rightarrow 0}(Z_{i-1}-4Z_{i}+3Z_{i+1}) = 0. \end{aligned}$$

Using the results in (19) and simplifying, we obtain

$$\begin{aligned} \frac{\sigma }{\rho }(e^{ b(0)\rho } + e^{-b(0)\rho } - 2) = \frac{b_0(t_{n})}{60}(-30e^{b(0)\rho }+30e^{-b(0)\rho }). \end{aligned}$$

Solving for \(\sigma\) in (22), the exponential fitting factor \(\sigma\) is obtained as

$$\begin{aligned} \sigma = \frac{\rho b_0(t_{n})}{2}\coth \left(\frac{\rho b(0)}{2}\right). \end{aligned}$$

The discrete scheme

Using the higher order finite difference scheme of (16) and inducing the exponential fitting factor (23) for \(1\le i \le N-1\) and \(0 \le n \le M-1\), then the fully discrete scheme is given as

$$\begin{aligned} L^{\Delta t,h}Z_{i}^{n}= \frac{1}{30}(R(s_{i-1},t_{n}) +28R(s_{i},t_{n})+R(s_{i+1},t_{n})), \end{aligned}$$

where \(L^{\Delta t,h}Z_{i}^{n}=-(\varepsilon \sigma -\frac{hb(s_{i-1},t_{n})}{30}+\frac{hb(s_{i+1},t_{n})}{30})(\frac{Z_{i-1}^{n} - 2Z_{i}^{n} + Z_{i+1}^{n}}{h^2}) + \frac{b(s_{i-1},t_{n})}{60h}(-3Z_{i-1}^{n}+4Z_{i}^{n}-Z_{i+1}^{n}) +\frac{28b(s_{i},t_{n})}{60h}(Z_{i+1}^{n}-Z_{i-1}^{n}) +\frac{b(s_{i+1},t_{n})}{60h}(Z_{i-1}^{n}-4Z_{i}^{n}+3Z_{i+1}^{n}) + \frac{P(s_{i-1},t_{n})}{30}Z_{i-1}^{n}+\frac{28P(x_{i},t_{n})}{30}Z_{i}^{n} +\frac{P(s_{i+1},t_{n})}{30}Z_{i+1}^{n}.\)

In the explicit form, we write

$$\begin{aligned} r_{i}^{-}Z_{i-1}^{n}+r_{i}^{0}Z_{i}^{n}+r_{i}^{+}Z_{i+1}^{n} = H_{i}^{n}, \end{aligned}$$


$$\begin{aligned} \begin{aligned} r_{i}^{-} =&-\frac{1}{h^2}\left(\varepsilon \sigma - \frac{hb(s_{i-1},t_{n})}{30} +\frac{hb(s_{i+1},t_{n})}{30}\right) - \frac{3b(s_{i-1},t_{n})}{60h}+\frac{P(s_{i-1},t_{n})}{30}\\&-\frac{28b(s_{i},t_{n})}{60h}+\frac{b(s_{i+1},t_{n})}{60h}, \\ r_{i}^{0} =&\frac{2}{h^2}\left(\varepsilon \sigma - \frac{hb(s_{i-1},t_{n})}{30} +\frac{hb(s_{i+1},t_{n})}{30}\right) + \frac{4b(s_{i-1},t_{n})}{60h}-\frac{4b(s_{i+1},t_{n})}{60h}+\frac{28P(s_{i},t_{n})}{30}, \\ r_{i}^{+} =&-\frac{1}{h^2}\left(\varepsilon \sigma - \frac{hb(s_{i-1},t_{n})}{30} +\frac{hb(s_{i+1},t_{n})}{30}\right) - \frac{b(s_{i-1},t_{n})}{60h}+\frac{28b(s_{i},t_{n})}{60h}\\ {}&+\frac{3b(s_{i+1},t_{n})}{60h} +\frac{P(s_{i+1},t_{n})}{30}, \\ H_{i}^{n} =&\frac{1}{30}(R(s_{i-1},t_{n}) +28R(s_{i},t_{n})+R(s_{i+1},t_{n})). \end{aligned} \end{aligned}$$

Stability and convergence analysis

In this part, for the developed scheme of (24) we take to prove the discrete comparison principle.

Lemma 3.3

There is a comparison function \(\nu _{i}^{n}\) such that \({\mathcal {L}}^{\Delta t,h}Z_{i}^{n} \le {\mathcal {L}}^{\Delta t,h}\nu _{i}^{n}\) for \(1 \le i \le N-1\) and if \(Z_{0}^{n}\le \nu _{0}^{n}\) and \(Z_{N}^{n}\le \nu _{N}^{n}\), then  \(Z_{i}^{n} \le \nu _{i}^{n}\) for \(1 \le i\le N\).


The discrete operator \({\mathcal {L}}^{\Delta t,h}Z_{i}^{n}\) is matrix of size \((N+1) \times (N+1)\) with its entries for \(1\le i \le N-1\) are \(r_{i}^{-}, \;r_{i}^{0},\) and \(r_{i}^{+}.\) We observe that \(\vert r_{i}^{-} |> 0, \ \vert r_{i}^{0} |> 0, \ \vert r_{i}^{+} |>0\) and \(\vert r_{i}^{0} |\ge \vert r_{i}^{-} |+ \vert r_{i}^{+} |\), giving that the matrix is diagonally dominant. Then, it satisfies the property of M matrix. Thus, the non-negative inverse of the matrix exists. So, it guarantees the existence of unique discrete solution [22, 24]. \(\square\)

Lemma 3.4

(Stability result) If the solution of (24) be \(Z_{i}^{n}\), then it satisfies

$$\begin{aligned} \big \vert Z_{i}^{n}\big |\le \frac{\Vert \ {\mathcal {L}}^{\Delta t,h}Z_{i}^{n}\Vert }{ \zeta } + \max \{\vert \psi _{l}(t_{n})\vert ,\vert \psi _{r}(t_{n})|\}, \end{aligned}$$

where \(P(s_{i},t_{n}) \ge \zeta > 0\).


let \(\Pi = \frac{\Vert \ {\mathcal {L}}^{\Delta t,h}Z_{i}^{n}\Vert }{ \zeta } + \max \{\vert \psi _{l}(t_{n})\vert ,\vert \psi _{r}(t_{n})|\}\) and set the barrier functions \(\vartheta _{i,n}^{\pm }\) by \(\vartheta _{i,n}^{\pm } = \Pi \pm Z_{i}^{n}\). On the boundaries, we obtain \(\vartheta _{0,n}^{\pm } = \Pi \pm Z_{0}^{n} = \frac{\Vert \ {\mathcal {L}}^{\Delta t,h}Z_{i}^{n}\Vert }{ \zeta } + \max \{\vert \psi _{l}(t_{n})\vert ,\vert \psi _{r}(t_{n})|\} \pm \psi _{l}(t_{n}) \ge 0\),

\(\vartheta _{N,n}^{\pm } = \Pi \pm Z_{N}^{n} = \frac{\Vert \ {\mathcal {L}}^{\Delta t,h}Z_{i}^{n}\Vert }{ \zeta } + \max \{\vert \psi _{l}(t_{n})\vert ,\vert \psi _{r}(t_{n})|\} \pm \psi _{r}(t_{n}) \ge 0\). On the discretized spatial domain \(s_{i} , 1< i < N-1,\) we have

$$\begin{aligned} \begin{aligned} {\mathcal {L}}^{\Delta t,h}\vartheta _{i,n} =&- \left(\varepsilon \sigma -\frac{hb(s_{i-1},t_{n})}{30}+\frac{hb(s_{i+1},t_{n})}{30}\right)\left(\frac{\Pi \pm Z_{i-1}^{n} - 2(\Pi \pm Z_{i}^{n}) + \Pi \pm Z_{i+1}^{n}}{h^2}\right)\\&+ \frac{ b(s_{i-1},t_{n})}{60h}(-3(\Pi \pm Z_{i-1}^{n})+4(\Pi \pm Z_{i}^{n})-(\Pi \pm Z_{i+1}^{n})) \\ {}&+\frac{28b(s_{i},t_{n})}{60h}(\Pi \pm Z_{i+1}^{n}-(\Pi \pm Z_{i-1}^{n})) \\&+\frac{b(s_{i+1},t_{n})}{60h}(\Pi \pm Z_{i-1}^{n}-4(\Pi \pm Z_{i}^{n})+3(\Pi \pm Z_{i+1}^{n})) \\&+ \frac{P(s_{i-1},t_{n})}{30}(\Pi \pm Z_{i-1}^{n}) +\frac{28P(s_{i},t_{n})}{30}(\Pi \pm Z_{i}^{n}) +\frac{P(s_{i+1},t_{n})}{30}(\Pi \pm Z_{i+1}^{n}) \\ =&\left(\frac{P(s_{i-1},t_{n})}{30} +\frac{28P(s_{i},t_{n})}{30} +\frac{P(s_{i+1},t_{n})}{30}\right)\Pi \pm {\mathcal {L}}^{\Delta t,h}Z_{i}^{n}\\ =&\left(\frac{P(s_{i-1},t_{n})}{30} +\frac{28P(s_{i},t_{n})}{30} +\frac{P(s_{i+1},t_{n})}{30}\right)\left(\frac{\Vert \ {\mathcal {L}}^{\Delta t,h}Z_{i}^{n}\Vert }{ \zeta } + \max \{\vert \psi _{l}(t_{n})\vert ,\vert \psi _{r}(t_{n})|\}\right) \pm \frac{1}{30}(R_{i-1}^{n} +28R_{i}^{n} +R_{i+1}^{n})\ge 0. \end{aligned} \end{aligned}$$

From the discrete comparison principle, we get \(\vartheta _{i,n}^{\pm } \ge 0,\ i = 0,1,2,...,N\). Hence, the necessary bound is satisfied. \(\square\)

Lemma 3.5

If \(\nu _{i}^{n}\) be any mesh function such that \(\nu _{0}^{n} = \nu _{N}^{n} = 0\). Then it satisfies

$$\begin{aligned} \big \vert \nu _{i}^{n}\big |\le \frac{1}{\zeta }\max _{m}\big \vert {\mathcal {L}}^{\Delta t, h}\nu _{m}^{n}\big |. \end{aligned}$$

Using Taylor’s series approximation, we have the bounds

$$\begin{aligned} \begin{aligned} \bigg \vert - (\frac{d}{ds^2} - \delta ^{2}_{s})Z^{n}(s_i)\bigg |\le&Ch^{2} \bigg \Vert \frac{d^{4}Z^{n}(s_i)}{ds^4}\bigg \Vert ,\\ \bigg \vert \frac{dZ^{n}(s_{i-1})}{ds} -(\frac{-Z_{i+1}^{n} +4Z_{i}^{n}- 3Z_{i-1}^{n}}{2h} + h\frac{d^{2}Z^{n}(s_{i})}{ds^{2}}) \bigg |\le&Ch^{2} \bigg \Vert \frac{d^{3}Z^{n}(s_i)}{ds^3}\bigg \Vert ,\\ \bigg \vert \frac{dZ^{n}(s_{i+1})}{ds} -(\frac{3Z_{i+1}^{n} -4Z_{i}^{n}+ Z_{i-1}^{n}}{2h} - h\frac{d^{2}Z^{n}(s_{i})}{ds^{2}}) \bigg |\le&Ch^{2}\bigg \Vert \frac{d^{3}Z^{n}(s_i)}{ds^3}\bigg \Vert ,\\ \bigg \vert (\frac{d}{ds}-\delta ^{0}_{s})Z^{n}(s_i)\bigg |\le Ch^{2} \bigg \Vert \frac{d^{3}Z^{n}(s_i)}{ds^3}\bigg \Vert , \bigg \vert \delta ^{2}_{s}Z^{n}(s_i) \bigg |\le&C\bigg \Vert \frac{d^{2}Z^{n}(s_i)}{ds^{2}}\bigg \Vert , \end{aligned} \end{aligned}$$

where \(\Vert Z^{(k)}(s_i)\Vert = \max _{s_{i}}\vert Z_{n}^{(k)}(s_i) |\ , k = 2, 3, 4\).

In the next theorem we bound the truncation error in space direction discretization.

Theorem 3.1

Consider the sufficiently smooth functions \(a(s), \ b(s,t_{n})\) and \(c(s,t_{n})\) of (67) so that \(Z^{n}(s)\in C^{4}[0,1]\). Then, the solution \(Z_{i}^{n}\) of (24) satisfies the bound

$$\begin{aligned} \big \vert {\mathcal {L}}^{\Delta t, h}(Z^{n}(s_{i})-Z_{i}^{n})\big |\le \frac{Ch^{2}}{h+\varepsilon }\left(1+\varepsilon ^{-3}\exp \left(-\frac{\beta }{\varepsilon }s_{i}\right)\right). \end{aligned}$$


In the spatial direction the local truncation error is given by

$$\begin{aligned} \begin{aligned} \bigg \vert {\mathcal {L}}^{\Delta t, h}(Z^{n}(s_{i})&-Z_{i}^{n})\bigg |= \bigg \vert - \varepsilon \left(\frac{d}{ds^2} -\sigma \delta ^{2}_{s}\right)Z^{n}(s_i) + \frac{ b(s_{i-1},t_{n})}{60h}\\ {}&\left(\frac{dZ^{n}(s_{i-1})}{ds} -\left(\frac{-Z_{i+1}^{n} +4Z_{i}^{n}- 3Z_{i-1}^{n}}{2h} + h\frac{d^{2}Z^{n}(s_{i})}{ds^{2}}\right)\right) \\&+\frac{28b(s_{i},t_{n})}{60h}\left(\frac{d}{ds}-\delta ^{0}_{s}\right)Z^{n}(s_i) +\frac{b(s_{i+1},t_{n})}{60h}\\ {}&+ \left(\frac{dZ^{n}(s_{i+1})}{ds}-\left(\frac{3Z_{i+1}^{n} -4Z_{i}^{n}+ Z_{i-1}^{n}}{2h} - h\frac{d^{2}Z^{n}(s_{i})}{ds^{2}}\right)\right)\bigg |\\ \le&\bigg \vert \left(\varepsilon (b(s_{i},t_{n})\frac{\rho }{2}\coth \left(b(0)\frac{\rho }{2}\right) -1\right)\delta ^{2}_{s}Z^{n}(s_i))\bigg |+\bigg \vert \varepsilon \left(\frac{d^2}{ds^2}-\delta ^{2}_{s}\right)Z^{n}(s_{i})\bigg |\\ {}& \bigg \vert \frac{ b(s_{i-1},t_{n})}{60h}(\frac{dZ^{n}(s_{i-1})}{d s}-(\frac{-Z_{i+1}^{n} +4Z_{i}^{n}- 3Z_{i-1}^{n}}{2h} + h\frac{d^{2}Z^{n}(s_{i})}{d s^{2}}))\bigg |\\&+\bigg \vert \frac{28b(s_{i},t_{n})}{60h}\left(\frac{d}{ds}-\delta ^{0}_{s}\right)Z^{n}(s_i) \bigg |+\bigg \vert \frac{b(s_{i+1},t_{n})}{60h}\\ {}&\left(\frac{d Z^{n}(s_{i+1})}{ds}-\left(\frac{3Z_{i+1}^{n} -4Z_{i}^{n}+ Z_{i-1}^{n}}{2h} - h\frac{d^{2}Z^{n}(s_{i})}{ds^{2}}\right)\right) \bigg |, \end{aligned} \end{aligned}$$

where \(\sigma = b(s_{i},t_{n})\frac{\rho }{2}\coth (b(0)\frac{\rho }{2}) \ and \ \rho = \frac{h}{\varepsilon }\).

For the constants \(C_{1}\) and \(C_{2}\) we have \(\big \vert \rho \coth (\rho )-1\big |\le C_{1}\rho ^{2}\), for \(\rho \le 1\). For \(\rho \rightarrow \infty\), since \(\lim _{\rho \rightarrow \infty }\coth (\rho )=1\) which gives \(\big \vert \rho \coth (\rho )-1\big |\le C_{1}\rho\).

Generally, \(\forall \rho > 0\), we put as

$$\begin{aligned} C_{1}\frac{\rho ^{2}}{\rho +1}\le \rho \coth (\rho )-1\le C_{2}\frac{\rho ^{2}}{\rho +1}, \end{aligned}$$

we obtain

$$\begin{aligned} \varepsilon [b(s_{i},t_{n})\frac{\rho }{2}\coth \left(b(0)\frac{\rho }{2}\right)-1]\le \varepsilon \frac{(h/\varepsilon )^{2}}{h/\varepsilon +1} = \frac{h^2}{h+\varepsilon }. \end{aligned}$$

Using the bound for the difference of the derivatives in (29) and (32), we obtain

$$\begin{aligned} \bigg \vert {\mathcal {L}}^{\Delta t, h}(Z^{n}(s_{i})-Z_{i}^{n})\bigg \vert \le \frac{Ch^2}{h+\varepsilon } \bigg \Vert \frac{d^{2}Z^{n}(s_i)}{ds^{2}}\bigg \Vert +Ch^{2}\bigg \Vert \frac{d^{3}Z^{n}(s_i)}{ds^3}\bigg \Vert + C\varepsilon h^{2}\bigg \Vert \frac{d^{4}Z^{n}(s_{i})}{ds^{4}}\bigg \Vert . \end{aligned}$$

From Lemma 3.2, we obtain the bound for the derivatives

$$\begin{aligned} \begin{aligned} \bigg \vert {\mathcal {L}}^{\Delta t, h}(Z^{n}(s_{i})-Z_{i}^{n})\bigg \vert \le&\frac{Ch^2}{h+\varepsilon }\left(1+\varepsilon ^{-2}\exp \left(-\frac{\beta }{\varepsilon }s_{i}\right)\right) +Ch^{2}\left[\left(1+\varepsilon ^{-3}\exp \left(-\frac{\beta }{\varepsilon }s_{i}\right)\right)+ \left(\varepsilon +\varepsilon ^{-3}\exp \left(-\frac{\beta }{\varepsilon }s_{i}\right)\right)\right]. \end{aligned} \end{aligned}$$

Evidently, \(\varepsilon ^{-2} \le \varepsilon ^{-3}\), then we obtain \(\big \vert {\mathcal {L}}^{\Delta t, h}(Z^{n}(s_{i})-Z_{i}^{n})\big |\le \frac{Ch^2}{h+\varepsilon }(1+\varepsilon ^{-3}\exp (-\frac{\beta }{\varepsilon }s_{i}))\) thus, it gives the desired bound. \(\square\)

Lemma 3.6

For a fixed mesh and as \(\varepsilon \rightarrow 0,\) it gives

$$\begin{aligned} \lim _{\varepsilon \rightarrow 0}\max _{i}\frac{\exp (-\beta s_{i}/\varepsilon )}{\varepsilon ^n} = 0, \ n = 1,2,3,..., \end{aligned}$$

where \(s_{i} = ih, \ 1 \le i \le N-1\).


The proof is considered in [9]. \(\square\)

Theorem 3.2

Let \(Z_{i}^{n}\) be the solution of (24), then we have the following uniform error bound

$$\begin{aligned} \sup _{\varepsilon \in (0, 1]}\max _{i}\big \vert Z^{n}(s_{i})-Z_{i}^{n}\big |\le Ch, \ i = 0,1,2,...N. \end{aligned}$$


Plugging the result in Lemma 3.6 into (30), we arrive at

$$\begin{aligned} \bigg \vert {\mathcal {L}}^{\Delta t, h}(Z^{n}(s_{i})-Z_{i}^{n})\bigg |\le \frac{Ch^2}{h+\varepsilon }. \end{aligned}$$

Hence, the result leads \(\big \vert Z^{n}(s_i)-Z_{i}^{n}\big |\le \frac{Ch^2}{h+\varepsilon }.\) Using the \(\sup\) over all \(\varepsilon \in (0,1]\), we get

$$\begin{aligned} \sup _{\varepsilon \in (0, 1]} \max _{i}\big \vert Z^{n}(s_{i})-Z_{i}^{n}\big |\le Ch. \end{aligned}$$

From the preceding theorem for the case when \(\varepsilon > h\), the obtained scheme gives second order uniformly convergent. For the case when \(\varepsilon \ll h\), the scheme is first order uniformly convergent in spatial direction. \(\square\)

Theorem 3.3

Let z and Z are the solutions of (1) and (24) respectively, then we have the following uniform error bound

$$\begin{aligned} \sup _{\varepsilon \in (0, 1]}\big \vert z-Z \big |\le C(h+(\Delta t)). \end{aligned}$$


The proof can be done by the combination of Lemma 3.1 and Theorem 3.2. \(\square\)

Numerical results and discussions

In this part, we are considering two model examples to validate the theoretical results obtained by the proposed method. If the exact solutions of the considered examples are not known the maximum pointwise error is estimated by using the double mesh principle. So, the maximum pointwise error is calculated by \(E_{\varepsilon }^{N,M} = \max _{i,n}\vert Z_{i,n}^{N,M}- Z_{i,n}^{2N,2M}|,\) and the \(\varepsilon\)-uniform error is estimated by \(E^{N,M} = \max _{i,n}(E_{\varepsilon }^{N,M}).\) The rate of convergence is calculated by \(r^{N,M}_{\varepsilon } = \log 2(E^{N,M}_{\varepsilon }/E^{2N,2M}_{\varepsilon }),\) and the \(\varepsilon\)-uniform rate of convergence is computed by \(r^{N,M} = \log 2(E^{N,M}/E^{2N,2M}).\)

Example 4.1

Consider the problem \(\frac{\partial u}{\partial t} - \varepsilon \frac{\partial ^2u}{\partial x^2} - \frac{\partial u}{\partial x}+(\frac{1+x^2}{2})u = t^3 - u(x,t-\tau ),~~(x,t)\in (0,1)\times (0,2]\) with interval condition u(x,t) = 0, on \((x,t) \in [0,1]\times [-\tau ,0]\) and the boundary conditions \(u(0,t) = 0~~ and ~~u(1,t) = 0,~~~t\in (0,2].\)

Example 4.2

Consider the problem \(\frac{\partial u}{\partial t} - \varepsilon \frac{\partial ^2u}{\partial x^2} -(2+x^{2}) \frac{\partial u}{\partial x}+xu(x,t) + u(x,t-\tau ) = 10t^{2}\exp (-t)x(1-x),\)

\((x,t)\in (0,1)\times (0,2]\)

with interval condition \(u(x,t) = 0\), on \((x,t) \in [0,1]\times [-\tau ,0]\) and the boundary conditions \(u(0,t) = 0~~ and ~~u(1,t) = 0,~~~t\in (0,2].\)

Fig. 1
figure 1

Numerical solution of Example 4.1 on a, \(\varepsilon = 2^{-6}\) and b, \(\varepsilon = 2^{-20}\)

Table 1 Example 4.1, maximum absolute errors for the case \(\tau = 0.5\varepsilon\)
Table 2 Example 4.2, maximum absolute errors for the case \(\tau = 0.5\varepsilon\)

For different values of \(\varepsilon\) and an equal number of mesh points the maximum pointwise error, \(\varepsilon\) -uniform error and \(\varepsilon\)-uniform rate of convergence of the proposed method are displayed in Tables 1 and 2. We observe from these Tables, as \(\varepsilon \rightarrow 0\), the maximum pointwise error after showing increment remains uniform. This shows that the scheme is stable and uniformly convergent irrespective of the values of \(\varepsilon\). The \(\varepsilon\)-uniform error and \(\varepsilon\)-uniform rate of convergence of the method are indicated in the last rows of each Tables and it confirms that the numerical results agree with the theoretical result.

In Fig. 1, we show the numerical solution of the scheme for Example 4.1 for different values of \(\varepsilon\). From these Figures, it can be seen that a strong boundary layer is created on the left side of the spatial domain for small \(\varepsilon\).

In the last section of Table 1, the comparison of the proposed scheme with the results of the existing published work of [14] is given. As one observes, the developed scheme gives more accurate result than the scheme in [14].


We developed a numerical scheme to solve a singularly perturbed parabolic convection-diffusion equation that exhibits a boundary layer. The proposed scheme consists of the backward-Euler method in the time direction and an exponentially fitted finite difference scheme for the spatial direction. Using the comparison principle, the stability of the discrete scheme is examined and analysed. The uniform convergence of the scheme is discussed theoretically. To validate the theoretical finding of the scheme, we considered two model examples and the numerical results are given by applying maximum point-wise absolute error, \(\varepsilon\)-uniform error and \(\varepsilon\)-uniform rate of convergence in Tables. The proposed method contributes more accurate, stable and \(\varepsilon\)-uniform numerical result with linear order of convergence.


  • The proposed scheme is not layer resolving method since there is no sufficient number of mesh points in the boundary layer region.

Availability of data and materials

No additional data is used for this research work.


  1. Baker CTH, Bocharov GA, Rihan FA. A report on the use of delay differential equations in numerical modelling in the biosciences. Citeseer. 1999.

  2. Kumar K, Gupta T, Chakravarthy PP, Rao RN. An adaptive mesh selection strategy for solving singularly perturbed parabolic partial differential equations with a small delay. Appl Math Sci Comput. 2019; 67–76.

  3. Das A, Natesan S. Uniformly convergent hybrid numerical scheme for singularly perturbed delay parabolic convection-diffusion problems on Shishkin mesh. Appl Math Comput. 2015;271:168–86.

    Google Scholar 

  4. Govindarao L, Mohapatra J. A second order numerical method for singularly perturbed delay parabolic partial differential equation. Comput Eng. 2019.

  5. Gowrisankar S, Natesan S. \(\varepsilon\)-Uniformly convergent numerical scheme for singularly perturbed delay parabolic partial differential equations. Int J Comput Math. 2017;94(5):902–21.

    Article  Google Scholar 

  6. Das A, Natesan S. Second-order uniformly convergent numerical method for singularly perturbed delay parabolic partial differential equations. Int J Comput Math. 2018;95(3):490–510.

    Article  Google Scholar 

  7. Podila PC, Kumar K. A new stable finite difference scheme and its convergence for time-delayed singularly perturbed parabolic PDES. Comput Appl Math. 2020;39(3):1–16.

    Article  Google Scholar 

  8. Negero NT, Duressa GF. An efficient numerical approach for singularly perturbed parabolic convection-diffusion problems with large time-lag. J Math Model. 2021; 1–18.

  9. Woldaregay MM, Aniley WT, Duressa GF. Novel numerical scheme for singularly perturbed time delay convection-diffusion equation. Adv Math Phys. 2021.

  10. Babu G, Bansal K. A high order robust numerical scheme for singularly perturbed delay parabolic convection diffusion problems. J Appl Math Comput. 2022;68(1):363–89.

    Article  Google Scholar 

  11. Kumar D, Kumari P. A parameter-uniform numerical scheme for the parabolic singularly perturbed initial boundary value problems with large time delay. J Appl Math Comput. 2019;59(1–2):179–206.

    Article  Google Scholar 

  12. Kumar D, Kumari P. A parameter-uniform scheme for singularly perturbed partial differential equations with a time lag. Numer Methods Partial Differ Equ. 2020;36(4):868–86.

    Article  Google Scholar 

  13. Negero NT, Duressa GF. A method of line with improved accuracy for singularly perturbed parabolic convection-diffusion problems with large temporal lag. Results Appl Math. 2021;11:100–74.

    Article  Google Scholar 

  14. Negero NT, Duressa GF. Uniform convergent solution of singularly perturbed parabolic differential equations with general temporal-lag. Iran. J Sci Technol Trans A Sci. 2022; 1–18.

  15. Salama A, Al-Amery D. A higher order uniformly convergent method for singularly perturbed delay parabolic partial differential equations. Int J Comput Math. 2017;94(12):2520–46.

    Article  Google Scholar 

  16. Selvi PA, Ramanujam N. A parameter uniform difference scheme for singularly perturbed parabolic delay differential equation with Robin type boundary condition. Appl Math Comput. 2017;296:101–15.

    Google Scholar 

  17. Das P. A higher order difference method for singularly perturbed parabolic partial differential equations. J Differ Equations Appl. 2018;24(3):452–77.

    Article  Google Scholar 

  18. Das P, Vigo-Aguiar J. Parameter uniform optimal order numerical approximation of a class of singularly perturbed system of reaction diffusion problems involving a small perturbation parameter. J Comput Appl Math. 2019;1(354):533–44.

    Article  Google Scholar 

  19. Das P, Mehrmann V. Numerical solution of singularly perturbed convection-diffusion-reaction problems with two small parameters. BIT Numer Math. 2016;56(1):51–76.

    Article  Google Scholar 

  20. Clavero C, Jorge J, Lisbona F. A uniformly convergent scheme on a nonuniform mesh for convection-diffusion parabolic problems. Comput Appl Math. 2003;154(2):415–29.

    Article  Google Scholar 

  21. Kumar K, Podila PC, Das P, Ramos H. A graded mesh refinement approach for boundary layer originated singularly perturbed time-delayed parabolic convection diffusion problems. Math Methods Appl Sci. 2021;44(16):12332–50.

    Article  Google Scholar 

  22. Kellogg RB, Tsan A. Analysis of some difference approximations for a singular perturbation problem without turning points. Math Comp. 1978;32(144):1025–39.

    Article  Google Scholar 

  23. O’malley RE. Singular perturbation methods for ordinary differential equations. vol. 89, Springer: USA; 1991

  24. Das P, Natesan S. Higher-order parameter uniform convergent schemes for Robin type reaction-diffusion problems using adaptively generated grid. 2012;9(04):1250052.

    Google Scholar 

Download references


The authors thanks the editor and reviewers for their constructive comments.


No funding organization for this research work.

Author information

Authors and Affiliations



SKT  and MMW carried out the numerical scheme development, algorithms writing, MATLAB code writing and the numerical experimentation. TGD and GFD formulated the problem, design, and write draft of the manuscript. All authors read, commented and approved the final version of the manuscript.

Corresponding author

Correspondence to Mesfin Mekuria Woldaregay.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.



Proof of Lemma 3.1 Since the function \(Z^{n}(s)\) satisfies

$$\begin{aligned}&(1-\tau c(s,t_{n}) + \Delta t {\mathcal {L}})Z^{n}(s) - \Delta g(s, t_{n}) = Z^{n-1}(s)\\ Z(t_{n-1}) =&(1-\tau c(t_n)+\Delta t {\mathcal {L}})Z(t_n) -\Delta g(t_n)+ \int _{t_{n-1}}^{t_n}(t_{n-1}-s){\frac{\partial ^{2}Z(s)}{\partial t^2}}ds\\ =&(1-\tau c(t_n)+\Delta t {\mathcal {L}})Z(t_n) -\Delta g(t_n) +O(\Delta t^2). \end{aligned}$$

Then \(e_{n}(s)\) is the solution of boundary value problem of type

\(((1-\tau c(s,t_n)) +\Delta t {\mathcal {L}})e_{n}(s) = O(\Delta t^2)\).

\(e_{n}(0) = 0 = e_{n}(1)\) by manimum principle, we obtain \(\Vert e_{n} \Vert \le C_{1}(\Delta t^2)\).

Taking the summation of all local error estimate up to \(n^{th}\) time step the global error estimate is given by

$$\begin{aligned} \begin{aligned} \Vert E_{n}\Vert =&\bigg \Vert \sum _{l=1}^{n}e_{l}\bigg \Vert \le \Vert e_{1}\Vert + \Vert e_{2}\Vert + \Vert e_{3} \Vert + . . . + \Vert e_{n}\Vert \\ \le&C_{1}T(\Delta t) , ~~since~~ (n)\Delta t \le T \\ =&C(\Delta t), ~~~ where~~ C_{1}T = C, \end{aligned} \end{aligned}$$

where C is a constant not depends on \(\varepsilon\) and \(\Delta t\).

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Tesfaye, S.K., Woldaregay, M.M., Dinka, T.G. et al. Fitted computational method for solving singularly perturbed small time lag problem. BMC Res Notes 15, 318 (2022).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • Accurate numerical method
  • Exponentially fitted method
  • Stability and uniform convergence