This work is licensed under the Creative Commons Attribution 4.0 International License.
Introduction
Consider an initial value problem for the semilinear second order singularly perturbed delay differential equation [12] in the interval Ī = [0,T]:
\varepsilon {u^{\prime\prime}}(t) + a(t){u^\prime}(t) + f(t,u(t),u(t - r)) = 0,\quad t \in I,u(t) = \varphi (t),\quad t \in {I_0},{u^\prime}(0) = A/\varepsilon ,
where
I = (0,T] = \mathop \cup \limits_{p = 1}^m {I_p}
, Ip = {t : rp−1 <t ≤ rp}, 1 ≤ p ≤ m and rs = sr, for 0 ≤ s ≤ m and I0 = (−r,0]. 0 < ɛ ≤ 1 is the perturbation parameter, a(t) ≥ α > 0, b(t), c(t), f (t) and ϕ(t) are given sufficiently smooth functions satisfying certain regularity conditions to be specified and r is a constant delay, which is independent of ɛ. Moreover
\left| {{{\partial f} \over {\partial u\;}}} \right| \le {b^ *}\,{\rm{and}}\;\left| {{{\partial f} \over {\partial v\;}}} \right| \le {c^ * }.
Delay differential equations play an important role in the mathematical modelling of various practical phenomena in the biosciences and control theory. Any system involving a feedback control will almost always involve time delays. These arise because a finite time is required to sense information and then react to it. A singularly perturbed delay differential equation is an ordinary differential equation in which the highest derivative is multiplied by a small parameter and involving at least one delay term [1,2,3,4]. Such problems arise frequently in the mathematical modelling of various practical phenomena, for example, in the modelling of several physical and biological phenomena like the optically bistable devices [5], description of the human pupil-light reflex [6], a variety of models for physiological processes or diseases and variational problems in control theory where they provide the best and in many cases the only realistic simulation of the observed phenomena [7]. An overview of numerical treatment for first and second order singularly perturbed delay differential equations, may be obtained in [8,9,10,11,12,13,14,15,16] (see, also references therein).
The numerical analysis of singular perturbation cases has always been far from trivial because of the boundary layer behavior of the solution. Such problems undergo rapid changes within very thin layers near the boundary or inside the problem domain. It is well known that standard numerical methods for solving singular perturbation problems do not give satisfactory result when the perturbation parameter is sufficiently small. Therefore, it is important to construct suitable numerical methods for these problems, whose accuracy does not depend on the perturbation parameter, i.e. methods that are uniformly convergent with respect to the perturbation parameter [17,18,19,20,21,22].
In a singularly perturbed delay differential equation, one encounters with two difficulties, one is because of occurrence of the delay term and another one is due to presence of perturbation parameter. To overcome the first difficulty, we employed the numerical method of steps [2] for the delay argument which converted the problem to a initial value problem for a singularly perturbed differential equation. Then we constructed a numerical scheme based on finite difference method on non uniform Shishkin mesh for the numerical solution.
In the present paper we discretize (1)–(2) using a numerical method, which is composed of an exponentially fitted difference scheme on piecewise uniform Shishkin mesh on each time-subinterval. In section 2, we state some important properties of the exact solution. In section 3, we describe the finite difference discretization and introduce the piecewise uniform mesh. In section 4, we present convergence analysis for approximate solution. Uniform convergence is proved in the discrete maximum norm. Some numerical results are being presented in section 5. The technique to construct discrete problem and error analysis for approximate solution is similar to those in [8,9,23,24]
Throughout the paper, C will denote a generic positive constant independent of ɛ and of the mesh parameter.
The Continuous Problem
Here we show some properties of the solution of (1)–(3), which are needed in later sections for the analysis of the appropriate numerical solution. For any continuous function g(t), we use ‖g‖∞ for the continuous maximum norm on the corresponding interval.
Lemma 2.1
Let δ (t) be nonnegative and continuous function such that
\delta (t) \le {\delta _0} + \int\limits_0^t \left\{ {{c_0}\delta (\tau ) + {c_1}\delta (\tau - r)} \right\}d\tau ,{\kern 1pt} {\kern 1pt} t > 0,\delta (t) = {\delta _0},\;{\kern 1pt} - r \le t \le 0,
where δ0, c0 and c1 are nonnegative constants. Then it holds that
\delta (t) \le {\gamma _0}\sum\limits_{s = 0}^p \left\{ {{{\gamma _1^s{{(t - {r_{p - 1}})}^s}} \over {s!}}} \right\},{\kern 1pt} \;{\kern 1pt} t \in {I_p}{\kern 1pt} \;{\kern 1pt} ,\,p = 1,2,...,m,
where γ0 = δ0exp(c0T) and γ1 = c1exp(c0T).
Now, differentiating the Eq.(1) we have
\varepsilon {u^{\prime\prime\prime}}(t) + a(t){u^{\prime\prime}}(t) = \Phi (t),\quad t \in {I_{k + 1}},
where
\Phi (t) = {f^\prime}(t,0,0) - ({a^\prime}(t) + b(t))(t){u^\prime}(t) - {b^\prime}(t)u(t) - {c^\prime}(t)u(t - r) - c(t){u^\prime}(t - r).
From (14) we have the following relation for u″(t)
{u^{\prime\prime}}(t) = {u^{\prime\prime}}({r_k})\exp \left( {\; - {1 \over \varepsilon }\int\limits_{{r_k}}^t a(s)ds} \right) + {1 \over \varepsilon }\int\limits_{{r_k}}^t \Phi (\tau )\exp \left( {\; - {1 \over \varepsilon }\int\limits_\tau ^t a(s)ds} \right)d\tau .
and it is easy to see that
\matrix{ {\left| {\Phi (t)} \right| \le \left| {{f^\prime}(t)} \right| + \left| {({a^\prime}(t) + b(t)){u^\prime}(t)} \right| + \left| {{b^\prime}(t)u(t)} \right| + \left| {{c^\prime}(t)u(t - r)} \right| + \left| {c(t){u^\prime}(t - r)} \right|} \cr { \le C\left( {1 + \left| {({a^\prime}(t) + b(t)){u^\prime}(t)} \right| + \left| {c(t){u^\prime}(t - r)} \right|} \right).}}
It follows from (1.1) that
\left| {{u^{\prime\prime}}(0)} \right| \le {C \over {{\varepsilon ^2}}}.
Using the estimates (8), (16) and (17) in (21) for k = 0, we have
\matrix{ {\left| {{u^{\prime\prime}}(t)} \right| \le {C \over {{\varepsilon ^2}}}\exp \left( {{{ - \alpha t} \over \varepsilon }} \right) + {1 \over \varepsilon }C\int\limits_0^t \left( {1 + {1 \over \varepsilon }exp \left( {{{ - \alpha \tau } \over \varepsilon }} \right)} \right)\exp \left( {{{ - \alpha (t - \tau )} \over \varepsilon }} \right)d\tau } \cr { \le C\left( {1 + {1 \over {{\varepsilon ^2}}}exp \left( {{{ - \alpha t} \over \varepsilon }} \right)} \right),{\kern 1pt} \;{\kern 1pt} t \in {I_1},}}
i.e., the inequality (9) is valid for p = 1.
Next, from the last inequality for t = r, we have
\left| {{u^{\prime\prime}}(r)} \right| \le C
and thereby from (15) it follows that
\matrix{ {\left| {{u^{\prime\prime}}(t)} \right| \le |{u^{\prime\prime}}(r)|\exp \left( {{{ - \alpha t} \over \varepsilon }} \right)} \cr { + {1 \over \varepsilon }C\int\limits_0^t \left( {1 + {1 \over \varepsilon }exp \left( {\;{{ - \alpha (\tau - r)} \over \varepsilon }} \right)} \right)\exp \left( {{{ - \alpha (t - \tau )} \over \varepsilon }} \right)d\tau } \cr { \le C\left( {1 + {1 \over {{\varepsilon ^2}}}(t - r)exp \left( {{{ - \alpha t} \over \varepsilon }} \right)} \right),{\kern 1pt} \;{\kern 1pt} t \in {I_2}.}}
i.e. the inequality (9) is also valid for p = 2.
Since
\left| {{\Phi ^\prime}(t)} \right| \le C,{\kern 1pt} \;{\kern 1pt} t \in {I_p}{\kern 1pt} \;{\kern 1pt} ,{\kern 1pt} \;{\kern 1pt} p = 3,...,m,
the estimate (10) follows immediately from (15).
Discretization and Mesh
In this section, we construct a numerical scheme for solving the initial value problem (1)–(2).
Let
{\bar \omega _{{N_0}}}
be any non-uniform mesh on Ī:
{\bar \omega _{{N_0}}} = \left\{ {0 = {t_0} < {t_1} < ... < {t_{{N_0}}} = T,{\kern 1pt} \;{\kern 1pt} {h_i} = {t_i} - {t_{i - 1}}} \right\}
which contains by N mesh point at each subinterval Ip(1 ≤ p ≤ m) :
{\omega _{N,p}} = \left\{ {{t_i}:(p - 1)N + 1 \le i \le pN} \right\},{\kern 1pt} \;{\kern 1pt} 1 \le p \le m,
and consequently
{\omega _{{N_0}}} = \bigcup\limits_{p = 1}^m {\omega _{N,p}}.
To simplify the notation we set wi = w(ti) for any function w(t), and moreover yi denotes an approximation of u(t) at ti. For any mesh function {wi} defined on
{\bar \omega _{{N_0}}}, we use
\matrix{ {{w_{\bar t,i}} = ({w_i} - {w_{i - 1}})/{h_i},{\kern 1pt} \;{\kern 1pt} {w_{t,i}} = ({w_{i + 1}} - {w_i})/{h_{i + 1}},} \cr {{w_{\mathop t\limits^0 ,i}} = ({w_{\bar t,i}} + {w_{t,i}})/2,{\kern 1pt} \;{\kern 1pt} {w_{\bar tt,i}} = ({w_{t,i}} - {w_{\bar t,i}})/{h_i},} \cr {{{\left\| w \right\|}_{\infty ,N,p}} = {{\left\| w \right\|}_{\infty ,{\omega _{N,p}}}}: = \mathop {\max }\limits_{1 \le i \le N} \left| {{w_i}} \right|.}}
For the difference approximation of (1), we are using the following identity
h_i^{ - 1}\int\limits_{{t_{i - 1}}}^{{t_{i + 1}}} Lu(t){\psi _i}(t)dt = 0,1 \le i \le {N_0} - 1,
where exponential basis functions
{\psi _i}(t) = \left\{ {\matrix{ {{\psi _{1i}}(t) \equiv {{{e^{{a_i}(t - {t_{i - 1}})/\varepsilon }} - 1} \over {{e^{{a_i}h/\varepsilon }} - 1}},{t_{i - 1}} < t < {t_i}} \hfill \cr {{\psi _{2i}}(t) \equiv {{1 - {e^{ - {a_i}({t_{i + 1}} - t)/\varepsilon }}} \over {1 - {e^{ - {a_i}h/\varepsilon }}}},{t_i} < {\kern 1pt} t < {t_{i + 1}},} \hfill \cr {0\quad \quad \quad \quad \quad \quad \quad \quad \quad ,t \notin ({t_{i - 1}},{t_{i + 1}}),} \hfill \cr } } \right.
and
h_i^{ - 1}\int_{{t_{i - 1}}}^{{t_{i + 1}}} {\psi _i}\left( t \right)dt = 1.
Using interpolating quadrature rules with the weight and remainder terms in integral form on subintervals [ti−1, ti+1], consistent with [23,24], after a simple calculation, we have the following relation:
\varepsilon {\theta _i}{u_{\bar tt,i}} + {a_i}{u_{\mathop t\limits^0 ,i}} + f({t_i},{u_i},{u_{i - N}}) + {R_i} = 0,\quad i = 1,2,...,{N_0} - 1,
where the factor
{\theta _i} = {{\rho {a_i}} \over 2}\coth ({a_i}\rho /2){\kern 1pt} {\kern 1pt} ,{\kern 1pt} {\kern 1pt} \rho = h/\varepsilon
and remainder term
\matrix{ {{R_i} = h_i^{ - 1}\int\limits_{{t_{i - 1}}}^{{t_{i + 1}}} [a(t) - a({t_i})]{u^\prime}(t){\psi _i}(t)dt} \cr { + h_i^{ - 1}\int\limits_{{t_{i - 1}}}^{{t_{i + 1}}} {d \over {dt}}f(\tau ,u\left( \tau \right),u\left( {\tau - r} \right))K_{0,i}^ * \left( {t,\tau } \right)d\tau } \cr {K_{0,i}^ * \left( {t,\tau } \right) = {T_0}\left( {{t_i} - \tau } \right) - {T_0}\left( {{t_i} - \tau } \right),{T_0}(t) = 1,\;t > 0;{T_0}(t) = 0,\;t < 0.}}
To define an approximation for the boundary condition (3), we proceed our discretization process by
\int\limits_{{t_0}}^{{t_1}} Lu(t){\varphi _0}(t)dt = 0
where the exponential basis function
{\psi _0}(t) = \left\{ {\matrix{ {{{1 - {e^{ - {a_0}({t_1} - t)/\varepsilon }}} \over {1 - {e^{ - {a_0}h/\varepsilon }}}},} \hfill & {{t_0} < {\kern 1pt} t < {t_1}} \hfill \cr {0\quad \quad \quad \quad\,\,\,\,\,,} \hfill & {t \notin ({t_0},{t_1})} \hfill \cr } } \right.
We note that functions ϕ0(t) is the solution of the following problem:
\eqalign{& {\varepsilon {\psi ^{\prime\prime}}(t) - {a_0}{\psi ^\prime}(t) = 0{\kern 1pt} {\kern 1pt} ,{\kern 1pt} {\kern 1pt} {t_0} < t < {t_1},} \cr & {\psi ({t_0}) = 0{\kern 1pt} {\kern 1pt} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \psi ({t_1}) = 1.}}
Whence, as similar above we can write the following difference relation:
\varepsilon \sigma {u_{t,0}} - A + {r^{(0)}} = 0
where the coefficient and the remainder term
\sigma = 1 + {a_0}{\varepsilon ^{ - 1}}\int\limits_{{t_0}}^{{t_1}} {\psi _0}(t)dt = {{{a_0}\rho } \over {1 - {e^{ - {a_0}\rho }}}}{r^{(0)}} = \int\limits_{{t_0}}^{{t_1}} \left[ {a(t) - {a_0}} \right]{u^\prime}(t){\psi _0}(t)dt + \int\limits_{{t_0}}^{{t_1}} f(t,u(t),u(t - r)){\psi _0}(t)dt.
Neglecting Ri and r(0) in (19) and (22), we propose the following difference scheme for approximation (1)–(2):
\varepsilon {\theta _i}{y_{\bar tt,i}} + {a_i}{y_{\mathop t\limits^0 ,i}} + f({t_i},{y_i},{y_{i - N}}) = 0,\quad i = 1,2,...,{N_0} - 1,{y_i} = {\varphi _i},\ - N \le i \le 0,\varepsilon \sigma {y_{t,0}} - A = 0,
where θi and σ are defined by (20),(23), respectively.
The difference scheme (25)–(26), in order to be ɛ- uniform convergent, we will use the non uniform Shishkin mesh. For the even number N, the piecewise uniform mesh ωN, p divides each of the interval [rp−1, σp] and [σp, rp] into N/2 equidistant subintervals, where the transition point σp, which separates the fine and coarse portions of the mesh is obtained by
{\sigma _p} = {r_{p - 1}} + \min \left\{ {r/2,{\kern 1pt} \;{\kern 1pt} {\alpha ^{ - 1}}{\beta _p}\varepsilon \ln N} \right\},
where β1 ≥ 1 and βp > 1 (2 ≤ p ≤ m) are some constants. Hence, if
h_p^{(1)}
and
h_p^{(2)}
denote the stepsizes in [rp−1, σp] and [σp, rp] respectively, we have
h_p^{(1)} = 2({\sigma _p} - {r_{p - 1}}){N^{ - 1}},{\kern 1pt} \;{\kern 1pt} h_p^{(2)} = 2({r_p} - {\sigma _p}){N^{ - 1}},{\kern 1pt} \;{\kern 1pt} 1 \le p \le m,
so
\matrix{ {{{\bar \omega }_{N,p}} = \left\{ {\matrix{ {{t_i} = {r_{p - 1}} + (i - (p - 1)N)h_p^{(1)},{\kern 1pt} \;{\kern 1pt} i = (p - 1)N,...,(p - 1/2)N} \hfill \cr {{t_i} = {\sigma _p} + (i - (p - 1/2)N)h_p^{(2)},{\kern 1pt} \;{\kern 1pt} i = (p - 1/2)N + 1,...,pN,} \hfill \cr } } \right.} \cr {1 \le p \le m.}}
In the rest of the paper we only consider this type mesh.
Convergence Analysis
We now estimate the approximate error zi = yi − ui, 0 ≤ i ≤ N0, which satisfies the discrete problem
\varepsilon {\theta _i}{z_{\bar tt,i}} + {a_i}{z_{\mathop t\limits^0 ,i}} + f({t_i},{y_i},{y_{i - N}}) - f({t_i},{u_i},{u_{i - N}}) = {R_i},i = 1,...,{N_0} - 1,{z_i} = 0,\ - N \le i \le 0,\varepsilon \sigma {z_{t,0}} - A + {r^{(0)}} = 0,
where the truncation error Ri and r(0) are given by (21) and (24).
Lemma 4.1
Let zi be the solution of the problem (27)–(28). Then the following estimate holds
{\left\| z \right\|_{\infty ,N,p}} \le \gamma (({\theta _ * }/\sigma )|{r^{(0)}}| + C\sum\limits_{k = 1}^p {\left\| R \right\|_{\infty ,{\omega _{N,k}}}},{\kern 1pt} \;{\kern 1pt} 1 \le p \le m.
where
{\theta _ * } = {{\rho \alpha } \over 2}\coth (\alpha \rho /2)
, γ = 4α−1 exp(4α−1 (b* + c*))
Proof
Let zt, i = vi. Then the relation (27) can be rewritten as
\varepsilon {\kern 1pt} {\theta _i}{v_{\bar t,i}} + {{{a_i}} \over 2}({v_i} + {v_{i - 1}}) = {F_i},
where
{F_i} = {R_i} - f({t_i},{y_i},{y_{i - N}}) + f({t_i},{u_i},{u_{i - N}}).
Solving the first order difference equation with respect to vi, we get
{v_i} = {v_0}{Q_i} + h\sum\limits_{k = 1}^i {{{F_\ell }} \over {\varepsilon {\theta _k} + h{a_k}/2}}{Q_{i - k}},
where
{Q_{i - k}} = \left\{ {\matrix{ {1{\kern 40pt},{\kern 1pt} {\kern 1pt} k = i,} \hfill \cr {\prod\limits_{j = k + 1}^i \left( {{{1 - {a_j}\rho /(2{\theta _j})} \over {1 + {a_j}\rho /(2{\theta _j})}}} \right),{\kern 1pt} {\kern 1pt} 0 \le k \le i - 1.} \hfill \cr } } \right.
Then, since
{z_p} = h\sum\limits_{i = 0}^{p - 1} {v_i} = h\sum\limits_{i = 1}^p {v_{i - 1}}
and taking into consideration (5), from (30) after some manipulations we have the following inequality
\left| {{z_p}} \right|{\kern 1pt} \le 4{\alpha ^{ - 1}}(\varepsilon {\theta _ * }\left| {{z_{t,0}}} \right| + h\sum\limits_{i = 1}^{p - 1} (|{R_i}| + {b^ * }|{z_i}| + {c^ * }|{z_{i - N}}|)),{\kern 1pt} 1 \le p \le {N_0} - 1.
Taking also into account (28), and using difference analogue of Gronwall’s inequality this leads to (29).
Lemma 4.2
Let a,b,c, f ∈ C1 (Ī), ψ ∈ C1 (Ī0). Then for the truncation errors Ri and r(0), the following estimates hold
\matrix{ {{{\left\| R \right\|}_{\infty ,{\omega _N},p}} \le C{N^{ - 1}}\ln N,{\kern 1pt} \;{\kern 1pt} 1 \le p \le m,} \cr {|{r^{(0)}}| \le C{N^{ - 1}}\ln N.}}
Proof
From explicit expression (21) for Ri on an arbitrary mesh, we have
\matrix{ {|{R_i}| = h_i^{ - 1}\int\limits_{{t_{i - 1}}}^{{t_{i + 1}}} |(a(t) - a({t_i}))||{u^\prime}(t)|{\psi _i}(t)dt} \hfill \cr { + h_i^{ - 1}\int\limits_{{t_{i - 1}}}^{{t_{i + 1}}} \left| {{d \over {dt}}f(t,u\left( t \right),u\left( {t - r} \right))} \right|dt,{\kern 1pt} \;{\kern 1pt} 1 \le i \le {N_0}.} \hfill}
This inequality together with (7), enable us to write
\left| {{R_i}} \right| \le C\left\{ {{h_i} + \int\limits_{{t_{i - 1}}}^{{t_{i + 1}}} \left( {\left| {{u^\prime}(t)} \right| + \left| {{u^\prime}(t - r)} \right|} \right)dt} \right\},{\kern 1pt} \;{\kern 1pt} 1 \le i \le {N_0}.
From here, in view of (8), it follows that
\left| {{R_i}} \right| \le C\left\{ {{h_i} + {1 \over \varepsilon }\int\limits_{{t_{i - 1}}}^{{t_{i + 1}}} {e^{ - {{\alpha t} \over \varepsilon }}}dt} \right\},{\kern 1pt} \;{\kern 1pt} 1 \le i \le {N_0}
in which
{h_i} = \left\{ {\matrix{ {h_p^{(1)},{\kern 1pt} \;{\kern 1pt} (p - 1)N \le i \le (p - 1/2)N} \cr {h_p^{(2)},{\kern 1pt} \;{\kern 1pt} (p - 1/2)N + 1 \le i \le pN.} \cr } } \right.
At the each submesh ωN, p, we estimate the truncation error R as follows. We consider first the case σp = rp−1 + r/2, and so r/2 ≤ α−1θpɛ lnN and
h_p^{(1)} = h_p^{(2)} = {\tau _p} = r/N
.
Hereby, since
\matrix{ {h_i^{ - 1}{\varepsilon ^{ - 1}}\int\limits_{{t_{i - 1}}}^{{t_{i + 1}}} {e^{ - {{\alpha t} \over \varepsilon }}}dt \le {\varepsilon ^{ - 1}}{h_p} \le {{2{\beta _p}\ln N} \over {\alpha r}}{r \over N} = 2{\beta ^{ - 1}}{\beta _p}{N^{ - 1}}\ln N,} \cr {(p - 1)N \le i \le pN,{\kern 1pt} \;{\kern 1pt} 1 \le p \le m - 1}}
it follows from (32) that
\left| {{R_i}} \right| \le C{N^{ - 1}}\ln N,{\kern 1pt} \;{\kern 1pt} (p - 1)N \le i \le pN,{\kern 1pt} \;{\kern 1pt} 1 \le p \le m - 1.
We now consider the case σp = rp−1 + α−1βpɛ lnN and estimate Ri on [rp−1, σp] and [σp, rp] separately. In the layer region [rp−1, σp], the inequality (32) reduces to
\matrix{ {\left| {{R_i}} \right| \le C(1 + {\varepsilon ^{ - 1}})h{p^{(1)}} = C(1 + {\varepsilon ^{ - 1}}){{{\alpha ^{ - 1}}{\beta _p}\varepsilon \ln N} \over {N/2}},(p - 1)N \le i \le (p - 1/2)N,} \cr {1 \le p \le m - 1.}}
Next we estimate Ri for (p−1/2)N +1 ≤ i ≤ pN. In this case, recalling that
{t_i} = {\sigma _p} + (i - (p - 1/2)N)h_p^{(2)} = {r_{p - 1}} + {\alpha ^{ - 1}}{\beta _p}\varepsilon \ln N + (i - (p - 1/2)N)h_p^{(2)}
, we obtain from (32)\left| {{R_i}} \right| \le C\left\{ {h_p^{(2)} + {\alpha ^{ - 1}}{\beta _p}({e^{ - {{\alpha ({t_{i - 1}})} \over {{\beta _p}\varepsilon }}}} - {e^{ - {{\alpha ({t_i})} \over {{\beta _p}\varepsilon }}}})} \right\}
and this implies that
\left| {{R_i}} \right| \le C{N^{ - 1}}.
Then, we estimate the remainder term r(0). From the explicit expression (32) and |ϕ0 (t)| ≤ 1, after similar calculation as above, it follows that
\left| {{r^{\left( 0 \right)}}} \right| \le C{N^{ - 1}}\ln N.
Combining the two previous lemmas gives us the following main convergence result.
Theorem 4.1
Let a ∈ C1 (Ī), ϕ ∈ C1 (I0). The continuously differentiable function f (t,u,v) satisfies the conditions (1.4) and the derivative
{\partial \over {\partial t}}f(t,u,v)
is bounded for 0 ≤ t ≤ T and |u|, |v| ≤ C. Then the following estimate holds
\left| {{y_i} - {u_i}} \right| \le C{N^{ - 1}}\ln N,\quad 0 \le i \le {N_0}.
Numerical Results
In this section, a simple numerical example is devised to verify the validity for the proposed method. Consider the test problem with
a = 1,f = - u(t - 1) + {t^2},T = 2,\varphi = 1 + t,A = - 1.
We use the double mesh principle to estimate the errors and compute the experimental rates of convergence in our computed solution, i.e. We compare the computed solution with the solution on a mesh that is twice as fine [17]. The error estimate
e_\varepsilon ^{N,p}
and the computed convergence rate
r_\varepsilon ^{N,p}
obtained in this way are denoted by
e_\varepsilon ^{N,p} = \mathop {\max }\limits_{{\omega _{N,p}}} \left| {{y^{\varepsilon ,N}} - {y^{\varepsilon ,2N}}} \right|,p = 1,2;r_\varepsilon ^{N,p} = \ln \left( {e_\varepsilon ^{N,p}/e_\varepsilon ^{2N,p}} \right)/\ln 2.
The resulting errors
e_\varepsilon ^{N,p}
and the corresponding numbers
r_\varepsilon ^{N,p}
for particular values of ɛ, N, for first and second subinterval are listed in the Tables 1 and 2.
Maximum errors and rates of convergence on ωN,1
ɛ
N = 64
N = 128
N = 256
N = 512
N = 1024
2−1
0.0085116
0.0042832
0.0021485
0.0010760
0.0005384
0.99
1.00
1.00
1.00
2−2
0.0140007
0.0070930
0.0035693
0.0017090
0.0008906
0.98
0.98
0.99
1.00
2−4
0.0241181
0.0143603
0.0083166
0.0047146
0.0026331
0.75
0.78
0.81
0.84
2−6
0.0230541
0.0137267
0.0079497
0.0045066
0.0025149
0.75
0.78
0.81
0.84
2−8
0.0227881
0.0135684
0.0078580
0.0044546
0.0024859
0.75
0.78
0.81
0.84
2−10
0.0227216
0.0135288
0.0078350
0.0044416
0.0024786
0.75
0.78
0.81
0.84
2−12
0.0227050
0.0135189
0.0078293
0.0044382
0.0024768
0.75
0.78
0.81
0.84
2−14
0.0227008
0.0135164
0.0078279
0.0044361
0.0024763
0.75
0.78
0.81
0.84
2−16
0.0226998
0.0135158
0.00782756
0.0044374
0.0024762
0.75
0.78
0.81
0.84
Maximum errors and rates of convergence on ωN,2
ɛ
N = 64
N = 128
N = 256
N = 512
N = 1024
2−1
0.0079814
0.0040103
0.0020101
0.0010063
0.0005033
0.98
0.99
1.00
1.00
2−2
0.0012138
0.0061392
0.0030875
0.0015482
0.0007775
0.98
0.99
1.00
1.00
2−4
0.0202600
0.0120754
0.0069885
0.0039609
0.0022110
0.75
0.79
0.81
0.84
2−6
0.0194243
0.0115426
0.0066802
0.0037862
0.0021126
0.75
0.78
0.82
0.84
2−8
0.0191427
0.0114094
0.0066031
0.0037425
0.0020882
0.75
0.78
0.82
0.84
2−10
0.0190868
0.0113761
0.0065838
0.0037315
0.0020821
0.75
0.78
0.82
0.84
2−12
0.0190729
0.0113678
0.0065790
0.0037288
0.0020806
0.75
0.78
0.82
0.84
2−14
0.0190696
0.0113657
0.0065778
0.0037281
0.0020668
0.75
0.78
0.82
0.84
2−16
0.0190685
0.0113652
0.0065754
0.0037280
0.0020801
0.75
0.78
0.82
0.84
The values of ɛ for which we solve the test problem are ɛ = 2−i,i = 1,2,4,...,16. These convergence rates are increasing as N increases for any fixed ɛ. Tables 1 and 2 verify the ɛ-uniform convergence of the numerical solution on both subintervals and computed rates are essentially in agreement with our theoretical analysis.
Conclusions
A numerical method is developed to solve singularly perturbed initial value problem for semilinear second-order delay differential equation. This method is based on exponentially fitted difference scheme on an equidistant mesh, which gives first order uniform convergence in the discrete maximum norm. The method is shown to be uniformly convergent to the continuous solution i.e., independent of mesh parameter and perturbation parameter. Efficiency of the present method is demonsrated by a numerical example and also by comparing the results with exact solution of the problem. Presented numerical results are essentially in agreement with our theoretical analysis.