\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}
\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2013 (2013), No. 30, pp. 1--10.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2013 Texas State University - San Marcos.}
\vspace{9mm}}
\begin{document}
\title[\hfilneg EJDE-2013/30\hfil A computational technique]
{A computational technique for solving boundary value problem with two
small parameters}
\author[Devendra Kumar \hfil EJDE-2013/30\hfilneg]
{Devendra Kumar}
\address{Devendra Kumar \newline
Department of Mathematics, BITS Pilani
Pilani - 333031, Rajasthan, India\newline
Phone +919413487072; fax: +911596244183}
\email{dkumar2845@gmail.com, dkumar@bits-pilani.ac.in}
\thanks{Submitted November 2, 2012. Published January 28, 2013.}
\subjclass[2000]{65L03, 65L10, 65L11, 65L20}
\keywords{Finite difference method; singular perturbation;
\hfill\break\indent delay differential equation; negative shift}
\begin{abstract}
In this article we study a singularly perturbed boundary-value problem
for a delay differential equation with a small delay parameter in the
first derivative term whose solution has a single boundary layer.
The proposed method is shown to be stable, and its performance is
confirmed with examples.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{example}[theorem]{Example}
\allowdisplaybreaks
\section{Introduction}
Delay differential equations arise in the mathematical modelling of various
practical phenomena, for instance, micro scale heat transfer \cite{tzou},
hydrodynamics of liquid helium \cite{joseph1}, second-sound theory \cite{joseph2},
thermo elasticity \cite{ezzat}, diffusion in polymers \cite{liu},
reaction-diffusion equations \cite{best}, stability \cite{burton},
control including control of chaotic systems \cite{liao}, a variety of
model for physiological processes or diseases \cite{waz,mack} etc.
As a consequence, they have received a lot of interest in recent years,
especially for linear problems, for which one can obtain analytical solutions
by means of, for example, the Laplace transform in time, separation of variables
in finite spatial domains, etc.
A delay differential equation is said to be of retarded type if the delay
argument does not occur in the highest order derivative term.
If we restrict it to a class in which the highest derivative term
is multiplied by a small parameter, then we obtain singularly perturbed delay
differential equations of the retarded type. Frequently, delay
differential equations have been reduced to differential equations with
coefficients that depend on the delay by means of first order accurate
Taylor's series expansions of the terms that involve delay and the
resulting differential equations have been solved either analytically
when the coefficients of these equations are constant or numerically,
when they are not \cite{ramos}.
It is well-known that the standard discretization methods for solving
singular perturbation problems are unstable and fail to give accurate
results when the perturbation parameter $\epsilon$ is small.
Therefore, it is important to develop suitable numerical methods for these
problems, whose accuracy does not depend on the parameter value $\epsilon$;
i.e., the methods that are convergent $\epsilon$-uniformly
\cite{doolan, farrell, roos}. There are essentially two strategies
to design schemes which have small truncation errors inside the layer region(s).
The first approach which forms the class of fitted mesh methods consists
in choosing a fine mesh in the layer region(s). The second approach is
in the context of the fitted operator methods in which the mesh remains
uniform and the difference schemes reflect the qualitative behavior of
the solution inside the layer region(s). A nice discussion using one or
both of the above strategies can be found in Miller {\em et al.} \cite{doolan}.
The work in this paper falls under the first category. We have derived a
finite difference scheme on a non uniform mesh for the boundary value problems
for a class of singularly perturbed delay differential equations.
\section{Statement of the problem}
Consider the following boundary-value problem for a singularly perturbed
delay differential equation with a small parameter multiplying to
the second derivative and containing negative shift in the first
derivative term
\begin{equation}
\label{s1}\epsilon y''(x)+p(x)y'(x-\delta)-q(x)y(x)=r(x),\quad x\in [0,1],
\end{equation}
under the interval and boundary conditions
\begin{equation}
\label{s2} y(x)=\phi(x),\quad -\delta\leq x\leq 0,\quad y(1) = \beta,
\end{equation}
where $\epsilon$ is a small parameter, $0<\epsilon\leq 1$, and $\delta$
is also a small shift parameter of $o(\epsilon)$. The functions
$p(x)$, $q(x)$, $r(x)$ and $\phi(x)$ are sufficiently smooth.
It is assumed that $p(x)\geq p^*>0$, $q(x)\geq q^* > 0$, for all $x\in [0,1]$
for some positive constants $p^*$ and $q^*$. For $\delta=0$
Equation \eqref{s1} reduces to a singularly perturbed ordinary differential
equation with only a single parameter $\epsilon$, which has a boundary
layer at one end or both ends depending on $p(x)$ and $q(x)$.
It is well-known that if $p(x)>0$ throughout the interval $[0,1]$,
the boundary layer exists near left end $x=0$ and if $p(x)<0$ throughout
the interval $[0,1]$, the boundary layer exists near right end $x=1$.
However, for $p(x)=0$ the layer exists at interior of the interval $[0,1]$
and such point is called the turning point. In this paper, we consider
the problems which have boundary layers (not interior layers) only.
This problem has already been considered by Kadalbajoo and Kumar
in \cite{kadal}. There we have used B-spline collocation method with
piecewise-uniform mesh and the method was shown to be parameter uniform
of order almost two. In this paper we have generated a geometric mesh
and used finite difference method with geometric mesh.
The proposed method is not very useful in the case when the delay
parameter is relatively large as the Taylor's series expansion is not
valid in that case.
Since $\epsilon$ is small and $\delta=o(\epsilon)$, so taking the
Taylor's series expansion for the term $y'(x-\delta)$,
from \eqref{s1}-\eqref{s2}, we obtain
\begin{equation}
\label{s5}\epsilon y''(x)+a(x)y'(x)-b(x)y(x)=f(x),
\end{equation}
with
\begin{equation}
\label{s6}y(0)=\phi_0=\phi(0),\quad y(1)=\beta,
\end{equation}
where
\[
a(x) = \frac{p(x)}{1-(\delta/\epsilon)p(x)},\quad
b(x) = \frac{q(x)}{1-(\delta/\epsilon)p(x)},\quad
f(x) = \frac{r(x)}{1-(\delta/\epsilon)p(x)}.
\]
Here we also assume that $\epsilon-\delta p(x)> 0$,
then we have $a(x)\geq a^*>0$ and $b(x)\geq b^*>0$ for some positive
constants $a^*$ and $b^*$.
The operator $L^{c}_{\epsilon}=\epsilon \frac{d^2}{dx^2}+a(x)\frac{d}{dx}-b(x)I$
in \eqref{s5} satisfies the following minimum principle:
\begin{lemma}\label{l1}
Suppose $\pi(x)$ be any sufficiently smooth function satisfying
$\pi(0)\geq 0$ and $\pi(1)\geq 0$. Then $L^{c}_{\epsilon}\pi(x)\leq 0$
for all $x\in(0,1)$ implies that $\pi(x)\geq 0$ for all $x\in [0,1]$.
\end{lemma}
\begin{proof}
Let $z\in[0,1]$ be such that $\pi(z)<0$ and $\pi(z)=\min_{0\leq x\leq 1} \pi(x)$.
Clearly $z \notin\{0,1\}$, therefore $\pi'(z)=0$ and $\pi''(z)\geq 0$.
Now we have from Eq.~\eqref{s5}
\[
L^{c}_{\epsilon} \pi(z)=\epsilon \pi''(z)+a(z)\pi'(z)-b(z)\pi(z)>0,
\]
which contradict our assumption, therefore we must have $\pi(z)\geq 0$
and thus $\pi(x)\geq 0,\;\forall\; x\in [0,1]$.
\end{proof}
Now we can show the boundedness of the solutions of the continuous
problem \eqref{s5}-\eqref{s6}.
\begin{lemma}\label{l2}
The solution $y(x)$ of \eqref{s5}-\eqref{s6} satisfies the inequality
\[
\|y\|\leq C \max\big\{|\phi_0|, |\beta|, \frac{1}{b^*}\|f\|\big\},
\]
where $\|\cdot\|$ is the $l_{\infty}$ norm given by
$\|y\|=\max_{0\leq x\leq 1}|y(x)|$.
\end{lemma}
\begin{proof}
Consider the barrier functions $\psi^{\pm}(x)$ defined by
\[
\psi^{\pm}(x)=\max\big\{| \phi_0|, | \beta|,\frac{1}{b^*}
\| f\|\big\}\pm y(x).
\]
Then we have
\begin{align*}
\psi^{\pm}(0)
&= \max\big\{|\phi_0|, | \beta|, \frac{1}{b^*}\| f\|\big\}\pm y(0)\\
&= \max\big\{| \phi_0|, | \beta|, \frac{1}{b^*}\| f\|\big\}\pm \phi_0
\geq 0,\\
\psi^{\pm}(1)
&= \max\big\{| \phi_0|, | \beta|,\frac{1}{b^*}\| f\|\big\}\pm y(1)\\
&= \max\big\{| \phi_0|, | \beta|, \frac{1}{b^*}\| f\|\big\}\pm \beta
\geq 0.
\end{align*}
Now we have for $x\in(0,1)$
\begin{align*}
L^{c}_{\epsilon}\psi^{\pm}(x)
&= \epsilon \frac{d^2}{dx^2}\psi^{\pm}(x)+a(x)\frac{d}{dx}\psi^{\pm}(x)-b(x)\psi^{\pm}(x)\\
&= \pm L^{c}_{\epsilon} y(x)-b(x)\max\big\{| \phi_0|, |\beta|, \frac{1}{b^*}
\| f\|\big\}\\
&\leq \pm f-b^* \max\big\{| \phi_0|, | \beta|,\frac{1}{b^*}
\| f\|\big\}
\leq 0.
\end{align*}
Using the minimum principle we obtain the required estimate.
\end{proof}
\begin{lemma}\label{l3}
The derivatives of the solution $y(x)$ of the problem \eqref{s5}-\eqref{s6}
satisfies
\[
\|y^{(k)}\|\leq C\epsilon^{-k},\quad k=1,2,3,
\]
where $C$ is a positive constant independent of $\epsilon$.
\end{lemma}
\begin{proof}
Let $x\in (0,1)$ and let $V=(c,c+\gamma)$ be a neighborhood of $x$,
where $\gamma >\epsilon$ is a constant, so that $V\subset (0,1)$,
then by mean value theorem there exists a point $\zeta\in V$ such that
\[
y'(\zeta)=\frac{y(c+\epsilon)-y(c)}{\epsilon},
\]
so
\begin{equation}
\label{s7}\epsilon \| y'(\zeta)\|\leq 2\|y\|.
\end{equation}
Now differentiating \eqref{s5} from $\zeta$ to $x$ and taking the
modulus from both sides, we obtain
\begin{equation}
\label{s8}\epsilon | y'(x)|\leq \epsilon| y'(\zeta)|+\|f\|| x-\zeta|
+\int_{\zeta}^{x}|a(t)y'(t)|\,dt+\| b\|\| y\| |x-\zeta|.
\end{equation}
Now since
\begin{equation}
\label{s9}\int_{\zeta}^{x}|a(t)y'(t)|\,dt\leq (2\|a\|+\|a'\||x-\zeta|)\| y\|,
\end{equation}
with inequalities \eqref{s7}, \eqref{s9} and using the fact $|x-\zeta|<\gamma$
and Lemma \ref{l2}, from \eqref{s8} we obtain
\[
\| y'\| \leq C\epsilon^{-1},
\]
where $C$ is a positive constant independent of $\epsilon$.
The bounds for $\|y''\|$ and $\|y'''\|$ can be obtained similarly.
\end{proof}
The solution $y(x)$ of the problem \eqref{s5}-\eqref{s6} can be decomposed
into a smooth and singular components as
\[
y=u+v,
\]
where $u$ and $v$ are smooth and singular components respectively.
The smooth component $u$ can be written in three term asymptotic expansion as
\[
u(x)=u_0(x)+u_{1}(x)\epsilon+u_{2}(x)\epsilon^{2},
\]
where $u_0, u_{1}$ and $u_{2}$ satisfies
\begin{gather*}
a(x)u'_0(x)+b(x)u_0(x)=f(x),\quad u_0(1)=u(1),\\
a(x)u'_{1}(x)+b(x)u_{1}(x)=-\epsilon u''_0(x),\quad u_{1}(1)=0,\\
u''_{2}(x)=0,\quad u_{2}(0)=0,\; u_{2}(1)=0.
\end{gather*}
The smooth component $u$ is the solution of
\begin{equation}
\label{s10} L^{c}_{\epsilon}u(x)=f(x),\quad u(0)=u_0(0)+\epsilon u_{1}(0),\quad
u(1)=y(1),
\end{equation}
and the singular component $v$ is the solution of the homogeneous problem
\begin{equation}
\label{s11} L^{c}_{\epsilon}v(x)=0,\quad v(0)=y(0)-u(0),\quad v(1)=0.
\end{equation}
Now we can state the following theorem on the bounds for the solutions
and derivatives of \eqref{s10} and \eqref{s11}
\begin{theorem}
The solutions and the derivatives of \eqref{s10} and \eqref{s11} satisfy
the following estimates
\begin{gather*}
\|u\| \leq C(1+\exp(-a^* x/\epsilon)),\\
\|v\| \leq C\exp(-a^* x/\epsilon),\\
\|u^{(k)}\| \leq C(1+\epsilon^{2-k}\exp(-a^* x/\epsilon)),\\
\|v^{(k)}\| \leq C\epsilon^{-k}\exp(-a^* x/\epsilon).
\end{gather*}
\end{theorem}
The proof of the above theorem can be found in \cite{miller}.
\section{The discrete problem}
Let $\Omega^N=\{x_0,x_{1},x_{2},\ldots,x_{N}\}$ be the partition of $[0,1]$
such that $x_0=0$, $x_i=\sum_{k=0}^{i-1}h_{k}$,
$i=1(1)N,\,h_{k}=x_{k+1}-x_{k},\,x_{N}=1$.
Let $r=\frac{h_i}{h_{i-1}},\,i=1(1)N$ be the common mesh ratio.
Taking the Taylor's series expansion and neglecting the term of third
and higher orders, we obtain the following expansions for $y_{i+1}$ and
$y_{i-1}$,
\begin{gather}
\label{d1} y_{i+1}\simeq y_i+h_iy_i'+\frac{h_i^{2}}{2}y_i'',\\
\label{d2} y_{i-1}\simeq y_i-h_{i-1}y_i'+\frac{h_{i-1}^{2}}{2}y_i''.
\end{gather}
If the boundary layer occurs at the left end then we choose $r>1$,
this gives more mesh points near the left boundary layer and if the
boundary layer occurs at the right end then we choose $r<1$,
this gives more mesh points near the right boundary layer.
Multiplying \eqref{d1} by $r$ and adding it to \eqref{d2}, we obtain
the approximation
\begin{equation}
\label{d3}y_i''\simeq\frac{2r}{h_i^{2}(1+r)}\left[y_{i+1}-(1+r)y_i+ry_{i-1}\right].
\end{equation}
Similarly we can get the two terms expression for $y_i'$ as
\begin{equation}
\label{d4} y'_i\simeq \frac{y_{i+1}-y_i}{rh_i}.
\end{equation}
Here we can use the central difference formula also in place of
forward difference formula. With the help of \eqref{d3} and \eqref{d4},
Equation \eqref{s5} can be discretized as
\begin{equation}
\label{d5} E_iy_{i-1}-F_iy_i+G_iy_{i+1}=H_i,
\end{equation}
with
\begin{equation}
\label{d6}y_0=\phi_0,\quad y_{N}=\beta,
\end{equation}
where
\begin{equation} \label{deq}
\begin{gathered}
E_i = \frac{2\epsilon r^2}{(1+r)h_i^2},\\
F_i = \frac{2\epsilon r}{h_i^2}+\frac{a_i}{rh_i}+b_i,\\
G_i = \frac{2\epsilon r}{(1+r)h_i^2}+\frac{a_i}{rh_i},\\
H_i = f_i,\quad i=1,2,\ldots,N-1.
\end{gathered}
\end{equation}
This can be written in matrix form as
\begin{equation}
\label{d7}AY=B,
\end{equation}
where $A=[c_{i,j}]$ is a tridiagonal matrix of order $N-1$ with entries
\begin{gather*}
c_{i,i-1} = E_i,\quad i=2(1)N,\\
c_{i,i} = -F_i,\quad i=1(1)N-1,\\
c_{i,i+1} = G_i,\quad i=1(1)N-1,
\end{gather*}
$Y=(y_{1},y_{2},\ldots y_{N-1})'$ and $B=(f_{1},f_{2},\ldots f_{N-1})'$
are the column vectors. Equation \eqref{d7} represents a system of linear
equations with $N-1$ equations in $N-1$ unknowns,
$y_{1},y_{2},\ldots y_{N-1}$. The system of equations can be easily solved
by discrete invariant imbedding algorithm given in \cite{angel}.
Note that one can use any other algorithm also such as Thomas algorithm.
The discrete problem \eqref{d5} satisfies the following discrete
minimum principle.
\begin{lemma}[Discrete minimum principle]\label{ld1}
Let $\psi_i$ be any mesh function such that $\psi_0,\,\psi_N\geq 0$,
then $L_{\epsilon}^d\psi_i\leq 0$ for $1\leq i \leq N-1$ implies that
$\psi_i\geq 0$ for all $0\leq i\leq N$.
\end{lemma}
\begin{proof}
Suppose there is a positive integer $k$ such that
$0>\psi_{k}=\min_{1\leq i\leq N-1}\psi_i$. Then we have
\begin{align*}
L_{\epsilon}^d\psi_{k}
&= E_k\psi_{k-1}-F_k\psi_{k}+G_{k}\psi_{k+1}\\
&= \frac{2\epsilon r^2}{(1+r)h_{k}^2}(\psi_{k-1}-\psi_k)
+\frac{2\epsilon r}{(1+r)h_k^2}(\psi_{k+1}-\psi_k)\\
&\quad +\frac{a_k}{h_{k+1}}(\psi_{k+1}-\psi_k)-b_k\psi_k
> 0,
\end{align*}
which contradict the hypothesis and hence $\psi_i\geq 0$ for all $0\leq i\leq N$.
\end{proof}
The existence, uniqueness and the stability of the solution of
problem \eqref{d5}-\eqref{d6} are given by the following theorem.
\begin{theorem}
The solution of the discrete problem \eqref{d5} together with the boundary
condition \eqref{d6} exists, is unique and satisfies
\[
|y_i| \leq C\max\{| y(0)|, | y(1)|, \frac{1}{b^*}
\max_{1\leq j\leq N-1} | L_{\epsilon}^d y_{j}|\}.
\]
\end{theorem}
\begin{proof}
Let $\psi_i$ be any mesh function satisfying $L_{\epsilon}^d(\psi_i)=f_i$.
Taking absolute value on both sides, using \eqref{d5}, we obtain
\[
|F_i||\psi_i|\leq |H_i|+|E_i||\psi_{i-1}|+|G_i||\psi_{i+1}|,\quad
i=1,2,\ldots,N-1.
\]
This gives
\begin{equation} \label{c2}
\begin{aligned}
&\frac{2\epsilon r^2}{(1+r)h_i^2}(|\psi_{i-1}|-|\psi_i|)
+\frac{2\epsilon r}{h_i^2(1+r)}(|\psi_{i+1}|-|\psi_i|)\\
&+\frac{a_i}{h_{i+1}}(|\psi_{i+1}|-|\psi_i|)-b_i|\psi_i|+|H_i|\geq 0.
\end{aligned}
\end{equation}
Let $u_i,\,v_i$ be two solutions of the difference equation \eqref{d5}
satisfying the boundary condition \eqref{d6}.
Then $w_i=u_i-v_i$ satisfies $L_{\epsilon}^d(w_i)=0$, with $w_0=w_{N}=0$.
Let $k$ be the integer such that $w_k=\max_{1\leq i\leq N-1}w_i$,
then from \eqref{c2} we have
\begin{equation}\label{c3}
\begin{aligned}
&\frac{2\epsilon r^2}{(1+r)h_{k}^2}(|w_{k-1}|-|w_k|)
+\frac{2\epsilon r}{h_k^2(1+r)}(|w_{k+1}|-|w_k|)\\
&+\frac{a_k}{h_{k+1}}(|w_{k+1}|-|w_k|)-b_k|w_k|\geq 0.
\end{aligned}
\end{equation}
Since $a_k> 0$, $b_{k}>0$, so the inequality \eqref{c3} gives $w_{k}=0$
and so $w_i\leq 0$ for $i=1,2,\ldots,N-1$.
Hence $u_i\leq v_i$ for $i=1,2,\ldots,N-1$.
Again if we set $z_i=v_i-u_i$, then $z_i$ is a mesh function satisfying
$z_0=z_N=0$. Continuing in the same way as above, we obtain
$u_i\geq v_i$ for $i=1,2,\ldots,N-1$.
Hence $u_i=v_i$ for $i=1,2,\ldots,N-1$, which shows the uniqueness
of the solution.
Now we define two mesh functions $\varphi_i^{\pm}$ such that
$$
\varphi_i^{\pm}=\max\{|y(0)|, | y(1)|,\max_{1\leq j\leq N-1}
| L_{\epsilon}^d y_{j}|\}\pm y_i.
$$
Then $\varphi_0^{\pm}\geq 0$ and $\varphi_{N}^{\pm}\geq 0$ and
for $1\leq i\leq N-1$ we have
\begin{align*}
L_{\epsilon}^d\varphi_i^{\pm}
&= -b_i\max\{| y(0)|, | y(1)|,
\frac{1}{b^*}\max_{1\leq j\leq N-1} | L_{\epsilon}^d y_{j}|\}
\pm L_{\epsilon}^d y_i \\
&\leq -b^*\max\{| y(0)|, | y(1)|,
\frac{1}{b^*}\max_{1\leq j\leq N-1} | L_{\epsilon}^d y_{j}|\}\pm
L_{\epsilon}^d y_i < 0.
\end{align*}
A consequence of Lemma \ref{ld1} gives the required estimate.
\end{proof}
\section{Stability Analysis}
Consider a difference relation
\begin{equation}
\label{sb1} y_i=S_iy_{i+1}+T_i,
\end{equation}
where $S_i=S(x_i)$ and $T_i=T(x_i)$ are unknowns which are
to be determined. From \eqref{sb1}, we have
\begin{equation}
\label{sb2}y_{i-1}=S_{i-1}y_i+T_{i-1}.
\end{equation}
Using \eqref{sb2} in \eqref{d5}, we obtain
\begin{equation}
\label{sb3}y_i=\frac{G_i}{F_i-E_iS_{i-1}}y_{i+1}
+\frac{E_iT_{i-1}-H_i}{F_i-E_iS_{i-1}}.
\end{equation}
On comparing \eqref{sb1} and \eqref{sb3}, we obtain the
recurrence relations
\begin{gather}
\label{sb4}S_i = \frac{G_i}{F_i-E_iS_{i-1}},\\
\label{sb5}T_i = \frac{E_iT_{i-1}-H_{i-1}}{F_i-E_iS_{i-1}}.
\end{gather}
To solve above recurrence relations for $i=1,2,\ldots N-1$, we need
$S_0$ and $T_0$.
Now it is given that $y_0=\phi_0$, therefore we have $S_0y_{1}+T_0=\phi_0$.
So we can choose $S_0=0$ and then $T_0=\phi_0$. Now by using
these initial conditions, we can compute $S_i$ and $T_i$ for
$i=1,2,\ldots,N-1$ and using these values of $S_i$ and $T_i$ in
\eqref{sb1}, we obtain $y_i$ for $i=1,2,\ldots,N-1$.
Now we give the proof of the stability of our scheme.
Suppose a small error $e_{i-1}$ has been made in the calculation
of $S_{i-1}$, then we have, $\bar{S}_{i-1}=S_{i-1}+e_{i-1}$, and we
are actually calculating
\begin{equation}
\label{sb6} \bar{S}_i=\frac{G_i}{F_i-E_i\bar{S}_{i-1}}.
\end{equation}
From \eqref{sb4} and \eqref{sb6}, we have
\begin{equation}\label{sb7}
\begin{aligned}
e_i&=\frac{G_i}{F_i-E_i(S_{i-1}+e_{i-1})}-\frac{G_i}{F_i-E_iS_{i-1}}\\
&=\frac{G_iE_ie_{i-1}}{(F_i-E_i(S_{i-1}+e_{i-1}))(F_i-E_iS_{i-1})}.
\end{aligned}
\end{equation}
Under the assumption, the error is small initially, from \eqref{sb7} we obtain
\begin{equation}
\label{sb8} e_i=\Big(\frac{S_i^2E_i}{G_i}\Big)e_{i-1}.
\end{equation}
Now we have
\[
G_{1}-F_{1}= -\frac{2\epsilon r^2}{(1+r)h_{1}^2}-b_{1}<0,
\]
so $\frac{G_{1}}{F_{1}}<1$. Therefore from \eqref{sb4}, we have
$S_{1}=\frac{G_{1}}{F_{1}} <1$.
Again from \eqref{sb4} we have
\[
S_{2}=\frac{G_{2}}{F_{2}-E_{2}S_{1}}<\frac{G_{2}}{F_{2}-E_{2}}
<\frac{G_{2}}{E_{2}+G_{2}-E_{2}}=1.
\]
Similarly we can show that $S_i<1$ for $1\leq i\leq N-1$. Now we have
\begin{align*}
|E_i|-|G_i|
&= \frac{2\epsilon r^2}{(1+r)h_i^2}-\frac{2\epsilon r}{(1+r)h_i^2}
-\frac{a_i}{rh_i},\\
&= \frac{2\epsilon r(r-1)}{(1+r)h_i^2}-\frac{a_i}{rh_i}\\
&< 0\quad \textrm{as $\epsilon$ is small and $r$ is close to $1$}.
\end{align*}
Thus $\frac{|E_i|}{|G_i|}<1$, it follows from \eqref{sb8} that
\[
|e_i|=|S_i|^2\big| \frac{E_i}{G_i}\big| |e_{i-1}|<|e_{i-1}|.
\]
Which confirm the stability of the recurrence relation \eqref{sb4}. Similarly
we can prove that the recurrence relation \eqref{sb5} is also stable.
\section{Numerical results and discussions}
To validate the theoretical results, we apply the proposed numerical
scheme to a test problem having a left boundary layer.
\begin{example}\label{ex:1} \rm
Consider the problem $\epsilon y''(x)+y'(x-\delta)-y(x)=0$,
under the interval conditions
$y(x)=1$, $-\delta \leq x \leq 0,\ y(1)=1$.
\end{example}
\begin{table}[ht]
\caption{Maximum absolute error for Example \ref{ex:1}
for $\delta=0.001\times\epsilon$ using $N=100$}\label{tab:1}
\begin{tabular*}{1.0\textwidth}{@{\extracolsep{\fill}}cccc}\hline
$\epsilon$ & $r=1.1$ & $r=1.01$ & $r=1.001$\\\hline
$10^{-2}$ & 2.42E-02 & 5.65E-02 & 7.96E-02\\
$10^{-4}$ & 2.95E-02 & 9.05E-03 & 7.97E-03\\
$10^{-6}$ & 6.39E-02 & 1.65E-03 & 1.53E-03\\
$10^{-8}$ & 2.57E-02 & 1.66E-03 & 1.47E-03\\
$10^{-10}$ & 2.57E-02 & 1.66E-03 & 1.47E-03\\
$10^{-12}$ & 2.57E-02 & 1.66E-03 & 1.46E-03\\\hline
\end{tabular*}
\end{table}
A numerical method for solving singularly perturbed boundary
value problem with a negative shift in the first derivative
term is considered. It is a practical method and can be easily
implemented on a computer to solve such problems. An example is
given to demonstrate the efficiency of the proposed method.
The maximum absolute errors $E_{\epsilon}^N=\max_i| y(x_i)-y_i|$
at the nodal points are tabulated in the table for different values of
perturbation parameter $\epsilon$ and different values of mesh ratio
$r$ by using $N=100$.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.9\textwidth]{fig1}
\end{center}
\caption{Numerical solution of Example \ref{ex:1} for
$\epsilon=0.1$}\label{fig1}
\end{figure}
The graph of the solution of the considered example for different values
of delay is plotted in Figure \ref{fig1} to examine the questions on
the effect of delay on the boundary layer behavior of the solution.
One can observe that if $\delta=o(\epsilon)$, the layer behavior is
maintained in the case of left boundary layer (the layer behavior is also
maintained in the case of right boundary layer). As the delay increases,
the thickness of the layer decreases in the case when the solution exhibits
layer behavior on the left side (as shown in Figure \ref{fig1}) while in the
case of the right side boundary layer, it increases.
The delay affects more in the case when the solution of the boundary
value problem exhibits layer behavior on the left side in comparison
to the right side boundary layer case.
\subsection*{Acknowledgements}
Author is very thankful to the reviewer for his/her valuable suggestions
and comments for the improvement of this paper.
\begin{thebibliography}{20}
\bibitem{angel} E. Angel, R. Bellman;
Dynamic Programming and Partial Differential Equations,
Academic, New York, (1972).
\bibitem{best} M. Bestehorn, E. V. Grigorieva;
Formation and propagation of localized states in extended systems,
\emph{Ann. Phys.}, (Leipzig) Vol. 13
(2004) 423-431.
\bibitem{burton} T. A. Burton;
Fixed points, stability, and exact linearization,
\emph{Nonlinear Anal.}, Vol. 61 (2005) 857-870.
\bibitem{waz} M. Wazewska-Czyzewska, A. Lasota;
Mathematical models of the red cell system, \emph{Mat. Stos.},
Vol. 6 (1976) 25-40.
\bibitem{doolan} E. P. Doolan, J. J. H. Miller, W. H. A. Schilders;
Uniform Numerical Methods for Problems with Initial and Boundary Layers,
Boole Press, Dublin, (1980).
\bibitem{ezzat} M. A. Ezzat, M. I. Othman, A. M. S. El-Karamany;
State space approach to two-dimensional generalized thermo-viscoelasticity
with two relaxation times, \emph{Int. J. Eng. Sci.}, Vol. 40 (2002) 1251-1274.
\bibitem{farrell} P. A. Farrell, A. F. Hegarty, J. J. H. Miller, R. E. O’Riordan,
G. I. Shishkin;
Robust Computational Techniques for Boundary Layers,
Chapman- Hall/CRC, New York, (2000).
\bibitem{joseph1} D. D. Joseph, L. Preziosi; Heat waves, \emph{Rev. Mod. Phys.}, Vol. 61 (1989) 41-73.
\bibitem{joseph2} D. D. Joseph, L. Preziosi;
Addendum to the paper ``Heat waves'', \emph{Rev. Mod. Phys.},
Vol. 62 (1990) 375-391.
\bibitem{kadal} M. K. Kadalbajoo, D. Kumar;
Fitted mesh B-spline collocation method for singularly perturbed
differential-difference equations with small delay,
\emph{Appl. Math. Comput.}, Vol. 204 (2008) 90-98.
\bibitem{liao} X. Liao;
Hopf and resonant codimension two bifurcation in van der Pol equation
with two time delays, \emph{Chaos, Solitons \& Fractals},
Vol. 23 (2005) 857-871.
\bibitem{liu} Q. Liu, X. Wang, D. De Kee;
Mass transport through swelling membranes, \emph{Int. J. Eng. Sci.},
Vol. 43 (2005) 1464-1470.
\bibitem{mack} M. C. Mackey, L. Glass;
Oscillations and chaos in physiological control systems, \emph{Science},
Vol. 197 (1977) 287-289.
\bibitem{miller} J. J. H. Miller, E. O'Riordan, G. I. Shishkin;
Fitted numerical methods for singular perturbation problems, World
Scientific, (1996).
\bibitem{ramos} J. I. Ramos;
On the numerical treatment of an ordinary differential equation arising
in one-dimensional non-Fickian diffusion problems,
\emph{Comput. Phys. Commun.}, Vol. 170 (2005) 231-238
\bibitem{roos} H. G. Roos, M. Stynes, L. Tobiska;
Numerical Methods for Singularly Perturbed Differential Equations,
Convection-Diffusion and Flow Problems, Springer Verlag, Berlin, (1996).
\bibitem{tzou} D. Y. Tzou;
Micro-to-macroscale heat transfer, Taylor \& Francis, Washington, DC, 1997.
\end{thebibliography}
\end{document}