Original article Biomath 3 (2014), 1411041, 1–8 B f Volume ░, Number ░, 20░░ BIOMATH ISSN 1314-684X Editor–in–Chief: Roumen Anguelov B f BIOMATH h t t p : / / w w w . b i o m a t h f o r u m . o r g / b i o m a t h / i n d e x . p h p / b i o m a t h / Biomath Forum A Parameter Uniform Almost First Order Convergent Numerical Method for a Semi-Linear System of Singularly Perturbed Delay Differential Equations N. Shivaranjani ∗, J.J.H. Miller† and S.Valarmathi∗ ∗Department of Mathematics, Bishop Heber College, Tiruchirappalli, TamilNadu, India. Email: shivaranjaninagarajan@gmail.com, valarmathi07@gmail.com †Institute for Numerical Computation and Analysis, Dublin, Ireland. Email: jm@incaireland.org Received: 11 August 2014, accepted: 4 November 2014, published: 24 November 2014 Abstract—In this paper an initial value problem for a semi-linear system of two singularly perturbed first order delay differential equations is considered on the interval (0,2]. The components of the solution of this system exhibit initial layers at 0 and interior layers at 1. A numerical method composed of a classical finite difference scheme on a piecewise uniform Shishkin mesh is suggested. This method is proved to be almost first order convergent in the maximum norm uniformly in the perturbation parameters. Keywords-Singular Perturbation problems, boundary layers, semi-linear delay differential equations, finite difference schemes, Shishkin mesh, parameter uniform convergence. I. INTRODUCTION Singularly perturbed delay differential equations play an important role in the modelling of sev- eral physical and biological phenomena like first exit time problems in modelling of activation of neuronal variability [3], bistable devices [8] and evolutionary biology [6] and in a variety of models for physiological processes or diseases [9],[10] and [11]. These systems also find applications in Belousov- Zhabotinskii reaction (BZ reaction) models and the modelling of biological oscillators [6]. A model of tumor growth that includes the im- mune system response and a cycle-phase-specific drug presented in [13] is cited here. The model considers three populations: immune system, pop- ulation of tumor cells during interphase and pop- ulation of tumor during mitosis. The governing equations of the system are Citation: N. Shivaranjani, J.J.H. Miller, S.Valarmathi, A Parameter Uniform Almost First Order Convergent Numerical Method for a Semi-Linear System of Singularly Perturbed Delay Differential Equations, Biomath 3 (2014), 1411041, http://dx.doi.org/10.11145/j.biomath.2014.11.041 Page 1 of 8 http://www.biomathforum.org/biomath/index.php/biomath http://dx.doi.org/10.11145/j.biomath.2014.11.041 N. Shivaranjani et al., A Parameter Uniform Almost First Order ... dTI dt =2a4TM − (c1I + d2)TI −a1TI(t− τ) dTM dt =a1TI(t− τ) −d3TM −a4TM − c3TMI −k1(1 −e−k2u)TM dI dt =k + ρI(TI +TM ) n α+(TI +TM )n − c2ITI − c4TMI −d1I − k3(1 −ek4u)I du dt = −γu with TI(t) =φ1(t) for t ∈ [−τ, 0] TM (t) =φ2(t) for t ∈ [−τ, 0] I(t) =φ3(t) for t ∈ [−τ, 0] u(0) =u0. Here, TI(t) - population of tumor cells during inter- phase at time t TM (t) -population to tumor cells during mitosis at time t I(t) -population of immune system at time t u(t) -amount of drug present at time t τ -the resident time of cells in interphase d2TI,d3TM,d1I - proportions of natural cell death or apoptosis a1,a4 - the rate at which cells cycle are repro- duce ci -losses from encounters of tumor cells with immune cells ρI(TI +TM ) n α+(TI +TM )n - non-linear growth of the immune population due to stimulus by tumor cells k -constant rate at which the immune cells grow, in the absence of tumor cells ρ,α,n -parameters depending on the type of tumor being considered and the health of the immune system. Thus, an initial value problem for a system of semilinear delay differential equations is used to model tumor growth. Here, the parameters may take large values, for instance the value of k is 1.3 × 104 in the paper cited. In these cases, the system becomes singularly perturbed. Motivated by this, in this paper, the following semilinear system of singularly perturbed delay differential equations is considered: ~T~u = E~u′(x) + ~f(x,u1,u2) + B(x)~u(x− 1) = ~0 on (0, 2], ~u = ~φ on [−1, 0]. (1) For all x ∈ [0, 2], ~u(x) = (u1(x),u2(x))T and ~f(x,u1,u2) = (f1(x,u1,u2),f2(x,u1,u2))T . E, B(x) are 2 × 2 matrices. E = diag(~ε), ~ε = (ε1,ε2) with 0 < ε1 ≤ ε2 ≤ 1, B(x) = diag(~b), ~b = (b1(x),b2(x)). It is assumed that the nonlinear terms satisfy ∂fk(x) ∂uk ≥ β > 0, ∂fk(x) ∂uj ≤ 0, k,j = 1, 2, k 6= j (2) min 1≤i≤2   2∑ j=1 ∂fi(x) ∂uj + bi(x)   ≥ α > 0, (3) bi(x) ≤ 0, i = 1, 2 (4) for x in [0, 2] × C2 where C = C0([−1, 2]) ∩ C1((0, 2]) ∩C2((0, 1) ∪ (1, 2)). These conditions and the implicit function theorem ensure that a unique solution ~u ∈ C2 exists for the problem (1). The solution ~u(x) has initial layers at x = 0 and interior layers at x = 1. Both the components u1 and u2 have layers of width O(ε2) and the component u1 has an additional sublayer of width O(ε1). For any vector-valued function ~y on [0, 2] the following norms are introduced: ‖ ~y(x) ‖= maxi |yi(x)|, i = 1, 2 and ‖ ~y ‖= sup{‖ ~y(x) ‖: x ∈ [0, 2]}. A mesh Ω̄N = {xi}Ni=0 is a set of points satisfying 0 = x0 < x1 < ... < xN = 2. A mesh function V = {V (xi)}Ni=0 is a real valued function defined on Ω̄N . The discrete maximum norm for the above function is defined by ‖ V ‖Ω̄N = maxi=0,1,...,N |V (xi)| and ‖ ~V ‖Ω̄N = max{‖ V1 ‖Ω̄N ,‖ V2 ‖Ω̄N} where the vector mesh functions ~V = (V1,V2) T = {V1(xi),V2(xi)}, i = 0, 1, ..,N. Biomath 3 (2014), 1411041, http://dx.doi.org/10.11145/j.biomath.2014.11.041 Page 2 of 8 http://dx.doi.org/10.11145/j.biomath.2014.11.041 N. Shivaranjani et al., A Parameter Uniform Almost First Order ... Throughout the paper C denotes a generic positive constant, which is independent of x and of all singular perturbation and discretization parameters. Furthermore, inequalities between vectors are understood in the componentwise sense. II. ANALYTICAL RESULTS The problem (1) can be rewritten in the form ε1u ′ 1(x) + f1(x,u1,u2) + b1(x)φ1(x− 1) = 0 ε2u ′ 2(x) + f2(x,u1,u2) + b2(x)φ2(x− 1) = 0, x ∈ (0, 1] ~u(0) = ~φ(0) (5) and ε1u ′ 1(x) + f1(x,u1,u2) + b1(x)u1(x− 1) = 0 ε2u ′ 2(x) + f2(x,u1,u2) + b2(x)u2(x− 1) = 0, x ∈ (1, 2] ~u(1) known from (5). (6) ~T1~u := E~u ′(x) + ~g(x,u1,u2) = ~0, x ∈ (0, 1] ~T2~u := E~u ′(x) + ~f(x,u1,u2) +B(x)~u(x− 1) = ~0, x ∈ (1, 2] where ~g(x,u1,u2) = ~f(x,u1,u2) + B(x)~φ(x− 1). (7) The reduced problem corresponding to (7) is given by ~g(x,r1,r2) = ~0, x ∈ (0, 1] (8) ~f(x,r1,r2) + B(x)~r(x− 1) = ~0, x ∈ (1, 2]. (9) The implicit function theorem and conditions (2),(3) and (4) ensure the existence of a unique solution for (8) and (9). This solution ~r has derivatives which are bounded independently of ε1 and ε2. Hence, |r(k)1 (x)| ≤ C; |r (k) 2 (x)| ≤ C; k = 0, 1, 2, 3; x ∈ [0, 2]. The following Shishkin decomposition [1], [2] of the solution ~u is considered: ~u = ~v + ~w, where the smooth component ~v(x) is the solution of the problem E~v′(x) + ~g(x,v1,v2) = ~0, x ∈ (0, 1] E~v′(x) + ~f(x,v1,v2) + B(x)~v(x− 1) = ~0, x ∈ (1, 2] ~v(0) = ~r(0) (10) and the singular component ~w(x) satisfies E~w′(x) + ~g(x,v1 + w1,v2 + w2) −~g(x,v1,v2) = ~0,x ∈ (0, 1] E~w′(x) + ~f(x,v1 + w1,v2 + w2) − ~f(x,v1,v2) +B(x)~w(x− 1) = ~0, x ∈ (1, 2] ~w(0) = ~u(0) −~v(0). (11) The bounds of the derivatives of the smooth com- ponent are contained in Lemma 1: The smooth component ~v(x) satis- fies |v(i)k (x)| ≤ C, k = 1, 2; i = 0, 1 and |v′′k(x)| ≤ Cε −1 k , k = 1, 2. Proof: The smooth component ~v is further decomposed as follows: ~v = ~̃q + ~̂q where ~̂q is the solution of g1(x, q̂1, q̂2) = 0 (12) ε2 dq̂2 dx + g2(x, q̂1, q̂2) = 0, x ∈ (0, 1] (13) q̂2(0) = v2(0); q̂1(0) = v1(0) (14) and f1(x, q̂1, q̂2) + b1(x)q̂1(x− 1) = 0 (15) ε2 dq̂2 dx + f2(x, q̂1, q̂2) + b2(x)q̂2(x− 1) = 0, x ∈ (1, 2] (16) q̂2(1) and q̂1(1) are known from (12) and (13). ~̃q is the solution of ε1 dq̃1 dx + g1(x, q̃1 + q̂1, q̃2 + q̂2) −g1(x, q̂1, q̂2) = −ε1 dq̂1 dx (17) Biomath 3 (2014), 1411041, http://dx.doi.org/10.11145/j.biomath.2014.11.041 Page 3 of 8 http://dx.doi.org/10.11145/j.biomath.2014.11.041 N. Shivaranjani et al., A Parameter Uniform Almost First Order ... ε2 dq̃2 dx + g2(x, q̃1 + q̂1, q̃2 + q̂2) −g2(x, q̂1, q̂2) = 0, x ∈ (0, 1] q̃1(0) = q̃2(0) = 0 (18) and ε1 dq̃1 dx + f1(x, q̃1 + q̂1, q̃2 + q̂2) −f1(x, q̂1, q̂2) +b1(x)q̃1(x− 1) = −ε1 dq̂1 dx (19) ε2 dq̃2 dx + f2(x, q̃1 + q̂1, q̃2 + q̂2) −f2(x, q̂1, q̂2) +b2(x)q̃2(x− 1) = 0, x ∈ (1, 2] (20) q̃1(1) and q̃2(1) are known from (17) and (18). Let x ∈ [0, 1]. Using (8), (12) and (13), a11(x)(q̂1 −r1) + a12(x)(q̂2 −r2) = 0 (21) ε2 d dx (q̂2 −r2) + a21(x)(q̂1 −r1) +a22(x)(q̂2 −r2) = −ε2 dr2 dx , (22) where, aij(x) = ∂gi ∂uj (x,ξi(x),ηi(x)), i,j = 1, 2; ξi(x),ηi(x) are intermediate values. Using (21) in (22), ε2 d dx (q̂2 −r2) + ( a22(x) − a12(x)a21(x) a11(x) ) ×(q̂2 −r2) = −ε2 dr2 dx Consider the linear operator, l1(z) := ε2z ′ + ( a22(x) − a12(x)a21(x) a11(x) ) z = −ε2 dr2 dx , (23) where, z = q̂2 −r2. This operator satisfies the maximum principle [1]. Thus, ‖ q̂2 −r2 ‖≤ Cε2 and ‖ d(q̂2 −r2) dx ‖≤ C. Using this in (21), ‖ q̂1 −r1 ‖≤ Cε2. Hence, ‖ q̂2 ‖≤ C, ‖ dq̂2 dx ‖≤ C and ‖ q̂1 ‖≤ C. Differentiating (22), ε2 d2 dx2 (q̂2 −r2) + a′21(x)(q̂2 −r2) + a21(x) d dx (q̂2 −r2) + a′22(x)(q̂1 −r1) + a22(x) d dx (q̂1 −r1) = −ε2 d2r2 dx2 . (24) Hence,‖ d2q̂2 dx2 ‖≤ Cε−12 . Differentiating (21) twice and using the above estimates of d2q̂2 dx2 , ‖ d2q̂1 dx2 ‖≤ Cε−12 . (25) From (17) and (18), ε1 dq̃1 dx + a∗11(x)q̃1 + a ∗ 12(x)q̃2 = −ε1 dq̂1 dx (26) ε2 dq̃2 dx + a∗21(x)q̃1 + a ∗ 22(x)q̃2 = 0 (27) q̃1(0) = q̃2(0) = 0 (28) where, a∗ij(x) = ∂gi ∂uj (x,ζi(x),χi(x)), i,j = 1, 2; ζi(x),χi(x) are intermediate values. From equations (26) and (27), ‖ q̃i ‖≤ C, i = 1, 2 (29) ‖ dq̃i dx ‖≤ C, i = 1, 2 (30) ‖ d2q̃i dx2 ‖≤ Cε−1i , i = 1, 2. (31) Hence from the bounds for ~̃q and ~̂q, the required bounds of ~v follow. Biomath 3 (2014), 1411041, http://dx.doi.org/10.11145/j.biomath.2014.11.041 Page 4 of 8 http://dx.doi.org/10.11145/j.biomath.2014.11.041 N. Shivaranjani et al., A Parameter Uniform Almost First Order ... Let x ∈ [1, 2]. Using (9), (15) and (16), p11(x)(q̂1 −r1) + p12(x)(q̂2 −r2)+ b1(x)(q̂1(x− 1) + r1(x− 1)) = 0 (32) ε2 d dx (q̂2 −r2) + p21(x)(q̂1 −r1) +p22(x)(q̂2 −r2) + b2(x)(q̂2(x− 1) −r2(x− 1)) = −ε2 dr2 dx (33) where, pij(x) = ∂fi ∂uj (x,κi(x),λi(x)), i,j = 1, 2; κi(x),λi(x) are intermediate values. Using (32) in (33), ε2 d dx (q̂2 −r2) + ( p22(x) − p12(x)p21(x) p11(x) ) ×(q̂2 −r2) − p21p11 (x)b1(x)(q̂1(x− 1) −r1(x− 1)) + b2(x)(q̂2(x− 1) −r2(x− 1)) = −ε2 dr2 dx Consider the linear operator, l2(z) := ε2z ′ + ( p22(x) − p12(x)p21(x) p11(x) ) z +b2(x)z(x− 1) = −ε2 dr2 dx − p21 p11 (x)b1(x)(q̂1(x− 1) −r1(x− 1)), (34) where, z = q̂2 −r2. This operator satisfies the maximum principle [12]. Hence using similar arguments as in the interval [0, 1] and the bounds of ~̂q and ~̃q in the interval [0, 1], the required bounds in the interval [1, 2] are derived. Lemma 2: The singular component ~w(x) satis- fies, for any x ∈ [0, 1], |wi(x)| ≤ Ce −αx ε2 ; i = 1, 2 |w′1(x)| ≤ C(ε −1 1 e −αx ε1 + ε−12 e −αx ε2 ) |w′2(x)| ≤ Cε −1 2 e −αx ε2 |w′′i (x)| ≤ Cε −1 i (ε −1 1 e −αx ε1 + ε−12 e −αx ε2 ), i = 1, 2 For x ∈ [1, 2], |wi(x)| ≤ Ce −α(x−1) ε2 ; i = 1, 2 |w′1(x)| ≤ C(ε −1 1 e −α(x−1) ε1 + ε−12 e −α(x−1) ε2 ) |w′2(x)| ≤ Cε −1 2 e −α(x−1) ε2 |w′′i (x)| ≤ Cε −1 i (ε −1 1 e −α(x−1) ε1 +ε−12 e −α(x−1) ε2 ), i = 1, 2 Proof: From equations (11), ε1w ′ 1(x) + s11(x)w1(x) + s12(x)w2(x) = 0 (35) ε2w ′ 2(x) + s21(x)w1(x) + s22(x)w2(x) = 0, x ∈ (0, 1] (36) w1(0) = u1(0) −v1(0); w2(0) = u2(0) −v2(0) and ε1w ′ 1(x) + s ∗ 11(x)w1(x) + s ∗ 12(x)w2(x) +b1(x)w1(x− 1) = 0 (37) ε2w ′ 2(x) + s ∗ 21(x)w1(x) + s ∗ 22(x)w2(x)+ b2(x)w2(x− 1) = 0, x ∈ (1, 2] (38) w1(1) = u1(1) −v1(1); w2(1) = u2(1) −v2(1) Here, sij(x) = ∂gi ∂uj (x,νi(x),υi(x)) and s∗ij(x) = ∂fi ∂uj (x,φi(x),φ ∗ i (x)); νi(x),υi(x),φi(x), φ∗i (x) are intermediate values. Biomath 3 (2014), 1411041, http://dx.doi.org/10.11145/j.biomath.2014.11.041 Page 5 of 8 http://dx.doi.org/10.11145/j.biomath.2014.11.041 N. Shivaranjani et al., A Parameter Uniform Almost First Order ... From equations (35),(36),(37) and (38), the bounds of the singular component ~w can be derived as in [5] in the domains [0, 1] and [1, 2]. III. SHISHKIN MESH A piecewise uniform Shishkin mesh Ω N = Ω− N ∪ Ω+N where Ω− N = {xj} N 2 0 and Ω +N = {xj}NN 2 +1 with N mesh-intervals is now con- structed on Ω = [0, 2],as follows, for the case ε1 < ε2. In the case ε1 = ε2 a simpler con- struction requiring just one parameter τ suffices. The interval [0, 1] is subdivided into 3 sub- intervals [0,τ1] ∪ (τ1,τ2] ∪ (τ2, 1]. The parameters τr, r = 1, 2, which determine the points separat- ing the uniform meshes, are defined by τ0 = 0, τ3 = 1 2 , τ2 = min { 1 2 , ε2 α ln N } and τ1 = min {τ2 2 , ε1 α ln N } . (39) Clearly 0 < τ1 < τ2 ≤ 12. Then, on the sub- interval (τ2, 1] a uniform mesh with N4 mesh points is placed and on each of the sub-intervals (0,τ1] and (τ1,τ2], a uniform mesh of N8 mesh points is placed. Similarly, the interval [1, 2] is also divided into 3 sub-intervals [1, 1 + τ1], (1 + τ1, 1 + τ2], (1 + τ2, 2] having the same number of mesh intervals as in [0, 1]. Note that, when both the parameters τr, r = 1, 2, take on their lefthand value, the Shishkin mesh becomes a classical uniform mesh on [0, 2]. IV. DISCRETE PROBLEM The initial value problems (5) and (6) are dis- cretised using the backward Euler scheme on the piecewise uniform fitted mesh Ω̄N. The discrete problem is TN ~U(xj) := ED −~U(xj) + ~g(xj,U1(xj),U2(xj)) = 0, j = 1(1) N 2 (40) T̃N ~U(xj) := ED −~U(xj) + ~f(xj,U1(xj),U2(xj)) = −B(xj)~U(xj − 1), j = N 2 + 1(1)N (41) ~U(0) = ~u(0) and D−~U(xj) = ~U(xj) − ~U(xj−1) xj −xj−1 , j = 1(1)N. Lemma 3: For any mesh functions ~Y and ~Z with ~Y (0) = ~Z(0), ‖ ~Y − ~Z ‖≤ C ‖ TN ~Y −TN ~Z ‖ Proof: TN ~Y −TN ~Z = ED−~Y (xj) + ~g(xj,Y1(xj),Y2(xj)) −ED−~Z(xj) −~g(xj,Z1(xj),Z2(xj)) = ED−(~Y − ~Z)(xj) + ∂~g ∂u1 (xj,~ξ(xj),~η(xj))(Y1 −Z1) + ∂~g ∂u2 (xj,~ξ(xj),~η(xj))(Y2 −Z2) = (T ′N )( ~Y − ~Z) where T ′N is the Frechet derivative of TN and the notation ∂~g ∂ui (xj,~ξ(xj),~η(xj)), i = 1, 2 is used to express the difference between the mid-values for the components g1 and g2. Since T ′N is linear, it satisfies the discrete maximum principle and discrete stability result [5].Hence ‖ ~Y−~Z ‖≤ C ‖ T ′N (~Y−~Z) ‖= C ‖ TN ~Y−TN ~Z ‖ and the lemma is proved. Parameter - uniform bounds for the error are given in the following theorem, which is the main result of this paper. Theorem 1: Let ~u be the solution of the prob- lem (1) and ~U be the solution of the discrete problem (40),(41). Then ‖ ~U −~u ‖≤ CN−1 ln N (42) Biomath 3 (2014), 1411041, http://dx.doi.org/10.11145/j.biomath.2014.11.041 Page 6 of 8 http://dx.doi.org/10.11145/j.biomath.2014.11.041 N. Shivaranjani et al., A Parameter Uniform Almost First Order ... TABLE I Values of DNε , D N, pN, p∗ and CNp∗ for ε1 = η 16 , ε2 = η 4 and α = 0.9. η Number of mesh points N 128 256 · · · 8192 16384 20 0.150E-01 0.806E-02 · · · 0.271E-03 0.136E-03 2−3 0.211E-01 0.121E-01 · · · 0.619E-03 0.336E-03 2−6 0.218E-01 0.125E-01 · · · 0.619E-03 0.336E-03 2−9 0.218E-01 0.125E-01 · · · 0.619E-03 0.336E-03 2−12 0.218E-01 0.125E-01 · · · 0.619E-03 0.336E-03 ... ... ... · · · ... ... 2−27 0.218E-01 0.125E-01 · · · 0.619E-03 0.336E-03 DN 0.218E-01 0.125E-01 · · · 0.619E-03 0.336E-03 pN 0.800E+00 0.854E+00 · · · 0.880E+00 CNp 0.249E+01 0.249E+01 · · · 0.196E+01 0.186E+01 Computed order of ~ε -uniform convergence, p∗ = 0.8 Computed ~ε -uniform error constant, CNp∗ = 2.48 Proof: Let x ∈ [0, 1]. From the above lemma, ‖ ~U −~u ‖≤ C ‖ TN ~U −TN~u ‖ Consider ‖ TN~u ‖=‖ TN~u−TN ~U ‖ Hence, ‖ TN~u−TN ~U ‖=‖ TN~u ‖ =‖ TN~u− ~T1~u ‖ = E|(D−~u−~u′)(x)| ≤ E|(D−~v −~v′)(x)| +E|(D−~w − ~w′)(x)| Since the bounds for ~v and ~w are the same as in [5] , the required result follows. Let x ∈ [1, 2]. From the above lemma, ‖ ~U −~u ‖ ≤ C ‖ T̃N ~U − T̃N~u ‖ ≤ C ‖ B(xj)(~U −~u)(xj − 1) ‖ ≤ C ‖ ~U −~u ‖ ≤ CN−1 ln N V. NUMERICAL RESULTS The numerical method proposed in this paper is illustrated through an example presented in this section. Example Consider the initial value problem ε1u ′ 1(x) + 3u1(x) − 1 4 exp(−u21)(x) −u2(x) −x2 + 1 −u1(x− 1) = 0 ε2u ′ 2(x) + 4u2(x) − cos(u2(x)) −u1(x)− ex −u2(x− 1) = 0; x ∈ (0, 1] ~u(x) = ~0; x ∈ [−1, 0]. The above quasi linear problem is solved using the numerical method suggested in this paper utilising the continuation method found in [2]. The maximum pointwise errors and the rate of convergence for this IVP are calculated using the two - mesh algorithm in [2] and are presented in Table 1. The notations DN,pN,CNp ,C N p∗ and p ∗ bear the same meaning as in [2] but the methods to arrive at them are modified for the vector solution. A graph of the numerical solution is presented in Figure 1 for N = 2048 and η = 2−15. The sharper initial layers at x = 0 and interior layers at x = 1 are evident. Biomath 3 (2014), 1411041, http://dx.doi.org/10.11145/j.biomath.2014.11.041 Page 7 of 8 http://dx.doi.org/10.11145/j.biomath.2014.11.041 N. Shivaranjani et al., A Parameter Uniform Almost First Order ... Fig. 1. Numerical solution -2 0 2 4 6 8 10 12 14 0 0.5 1 1.5 2 u1 u2 ACKNOWLEDGMENT The first author wishes to acknowledge the financial assistance extended through INSPIRE fellowship by the Department of Science and Tech- nology, Government of India. REFERENCES [1] J. J. H. Miller, E. O’Riordan, G.I. Shishkin, Fitted Numer- ical Methods for Singular Perturbation Problems, World Scientific, Revised edition (2012). [2] P.A. Farrell, A. Hegarty, J. J. H. Miller, E. O’Riordan, G. I. Shishkin, Robust Computational Techniques for Boundary Layers, Applied Mathematics & Mathemati- cal Computation (Eds. R. J. Knops & K. W. Morton), Chapman & Hall/CRC Press (2000). [3] C.G.Lange, R.M.Miura, Singular perturbation analysis of boundary-value problems for differential-difference equations. SIAM J.Appl.Math.42(3),502-530(1982). http://dx.doi.org/10.1137/0142036 [4] Zhongdi Cen, A second-order hybrid finite difference scheme for system of singularly perturbed initial value problems, Journal of Computational and Applied Mathe- matics 234, 3445-3457 (2010). http://dx.doi.org/10.1016/j.cam.2010.05.006 [5] S.Valarmathi, J.J.H.Miller, A parameter uniform finite difference method for a singularly perturbed linear dy- namical systems,International Journal of Numerical Anal- ysis and Modelling, 7(3), 535-548 (2010). [6] J.D.Murray Mathematical Biology: An Introduction (Third Edition),Springer (2002) . [7] T.Linss and N.Madden, Accurate solution of a system of coupled singularly perturbed reaction-diffusion equations, Computing, vol 73, 121-133,(2004). [8] M.W.Derstine, H.M.Gibbs, F.A.Hopf and D.L.Kaplan, Bifurcation gap in a hybrid optically bistable system. Physical Review 26(6),3720-3722(1982). http://dx.doi.org/10.1103/PhysRevA.26.3720 [9] Rebecea V.Culshaw, Shigui Ruan A delay differential equation model of HIV infection of CD4+ T-Cells. Math- ematical Biosciences 165,27-39(2000). http://dx.doi.org/10.1016/S0025-5564(00)00006-7 [10] A.Longtin, J.G.Milton Complex oscillations in the hu- man pupil light reflex with mixed and delayed feedback. Mathematical Biosciences.90(1-2),183-199,(1988). http://dx.doi.org/10.1016/0025-5564(88)90064-8 [11] Patrick W.Nelson, Alan Perelson Mathematical analysis of delay differential equation models of HIV - 1 infection. Mathematical Biosciences.179,73-94,(2002). http://dx.doi.org/10.1016/S0025-5564(02)00099-8 [12] Zhongdi Cen A hybrid finite difference scheme for a class of singularly perturbed delay differential equations. Neural, Parallel and Scientific Computations .16,303- 308(2008). [13] Minaya Villasana, Ami Radunskaya A delay differential equation model for tumor growth. J.Math. Biol. 47,270- 294(2003)http://dx.doi.org/10.1007/s00285-003-0211-0 Biomath 3 (2014), 1411041, http://dx.doi.org/10.11145/j.biomath.2014.11.041 Page 8 of 8 http://dx.doi.org/10.1137/0142036 http://dx.doi.org/10.1016/j.cam.2010.05.006 http://dx.doi.org/10.1103/PhysRevA.26.3720 http://dx.doi.org/10.1016/S0025-5564(00)00006-7 http://dx.doi.org/10.1016/0025-5564(88)90064-8 http://dx.doi.org/10.1016/S0025-5564(02)00099-8 http://dx.doi.org/10.1007/s00285-003-0211-0 http://dx.doi.org/10.11145/j.biomath.2014.11.041 Introduction Analytical Results Shishkin Mesh Discrete Problem Numerical Results References